#---------------------------------------
# circuito RCL sinusoidale
#
# GdA 18/5/2011
#--------------------------------------


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

# funzioni per impedenze
zRCL <- function(r0=50, r, c, l, rl=0, nu) r0 + r + rl +zC(c, nu) +zL(l,nu)
zC <-  function(c, nu) -1i/(2*pi*nu*c)
zL <-  function(l, nu) 1i*2*pi*nu*l

# 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

nu0 <- 1/(2*pi*sqrt(l*c))  # frequenza di risonanza

nu <- seq(10000, 50000, len=200)  # valori di frequenze

# grandezze complesse:
i<-v0/zRCL(r0,r,c,l,rl, nu)
vr<-i*r
vc<-i*zC(c, nu)
vl<-i*zL(l,nu)

# moduli e fasi
Vr<-abs(vr)
phi.r<-Arg(vr)
Vc<-abs(vc)
phi.c<-Arg(vc)
Vl<-abs(vl)
phi.l<-Arg(vl)

# plot vari 
layout(matrix(1:3, 3,1))   # tre plot sovrapposti
plot(nu,Vr,type="l")
plot(nu,Vc,type="l", col='red')
plot(nu,Vl,type="l", col='blue')
layout(1)                  # rimettiamo a posto il layout

pausa()
vmax = max(Vr, Vc, Vl)*1.01
plot(nu,Vr,type="l", ylim=c(0,vmax))
points(nu,Vc,type="l", col='red')
points(nu,Vl,type="l", col='blue')
points(c(nu0,nu0),c(0,vmax), ty='l', col='gray')
legend("topright", lty=1, c("R", "C", "L"), col=c("black","red","blue"))

pausa()
layout(matrix(1:3, 3,1))   # tre plot sovrapposti
plot(nu,phi.r,type="l")
plot(nu,phi.c,type="l", col='red')
plot(nu,phi.l,type="l", col='blue')
layout(1)                  # rimettiamo a posto il layout

pausa()
plot(nu,phi.r,type="l",ylim=c(-pi,pi))
points(nu,phi.c,type="l", col='red')
points(nu,phi.l,type="l", col='blue')
points(c(min(nu),max(nu)), c(0,0), ty='l', col='gray')
points(c(nu0,nu0),c(-pi,pi), ty='l',col='gray')
legend("topright", lty=1, c("R", "C", "L"), col=c("black","red","blue"))

# esempio di come si salva un plot:
salva=FALSE  # sostituire con TRUE (o !FALSE) se si vuole salvare su file
if(salva) {
  png("RCL_ampiezze.png", width=800, height=600)
  plot(nu,Vr,type="l", ylim=c(0,vmax))
  points(nu,Vc,type="l", col='red')
  points(nu,Vl,type="l", col='blue')
  points(c(nu0,nu0),c(0,vmax), ty='l', lty='dotted',col='gray')
  legend("topright", lty=1, c("R", "C", "L"), col=c("black","red","blue"))
  dev.off()

  png("RCL_fasi.png", width=800, height=600)
  plot(nu,phi.r,type="l",ylim=c(-pi,pi))
  points(nu,phi.c,type="l", col='red')
  points(nu,phi.l,type="l", col='blue')
  points(c(min(nu),max(nu)), c(0,0), ty='l', col='gray')
  points(c(nu0,nu0),c(-pi,pi), ty='l',col='gray')
  legend("topright", lty=1, c("R", "C", "L"), col=c("black","red","blue"))
  dev.off()
}