#--------------------------------------
#  oscillatore forzato 
#   (per la teoria si veda oscillatore_forzato.pdf)
#
# GdA, maggio 2025
#--------------------------------------

pausa <- function() { cat ("\n >> press Enter to continue\n"); scan() }

# per la molla usiamo gli stessi dati di 'oscillazioni_smorzate_molla.R'
# (che poi saranno cambiati)
k    = 50      # (N/m)    costante della molla
m    = 0.5     # kg       massa sospesa
beta = 1     # N/(m/s ) coefficiente forza di resistenza aria

# parametri equazione differnzale
om0   <- sqrt( k/m )
T0    <- 2*pi/om0
gamma <- beta / m

# parametro termine forzante
f0  <- 1           # N 
g0  <- f0 / m

# om <- seq(0.5, 1.5, len=201) * om0
# om <- seq(0.01, 10, len=501) * om0
om <- seq(0.1, 10.0, len=501) * om0


# valori di omega equispazioati spaziati logaritmicamente
#om <- 10^seq(-1, 1, len=101) * om0


par(mfrow= c(2,1) )
old.mar = par("mar")
par(mar=c(1.5,4.2,0.5,0.5))
kappa <- g0 / ( om0^2 - om^2 + 1i * gamma * om  )
plot(kappa, col='blue', asp='1')
text(0.07, -0.097, 'Re( kappa)', cex=1.5 )
grid()
pausa()

par(mar=c(2.5,4.2,0.5,0.5))
plot(om/om0, Re(kappa), ty='l', col='blue')
grid()
text(8.5, -0.043, 'omega / omega0', cex=1.5 )

pausa()


plot(om/om0, abs(kappa), ty='l', col='blue', log='x')
grid()
text(5, 0.01, 'omega / omega0', cex=1.5 )
pausa()

plot(om/om0, Arg(kappa)/(pi/2), ty='l',  col='blue', log='x')
grid()
text(5, -1.8, 'omega / omega0', cex=1.5 )

par(mar=old.mar)
par(mfrow= c(1,1) )