#----------------------------------------------------------
# Integratore e derivatore
#
# GdA 22/4/2015
#----------------------------------------------------------

library(deSolve)

# riscriviamo Vgen 
Vgen <- function(t, V0, t.max)  {
   V <- rep(0, length(t) ) 
   for (i in 1:length(t) ) {  # deve gestire vettori negli if!!
     if (t[i] < t.max/4) {
       V[i] <- V0*(t[i]/t.max)  
     }  else if (t[i] < t.max/2) {
       V[i] <- V0/4      
     }  else if(t[i] < t.max*3/4) {
       V[i] <- V0/4 + V0*3/4* ( (t[i]-t.max/2)/(t.max/4) )^2  
     } else {
       V[i] <- V0 - V0/(t.max/4)*(t[i]-t.max*3/4)   
     }
   }
   return(V)
}

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

#--- parametri generatore di onde
# Integratore: nu.su.nuT <- 20    (frequenze elevate!)
# Derivatore: nu.su.nuT <- 1/20   (frequenze basse) 
nu.su.nuT <- 20     # frequenza in unità di frequenza di taglio
V0 <- 1              # Ampiezza (V)

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

#--- grandezze di  interesse  
tau <- R*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, T, by=T/1000)
t.max <- max(t)
par <- c(tau=tau, V0=V0, t.max=t.max)
state <- c(Vc=Vc0)    

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

# aggiustiamo scale a seconda di integratore o derivatore
leg.C <- 'Vc * 10'
leg.R <- 'Vr'
scale.C <- 10
scale.R <- 1
Vlim <- c(0, V0)
main <- 'RC integratore'
if (nu.su.nuT < 1) {
  leg.C <- 'Vc'
  leg.R <- 'Vr * 10'
  scale.C <- 1
  scale.R <- 10
  Vlim <- c(-V0/3, V0) # (andrebbe fatto un po' più professionale...)
  main <- 'CR derivatore'
}


plot(t, Vc*scale.C, ty='l', col='blue', ylim=Vlim,
     xlab='t (s)', ylab='tensioni', main=main )  # Vc
legend("topleft", legend=c('Vg', leg.C, leg.R),
       lty=c(3,1,1), 
       col=c('black', 'blue', 'magenta'))
abline(h=0, lty=2)
abline(v=0, lty=2)
abline(v=t.max/4, lty=2)
abline(v=t.max/2, lty=2)
abline(v=t.max*3/4, lty=2)
abline(v=t.max, lty=2)
V <- Vgen(t, V0, t.max)
points(t, V, ty='l', lty=3)                      # Vgen
Vc <- out[,2] 
I = (V-Vc)/R                   # corrente
points(t, I*R*scale.R, ty='l', col='magenta')     # Vr