#include<iostream>
#include<fstream>
#include<math.h>
#include<string.h>
#include<stdio.h>

#include "../def.h"
#include "../model.h"
#include "../nr.h"

using namespace std;

//modules

//Begin definition of module prototypes

void tabulateValues (Glob *globs,GlobExp *ex,double t0, double t1, double *t, 
		     long n, double *state, double **y, double ***dmds, 
		     double ***dmdp, double **dyds,double **dydp);

void simInit(GlobExp *ex,Glob *globs)
{
  long i,j;
  long nvar=ex->nvar;
  long nPoints=ex->nPoints;
  long idum=-42;
  double *t=dvector(1,2);
  double **y=dmatrix(1,2,1,nvar);
  
   if (ex->y0) 
    {
      for (i=1; i<=nvar; ++i) 
	ex->yTry[1][i]=ex->y0[i];
    } 
   else
     {
        for (i=1; i<=nvar; ++i) 
	  {
	    ex->yTry[1][i]=DefYValues[i-1];
	  }
     }

   //integrating 

   for(i=2;i<=nPoints;i++)
     {
       for(j=1;j<=nvar;j++)
	 ex->yTry[i][j]=ex->yTry[i-1][j];
       t[1]=ex->mesh[i-1];
       t[2]=ex->mesh[i];
       tabulateValues(globs,ex,t[1],t[2],t,2,ex->yTry[i],y,NULL,NULL,NULL,NULL);
      for(j=1;j<=nvar;j++)
	ex->yTry[i][j]+=globs->pert*ex->yTry[i][j]*gasdev(&idum);
     }
   


   free_dvector(t,1,2);
   free_dmatrix(y,1,2,1,nvar);
}
