#---------------------------------------------------------- 
# Transiente nell'RCL risolvendo numericamente l'equazione
# differenziale mediante deSolve
#
# GdA ottobre 2013
#---------------------------------------------------------- 


pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

library(deSolve)

tr.RCL <- function(t,y,params) {
  retv <- {c(y[2],
             -params[1]*y[2] - params[2]*y[1] +
             params[2]*params[3]*ifelse(t<=params[5], cos(params[4]*t), 0)
             )}
  list(dy=retv,global=c())                                        
}

R=160
L=0.01
C=10e-9
tau=R*C
gamma <- R/L
omega0 <- 1/sqrt(L*C)
nu0=omega0/(2*pi)

V0 = 1
nu = nu0
omega <- 2*pi*nu
T <- 2*pi/omega

t <- seq(0, 30*T, len=1000)
t.max = 20*T
phi.gen <- -pi/2                    # fase del gen per passare da cos a sin 
V.gen <- V0*cos(omega*t + phi.gen)   
V.gen[t>t.max] <- 0
risposta <- lsoda(c(0, 0), t, tr.RCL, c(gamma, omega0^2, V0, omega, t.max))

#plot(risposta, main='transiente RCL ', xlab='t', col='red')
#points(t, V0*cos(omega*t), ty='l', col='blue')

#pausa()
layout(matrix(1:2, c(2,1)))
plot(risposta[,1], risposta[,2], ty='l', col='red', lwd=1.5,
     xlab='t (s)', ylab='Vc (V)')
#      ylim=c(min(min(V.gen),min(risposta[,2])), max(max(V.gen),max(risposta[,2])) ))
points(t, V.gen, ty='l', col='blue')
plot(risposta[,1], risposta[,3]*tau, ty='l', col='green', lwd=1.5,
      xlab='t (s)', ylab='Vr (V)',ylim=range(V.gen))
points(t, V.gen, ty='l', col='blue')
layout(1)