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


void Einstein_de_Sitter (double) ;
void FriedmannLemaitre (double, double, double, double, double) ;


void main (int argc, char *argv[])
{
  if (argc < 5) {
    puts ("") ;
    puts ("Ce programme prend cinq arguments :") ;
    puts ("cosmod <H0> <Omega mat> <Omega vide> <dt>") ;
    puts ("où :\tH0 : constante de Hubble en (km/s)/Mpc") ;
    puts ("\tOmega mat : densité de matière") ;
    puts ("\tOmega vide : densité du vide") ;
    puts ("\tTo : température du RFC (2.736 K) ou 0 pour négliger le rayonnement") ;
    puts ("\tdt : intervalle d'intégration, en Gyr (p. ex. 0.1)") ;
    puts ("") ;
    exit (1) ;
  }

  // Lecture des arguments passés au programme
  double H0 = atof (argv[1]) ;
  double omega_mat = atof (argv[2]) ;
  double omega_vide = atof (argv[3]) ;
  double T0 = atof (argv[4]) ;
  double dt = atof (argv[5]) ;

  // Calcul du modèle par appel de la fonction FriedmannLemaitre ()
  // (voir à la fin)
  FriedmannLemaitre (H0, omega_mat, omega_vide, T0, dt) ;
}


void Einstein_de_Sitter (double H0)
{
  double A, step, R, dR, t ;

  A = 4.437e-36 ;
  t = 0.0 ;
  R = 1.0 ;

  double yr = 365.25 * 86400.0 ;
  double Gyr = 1e9 * yr ;
  step = -0.1 * Gyr ; // Pas d'intégration

  while (R > 0.0) {
    printf ("%+6.1lf  %5.3lf\n", t / Gyr, R) ;
    dR = sqrt (A / R) * step ;
    R += dR ;
    t += step ;
  }
}


void FriedmannLemaitre (double H0, double omega_m, double omega_v, double T0, double step)
{
  // Intégration numérique des équations de Friedmann-Lemaître générales
  // (constante cosmologique (densité du vide) nulle ou non nulle)

  // Quelques constantes physiques et mathématiques
  double pi = 4.0 * atan (1.0) ;           // pi
  double c = 299792458.0 ;                 // vitesse de la lumière (m/s)
  double G = 6.672 / (1.e11) ;             // constante de la gravitation universelle (m^3/(s^2kg))
  double sigma = 5.67e-8 ;                 // constante de Stefan-Boltzmann (W/(m^2K^4))
  double pc = 149.6e9 * 180 * 3600 / pi ;  // parsec (m)
  double Mpc = 1.0e6 * pc ;                // mégaparsec (m)
  double yr = 365.25 * 24 * 3600 ;         // année (s)
  double Gyr = 1.0e9 * yr ;                // milliard d'années (s)
  double km = 1.0e3 ;                      // kilomètre (m)

  // Paramètres du modèle cosmologique
  H0 *= km / Mpc ;                         // constante de Hubble actuelle (s^-1)
  double H = H0 ;                          // constante de Hubble (s^-1)
  double T = T0 ;                          // température absolue (2.736 K)
  double rho_c, rho_c0 = 3 * H0 * H0 / (8 * pi * G) ;
  double rho_m, rho_m0 = omega_m * rho_c0 ;
  double rho_r, rho_r0 = sigma * T0 * T0 * T0 * T0 / c / c / c ;
  double rho_v = omega_v * rho_c0 ;
  double k, k0 = ((rho_m0 + rho_r0 + rho_v) / rho_c0 - 1) * H0 * H0 ;
  double q, z, omega ;

  printf ("rho_c0 = %lg\nrho_m0 = %lg\nrho_r0 = %lg\nrho_v0 = %lg\n",
	  rho_c0, rho_m0, rho_r0, rho_v) ;

  double R = 1.0, dR ;
  double t = 0.0, dt = step * Gyr ;

  double printinterval = 0.5 ;  // en Gyr
  int i = 0, n ;
  if (step != 0.0)
    n = fabs (step) > printinterval ? 1 : (int) (printinterval / fabs (step) + 0.5) ;
  else
    return ;
  // Boucle ; s'interrompt dès que R <= 0 ou que t dépasse 15 Gyr
  while (1.0e-6 < R && t <= 15.0 * Gyr) {
    rho_m = rho_m0 / (R * R * R) ;
    rho_r = rho_r0 / (R * R * R * R) ;
    rho_c = 3 * H * H / (8 * pi * G) ;
    omega = (rho_m + rho_r + rho_v) / rho_c ;
    T = T0 / R ;
    q = (rho_m / 2.0 + rho_r - rho_v) / rho_c ;
    k = ((rho_m + rho_r + rho_v) / rho_c - 1) ;
    z = 1.0 / R - 1.0 ;    // z = R(tr) / R(te) - 1 avec R(tr) = 1 et R(te) = R
    if ((i % n) == 0) {
      // Affichage de résultats en unités courantes (t en Gyr, H en (km/s)/Mpc, etc.
      printf ("t=%+7.3lf | ", t / Gyr) ;
      printf ("R=%6.3lf | ", R) ;
      printf ("H=%9.4lg | ", H / (km / Mpc)) ;
      printf ("q=%+6.3lf | ", q) ;
      printf ("T=%9.3lf | ", T) ;
      printf ("W=%+9.3lf | ", omega) ;
      printf ("k=%+9.3lf | ", k) ;
      printf ("z=%8.3lf\n", z) ;
      // Eventuellement : adaptation du pas d'intégration
      // if (fabs (dt) > 0.01 * (1 / H)) { dt *= 0.01 , n /= n >= 100 ? 100 : 1 , i = 0 ; }
    }
    // Nouvelle valeur du facteur d'échelle R
    dR = sqrt (8.0 * pi / 3.0 * G * (rho_m0 / R + rho_r0 / R / R + rho_v * R * R) - k0) * dt ;
    H = dR / dt / R ;
    R += dR ;
    t += dt ;
    ++ i ;
  }
}