//RodentPreBotCNeuron_main.c
//Jessica Parker, October 5, 2024
//
//This is the main function, which begins your program and guides the actions in your program.

#include "RodentPreBotCNeuron.h" 
#include "time.h" 

int main (void) {  

  int run1 = 5;  //These numbers gets written into the names of the data files. You can use any 
  int run2 = 0;  //combination of numbers according to whatever organizational system you prefer.
  int run3 = 0;  
  int run4 = 2;  
  int run5 = 0;  
  int run6 = 0;  

  //parameters defined in header file
  
  double tf1 = 4000.0;  //The amount of time the code integrates before beginning to save the data
  double tf2 = 30.0;  //The amount of time integrated in each data file.
  double tf3 = 50.0;
  double tint = 0.0001;  //Time step between data points.
 
  clock_t start;
  clock_t end; 
  float seconds; 

  int nvar = 4;  //Number of variables or number of equations
  double yy[nvar]; 
  double yy0[nvar];  //initial condition variable
  int numIP=0;   //used to read initial condition file
  double Vclamp0 = -80.0;
  
  int a, b, c, d, g; //used for loops

  FILE * f = fopen("ipVoltageClamp4D.txt", "r");  //initial condition file
  double number = 0; 
  
  while( fscanf(f, "%lf,", &number) > 0 )  //reading initial conditions 
    { 
      yy0[numIP]= number; 
      numIP++; 
    }  
  fclose(f);

  for (a=0; a<nvar; a++)  //Setting time zero to initial conditions
    { 
      yy[a] = yy0[a]; 
    } 

  for (b=0; b<1; b++)  //Change maximum limit if you want to sweep a parameter
    {
      //kmK = -8.0-1.0*b;
      //run3 = b+1;    //Uncomment if you are doing a 3D parameter sweep, where you are sweeping 3 parameters independently. This makes sure each simulation has a unique file name.
      
      for (g=0; g<1; g++)
	{
	  //gNaP = 1.6+0.2*g;
	  //run4 = g+1;  //Uncomment if you are doing a 2D or 3D parameter sweep, where you are sweeping 2 or 3 parameters independently.
      
	  for (d=0; d<3; d++)  //Change maximum limit if you want to run a 2D sweep
	    {
	      Ipmpmx = 40.0 + 40.0*d;
	      run5 = d+1;    //Uncomment if you are sweeping a parameter. This also makes sure each simulation has a unique file name.
	      
	      printf("run3 = %i, run4 = %i, run5 = %i \n",run3,run4,run5);

	      for (a=0; a<nvar; a++)  //Uncomment if you want to use the same initial conditions for each integration   
		{ 
		  yy[a] = yy0[a]; 
		}
	      //Vclamp = -80.0;
	      
	      start = clock();  
	      //integrateNW(run1,run2,run3,run4,run5,0,nvar,tf1,tint,yy);   //Call to function defined in EpisodicLocomotion_integrateNW.c
	
	      Vclamp = -90.0;
	      integrate(run1,run2,run3,run4,run5,1,nvar,tf3,tint,yy);
	      Vclamp = -80.0;
	      integrate(run1,run2,run3,run4,run5,2,nvar,tf2,tint,yy);
	      Vclamp = -40.0;
	      integrate(run1,run2,run3,run4,run5,3,nvar,tf3,tint,yy);
	      	      
	      end = clock();
	      seconds = (float)(end-start)/CLOCKS_PER_SEC;
	      
	      printf("running time: %6.6g\n", seconds); //the code as it is should take about 15 min.
	    }
	}
    }
  return(0); 
}