/* henon-heiles-sez-poin.c : file sorgente di un programma che
                             determina numericamente le sezioni di
                             Poincare' per il modello di Henon e
                             Heiles. */

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

/* Definizione dei prototipi delle funzioni */
void corrector_motion(double x[], double y[], double omega[], double dt);
void da_cartesiane_a_polari(double x, double y, double *rho, double *theta);
void da_polari_a_cartesiane(double rho, double theta, double *x, double *y);
double energia(double x[], double y[], double omega[]);
void harm_osc_motion(double x[], double y[], double omega[], double dt);
void perturb_motion(double x[], double y[], double dt);
void sbab3(double x[], double y[], double omega[],
	   double c_dt[], double d_dt[]);
void sbab3c(double x[], double y[], double omega[],
	    double c_dt[], double d_dt[], double coef_corr_dt);
void sez_poin(double xt1[], double yt1[], double xt2[], double yt2[],
	      double omega[], double c_dt[], double d_dt[],
	      double coef_corr_dt, double err,
	      double x_sez_poin[], double y_sez_poin[]);


/* ### */
int main()
{
  int i, j, numintervalli;
  double x[2], y[2], xiniz[2], yiniz[2], omega[2];
  double xprec[2], yprec[2], xstop[2], ystop[2], xtmp[2], ytmp[2];
  double energia, energiainiz;
  double c[2], d[2], coef_corr, Tf, dt, c_dt[2], d_dt[2], coef_corr_dt;
  char riga[200];
  FILE *ifp, *ofp;

  printf("  Leggo dal file di input i parametri e le condizioni iniziali\n\n");
  /* Inizio delle operazioni di lettura dal file di input */
  ifp = fopen("henon-heiles-sez-poin.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", &energiainiz);
  /* 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 i valori iniziali della seconda coppia di x e y */
  sscanf(riga, "%lf %lf", &xiniz[1], &yiniz[1]);
  /* Pongo il valore iniziale della prima coordinata x[0] uguale a 0 */
  xiniz[0] = 0;
  /* Calcolo il valore del primo momento y[0] a partire dall'energia */
  energia = energiainiz - omega[1] * ( pow(xiniz[1],2) + pow(yiniz[1],2) ) / 2
            + pow(xiniz[1],3) / 3;
  yiniz[0] = sqrt( 2 * energia / omega[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 tempo finale di integrazione */
  sscanf(riga, "%lf", &Tf);
  /* 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 valore del passo di integrazione */
  sscanf(riga, "%lf", &dt);
  /* Fine delle operazioni di lettura dal file di input */
  fclose(ifp);

  /* Pongo i valori degli elementi dei vettori c e d uguali alle
     costanti relative al metodo SBAB3 (vedasi articolo di Laskar e
     Robutel del 2001 su Celestial Mechanics and Dynamical
     Astronomy) */
  c[0] = 0.5 - sqrt(5.) / 10.;  c[1] = sqrt(5.) / 5.;
  d[0] = 1. / 12.;  d[1] = 5. / 12.;
  /* Pongo il valore del coefficiente di correzione uguale
     all'opportuna costante del metodo SBAB3C (vedasi il suddetto
     articolo di Laskar e Robutel del 2001) */
  coef_corr = ( 13. - 5. * sqrt(5) ) / 288.;

  /* Pongo numintervalli uguale al rapporto tra Tf e il passo di
     integrazione */
  numintervalli = Tf / dt;
  /* 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
     SBAB3 */
  for(i=0; i<2; i++) {
    c_dt[i] = c[i] * dt;
    d_dt[i] = d[i] * dt;
  }
  /* Calcolo il valore di coef_corr_dt in modo che sia uguale al
     valore dei passi di integrazione che vengono effettuati lungo
     il flusso del termine di correzione, durante l'applicazione del
     metodo SBAB3C */
  coef_corr_dt = - coef_corr * 0.5 * dt * dt * dt;

  /* 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];
  }    
  /* Pongo i valori delle coordinate canoniche precedenti a una
     sezione di Poincare' uguali ai corrispondenti valori iniziali
     letti dal file di input */
  for(i=0; i<2; i++) {
    xprec[i] = xiniz[i];
    yprec[i] = yiniz[i];
  }    
  /* Inoltre, per far funzionare correttamente il programma, se
     xprec[0] e' proprio uguale a zero allora lo modifico in modo che
     sia molto piccolo e abbia lo stesso segno della derivata
     temporale di x[0] (per quanto riguarda il modello di Henon e
     Heiles, una sezione di Poincare' con tangente orizzontale e'
     possibile solo nel caso della curva limite) */
  if( xprec[0] == 0. )
    xprec[0] += omega[0] * y[0] * dt * 1.e-16;
 
  /* Apro il file di output */
  ofp = fopen("henon-heiles-sez-poin.out" , "w");
  /* Se le coordinate canoniche corrispondenti al tempo iniziale sono
     proprio in corrispondenza a una sezione di Poincare' ... */
  if( xiniz[0] == 0. ) {
    /* ... allora stampo su file i valori della seconda coppia di
       coordinate canoniche */
    fprintf(ofp, "%13.6le  %13.6le\n", xiniz[1], yiniz[1]);
  }

  /* Calcolo i valori delle coordinate canoniche fino al tempo Tf
     applicando numintervalli volte il metodo SBAB3C */
  for(j=0; j<numintervalli; j++) {
    sbab3c(x, y, omega, c_dt, d_dt, coef_corr_dt);
    /* Se le nuove coordinate canoniche sono proprio in corrispondenza
       a una sezione di Poincare' ... */
    if( x[0] == 0. ) {
      /* ... allora stampo su file i valori della seconda coppia di
	 coordinate canoniche ... */
      fprintf(ofp, "%13.6le  %13.6le\n", x[1], y[1]);
   }
    /* altrimenti, se il valore del primo momento y[0] ha cambiato
       segno rispetto al suo valore precedente, ... */
    else if( xprec[0] * x[0] < 0 ) {
      /* ... allora esiste una sezione di Poincare' all'interno
         dell'ultimo intervallo di integrazione (di passo dt) che e'
         stato appena effettuato; quindi, memorizzo i nuovi valori
         delle coordinate */
      for(i=0; i<2; i++) {
	xstop[i] = x[i];  xtmp[i] = x[i];
	ystop[i] = y[i];  ytmp[i] = y[i];
      }
      /* Determino la sezione di Poincare', con una chiamata ad
	 un'opportuna function, che effettuera' la ricerca della
	 sezione di Poincare' con un metodo di bisezione */
      sez_poin(xprec, yprec, xtmp, ytmp, omega, c_dt, d_dt, coef_corr_dt, 1.e-7,
	       x, y);
      /* Stampo su file i valori della seconda coppia di coordinate
	 canoniche (che individuano la sezione di Poincare') */
      fprintf(ofp, "%12.6le  %12.6le\n", x[1], y[1]);
      /* Ripristino i valori delle coordinate canoniche in corrispondenza
	 a quelli dove ho precedentemente arrestato l'integrazione numerica */
      for(i=0; i<2; i++) {
	x[i] = xstop[i];
	y[i] = ystop[i];
      }
      /* Ripristino i valori corretti 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 SBAB3 */
      for(i=0; i<2; i++) {
	c_dt[i] = c[i] * dt;
	d_dt[i] = d[i] * dt;
      }
      /* Ripristino (anche se non ce ne sarebbe bisogno) il valore di
	 coef_corr_dt in modo che sia uguale al valore dei passi di
	 integrazione che vengono effettuati lungo il flusso del
	 termine di correzione, durante l'applicazione del metodo
	 SBAB3C */
      coef_corr_dt = - coef_corr * 0.5 * dt * dt * dt;
    }

    /* Aggiorno i valori delle coordinate canoniche precedenti a una
       sezione di Poincare'  */
    for(i=0; i<2; i++) {
      xprec[i] = x[i];
      yprec[i] = y[i];
    }
    /* Inoltre, per far funzionare correttamente il programma, se
       xprec[0] e' proprio uguale a zero allora lo modifico in modo
       che sia molto piccolo e abbia lo stesso segno della derivata
       temporale di x[0] (per quanto riguarda il modello di Henon e
       Heiles, una sezione di Poincare' con tangente orizzontale e'
       possibile solo nel caso della curva limite) */
    if( xprec[0] == 0. )
      xprec[0] += omega[0] * y[0] * dt * 1.e-16;

  }


  /* Chiudo il file di output */
  fclose(ofp);
  printf("  La tabella delle sezioni di Poincare' e' state memorizzata\n");
  printf("  nel file henon-heiles-sez-poin.out, che puo' essere messo in\n");
  printf("  grafico, utilizzando gnuplot \n\n");

  /* A mo' di test finale, calcolo il valore dell'energia corrispondente ai
     valori assunti dalle coordinate canoniche al tempo Tf */
  energia = 0;
  for(i=0; i<2; i++)
    energia += omega[i] * ( x[i] * x[i] +  y[i] * y[i] ) / 2.;
  energia += x[0] * x[0] * x[1] - x[1] * x[1] * x[1] / 3.;
  /* Stampo la differenza tra l'energia calcolata al tempo Tf e quella
     iniziale */
  printf("  A mo' di test finale sulla correttezza dell'integrazione\n");
  printf("  numerica, stampo la differenza tra energia finale ed iniziale :\n");
  printf("    %8.2le\n", fabs( energia - energiainiz ) );

  /* Stampo su video due righe "vuote" per distinguere meglio la fine
     del programma dai messaggi di "prompt" successivi, all'interno
     della finestra di terminale */
  printf("\n\n");

}


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

/* ### */
void da_cartesiane_a_polari(x, y, rho, theta)
     double x, y, *rho, *theta;
{
  *rho = sqrt(x*x + y*y);
  *theta = atan2(y, x);
}

/* ### */
void da_polari_a_cartesiane(rho, theta, x, y)
     double rho, theta, *x, *y;
{
  *x = rho * cos(theta);
  *y = rho * sin(theta);
}

/* ### */
double energia(x,y,omega)
     double x[2], y[2], omega[2];
{
  int i;
  double ris = 0;
  for(i=0; i<2; i++) 
    ris += omega[i] * ( x[i] * x[i] + y[i] * y[i] ) / 2;
  ris += x[1] * ( x[0] * x[0] - x[1] * x[1] / 3 );
  return ris;
}

/* ### */
void harm_osc_motion(x,y,omega,dt)
     double x[2], y[2], omega[2], dt;
{
  int i;
  double rho, theta;
  for(i=0; i<2; i++) {
    da_cartesiane_a_polari(x[i], y[i], &rho, &theta);
    theta -= omega[i] * dt;
    da_polari_a_cartesiane(rho, theta, x+i, y+i);
  }
}

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

/* ### */
void sbab3(x,y,omega,c_dt,d_dt)
     double x[2], y[2], omega[2], c_dt[2], d_dt[2];
{
  perturb_motion(x,y,d_dt[0]);
  harm_osc_motion(x,y,omega,c_dt[0]);
  perturb_motion(x,y,d_dt[1]);
  harm_osc_motion(x,y,omega,c_dt[1]);
  perturb_motion(x,y,d_dt[1]);
  harm_osc_motion(x,y,omega,c_dt[0]);
  perturb_motion(x,y,d_dt[0]);
}

/* ### */
void sbab3c(x,y,omega,c_dt,d_dt, coef_corr_dt)
     double x[2], y[2], omega[2], c_dt[2], d_dt[2], coef_corr_dt;
{
  corrector_motion(x, y, omega, coef_corr_dt);
  sbab3(x,y,omega,c_dt,d_dt);
  corrector_motion(x, y, omega, coef_corr_dt);
}

/* ### */
void sez_poin(xt1, yt1, xt2, yt2, omega, c_dt, d_dt, coef_corr_dt, err,
	      x_sez_poin, y_sez_poin)
     double xt1[2], yt1[2], xt2[2], yt2[2];
     double omega[2], c_dt[2], d_dt[2], coef_corr_dt, err;
     double x_sez_poin[2], y_sez_poin[2];
{
  int i;

  /* Ricerca della sezione di Poincare' con il metodo di bisezione */
  do {
    /* Dimezzo l'intervallo di tempo corrispondente al passo di
       integrazione, quindi ... */ 
    for(i=0; i<2; i++) {
      c_dt[i] /= 2;
      d_dt[i] /= 2;
    }
    coef_corr_dt /= 8;
    /* Pongo momentaneamente la nuova approssimazione della sezione
       di Poincare' uguale alla soluzione corrispondente al valore
       dell'estremo sinistro dell'intervallo di tempo */
    for(i=0; i<2; i++) {
      x_sez_poin[i] = xt1[i];
      y_sez_poin[i] = yt1[i];
    }
    /* Calcolo la nuova approssimazione della sezione di Poincare' in
       corrispondenza al centro dell'intervallo di tempo di
       integrazione, applicando il metodo SBAB3C */
    sbab3c(x_sez_poin, y_sez_poin, omega, c_dt, d_dt, coef_corr_dt);
    /* Se abbiamo "preso in pieno" la sezione di Poincare' esco dal
       ciclo */
    if( x_sez_poin[0] == 0. )
      break;
    /* Se la sezione di Polincare' e' compreso tra l'estremo sinistro
       dell'intervallo di integrazione e la nuova approssimazione ... */
    if( xt1[0] * x_sez_poin[0] < 0 ) {
      /* ... allora sposto l'estremo destro in corrispondenza alla
	 nuova approssimazione ... */
      for(i=0; i<2; i++) {
	xt2[i] = x_sez_poin[i];
	yt2[i] = y_sez_poin[i];
      }
    }
    /* altrimenti ... */
    else {
      /* ... sposto l'estremo sinistro in corrispondenza alla nuova
	 approssimazione ... */
      for(i=0; i<2; i++) {
	xt1[i] = x_sez_poin[i];
	yt1[i] = y_sez_poin[i];
      }
    }
  } while( fabs( x_sez_poin[0] ) < err );

}