/* Soluzione numerica del problema del       
   punto materiale lanciato con un certo angolo   

   Variante (rispetto a lancio_2D_file.c):
     -> resistenza dell'aria modellizzata realisticamente
        dipendente dal quadrato della velocità 
        (e non semplicemente dalla velocità)
           https://en.wikipedia.org/wiki/Drag_(physics)
         
      -> |F_A| = eta * v^2 

   Per la logica vedere la soluzione numerica della molla:
     -> molla_xva_smorz.R
        (http://www.roma1.infn.it/~dagos/LabC/molla_1.html)

   Legge da riga di comando
      - i parametri fisici (m dell'oggetto e eta di resistenza) 
      - i parametri del lancio (v0, theta)
      - il parametro di discretizzazione del tempo (dt)     

        ./lancio_2D_aria_file m eta v0 theta ylim dt

   Es:  
        ./lancio_2D  1.0 0.1 10 45 -1. 0.05

         (Att: ymin deve essere <= 0, altrimenti viene posto a 0)

*/
 

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

// #define TO_FILE 1
#define HEADER  "t x y vx vy"     // mette i nomi all'inizio delle variabili

int main(int argc, char *argv[]) {
  int step=0;
  float v0, theta, vx, vy, vx0, vy0, vxm0, vym0, ax, ay, amx, amy, dt;
  float ymin; 
  float Fx, Fy, m, v2, Fmx, Fmy, vm2;
  float t=0, x = 0., y = 0.;
  const float grad2rad = acos(-1) / 180.;
  const float g = 9.8;     // m/s^2
  float eta;               // N / (m/s)^2     -> da riga di comando
  int   toFile=0;          // 'flag' per scrittura su file
  FILE * fp;               // puntatore al file
  char nomeFile[50];       // vettore di caratteri per il nome del file
  char *p;                 // puntatore a carattere (-> argv[6])
  

  if (argc < 7) {
    puts("Devi dare anche: m eta v0 theta ymin dt");
    return 0;
  } 

  m     = atof(argv[1]);
  eta   = atof(argv[2]);
  v0    = atof(argv[3]);
  theta = atof(argv[4]);
  ymin  = atof(argv[5]);
  if(ymin > 0) ymin = 0. ; 
  dt    = atof(argv[6]);
  printf("Parametri fisici: m = %.2f kg,  eta = %.3f N / (m/s)^2\n", m, eta);
  printf("Parametri lancio: v0=%.2f m/s, theta=%.0f gradi\n", v0, theta); 
  printf("Discretizzazione: dt = %.3f s\n", dt);
  
  vx = v0 * cos(theta * grad2rad);  
  vy = v0 * sin(theta * grad2rad);
  printf("\nVelocità iniziale: (%.2f, %.2f) m/s\n", vx, vy);

  if (argc > 7) {
    p = argv[7];
    if (*p == 'f' || *p == 'F') {
      toFile=1;
    }
  } 
  
  // apertura file
  if (toFile) {
    sprintf(nomeFile, "lancio_%s_%s_%s_%s_%s_%s.txt",
	    argv[1], argv[2], argv[3], argv[4], argv[5], argv[6]);
    printf("I dati andranno sul file %s\n", nomeFile);
    fp = fopen(nomeFile, "w");
    #ifdef HEADER
    printf("Header definito: %s\n", HEADER);
    fprintf(fp, "%s\n", HEADER);
    #endif
    fprintf(fp, "%.5f %.5f %.5f %.5f %.5f\n", t, x, y, vx, vy);
  } else {
    printf("%.5f %.5f %.5f %.5f %.5f\n", t, x, y, vx, vy);
  }
  
  while (y >= ymin) {
    // salviamo la velocità in t (all'inizio dell'intervallino dt)
    vx0 = vx;
    vy0 = vy;
    
    // quadrato della velocità e forza calcolata in t
    v2 = vx*vx + vy*vy;
    Fx = 0    - eta * v2 * vx/sqrt(v2);  // vx/sqrt(v2) : versore x della velocità
    Fy = -m*g - eta * v2 * vy/sqrt(v2);  // vy/sqrt(v2) : versore y della velocità

    // accelerazione calcolata in t
    ax = Fx / m;
    ay = Fy / m;

    // velocità media (di prima approssimazione) nell'intervallino dt
    // (calcolata dall'accelerazione in t)
    vxm0 = vx + ax * dt/2;  // velocita' media: (v + (v+a*dt))/2
    vym0 = vy + ay * dt/2;  

    // riaggiorniamo la velocità quadra con quella media nell'intervallino
    // e ricalcoliamo forze e accelerazioni
    vm2  = vxm0*vxm0 + vym0*vym0;
    Fmx  = 0    - eta * vm2 * vxm0/sqrt(vm2); 
    Fmy  = -m*g - eta * vm2 * vym0/sqrt(vm2); 
    amx  = Fmx / m;                             // ax in t + dt/2
    amy  = Fmy / m;                             // ay in t + dt/2 
    
    // velocità dopo l'intervallo dt 
    vx += amx * dt;
    vy += amy * dt;

    // posizione dopo l'intervallino dt
    x  += (vx0+vx)/2. * dt;            /* coordinata x dopo l'intervallino dt,
                                          durante il quale il corpo si e' mosso 
                                          con v media (vx0+vx)/2                */
				       				 
    y  += (vy0+vy)/2 * dt;             // idem per la coordinata y

    t  += dt;                          // tempo dopo intervallo dt

    if (toFile) {
      fprintf(fp, "%.5f %.5f %.5f %.5f %.5f\n", t, x, y, vx, vy);
    } else {
      printf("%.5f %.5f %.5f %.5f %.5f\n", t, x, y, vx, vy);
    }
  }

  if (toFile) {
    fclose(fp);
  } 

  return 0;
}