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

# 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 )
gamma <- beta / m

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

om  <- 1 * om0   # espressa in unità di omega0

T   <- 2*pi/om

t.max <-  2*T       # limiti per plot

t <- seq(0, t.max, len=101)
F <- f0 * cos(om*t)

par(mfrow= c(2,1) )
old.mar = par("mar")
par(mar=c(0.5,4.2,2.5,0.5))

titolo = sprintf("omega / omega0  =  %.2f", om / om0)
plot(t/T, F, ty='l', xlab='t / T', col='red',  xaxs = 'i',
     ylab = "F (N)", main=titolo)
text(1.8, -0.95, "t / T", cex=1.5)
grid()

#pausa()

kappa <- g0 / ( om0^2 - om^2 + 1i * gamma * om  )
zeta  <- kappa * exp(1i * om * t)
x     <-  Re(zeta)
par(mar=c(2.5,4.2,2.5,0.5))
plot(t/T, x / max(x), ty='l', col='blue',
     ylab = " x / xM;   v / vM", xaxs = "i", ylim=c(-1,1.1))
testo.x <- sprintf("x / xM  [con  xM = %.1e m ]", max(x))
text(0.5, 1.1, testo.x, col='blue', cex=1.5)
grid()


#pausa() 
zeta.dot <- 1i * om * zeta
v <- Re(zeta.dot)
points(t/T, v / max(v), ty='l', col='green')
testo.v <- sprintf("v / vM  [con  vM = %.1e m/s ]", max(v))
text(1.5, 1.1, testo.v, col='green', cex=1.5)


# risistemiamo i parametri per i plot
par(mfrow= c(1,1) )
par(mar=old.mar)