#----------------------------------------------------------
# nuova versione, ma con resistenza del generatore
# e codice scritto in modo meno criptico
#
# GdA 15/4/2015
#----------------------------------------------------------

library(deSolve)

Vgen <- function(t, V0, nu, phi=0, off=0, oq=FALSE) {
   V <-  V0 * cos((2*pi*nu)*t + phi)
   if(oq) V <- sign(V)
   V <- V + off
   return(V)
}

RC <- function(t, state, par) {
  with( as.list(c(state, par)),
         { V <-  Vgen(t, V0, nu, phi, off, oq)
           dVc.dt <- ( V - Vc ) / tau 
           return( list( c(dVc.dt) ) )
         }
      )
}

#--- parametri generatore di onde
nu.su.nuT <- 1     # frequenza in unità di frequenza di taglio
R0 <- 50             # resistenza interna (Ohm)
V0 <- 1              # Ampiezza (V)
phi0 <- -pi/2        # -pi/4; 0; pi/2; -pi/2
off <- 0             #   offset (V)
oq <- TRUE           # onda quadra?

#--- parametri circuito -----
R  <- 300    # Ohm
C  <- 10e-9  # F
Vc0 <- 0     # Tensione iniziale del condensatore (V)

#--- grandezze di  interesse  
Rt <- R0 + R
tau <- Rt*C
nuT <- 1/(2*pi*tau)
nu <- nu.su.nuT*nuT 
T <- 1/nu

# inputs per ode()  [par viene passato da ode() a RC() ]
t <- seq(0, 3*T, by=T/500)
par <- c(tau=tau, V0=V0, nu=nu, phi=phi0, off=off, oq=oq)
state <- c(Vc=Vc0)    

out <- ode(y=state, times=t, func=RC, parms=par)

Vlim <- c(min(-(1+oq)*V0, -V0 + off), max((1+oq)*V0, V0+off) ) # range verticale
plot(out, col='blue', ylim=Vlim, xlab='t (s)', ylab='tensioni',
      main=sprintf("nu = %.1f nuT", nu.su.nuT) )  # Vc
legend("topright", legend=c('Vg', 'Vout', 'Vc', 'Vr', 'V(r+r0)'),
       lty=c(3,1,1,1,1), 
       col=c('black', 'cyan', 'blue', 'magenta', 'brown'))
#legend("topleft", legend=c('Vg', 'Vout', 'Vc', 'Vr', 'V(r+r0)'),
#       lty=c(3,1,1,1,1), 
#       col=c('black', 'cyan', 'blue', 'magenta', 'brown'))
abline(h=0, lty=2)
V <- Vgen(t, V0, nu, phi0, off, oq)
points(t, V, ty='l', lty=3)           # Vgen
Vc <- out[,2] 
I = (V-Vc)/Rt                    # corrente
points(t, I*Rt, ty='l', col='brown')  # somma cadute di tensioni su resistenze
points(t, I*R, ty='l', col='magenta') # Vr
points(t, Vc+I*R, ty='l', col='cyan')    # Vout