
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 + off
}

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

#--- parametri circuito -----
R  <- 10000  # Ohm
C  <- 10e-9  # F
nu0 <- 1/(2*pi*R*C)

#--- parametri generatore di onde
nu <- 10*nu0   # Hz
T <- 1/nu
V0 <- 1
phi0 <- 0 # -pi/2 # -pi/4; 0; pi/2; -pi/2
off <- 1
oq <- !TRUE   # onda quadra?

par <- c(tau=tau, V0=V0, nu=nu, phi=phi0, off=off, oq=oq)
state <- c(Vc=0)
t <- seq(0, 5*T, by=T/1000)
out <- ode(y=state, times=t, func=RC, parms=par)

Vlim <- c(min(-(1+oq)*V0, -V0 + off), max((1+oq)*V0, V0+off) )
plot(out, col='red', ylim=Vlim, xlab='t (s)', ylab='V, Vc, Vr',
     main=sprintf("nu = %.1f nu0", nu/nu0))
legend("topright", legend=c('Vg', 'Vc', 'Vr'), lty=1, 
       col=c('blue', 'red', 'green'))
legend("topleft", legend=c('Vg', 'Vc', 'Vr'), lty=1, 
       col=c('blue', 'red', 'green')) 
abline(h=0, lty=2)
V <- Vgen(t, V0, nu, phi0, off, oq)
points(t, V, ty='l', col='blue')
points(t, V-out[,2], ty='l', col='green')

