#--------------------------------------- # 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)")