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

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

Vgen <- function(t, V0, nu) V0 *sin( (2*pi*nu)*t)

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 = 0, y2 = 0)
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, Q*V0*(1-exp(-t/(2*L/R))), ty='l', lty=2)

pausa()
E.C <- 1/2 * C * V.C^2
E.L <- 1/2 * L * (V.R/R)^2
plot(t, E.C, cex=0.1, xlab='t (s)', ylab='Energia (J)', col='red')
grid()
points(t, E.L, cex=0.1, col='green')
points(t, E.C+E.L, cex=0.1, col='gray')

pausa()
P.gen <- Vgen(t, V0, nu) * V.R/R
P.R   <- V.R^2/R
plot(t, P.gen, cex=0.1, xlab='t (s)', ylab='P (W)')
grid()
points(t, P.R, cex=0.1, col='blue')