#---------------------------------------
# circuito RCL impulsato
#
# GdA 21/5/2011
#--------------------------------------


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


tensioni <- function(t, gamma, omega0, Vc0, f) {
  v <- NULL
  hgamma=gamma/2 + 01
  discr = (hgamma)^2 - omega0^2 +0i
  alpha1 = -hgamma - sqrt(discr)
  alpha2 = -hgamma + sqrt(discr)
  k1 = alpha2/(alpha2-alpha1)*(Vc0-f)
  k2 = alpha1/(alpha1-alpha2)*(Vc0-f)
  v$vc = Re( k1*exp(alpha1*t) + k2*exp(alpha2*t) ) + f
  v$vr = r*c * Re( alpha1*k1*exp(alpha1*t) + alpha2*k2*exp(alpha2*t) )
  return(v)
}

#--  valori dei componenti -------------------------
r <- 100
c  <- 4.7e-9   # Farad
l  <- 0.010    # Henry 
v0 <- 1      # tensione unitaria

omega0 = 1/sqrt(l*c)
gamma = r/l
hgamma=gamma/2 
tau = 2/gamma

tmax = 6*tau

t <- seq(0, tmax, len=2000)  # valori di tempi

# da su a giu'
v0 <- tensioni(t, gamma, omega0, 1, 0)
i <- v0$vr/r
p <- r*i^2
E <-  0.5*l*i^2 + 0.5*c*(v0$vc^2)

# da giu' a su
v1 <- tensioni(t, gamma, omega0, 0, 1)
i1 <- v1$vr/r
p1 <- r*i1^2
E1 <-  0.5*l*i1^2 + 0.5*c*(v1$vc^2)


# plot ---------------------------------------------------------------------
plot(t, v0$vc, ty='l', xlim=c(0,2*tmax), ylim=c(min(min(c(v0$vc,v1$vc)),0), max(c(v0$vc,v1$vc))),
     xlab="t (s)", ylab="Vc")
points(t+tmax, v1$vc, ty='l')
pausa()
plot(t, v0$vr, ty='l', xlim=c(0,2*tmax), ylim=c(min(min(c(v0$vr,v1$vr)),0), max(c(v0$vr,v1$vr))),
     xlab="t (s)", ylab="Vr" )
points(t+tmax, v1$vr, ty='l')
pausa()
plot(t, p, ty='l', xlim=c(0,2*tmax), ylim=c(min(c(p,p1)), max(c(p,p1))), 
     xlab="t (s)", ylab="Potenza dissipata (W)")
points(t+tmax, p1, ty='l')
pausa()
plot(t, E, ty='l',xlim=c(0,2*tmax), ylim=c(min(c(E,E1)), max(c(E,E1))),
     xlab="t (s)", ylab="Energia totale (J)")
points(t+tmax, E1, ty='l')

pausa()
old.mar = par("mar")
par(mar=c(2.5,4,1,1))  # cambia i margini dei plot
layout(matrix(c(1:2),c(2,1)))
plot(t, v0$vc, ty='l', xlim=c(0,2*tmax), ylim=c(min(min(c(v0$vc,v1$vc)),0), max(c(v0$vc,v1$vc))),
     xlab="t (s)", ylab="Vc")
points(t+tmax, v1$vc, ty='l')
plot(t, v0$vr, ty='l', xlim=c(0,2*tmax), ylim=c(min(min(c(v0$vr,v1$vr)),0), max(c(v0$vr,v1$vr))),
     xlab="t (s)", ylab="Vr" )
points(t+tmax, v1$vr, ty='l')
layout(1)
par(mar=old.mar)  # rimette i margini standard