/* scrivisegnale_HH.c : file sorgente di un programma che consente di
                        scrivere su file esterni i segnali associati
                        alle leggi del moto delle coppie di variabili
                        canoniche. L'output e' semplicissimo (ciascun
                        file e' costituito da una potenza intera di 2
                        piu' una righe e 2 colonne di numeri in double
                        precision) ed e' organizzato in modo che sia
                        possibile studiarlo con i programmi che
                        implementano l'analisi in frequenza.
                        L'integrazione numerica viene effettuata con
                        il metodo simplettico di tipo leap-frog (o
                        SABA1). Al termine del procedimento di
                        integrazione numerica, viene comunque
                        riportato il risultato del test sulla
                        conservazione dell'energia per il modello di
                        Henon-Heiles. */

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

/* Definizione dei prototipi delle funzioni */
double energia(double x[2], double y[2], double omega[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 j;
  long int i, numintervalli, i_graph;
  double x[2], y[2], xiniz[2], yiniz[2], omega[2];
  double en_t, en_0, delta_en, errmax_en;
  double c[2], d[2], Tf, dt, dt_graph, c_dt[2], d_dt[2];
  char riga[200], nomefile[20];
  FILE *ifp, *ofp[2];
  
  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("scrivisegnale_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 i valori iniziali della prima coppia di x e y */
  sscanf(riga, "%lf %lf", &xiniz[0], &yiniz[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 di dati i valori iniziali della seconda coppia di x e y */
  sscanf(riga, "%lf %lf", &xiniz[1], &yiniz[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 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 passo di integrazione */
  sscanf(riga, "%lf", &dt);
  /* 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 l'intervallo di tempo di scansione del segnale */
  sscanf(riga, "%lf", &dt_graph);
  /* Fine delle operazioni di lettura - chiusura del file di input */
  fclose(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.;
  
  /* Pongo i valori delle coordinate canoniche uguali ai
     corrispondenti valori iniziali letti dal file di input */
  for(j=0; j<2; j++) {
    x[j] = xiniz[j];
    y[j] = yiniz[j];
  }
  
  /* Pongo numintervalli uguale al rapporto tra il tempo finale e il
     passo di integrazione */
  numintervalli = Tf / dt;
  /* Pongo i_graph uguale al numero di passi che devono intercorrere
     tra il tracciamento di un punto per il segnale e il successivo */
  i_graph = dt_graph / 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
     leap-frog */
  for(j=0; j<1; j++) {
    c_dt[j] = c[j] * dt;
    d_dt[j] = d[j] * dt;
  }
  
  /* Apertura dei file di output */
  for(j=0; j<2; j++) {
    sprintf(nomefile, "segnale_HH_%d.out", j+1);
    ofp[j] = fopen(nomefile , "w");
  }
  /* Stampo sui file esterni i valori delle condizioni iniziali */
  for(j=0; j<2; j++)  
    fprintf(ofp[j], "     %24.16le      %24.16le\n", y[j], x[j]);
  /* Calcolo il valore dell'energia corrispondente ai valori assunti
     dalle coordinate canoniche al tempo iniziale, cioe' 0  */
  en_0 = energia(x, y, omega);
  /* Coerentemente, pongo inizialmente uguale a zero il valore
     dell'errore massimo sull'energia */
  errmax_en = 0.;
  
  /* Inizio l'integrazione numerica vera e propria delle equazioni del
     moto */
  printf("  Effettuo l'integrazione numerica delle equazioni del moto\n\n");
  for(i=1; i<=numintervalli; i++) {
    /* Calcolo i valori delle coordinate canoniche al tempo i * dt
       applicando numintervalli volte il metodo leap-frog */
    saba1(x, y, omega, c_dt, d_dt);
    if(i % i_graph == 0) { 
      /* Stampo il valori assunti all'istante di tempo (i * dt)
	 dalle coordinate canoniche */
      for(j=0; j<2; j++)  
	fprintf(ofp[j], "     %24.16le      %24.16le\n", y[j], x[j]);
      /* Calcolo il valore dell'energia corrispondente ai valori
	 assunti dalle coordinate canoniche al tempo i * dt e lo
	 scarto rispetto all'energia iniziale */
      en_t = energia(x, y, omega);
      delta_en = en_t - en_0;
      /* Eventualmente, aggiorno il valore dell'errore massimo
	 sull'energia */
      if(errmax_en < fabs(delta_en))
	errmax_en = fabs(delta_en);
    }
  }
  
  /* Fine delle operazioni di scrittura - chiusura dei file di output */
  for(j=0; j<2; j++)
    fclose(ofp[j]);
  /* Stampa finale a proposito della conservazione dell'energia */
  printf(" En. iniz.    max_t | E(t) - E(0) |   max_t | (E(t)-E(0))/E(0) |\n"); 
  printf("%le       %le              %le\n\n",
	 en_0, errmax_en, fabs(errmax_en/en_0) );
}


/* ### */
double energia(double x[2], double y[2], double 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 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]);
}