# soluzione numerica del problema della gittata
#
# per eseguire: da console R:
# > source("gittata1.R")
#
# rispetto a gittata.R:
#   - accelerazione data da funzione che dipende da v
#   - possibilita' di fissare il punto di partenza ( 'r0' )
#   - traiettoria a punti equidistanziati in tempo
#     ('visione stroboscopica')
#   - plot di posizione, velocita' e accelerazione
#   - plot di spazio percorso in funzione del tempo
#   - griglia per meglio rileggere i punti sul plot
#   - piu' plot sulla stessa schermata
#   - possibilita' di inserire valori da tastiera
#
#                          G. D'Agostini ottobre 09 (corso PED)
#--------------------------------------------------------------

# function accelerazione
a <- function(v)  c(0,-g) - beta.su.m * v

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}


#------------- Main ----------------
# condizioni iniziali e parametri
r0 <- c(0,0)     # posizione iniziale
beta.su.m <- 5   # beta/m
g <- 9.8
v0    <- 10            # ms
alpha <- 45 * (pi/180) # angolo, trasformato in radianti
dt    <- 0.005         # intervalo di campionamento in s 
nmax  <- 100000        # lunghezza max vettori (att. uso memoria!)

## in alternativa si possono dare dei PARAMETRI DA TASTIERA:
#  (per togliere l'opzione, commentare le seguenti linee)
cat ("\n >> dai velocita' (m/s)\n")
v0 <- scan(nmax=1)
cat ("\n >> dai angolo (gradi)\n")
alpha <- scan(nmax=1) * (pi/180)

# moto in x-z
v <- c(v0*cos(alpha), v0*sin(alpha))  # velocita' iniziale
r <- r0                               # raggio vettore
t <- 0
n <- 1
s <- numeric()                           # spazio percorso (per plot)
tv <- numeric()                          # tempo (per plot)
x  <- numeric()                          # coordinata x (per plot)
z  <- numeric()                          # coordinata z (per plot)
vx <- numeric()                          # velocita' x (per plot)
vz <- numeric()                          # velocita' z (per plot)
ax <- numeric()                          # accelerazione x (per plot)
az <- numeric()                          # accelerazione z (per plot)

s[n]  = 0
x[n]  <- r[1]
z[n]  <- r[2]
vx[n] <- v[1]
vz[n] <- v[2]
ax[n] <- a(v)[1]
az[n] <- a(v)[2]
tv[n] <- t

# loop principale
while ( r[2] >= 0 && n <= nmax) {
  n <- n+1

  # tutta la fisica e' nelle seguenti quattro istruzioni:
  vm <- v + a(v) * dt/2                       # velocita' media in dt
  v  <- v + a(v) * dt                         # velocita' dopo dt  
  r  <- r + vm*dt                             # posizione dopo dt
  s[n] <- s[n-1] + dt*sqrt(vm[1]^2 + vm[2]^2) # spazio percorso
  
  t  <- t + dt                                # nuovo tempo

  # memorizziamo le variabili cinematiche da analizzare/plottare:
  tv[n] <- t
  x[n]  <- r[1]
  z[n]  <- r[2]
  vx[n]  <- v[1]
  vz[n]  <- v[2]
  ax[n]  <- a(v)[1]
  az[n]  <- a(v)[2]
}

cat( sprintf("numero steps %d\n",n))
cat( sprintf("max x (gittata): %f m,   max z (altezza max) %f m\n",
             max(x),max(z)))
cat( sprintf("spazio percorso: %f m\n", max(s)))

plot(x,z,ty='l', main="traiettoria"); grid()
pausa()

# la vista a punti da' una immagine stroboscopica
# della traiettoria, dalla quale si ha una percezione
# della velocita' di percorrenza
# (il parametro cex influenza la grandezza dei puntini,
#  qui fatti piccoli per veder meglio la distanza fra di loro
# in funzione del tempo)
plot(x,z, main="traiettoria - plot stroboscopico", cex=0.2); grid()
pausa()

plot(tv,x,ty='l', main="eq. oraria x(t)", xlab="t"); grid()
pausa()

plot(tv,z,ty='l', main="eq. oraria z(t)", xlab="t"); grid()
pausa()

plot(tv, s, ty='l',  main="eq. oraria s(t)", xlab="t"); grid()
pausa()


# tre plot insieme
layout(matrix(1:3,3,1)) 

grid()
plot(tv,x,ty='l', main="eq. oraria x(t)", xlab="t"); grid()
plot(tv,vx,ty='l', main="eq. oraria vx(t)", xlab="t"); grid()
plot(tv,ax,ty='l', main="eq. oraria ax(t)", xlab="t"); grid()
pausa()

plot(tv,z,ty='l', main="eq. oraria z(t)", xlab="t"); grid()
plot(tv,vz,ty='l', main="eq. oraria vz(t)", xlab="t"); grid()
plot(tv,az,ty='l', main="eq. oraria az(t)", xlab="t"); grid()
pausa()

layout(matrix(1:2,2,1)) 
plot(tv, s, ty='l',  main="eq. oraria s(t)", xlab="t"); grid()
plot(tv, sqrt(vx^2+vz^2), ty='l',  main="eq. oraria v(t)", xlab="t"); grid()

#layout(matrix(1:3,3,1))
#plot(tv, s, ty='l',  main="eq. oraria s(t)", xlab="t"); grid()
#plot(tv, sqrt(vx^2+vz^2), ty='l',  main="eq. oraria v(t)", xlab="t"); grid()
#plot(tv, sqrt(ax^2+az^2), ty='l',  main="eq. oraria a(t)", xlab="t"); grid()

layout(1) # rimettiamo a posto il layout 

cat ("\n >> --- Fine ---\n")
cat ("\n >>>> cambia i parametri e riesegui lo script\n")