#---------------------------------------------
# moto dell'oscillatore RLC smorzato nello
# spazio delle fasi I Vs Vc
# 
# GdA 29/4/2015
#----------------------------------------------

osc.smorz.ode <- function(t, state, par) {
  with( as.list(c(state, par)),
        { dVc.dt <- vVc                       # derivata di Vc
          dvVc.dt <- -gamma*vVc - omega2*Vc   # derivata della derivata di Vc
          return( list(c(dVc.dt, dvVc.dt)) )
        }
      )
}

# parametri del circuito
L      <- 0.01   # H
C      <- 4.7e-9 # F
R      <- 100    # Ohm

# parametri dell'oscillatore smorzato
beta   <- 0.9
gamma  <- R/L
omega2 <- 1/(L*C)
omega  <- sqrt(omega2)

# range e intervallo di tempi per la soluzione numerica
T0 <- 2*pi/omega
t.max <- 20 * T0
dt    <- T0/500
t <- seq(0, t.max, by=dt)

# parametri da passare a osc.smorz.ode()
par   <- c(gamma=gamma, omega2=omega2)

# Stato iniziale: condensatore carico, corrente nulla
# le variabili sono Vc e la sua "velocità" vVc (-> derivata)
Vc0=1
state <- c(Vc=Vc0, vVc=0)  

library(deSolve)
out <- ode(y = state, times = t, func = osc.smorz.ode, parms = par)

# estraiamo i vettori di interesse
Vc <- out[,'Vc']     # tensione del condensatore
I  <- out[,'vVc']*C  # corrente:  C * dVc/dt
I.max.ideale <- omega*C*Vc0  # ampiezza di oscillazione di I se R=0

plot(Vc/Vc0, I/I.max.ideale, ty='l',
     xlim=c(-1.05,1.1), ylim=c(-1.0,1.0),
     col='blue', asp=1, xlab='Vc/Vc0', ylab='I/I_Max-Ideale')
abline(h=0, lty=2)
abline(v=0, lty=2)
points(1,0, pch=19, col='blue', cex=0.8) 
text(1.0, 0.1, "inizio")
x <- seq(-1,1,by=0.01)
points(x,sqrt(1-x^2), ty='l', lty=2, col='green')
points(x,-sqrt(1-x^2), ty='l', lty=2, col='green')
text(0.7, 0.9, "caso ideale", col='green')