/* introduciamo le strutture, ripassiamo l'I/O, le funzioni */
/* impariamo ad evitare gli effetti collaterali nell'uso di funzioni */

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

#define SUCCESSO 0

/* Notate che il dannato ; alla fine della prossima riga e' cruciale! */

struct configurazione {double x; double v;};


/* I PROTOTIPI */

/* m_ nelle righe sotto ci ricorda che si tratta di variabili mute, */
/* li' solo per memoria (potremmo non metterle nel prototipo, */
/* ma indicare solo il tipo) */

struct configurazione *oa_inizializza(double m_x0, double m_v0);
struct configurazione *eulero(double m_dt, double m_omega2, 
			      struct configurazione *m_p_x_e_v);
struct configurazione *eulero_cromer(double m_dt, double m_omega2, 
			      struct configurazione *m_p_x_e_v);

/* notate che energia() e' generica, */
/* mentre oa_forza() e' proprio per l'oscillatore armonico (da cui il nome) */

double energia(double m_omega2, struct configurazione *m_p_x_e_v);
double oa_forza(double m_omega2, double m_x);
double oa_energia_potenziale(double m_omega2, double m_x);
double my_read_double(char *m_messaggio_print);
long my_read_long(char *m_messaggio_print);

/*****************************************************************************/
int main(void)
{
  double  x0, v0;
  double dt, omega2, tempo_totale, energia0, energia1;
  long n_passi, algoritmo;

  long i;

  struct configurazione x_e_v;

  /* I # sono per far piacere a gnuplot, per esempio, o ad aiutarci */
  /* ad usare un filtro per eliminare le righe che non ci */
  /* serviranno nella grafica */

  printf("# Integrazione Oscillatore Armonico v1.1.1 5 ottobre 2002\n");

  /* queste my_read sono li' per evitarci un po' di codice... */

  tempo_totale = 
    my_read_double(
       "Inserire il tempo totale di integrazione (variabile reale double)");
  algoritmo = 
    my_read_long("Inserire 0 per Eulero, 1 per Eulero-Cromer");
  dt = my_read_double("Inserire dt (variabile reale double)");
  x0 = my_read_double("Inserire x0 (variabile reale double)");
  v0 = my_read_double("Inserire v0 (variabile reale double)");
  omega2 = my_read_double("Inserire omega2 (variabile reale double)");

  /* Il cast lo facciamo noi esplicitamente, per tenere */
  /* le cose sotto controllo!! */

  n_passi = (long)(tempo_totale/dt);

  printf("# dt = %g tempo_totale = %g omega2 = %g numero passi = %d\n",
	 dt, tempo_totale, omega2, n_passi);

  /* una funzione senza effetti collaterali. quel che cambia cambia con = */

  x_e_v = *oa_inizializza(x0, v0);

  printf("# Abbiamo scelto condizioni iniziali x0 = %g v0 = %g\n", 
	 x_e_v.x, x_e_v.v);

  energia0 = energia(omega2, &x_e_v);
  printf("# L'energia al tempo t = 0 vale %g\n", energia0);

  if(algoritmo==0){
    printf("# Usiamo l'algoritmo di Eulero\n");
    for(i=0; i<n_passi;i++){

      /* un'altra funzione senza effetti collaterali... */

      x_e_v = *eulero(dt, omega2, &x_e_v);
      energia1 = energia(omega2, &x_e_v);
      printf("%g %g %g %g\n", 
	     (double)i*dt, x_e_v.x, x_e_v.v,energia1-energia0);
    }
  }else if(algoritmo==1){
    printf("# Usiamo l'algoritmo di Eulero-Cromer\n");
    for(i=0; i<n_passi;i++){

      /* un'altra funzione senza effetti collaterali... */

      x_e_v = *eulero_cromer(dt, omega2, &x_e_v);
      energia1 = energia(omega2, &x_e_v);
      printf("%g %g %g %g\n", 
	     (double)i*dt, x_e_v.x, x_e_v.v,energia1-energia0);
    }
  }else{
    printf("# TUONI E FULMINI:");
    printf("L'algoritmo n. %ld non e' ancora stato implementato. Abort.\n",
	   algoritmo);
    exit(-9);
  }

  energia1 = energia(omega2, &x_e_v);
  printf("# L'energia al tempo t = %d vale %g\n", n_passi, energia1);

  return SUCCESSO;
}

/*****************************************************************************/
struct configurazione  *oa_inizializza(double x0, double v0)
{
  static struct configurazione p_x_e_v;
  p_x_e_v.x = x0;
  p_x_e_v.v = v0;
  return &p_x_e_v;
}

/*****************************************************************************/
struct configurazione *eulero(double dt, double omega2, 
			      struct configurazione *p_x_e_v_old)
{
  static struct configurazione x_e_v_new;

  /* che succede se omettete le prime parentesi nella riga sotto ?*/

  x_e_v_new.v = (*p_x_e_v_old).v + oa_forza(omega2, (*p_x_e_v_old).x) * dt;

  /* perche' qui usiamo -> e non . ? */

  x_e_v_new.x = p_x_e_v_old->x + p_x_e_v_old->v * dt;
  return &x_e_v_new;
}

/*****************************************************************************/
/* Qual'e' la differenza fra eulero ed eulero-cromer? */
struct configurazione *eulero_cromer(double dt, double omega2, 
			      struct configurazione *p_x_e_v_old)
{
  static struct configurazione x_e_v_new;
  x_e_v_new.v = (*p_x_e_v_old).v + oa_forza(omega2, (*p_x_e_v_old).x) * dt;
  x_e_v_new.x = p_x_e_v_old->x + x_e_v_new.v * dt;
  return &x_e_v_new;
}

/*****************************************************************************/
double energia(double omega2, struct configurazione *p_x_e_v)
{

/* Calcoliamo in realta' 2 E / m */

  double l_energia;
  l_energia = (*p_x_e_v).v * (*p_x_e_v).v 
    + oa_energia_potenziale(omega2, (*p_x_e_v).x);
  return l_energia;
}

/*****************************************************************************/
double oa_forza(double omega2, double x)
{
  return - omega2 * x;
}

/*****************************************************************************/
double oa_energia_potenziale(double omega2, double x)
{
  return omega2 * x * x;
}

/*****************************************************************************/
long my_read_long(char *messaggio_print)
{
  long parola_letta;

  printf("# %s\n",messaggio_print); fflush(stdout);
  scanf("%ld",&parola_letta);

  return parola_letta;
}

/*****************************************************************************/
double my_read_double(char *messaggio_print)
{
  double parola_letta;

  printf("# %s\n",messaggio_print); fflush(stdout);
  scanf("%lg", &parola_letta);

  return parola_letta;
}