#include "stdio.h"
#include "stdlib.h"
#include <fstream>
#include <iostream>
#include "math.h"
#include "string.h"
#include "float.h" //new!

using namespace std;
double firstDerivative ( double H0, double Om, double Or, double OL, double power, double Ok , double a );
int sign = 1.; //expansion or contraction?

int main ( int argc, char** argv ) {
  
  double scaleFactor[2], m, L, r, k, w;
  string line; double temp, Hubble, step, aMax;
  
  int i = 0;
  ifstream paramFile; paramFile.open("omegas.txt",ios::in);
  while ( paramFile >> line >> temp ) {
    switch ( i ) {
    case 0:
      L = temp;
      break;
    case 1:
      m = temp;
      break;
    case 2:
      r = temp;
      break;
    case 3:
      w = temp;
      break;
    case 4:
      Hubble = temp;
      break;
    case 5:
      aMax = temp;
      break;
    case 6:
      step = temp;
      break;
    }
    i++;
  } paramFile.close();
  
  Hubble *= 3.2407792700053954e-20 * 31536000.; // km/(sec.*Mpc) into yr^-1
  k = 1. - ( m + r + L );
  double p = -3. * ( 1. + w ); //for dark energy
  
  scaleFactor[0] = 1e-3; //arbitrary, hand-tuned for speed vs. accuracy
  double time = 0.;//1.709494926e-51; //Planck time in years
  double k1, k2, k3, k4;
  while ( scaleFactor[0] < (aMax+0.01) ) {

    k1 = sign * firstDerivative ( Hubble, m, r, L, p, k, scaleFactor[0] );
    k2 = sign * firstDerivative ( Hubble, m, r, L, p, k, scaleFactor[0] + (step/2.)*k1 );
    k3 = sign * firstDerivative ( Hubble, m, r, L, p, k, scaleFactor[0] + (step/2.)*k2 );
    k4 = sign * firstDerivative ( Hubble, m, r, L, p, k, scaleFactor[0] + step*k3 );
    
    scaleFactor[1] = scaleFactor[0] + (step/6.)*(k1+2.*k2+2.*k3+k4);
    
    printf("%e\t%.20f\n",time,scaleFactor[0]);
    
    time += step;
    if ( sign ) step *= 1.001; else step /= 1.001;
    if ( k1 < 2e-13 && sign && time > 1e9 ) sign = -1.;
    scaleFactor[0] = scaleFactor[1];
    
  }
  
  return 1; //as usual end of program with a return of an integer

}

double firstDerivative ( double H0, double Om, double Or, double OL, double power, double Ok , double a ) {

  return H0 * sqrt ( Om / a + Or / pow(a,2.) + OL * pow(a,2.+power) + Ok );

} // solves a-dot or da/dt for you for a given 'a,' H_0, and Omega parameters