#-------------------------------------------------
#  Oscillazioni smorzate: 
#  Riferimenti: circuiti_253-267.pdf , par 10.9
#
#   GdA maggio 2025
#--------------------------------------------------

T0  <-  1        # s -> periodo senza smorzamento)
om0 <-  2*pi/T0  #   -> pulsazione
z0  <-  1        # ampiezz<a arbitraria
t.max <- 10*T0   # max scala temporale dei plot

#  gamma in unità di 2*om0
#  
gamma <- ( 2* om0 ) / 10.00001
# ATT! se il denomiatore vale 1  DIVERGE -> (10.90)-10.91
#      => basta mettere 1.001 (vedi par.10.9.6)

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 

print(alpha1)
print(alpha2)

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

print(k1)
print(k2)

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)', main=titolo)
grid()
text(0.9*t.max, 0.9*min(x), 't (s)', cex=1.5)

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

# plot di x.dot  (-> 'velocità' o 'corrente')
par(mar=c(2,4.2,0.5,0.5))    # parametri leggermente diversi per il secondo plòot
plot(t, x.dot, ty='l', ylab='x.dot (t)', col='red')
grid()
text(0.9*t.max, 0.9*min(x.dot), 't (s)', cex=1.5)

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