#include<iostream>
#include<math.h>
#include<stdlib.h>
#include<fstream>
#include<bibli_fonctions.h>
using namespace std;

double m1 = 1.e-5, m2 = 2.5e-5;

double H2(double* q){
  return (q[2]*q[2]+q[3]*q[3]+m1*m1*q[0]*q[0]+m2*m2*q[1]*q[1])/6.;
}

void systeme(double* q, double t, double* qp, int n){
  double H = sqrt(H2(q));
  qp[0] = q[2];
  qp[1] = q[3];
  qp[2] = -3.*H*q[2] - m1*m1*q[0];
  qp[3] = -3.*H*q[3] - m2*m2*q[1];
}

void euler(void (*syst)(double*,double,double*,int), double* q, double t,
	   double dt, int n){
  int i;
  double* qp = (double*)malloc(n * sizeof(double));
  syst(q,t,qp,n);
  for(i = 0; i < n; i++)
    q[i] = q[i] + dt * qp[i];
  free(qp);
}

int main(){
  int i, n = 4, Nt = 1000;
  double f0 = 30., t = 0., tfin = 1.5*f0/m2, dt = (tfin-t)/(Nt-1);
  double* q = (double*)malloc(n * sizeof(double));
  fstream fich;

  q[0] = f0; q[1] = f0; q[2] = 0.; q[3] = 0.;
  fich.open("inflation.res", ios::out);
  for(i = 0; i < Nt; i++){
    fich << t << " " << q[0] << " " << q[1] << endl;
    euler(systeme,q,t,dt,n);
    t += dt;
  }
  fich.close();
  free(q);
  ostringstream pyth;
  pyth << "A = loadtxt('inflation.res')\n"
       << "plot(A[:,0],A[:,1])\n"
       << "plot(A[:,0],A[:,2])\n";
  make_plot_py(pyth);
  return 0;
}