/* sez_Poi_LF_HH.c : file sorgente di un programma che consente di
                     localizzare le sezioni di Poincare' (per diverse
                     orbite, ma corrispondenti allo stesso livello di
                     energia) per il modello di Henon e
                     Heiles. L'integrazione numerica delle equazioni
                     del moto viene effettuata con il metodo
                     simplettico di integrazione numerica di tipo
                     leap-frog (o SABA1) */

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

/* Definizione dei prototipi delle funzioni */
void interpola_linearmente(double xprec[2], double yprec[2],
			   double x[2], double y[2], double xy_sezPoi[2]);
void kin_flow(double x[2], double y[2], double omega[2], double dt);
void pot_flow(double x[2], double y[2], double omega[2], double dt);
void saba1(double x[], double y[], double omega[],
	   double c_dt[], double d_dt[]);


/* ### */
int main()
{
  int i, numfin_sezPoi, conta_sezPoi;
  double x[2], y[2], xiniz[2], yiniz[2], xprec[2], yprec[2], omega[2];
  double xy_sezPoi[2], en_0, delta_en;
  double c[2], d[2], Tf, dt, dt_graph, c_dt[2], d_dt[2];
  char riga[200];
  FILE *ifp, *ofp;
  
  printf("  Lettura dal file di input: parametri e condizioni iniziali\n\n");
  /* Inizio delle operazioni di lettura dal file di input */
  ifp = fopen("sez_Poi_LF_HH.inp" , "r");
  /* Leggo "a vuoto" dal file una riga di commenti */
  fgets(riga, 200, ifp);
  /* Leggo dal file una riga di dati */
  fgets(riga, 200, ifp);
  /* Leggo dalla riga le frequenze omega */
  sscanf(riga, "%lf %lf", &omega[0], &omega[1]);
  /* Leggo "a vuoto" dal file una riga di commenti */
  fgets(riga, 200, ifp);
  /* Leggo dal file una riga di dati */
  fgets(riga, 200, ifp);
  /* Leggo dalla riga di dati il valore dell'energia */
  sscanf(riga, "%lf", &en_0);
  /* Leggo "a vuoto" dal file una riga di commenti */
  fgets(riga, 200, ifp);
  /* Leggo dal file una riga di dati */
  fgets(riga, 200, ifp);
  /* Leggo dalla riga il numero delle sezioni di Poincare' da
     determinare per ciascuna orbita */
  sscanf(riga, "%d", &numfin_sezPoi);
  /* Leggo "a vuoto" dal file una riga di commenti */
  fgets(riga, 200, ifp);
  /* Leggo dal file una riga di dati */
  fgets(riga, 200, ifp);
  /* Leggo dalla riga il passo di integrazione */
  sscanf(riga, "%lf", &dt);
  /* Leggo "a vuoto" dal file una riga di commenti */
  fgets(riga, 200, ifp);
  
  /* Pongo i valori degli elementi dei vettori c e d uguali alle
     costanti relative al metodo leap-frog, che e' anche il SABA1,
     utilizzando la notazione in accordo con l'articolo di Laskar e
     Robutel del 2001 su Celestial Mechanics and Dynamical
     Astronomy) */
  c[0] = 0.5;
  d[0] = 1.;
  
  /* Calcolo i valori degli elementi dei vettori c_dt e d_dt in modo
     che corrispondano ai valori dei passi di integrazione che
     vengono effettuati (lungo il flusso delle due parti in cui
     viene divisa la Hamiltoniana) durante l'applicazione del metodo
     leap-frog */
  for(i=0; i<1; i++) {
    c_dt[i] = c[i] * dt;
    d_dt[i] = d[i] * dt;
  }
  
  /* Apertura del file di output */
  ofp = fopen("sez_Poi_LF_HH.out" , "w");

  printf("  Integrazione numerica a partire dalle seguenti cond. iniz.:\n");
  /* Eseguo il ciclo di istruzioni racchiuse tra il prossimo
     "while(...){" e la corrispondente "}", ogni qual volta che nel
     file esterno c'e' una riga non vuota di dati */
  while(fgets(riga, 200, ifp) != NULL) { 
    /* Leggo dalla riga di dati i valori iniziali della seconda coppia
       di x e y */
    printf("%s", riga);
    sscanf(riga, "%lf %lf", &xiniz[1], &yiniz[1]);
    /* Pongo il valore iniziale della prima coordinata x uguale a 0 */
    xiniz[0] = 0;
    /* Calcolo il primo momento y a partire dall'energia */
    delta_en = en_0 - omega[1] * ( pow(xiniz[1],2) + pow(yiniz[1],2) ) / 2
               + pow(xiniz[1],3) / 3;
    yiniz[0] = sqrt( 2 * delta_en / omega[0] );
    /* Stampo su file esterno i valori iniziali della seconda coppia
       di x e y poiche' sono tali da essere in corrispondenza alla
       sezione di Poincare' */
    fprintf(ofp, "     %le      %le\n", xiniz[1], yiniz[1] );
    /* Pongo i valori delle coordinate canoniche uguali ai
       corrispondenti valori iniziali letti dal file di input */
    for(i=0; i<2; i++) {
      x[i] = xiniz[i];
      y[i] = yiniz[i];
    }
    /* Aggiorno i valori precedenti delle coordinate canoniche in modo
       che siano uguali a quelli attuali */
    for(i=0; i<2; i++) {
      xprec[i] = x[i];
      yprec[i] = y[i];
    }
    /* Azzero il contatore delle sezioni di Poincare' */
    conta_sezPoi = 0;
    
    /* Eseguo il ciclo di istruzioni racchiuse tra il prossimo
       "while(...){" e la corrispondente "}", ogni qual volta che il
       valore del contatore delle sezioni di Poincare' e' inferiore al
       numero finale di sezioni di Poincare' per ciascuna orbita */
    while(conta_sezPoi < numfin_sezPoi) {
      /* Calcolo i nuovi valori delle coordinate canoniche applicando il
	 metodo leap-frog */
      saba1(x, y, omega, c_dt, d_dt);
      /* Effettuo un test per comprendere se il valore della x della
	 prima coppia di coordinate canoniche ha cambiato segno,
	 rispetto a prima dell'ultima applicazione del metodo
	 leap-frog... */
      if( xprec[0] * x[0] < 0 ) {
	/* ... allora esiste una sezione di Poincare' all'interno
	   dell'ultimo intervallo di integrazione (di passo dt),
	   quindi la determino, con una chiamata ad un'opportuna
	   function, che effettuera' la ricerca della sezione di
	   Poincare' per interpolazione */
	interpola_linearmente(xprec, yprec, x, y, xy_sezPoi);
	/* Stampo su file i valori della seconda coppia di coordinate
	   canoniche, in corrispondenza alla sezione di Poincare') */
	fprintf(ofp, "     %le      %le\n", xy_sezPoi[0], xy_sezPoi[1]);
	/* Incremento di uno il valore del contatore delle sezioni di
	   Poincare' */
	conta_sezPoi++;
      }
      /* Aggiorno i valori precedenti delle coordinate canoniche in
	 modo che siano uguali a quelli attuali */
      for(i=0; i<2; i++) {
	xprec[i] = x[i];
	yprec[i] = y[i];
      }
    }
  }
  
  /* Fine delle operazioni di scrittura - chiusura dei file di input e
     di output */
  fclose(ifp);
  fclose(ofp);
}


