#-------------------------------------------------
#  Oscillazioni smorzate: 
#  Applicazione di oscillazioni_smorzate.R
#  al caso della molla (lezione 11, probl. 10.4)
#
#   GdA maggio 2025
#--------------------------------------------------

# dati del sistema fisico
k    = 50      # (N/m)    costante della molla
m    = 0.5     # kg       massa sospesa
beta = 1.0     # N/(m/s ) coefficiente forza di resistenza aria
xM   = 0.05    # m        allungamento massimo (a t=0)

t.max <- 10*T0   # max scala temporale dei plot

# parametri equazione differnzale
om0   <- sqrt( k/m )
T0    <- 2*pi/om0
gamma <- beta / m
if(abs(gamma/2 - om0) < 0.0001) {
    gamma <- gamma * 1.0001        # protezione contro divergenze 
}

cat(sprintf("\nomega_0 = %.1f s^-1,  T0 = %.3f s,  gamma = %.2f s^-1\n",
            om0, T0, gamma))



alpha1 <- -gamma/2 + sqrt((gamma/2)^2 - om0^2 + 0i)   #  (10.84) 
alpha2 <- -gamma/2 - sqrt((gamma/2)^2 - om0^2 + 0i)   
# notare  '+0i': => variabili complesse 

cat("\nalpha1 e alpha2: \n")
print(c(alpha1, alpha2))

k1 <- alpha2 / (alpha2-alpha1) * xM                   #  (10.90)
k2 <- alpha1 / (alpha1-alpha2) * xM                   #  (10.91)

cat("\nk1 e k2: \n")
print(c(k1, k2))

# soluzione complessa e reale per 0 <= t <= t.max
t <- seq(0, t.max, len=501)                      
z <- k1*exp(alpha1*t) + k2*exp(alpha2*t)
x <- Re(z)

# parametri per i plot 
par(mfrow= c(2,1) )
old.mar = par("mar")
par(mar=c(2,4.2,2.5,0.5))

# plot di x(t)
titolo <- sprintf("omega_0 / (gamma/2): %.3f", om0/ (gamma/2))
plot(t, x, ty='l', col='blue', ylab='x (t)',  xaxs = "i", main=titolo)
grid()
text(0.9*t.max, 0.9*min(x), 't (s)', cex=1.5)

# dz/dt e sua parte reale  -> v(t)
z.dot <- alpha1*k1*exp(alpha1*t) + alpha2*k2*exp(alpha2*t)
v <- Re(z.dot)

# plot di v 
par(mar=c(2,4.2,0.5,0.5))    # parametri leggermente diversi per il secondo plòot
plot(t, v, ty='l', ylab='v (t)', col='red',  xaxs = "i")
grid()
text(0.9*t.max, 0.9*min(v), 't (s)', cex=1.5)

# ora plottiamo l'energia
pausa()

Ep <- 1/2 * k *x^2
Ec <- 1/2 * m * v^2

plot(t, Ep,  ty='l', ylab='Ep, Ec  (J)', col='brown', lwd=1.5,  xaxs = "i", yaxs = "i")
points(t, Ec,  ty='l', ylab='Ec (J)', col='blue', lwd=1.5)
text(0.9*t.max, 0.1*max(Ep), 't (s)', cex=1.5)
grid()
# plot(t, Ec,  ty='l', ylab='Ec (J)', col='blue', lwd=1.5)
# text(0.9*t.max, 0.1*max(Ec), 't (s)', cex=1.5)
# grid()

plot(t, Ep + Ec,  ty='l', ylab='Etot (J)', col='blue', lwd=1.5,  xaxs = "i", yaxs = "i")
text(0.9*t.max, 0.1*max(Ep+Ec), 't (s)', cex=1.5)
grid() 

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

pausa()
plot(t, Ep + Ec,  ty='l', xlab='t (s)', ylab='Etot (J)', col='blue', lwd=1.5,  log='y')
grid() 

# xaxs = "i", yaxs = "i"