# soluzione numerica del problema della gittata
#
# per eseguire: da console R:
# > source("gittata.R")
#
#                          G. D'Agostini ottobre 09
#-----------------------------------------------------
#
#
g = 9.8 
a = c(0, -g)           # vettore accelerazione

#------------- Main ----------------
# condizioni iniziali
v0    <- 10            # ms
alpha <- 45 * (pi/180) # angolo, trasformato in radianti
dt    <- 0.01          # intervalo di campionamento in s 
nmax  <- 100000        # lunghezza max vettori (att. uso memoria!)

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

x[n]  <- r[1]
z[n]  <- r[2]
tv[n] <- t
# loop principale
while ( r[2] >= 0 && n <= nmax) {
  n <- n+1
  vm <- v + a * dt/2                  # velocita' media in dt
  v  <- v + a * dt                    # velocita' dopo dt
  r  <- r + vm*dt                     # posizione dopo dt
  t  <- t + dt                        # nuovo tempo
  s <- s + dt*sqrt(vm[1]^2 + vm[2]^2) # spazio percorso
  x[n]  <- r[1]
  z[n]  <- r[2]
  tv[n] <- t
}
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", s))

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

cat ("\n >> Guarda il plot e dai enter per continuare\n")
scan()

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

cat ("\n >> Guarda il plot e dai enter per continuare\n")
scan()

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

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