/* ### */
void interpola_linearmente(double xprec[2], double yprec[2],
			   double x[2], double y[2], double xy_sezPoi[2])
{
  double delta_x0, coef_prop;
  delta_x0 = x[0] - xprec[0];
  coef_prop = - xprec[0] / delta_x0;
  xy_sezPoi[0] = xprec[1] + coef_prop * (x[1] - xprec[1]);
  xy_sezPoi[1] = yprec[1] + coef_prop * (y[1] - yprec[1]);
}


/* ### */
void kin_flow(double x[2], double y[2], double omega[2], double dt)
{
  x[0] += omega[0] * y[0] * dt;
  x[1] += omega[1] * y[1] * dt;
}


/* ### */
void pot_flow(double x[2], double y[2], double omega[2], double dt)
{
  y[0] -= ( omega[0] * x[0] + 2 * x[0] * x[1] ) * dt;
  y[1] -= ( omega[1] * x[1] + x[0] * x[0] - x[1] * x[1] ) * dt;
}


/* ### */
void saba1(double x[2], double y[2], double omega[2],
	   double c_dt[2], double d_dt[2])
{
  kin_flow(x, y, omega, c_dt[0]);
  pot_flow(x, y, omega, d_dt[0]);
  kin_flow(x, y, omega, c_dt[0]);
}