#----------------------------------------------
# potenza in regime sinusoidale su Z arbitraria
#
# GdA dicembre 2021
#----------------------------------------------


nu    = 50       # Hertz (tanto per metterci un valore)
T     = 1/nu     # periodo
omega = 2*pi*nu  # pulsazione
V0    = 1.       # tensione unitaria

t = seq(-0.5*T,2.5*T, len=300)

while (1) {
  # loop infinito, finche'  nusunut > 0
  cat ("\n >> dai modulo di Z (in Ohm; niente per terminare) e fase (gradi)\n")
  z.pol = scan("",nmax=2)
  if (length(z.pol) < 2) break
#  if (length(z.pol) == 1) z.pol[2] = 0  
  Z = z.pol[1]*cos(z.pol[2]*pi/180) + 1i*z.pol[1]*sin(z.pol[2]*pi/180)
  v = V0*exp(1i*omega*t)
  i = v/Z
  p = v*i
  V = Re(v)
  I = Re(i)
  P = V*I
#  P.mean = sum(P[t<=T])*(t[2]-t[1])/T * (length(t)-1)/length(t) # valore approsimativo
  P.mean.th = (V0^2/2)/abs(Z) *Re(Z)/abs(Z)
  cat( sprintf("Potenza media %.2e W\n", round(P.mean.th, 13) ) )
    
  layout(matrix(c(1:3),c(3,1))) # tre plot sovrapposti
# plot tensione ---------------------------------------------------------    
  titolo = sprintf('Tensione [ %.0f Hz,  T = %.3f s]', nu, T)
  plot(t, V,  ty='l', xlab='t (s)', ylab='f (V)', main=titolo)
  abline(v=0, lty=2);    abline(v=T, lty=2);  abline(v=2*T, lty=2)
  abline(h=0, lty=2)
    
# plot corrente ---------------------------------------------------------
  titolo = sprintf("Corrente,    Z = %.0f Ohm, %.0f gradi      [ ( %.1f + %.1f j ) Ohm ]",
                   z.pol[1], z.pol[2], Re(Z), Im(Z))
  plot(t, I,  ty='l', xlab='t (s)', ylab='I (A)', col='blue', main=titolo) 
  abline(v=0, lty=2, col='blue');
  abline(v=T, lty=2, col='blue');
  abline(v=2*T, lty=2, col='blue')
  abline(h=0, lty=2, col='blue')
    
# plot potenza ---------------------------------------------------------
  titolo = sprintf("Potenza,    Z = %.0f Ohm, %.0f gradi      [ ( %.1f + %.1f j ) Ohm ]",
           z.pol[1], z.pol[2], Re(Z), Im(Z)) 
  plot(t, P,  ty='l', xlab='t (s)', ylab='P (W)', col='red', main=titolo)
  abline(v=0, lty=2, col='red');
  abline(v=T, lty=2, col='red');
  abline(v=2*T, lty=2, col='red')   
  abline(c(0,0), lty='dashed')
  abline(P.mean.th, 0, col='red')
  layout(1)
}