# variante di RLC_transiente.R per fare la scarica del condensatore
# -> modifiche indicate con '### Mod 1', etc

library(deSolve)     # install.packages("deSolve")

Vgen <- function(t, V0, nu)  rep(0, length(t))  ### Mod 1:  generatore spento (ovvero ddp=0) 


RLC <- function(t, y, par) {
  list(c(
     y[2],
     (Vgen(t,V0,nu) - tau*y[2]-y[1])* om2   
  ))  
}

V0 <- 1
R <- 50     # Ohm
C <- 20e-9  # F
L <- 0.01   # H
tau <- R*C
om2 <- 1/(L*C)
nu0 <- sqrt(om2)/(2*pi)
nu <- 1*nu0  # 0.1*mu0 # 10*nu0

Q = sqrt(L/C)/R
cat(sprintf("Q = %.3f\n", Q) )

T <- 1/nu
t.max = 20*T
t <- seq(0, t.max, by=T/1000)
yini <- c(y1 = V0, y2 = 0)                 ### Mod 2: V.C posto inizialmente a V0
par <- c(tau=tau, om2=om2, V0=V0, nu=nu)
out <-  ode(y = yini, func = RLC, times = t, parms = par)

V.C <- out[,2]
V.R <- out[,3]*tau
V.C <- out[,2]
V.R <- out[,3]*tau
V.L <-  Vgen(t,V0,nu) - V.R - V.C
V.max <- max(c(V.C,V.R,V.L))
V.min <- min(c(V.C,V.R,V.L))
plot(t, V.C, cex=0.1, xlab='t (s)', ylab='tensioni (V)',
     ylim= c(V.min, V.max), col='red')
abline(h=0)
grid()
points(t, V.R, cex=0.1, col='blue')
points(t, Vgen(t, V0, nu), cex=0.1, col='black')
points(t, V.L, cex=0.1, col='green')
points(t, V0*exp(-t/(2*L/R)), ty='l', lty=2)