#-----------------------------------------
# forme d'onda circuito RC 
#
#  GdA dicembre 2021
#-----------------------------------------

Vc <- function(nu, nut, V0, t) {
  return( V0/sqrt(1+(nu/nut)^2 ) * cos( t*2*pi*nu - atan(nu/nut) ) )
}
Vr <- function(nu, nut, V0, t) {
  return( V0/sqrt(1+(nu/nut)^(-2) ) * cos( t*2*pi*nu - atan(nu/nut) + pi/2) ) 
}

r  = 10000 # Ohm
c  = 10e-9 # Farad
V0 = 1     # Volt   (valore numerico 'fisicamente irrilevante', in quanto gioca
           #          il ruolo di 'fattore di scala' delle ampiezze di tensione')         

Deltat.su.T = 3 # finestra temporale da guardare
nut = 1/(2*pi*r*c)

cat(sprintf("\n Valori elementi circuitali: R = %.0f Ohm, C = %.0e F\n", r, c))
while (1) {
  # loop infinito, finche'  nusunut > 0
  cat ("\n >> dai valore frequenza in unità di frequenza di taglio ( < 0 per terminare )\n")
  nur = scan("",nmax=1)
  if (nur < 0) break
  nu = nur*nut
  # stampe 
    cat(sprintf("V_R:  Ampiezza/V0 %.3f,  Fase/pi %.3f \n",                
                V0/sqrt(1+(nu/nut)^(-2) ) , ( - atan(nu/nut) + pi/2) / pi ))
    cat(sprintf("V_C:  Ampiezza/V0 %.3f,  Fase/pi %.3f \n",                
                 V0/sqrt(1+(nu/nut)^2 ) ,  - atan(nu/nut) / pi ))
  # plot ------------------------------  
  T <- 1/nu
  t <- seq(0, T*Deltat.su.T, len=500)
  titolo = sprintf("%.0f Hz  (T = %.2e s)", nu, T)
  plot (t, V0*cos(t*2*pi*nu), ty='l', xlab='t (s)', ylab='V (V)', main=titolo)
  points(t, Vc(nu, nut, V0, t), ty='l', col='red')
  points(t, Vr(nu, nut, V0, t), ty='l', col='blue')       
}