#-------------------------------------------------------------------- 
#
#  Sequenza di urti perfettamente elastici
#  al fine di riprodurre https://www.youtube.com/watch?v=HEfHFsfGXjs
#
#   GdA, maggio 2023 
#--------------------------------------------------------------------


urto.elastico <- function(m1, v1, m2, v2) {
    v1p <- ( 2*m2*v2 + (m1-m2)*v1 ) / (m1+m2)  # v.d. (242)
    v2p <- v1 + v1p -v2                        # v.d. (283) 
    return( list(v1=v1p, v2=v2p))
}

m1    = 1      # kg
m2    = 100  # kg
v1.in = 0      # m/s
v2.in = -10    # m/s

N = 100   # limite di stampa

v1 <- 0   # 
v2 <- -10 # m/s

v1 <- numeric()
v2 <- numeric()
i = 1
v1[i] <- v1.in
v2[i] <- v2.in
cat(sprintf("000  v1, v2:  %.3f, %.3f \n", v1[i], v2[i]))
n.urti <- 0 
while(1) {
    i <- i+1     
    # urto fra i due oggetti
    ris <- urto.elastico(m1, v1[i-1], m2, v2[i-1])
    n.urti <- n.urti + 1 
    v1[i] <- ris$v1    
    v2[i] <- ris$v2
    if(i<=N) cat(sprintf("%.3d  v1, v2:  %.3f, %.3f \n", i-1, v1[i], v2[i]))
    if(v1[i] >= 0) break
    # urto con la parete a sinistra   
    n.urti <- n.urti + 1
    v1[i] <- -v1[i]
    if(v1[i] <= v2[i]) break
    if(i<=N) cat(sprintf("%.3d  v1, v2:  %.3f, %.3f \n", i-1, v1[i], v2[i]))
}
fmt <- sprintf(" Nr di urti: %%d;  nr. di urti * sqrt(m1/m2): %%.%df \n",
               ceiling(log10(n.urti)) )
cat(sprintf(fmt,  n.urti, n.urti*sqrt(m1/m2))) 

titolo = sprintf("m2/m1 = %d/%d;    nr urti _totali_: %d", m2, m1, n.urti)
par(mfrow= c(2,1) )
old.mar = par("mar")
par(mar=c(2,4.2,4.1,0.3))
plot(v1, col='red', cex=1, pch=19, xaxs = "i", yaxs = "i",
     ylim=range(v1)*1.05, xlab='i', ylab='v1 (m/s)', main=titolo)
par(mar=c(4,4.2,0.1,0.3))
plot(v2, col='blue', cex=1, pch=19,  xaxs = "i", yaxs = "i",
     ylim=range(v2)*1.05, xlab='i', ylab='v2 (m/s)')
abline(h=0)
par(mar=old.mar)
par(mfrow= c(1,1) )