/* Trajectoire fusee */
/* (examen physique numerique janvier 2007) */
/* les variables x,y,vx,vy du modele sont stockees dans les tableaux selon
   la correspondance:
     x <-> S[0]
     y <-> S[1]
     vx <-> S[2]
     vy <-> S[3]
     m <-> S[4]
*/

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

// nombre de variables d'etat (=nombre d'equations d'evolution)
#define NEQ 5

double alpha = 60/*deg*/ *(M_PI/180); /* inclinaison fusee par rapport
                                         à la verticale */
double carburant = 0.97; /* fraction initiale carburant */
double poussee = 1.6; /* poussee moteur initiale en g */
double combust = 4; /* taux de combustion (-dm/dt)*/
int vollibre = 0; /* 0: avec moteur, >0: vol libre */

/* conditions initiales de notre probleme */
void init(double *S) {
  S[0]=1;// surface de la terre
  S[1]=0;
  S[2]=0;
  S[3]=0.058;// vitesse tangentielle de l'equateur
  S[4]=1;
}

/* Cette fonction contient toute la physique: les equations
   d'evolution sont utilisees pour calculer le changement d'etat dS a
   partir de l'etat S */
void derive(double *S, double* dS) {
  double r2,phi; //coordonnees cylindriques
  r2 = S[0]*S[0]+S[1]*S[1];
  phi = atan2(S[1],S[0]);
  //printf("%f %f\n",r2,phi);
  dS[0] = S[2];
  dS[1] = S[3];
  if (vollibre == 0) {// encore du carburant ?
    dS[2] = - cos(phi)/r2 + cos(phi+alpha)*poussee/S[4];
    dS[3] = - sin(phi)/r2 + sin(phi+alpha)*poussee/S[4];
    dS[4] = -combust;
  } else {// sinon vol libre
    dS[2] = - cos(phi)/r2;
    dS[3] = - sin(phi)/r2;
    dS[4] = 0;
  }
}

/* Cette fonction integre les equations d'evolution sur un pas de
   temps h, selon un schema de Runge-Kutta d'ordre 4, en partant de
   l'etat contenu dans le tableau 'etat'. Le nouvel etat du systeme
   est finalement stocke dans le meme tableau. */
void rk4(double *etat, double h) {
  double etatTest[NEQ], k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ];
  int i;

  // k1 = f(etat)
  derive(etat,k1);

  // k2 = f(etat+0.5*dt*k1)
  for (i=0; i<NEQ; i++)
    etatTest[i] = etat[i] + 0.5*h*k1[i];
  derive(etatTest,k2);

  // k3 = f(etat+0.5*dt*k2)
  for (i=0; i<NEQ; i++)
    etatTest[i] = etat[i] + 0.5*h*k2[i];
  derive(etatTest,k3);

  // k4 = f(etat+dt*k3)
  for (i=0; i<NEQ; i++)
    etatTest[i] = etat[i] + h*k3[i];
  derive(etatTest,k4);

  // nouvel etat = ancien etat + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
  for (i=0; i<NEQ; i++)
    etat[i] = etat[i] + h/6*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
}

/* Point d'entree du programme. Cette fonction contient
   essentiellement la boucle d'integration et l'affichage de l'etat a
   chaque pas de temps */
int main(void) {
  // tableau contenant l'etat du systeme (voir premieres lignes du
  // fichier pour la description des elements)
  double S[NEQ];
  // intervalle d'integration [t,tmax] et pas d'integration
  double t=0, tmax=100, dt=0.001;
  // compteur de pas, nombre de pas d'integration
  int n,nmax;
  nmax = (tmax-t)/dt;
  // initialisation
  init(S);
  // integration
  for (n=0; n<nmax ; n++) {
    t = n*dt;
    // afficher l'état actuel
    printf("%lf %lf %lf %lf %lf %lf\n",t,S[0],S[1],S[2],S[3],S[4]);
    // integrer d'un pas de temps dt
    rk4(S,dt);
    // verifier si nous avons encore du carburant
    if ((vollibre==0) && (S[4] < 1-carburant)) {
      vollibre = 1;
      dt = 100*dt;
    }
    // verifier si on est retombe sur la terre
    if ((vollibre>0) && (S[0]*S[0]+S[1]*S[1] <= 1)) {
      printf("# crash!");
      break;
    }
  }
  // affichage de l'état final
  t = n*dt;
  printf("%lf %lf %lf %lf %lf %lf\n",t,S[0],S[1],S[2],S[3],S[4]);
  return EXIT_SUCCESS;
}