#---------------------------------------
# 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) {
  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
  k2 = alpha1/(alpha1-alpha2)*Vc0
  v$vc = Re( k1*exp(alpha1*t) + k2*exp(alpha2*t) )
  v$vr = r*c * Re( alpha1*k1*exp(alpha1*t) + alpha2*k2*exp(alpha2*t) )
  return(v)
}

#------------ leggi valore resistenza ---------------------------------
#  Nota: se il valore e' negativo, il suo opposto e'
#        pari al rapporto fra resistenza totale e resistenza critica
#        Es: -0.1 produce oscillazioni sottosmorzate;
#            -2   produce oscillazioni sovrasmorzata
#
if(!exists("rr")) r.in0 = 10
cat ("\n >> Dai R (o frazione di 'R critica' se valore < 0) \n
      \n (default=usa il valore precedente)\n")
r.in = scan("",nmax=1)
if (length(r.in) == 0) r.in = r.in0

#------------- leggi fattore scala temporale ---------------------------------
cat ("\n >> Dai fattore di scala temporale (default: 1)\n")
tsf = scan("",nmax=1)
if (length(tsf) == 0) tsf = 1
#-----------------------------------------------------------------------

# parametri del circuito
r0 <- 50    # Ohm
#r  <- 10      # Ohm
c  <- 4.7e-9   # Farad
l  <- 0.010    # Henry 
rl <- 100    # Ohm
v0 <- 1      # tensione unitaria

r.c <- 2*sqrt(l/c)   # resistenza critica (discr=0)

if (r.in >= 0) {
  r = r.in
  rt <- r0 + r + rl
} else {
  rt  <- r.c * abs(r.in) 
  r   <-  max(rt - r0 - rl, 1) # vogliamo almeno 1 Ohmper avere Vr!=0
}
cat( sprintf("Resistenza: %f  (Res. tot. %f)\n",r,rt) )
cat( sprintf("[Resistenza critica: %f]\n",r.c) )
#rt <- r0 + r + rl

#r.c = 2*sqrt(l/c) - r0 - rl
#rt <- r0 + r.c + rl

omega0 = 1/sqrt(l*c)
gamma = rt/l
hgamma=gamma/2 
discr = hgamma^2 - omega0^2

tmax = 6*(2/gamma) * tsf

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

v <- tensioni(t, gamma, omega0, 1)
i <- v$vr/r

plot(t, v$vc, ty='l', ylim=c(min(min(v$vc),0), max(v$vc)),
     xlab="t (s)", ylab="Vc")
pausa()
plot(t, v$vr, ty='l', ylim=c(min(min(v$vr),0), max(v$vr)),
     xlab="t (s)", ylab="Vr" )
pausa()
plot(t, rt*i^2, ty='l', xlab="t (s)", ylab="Potenza dissipata (W)")
pausa()
plot(t, 0.5*l*i^2 + 0.5*c*(v$vc^2), ty='l',
     xlab="t (s)", ylab="Energia totale (J)")