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


using namespace std;

//parameter values for Gingrich and Byrne model in the journal of neurophysiology

// 1)  Fig.2 Simulations of the maintenance of short-term facilitation (STF) of nondepressed synapses
//             Control     ::  total_t=592, ISI =20, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=false, PKA_ON=false, PKC_ON=false  PKC_OnSetTime=430, PKA_OnSetTime=430
//             5-HT        ::  total_t=592, ISI =20, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=false, PKA_ON=true,  PKC_ON=true,  PKC_OnSetTime=430, PKA_OnSetTime=430
//             5-HT+Met    ::  total_t=592, ISI =20, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=false, PKA_ON=false, PKC_ON=true,  PKC_OnSetTime=430, PKA_OnSetTime=430
//             5-HT+Che    ::  total_t=592, ISI =20, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=false, PKA_ON=true,  PKC_ON=false, PKC_OnSetTime=430, PKA_OnSetTime=430

//2)   Fig. 3. Empirical analysis and simulation of the effects of PKC on STF of depressed synapses
//             5-HT+Met     ::  total_t=452, ISI =20, Time_FirstSpike=150, Time_SencondSpike=170, sucrose=false, PKA_ON=false, PKC_ON=true,  PKC_OnSetTime=290, PKA_OnSetTime=290

//3)   Fig. 4. Dynamics of PKC- and PKA- dependent processes in the maintenance of STF of nondepressed synapse
//             Control     ::  total_t=292, ISI =20, Time_FirstSpike=150, Time_SencondSpike=170, sucrose=false, PKA_ON=false, PKC_ON=false  PKC_OnSetTime=130, PKA_OnSetTime=130
//             5-HT        ::  total_t=292, ISI =20, Time_FirstSpike=150, Time_SencondSpike=170, sucrose=false, PKA_ON=true,  PKC_ON=true,  PKC_OnSetTime=130, PKA_OnSetTime=130
//             5-HT+Met    ::  total_t=292, ISI =20, Time_FirstSpike=150, Time_SencondSpike=170, sucrose=false, PKA_ON=false, PKC_ON=true,  PKC_OnSetTime=130, PKA_OnSetTime=130
//             5-HT+Che    ::  total_t=292, ISI =20, Time_FirstSpike=150, Time_SencondSpike=170, sucrose=false, PKA_ON=true,  PKC_ON=false, PKC_OnSetTime=130, PKA_OnSetTime=130

//4)   Fig. 5. PKC-dependent process enhanced the recovery of the RRP (sucrose modeling:: first time of application of sucrose happen at 150s , second time of applicatoin of sucrose happened at 450s. Note Calcium influx (Ica) was set to 0)
//             Control     ::  total_t=450, ISI =300, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=true, PKA_ON=false, PKC_ON=false  PKC_OnSetTime=390, PKA_OnSetTime=390
//             5-HT        ::  total_t=450, ISI =300, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=true, PKA_ON=true,  PKC_ON=true,  PKC_OnSetTime=390, PKA_OnSetTime=390
//             5-HT+Che    ::  total_t=292, ISI =300, Time_FirstSpike=150, Time_SencondSpike=450, sucrose=true, PKA_ON=true,  PKC_ON=false, PKC_OnSetTime=390, PKA_OnSetTime=390

//5)   Fig. 6. Contribution of PKC-dependent processes to the maintenance of STF was dependent on stimulus frequency. ( Fig. 6   recalibration of model : Kr=4.266 )
//             Control     ::  total_t=172, ISI =50, Time_FirstSpike=150, Time_SencondSpike=(151 (ISIs=1s at the paper); 155(ISIs=5s ISIs=5s at the paper); 10s(ISIs=10s); 20s(ISIs=20s)), sucrose=false, PKA_ON=false, PKC_ON=false  PKC_OnSetTime=130, PKA_OnSetTime=130
//             5-HT        ::  total_t=172, ISI =50, Time_FirstSpike=150, Time_SencondSpike=(151 (ISIs=1s at the paper); 155(ISIs=5s ISIs=5s at the paper); 10s(ISIs=10s); 20s(ISIs=20s)), sucrose=false, PKA_ON=true,  PKC_ON=true,  PKC_OnSetTime=130, PKA_OnSetTime=130
//             5-HT+Met    ::  total_t=172, ISI =50, Time_FirstSpike=150, Time_SencondSpike=(151 (ISIs=1s at the paper); 155(ISIs=5s ISIs=5s at the paper); 10s(ISIs=10s); 20s(ISIs=20s)), sucrose=false, PKA_ON=false, PKC_ON=true,  PKC_OnSetTime=130, PKA_OnSetTime=130
//             5-HT+Che    ::  total_t=172, ISI =50, Time_FirstSpike=150, Time_SencondSpike=(151 (ISIs=1s at the paper); 155(ISIs=5s ISIs=5s at the paper); 10s(ISIs=10s); 20s(ISIs=20s)), sucrose=false, PKA_ON=true,  PKC_ON=false, PKC_OnSetTime=130, PKA_OnSetTime=130

//6)   Fig. 7.  5-HT-induced STF of depressed synapses was dependent on the stimulus frequency  
//             5-HT (ISI=1s)   ::  total_t=199.5, ISI =1, Time_FirstSpike=150, Time_SencondSpike=151, sucrose=false, PKA_ON=true,  PKC_ON=true,  PKC_OnSetTime=169, PKA_OnSetTime=169
//             5-HT (ISI=20s)  ::  total_t=1132, ISI =20, Time_FirstSpike=150, Time_SencondSpike=170, sucrose=false, PKA_ON=true,  PKC_ON=true,  PKC_OnSetTime=530, PKA_OnSetTime=530
//orignal G-B model (1985)

/*double total_t=1250.5; // total time (seconds)

double ISI=100; // ITI=interspike interval (seconds)

double Time_FirstSpike=150; // Timing of first spike (EPSP0)
double Time_SecondSpike=250; // Timing of Second spike (EPSP1)
*/

/*double total_t=592; // total time (seconds)

double ISI=20; // ITI=interspike interval (seconds)

double Time_FirstSpike=150; // Timing of first spike (EPSP0)
double Time_SecondSpike=450; // Timing of Second spike (EPSP1)
*/

/*double total_t=530.5; // total time (seconds)

double ISI=20; // ITI=interspike interval (seconds)

double Time_FirstSpike=150; // Timing of first spike (EPSP0)
double Time_SecondSpike=170; // Timing of Second spike (EPSP1)
*/

 /*double Time_FirstSpike = 150;
 double Time_SecondSpike = 151;
 double ISI = 1;
 double total_t = 169.5;
 */

double total_t=292; // total time (seconds)

double ISI=20; // ITI=interspike interval (seconds)

double Time_FirstSpike=150; // Timing of first spike (EPSP0)
double Time_SecondSpike=170; // Timing of Second spike (EPSP1)


bool Sucrose=false;   //timing of first applicatoin of sucrose 150:  Time_FirstSpike=150;    timing of 2nd applicatoin of sucrose 450: ISI=300;  PKC_OnSetTime=390;PKA_OnSetTime=390;

bool PKC_ON=false;
bool PKA_ON=true;
//double PKC_OnSetTime=290;

double PKC_OnSetTime=130;
double PKA_OnSetTime=130;

double spikeWidth=0.0;    // spikeWidth= 0
double spikeWidth1=0.003;  //spikeWidth before 5HT induced activation of PKA
//double spikeWidth2=0.009;   //spikeWidth after 5HT induced activation of PKA
double spikeWidth2=0.0044;   //spikeWidth after 5HT induced activation of PKA 5 parameters
//double spikeWidth2=0.0073;   //spikeWidth after 5HT induced activation of PKA



double ti[100]={0};
double Release[100]={0}, Replenishment[100]={0};



double Mins=0.002, Mins1=150,Mins2=0.0;

//paramenter of Ca regulation in the G-B model
static double Kc=103.0;
static double Ta=0.001;
static double Ti=0.44;
static double Tr=15.6;

static double Ku=2907.0; //(orignal value from G-B model)
//static double Ku=407.0;  //new value 

//static double Ku=599.9713;  //new value 11 parameters 

static double Mu=28.107;  // Mu^2=709 (format different)
static double Kd=0.338;

static double Vc=2.15; //(orignal value from G-B model)
//static double Vc=0.22; //new value

//static double Vc=0.2142; //new value 11 parameters

//parameters for vesicle storage and release in the G-B model
//static double Vr=1.0; //(orignal value from G-B model)

static double Vr=1.0134;  // new value from PSO 5 parameters

//static double Kr=3.16; (orignal value from G-B model)
static double Kr=2.266; //new value from PSO 5 paramerters
//static double Kr=3.6235; //new value from PSO 11 parameters

static double Kvd=0.001; //(orignal value from G-B model)
//static double Kvd=0.003; //new value
//static double Kvd=0.0077; //new value 11 parameters

//static double Kf=21.0;  //(orignal value from G-B model)
//static double Kf=148.0;  //new value
static double Kf=34.7732; //new value from PSO 5 paramters
//static double Kf=198.7943; //new value from PSO

static double Mf=0.0008;  //(orignal value from G-B model)
//static double Mf=0.008;  //new value
//static double Mf=0.0008;  //new value 11 paramters

//static double Nf=2.83; //(orignal value from G-B model)
static double Nf=3.7970;  //new value from PSO 5 parameters
//static double Nf=24.1374;  //new value from PSO 11 paramters

static double Ks=35.0;
static double Ms=0.2276; // format issue
static double Ns=1.75;

static double Ts=213.0; //(orignal value from G-B model)
//static double Ts=3400.0; //new value
//static double Ts=340.0; //new value 11 parameters

//static double Vf=6.08; (orignal value from G-B model)
//static double Vf=3.04; //new value
static double Vf=2.4397; //new value from PSO 5 parameters
//static double Vf=2.4898; //new value from PSO 11 parameters

static double Kn=0.005;// (orignal value from G-B model)
//static double Kn=0.015; //new value





//initial conditions
double Ai,Ica, Iu, Id, If, Is, Ivd, In, Ir;
static double Ci=1.0;
static double Bi=1.0;

//static double Cf1=100.0;
static double Cf1=100.0;


double Cc=0.0;
double PVM=0.0;
double Cf=100.0;
double Cr=100.0;
double ePsp=0.0;

double y_temp[5]={0.0,0.0,100.0,100.0,0.0};   //y[0]: Ca concentration in the interior department. y[1]:PVM.    y[2]: vesicle concentration in feeding pool. y[3]: vesicle concentration in RRP. y[4]: accumulated vesicle released.
double y[5]={0.0,0.0,100.0,100.0,0.0};

//double spikeWidth=0.003;
//double spikeWidth=0.005;



//paramters for PKC dependent vesicle mobilization


double Ireplenish_temp=0, If_temp=0, Is_temp=0, Cc_temp=0,Cf_temp=0,Cr_temp=0;
double Accumulated_release=0, Accumulated_replenishment=0;

//double Gpkc=6.2;  // maximum contribution of PKC : Gpkc=6.0  
//double Tau1=-4.0, Tau2=-110;

double Gpkc=14.2031;  // maximum contribution of PKC : Gpkc=6.0   from PSO 5 paramters
//double Gpkc=11.784;  // maximum contribution of PKC from PSO 11 parameters


double Tau1=-0.4, Tau2=-67.8612; //5 parameters
//double Tau1=-0.4, Tau2=-78.4570; // 11 parameters

double Scale_Replenish=1;
//double Scale_Replenish=4.2820;

double Ireplenish=0; 
double Replenish=0;


double m_activation=0, h_inactivation=1.0;



int nu=0; // number of spike





double Cai_s(double);  //calcium function





double Diff=0,Diff_euler=0;

long double  t3=0.0001;  // numeric stepsize. t3=0.0001  denote stepsize = 0.1 ms.

double t2=0.1; // output stepsize: t2= 1 denote 1 output data per seconds.



void main()
{
	ofstream ofile("release.txt",ios::out);   //output file :: synaptic release
	ofstream ofile1("If.txt",ios::out);
	ofstream ofile2("supply_rate.txt",ios::out);  //output file::  calcium dynamics and vesicle storage and mobillization
	ofstream ofile3("supply_rate_1.txt",ios::out);
	ofstream ofile4("If_test.txt",ios::out);


	//for ( double i=2000;i<=51000;i=i+ISI)
	//for ( double i=Time_FirstSpike;i<=total_t-1000;i=i+ISI)

	//after PSO fitness
	/*
Ku=452.37;
Vc=0.1879;
Kr=3.4813;
Kvd=0.006;
Kf=274.2275;
Mf=0.0059;
Nf=35.2744;
Ts=3674.93;
Kn=0.0267;
Vf=1.6322;
*/

//	2 fold boundary . fig 1 control and Fig. 7
/*Ku=203.5    	 ;
Vc=0.2968  	;
Kr=3.1164   	;
Kvd=0.006   	;
Kf=177.1591      ;
Mf=0.004  ;
Nf=14.15;
Ts=6800;
Kn=0.0075;
Vf=4.1614 ;
*/

	// 10 fold boundary Fig 1 control and fig. 7
/*Ku=43.6589   	  ;
Vc=0.2870  	;
Kr=3.6251   	;
Kvd=0.0091   	;
Kf=117.7265         ;
Mf=0.0091   ;
Nf=12.6584    ;
Ts=34000   ;
Kn=0.0015   ;
Vf=4.2644 ;
*/

	ti[nu]=Time_FirstSpike;
	cout<<"ti["<<nu<<"]="<<Time_FirstSpike<<endl;
	nu++;

	//for ( double i=Time_SecondSpike;i<=total_t;i=i+ISI)
	for ( double i=Time_SecondSpike;i<=total_t;i=i+ISI) // original G-B model (1985)
	{

		ti[nu]=i;
		cout<<"ti["<<nu<<"]="<<i<<endl;
		nu++;
	}

	//ti[nu]=ti[nu-1]+100;
	//cout<<"ti["<<nu<<"]="<<ti[nu]<<endl;
	//nu++;
	cout<<"nu="<<nu<<endl;

	int ii=0;


	double Cai_spike=0;




	cout<< ti[0]<<endl;
	cout<<nu<<endl;
	ofile<<"Stim.No."<<" "<<"ReleaseNormToEPSP0"<<" "<<"Accumulated_release"<<" "<<"Accumulated_Replenishment"<<" "<<"ReleaseNormToEPSP1"<<endl;
	ofile1<<"t"<<" "<<"Cai"<<" "<<"N"<<" "<<"NF"<<" "<<"Rn"<<" "<<"R"<<" "<<"R+NF+N"<<" "<<"N1"<<" "<<"Rn1"<<" "<<"Pn1"<<endl;
	ofile2<<"t"<<" "<<"Ica"<<" "<<"supply_rate"<<" "<<"If"<<" "<<"Is"<<" "<<"Cc"<<" "<<"Cf"<<" "<<"Cr"<<" "<<"PKC"<<endl;
	ofile3<<"t"<<" "<<"Ica"<<" "<<"supply_rate"<<" "<<"If"<<" "<<"Is"<<" "<<"Cc"<<" "<<"Cf"<<" "<<"Cr"<<" "<<"PKC"<<endl;
	ofile4<<"t"<<" "<<"Ica"<<" "<<"Ireplenish"<<" "<<"If"<<" "<<"Is"<<endl;

	double tt=0;
	for (int i=0;i<100;i++)
	{   
		tt=tt+1;
		//ofile1<<tt<<" "<<0<<" "<<N<<" "<<Rn<<" "<<R<<" "<<Vpre1<<" "<<N1<<" "<<N2<<" "<<N3<<" "<<Rn1<<" "<<Rn2<<endl;
	}

	int Ca_Int=0;

	Mins1=Time_SecondSpike;


	for (long double t=0; t<total_t;t=t+t3)
	{
		if (t>=PKC_OnSetTime && PKC_ON==true)
		{


			m_activation=1.0-exp((t-PKC_OnSetTime)/Tau1);
			h_inactivation=exp((t-PKC_OnSetTime)/Tau2);
			
			Kn=0.015*m_activation*h_inactivation*Gpkc;
		
		}

		if (t>=PKA_OnSetTime && PKA_ON==true)
		{

            
 
			spikeWidth=spikeWidth2;

			Cai_spike=Cai_s(t);

			Ica=Cai_spike;
		}
		else if (t>=PKA_OnSetTime && PKA_ON==false && PKC_ON==true)
		{

			
				spikeWidth=spikeWidth1;
				Cai_spike=Cai_s(t);
				Ica=Cai_spike;
		}
		else
		{
				spikeWidth=spikeWidth1;
							
				Cai_spike=Cai_s(t);
				Ica=Cai_spike;
		}


			
	
	


       if (Sucrose==true )
	   {
		   Ica=0;

		   if ( abs((t-150.0))<t3)
		   {
			   y[3]=0;
			   y_temp[3]=0;
		   }

		  
	   }


		//Ica = Ai * Bi * Kc;
		Iu = Ku / (1.0 + pow(Mu / y_temp[0], 2.0));
		Id = y_temp[0] * Kd;
		If = (y_temp[2] * Kf / Cf1) * (1.0 / (1.0 + Mf / pow(y_temp[0], Nf)));
		Is = y_temp[2] * y_temp[1] / Cf1;
		Ivd = (y_temp[2] - y_temp[3]) * Kvd;

		In = (Cf1 - y_temp[2]) * Kn;
		if (In<=0)
			In=0;


		Ir = (y_temp[3] * Vr) * (Ica * Kr);

		//Is=0;

		//Ireplenish=Scale_Replenish*(If+Is+Ivd)*(1.0-y_temp[3]/100.0)*(1.0+7.9*m_activation*h_inactivation);
		Ireplenish=Scale_Replenish*(If+Is+Ivd)*(1.0-y_temp[3]/100.0)+Gpkc*m_activation*h_inactivation*(1.0-y_temp[3]/100.0);



		// 
		//if (Ireplenish<0)
		//	Ireplenish=0;


		 // sumulating sucrose experiments


		y[0] = y_temp[0]+t3*(Ica - Iu - Id) / Vc;
		y[1] = y_temp[1]+t3*((Ks / (1.0 + pow(Ms / y_temp[0], Ns))) - y_temp[1]) / Ts;
		y[2] = y_temp[2]+t3*(In - Ireplenish) / Vf;
		y[3] = y_temp[3]+t3*(Ireplenish- Ir) / Vr;
		y[4] = y_temp[4]+t3*Ir;

		y_temp[0]=y[0];
		y_temp[1]=y[1];
		y_temp[2]=y[2];
		y_temp[3]=y[3];
		y_temp[4]=y[4];

		Ireplenish_temp=Ireplenish_temp+t3*Ireplenish;
		//Accumulated_replenishment=Accumulated_replenishment+t3*Ireplenish;
		If_temp=If_temp+t3*Scale_Replenish*(If)*(1.0-y_temp[3]/100.0);
		Is_temp=Is_temp+t3*Scale_Replenish*(Is)*(1.0-y_temp[3]/100.0);
		Cc_temp=Cc_temp+t3*y[0];
		Cr_temp=Cr_temp+t3*y[3];
		Cf_temp=Cf_temp+t3*y[2];


		if (t>=149 && t<170)
			Replenish=Replenish+t3*Ireplenish;
      

		if (t>Time_FirstSpike+spikeWidth)
		{
			cout<<endl;
			cout<<"t="<<t<<"                "<<"y[0]="<<y[0]<<endl;
			cout<<"t="<<t<<"                "<<"y[1]="<<y[1]<<endl;
			cout<<"t="<<t<<"                "<<"y[2]="<<y[2]<<endl;
			cout<<"t="<<t<<"                "<<"y[3]="<<y[3]<<endl;
			cout<<"t="<<t<<"                "<<"y[4]="<<y[4]<<endl;
	   // ofile3<<t<<" "<<Ireplenish<<endl;
		//ofile4<<t/1000.0<<" "<<Cai<<" "<<N<<" "<<NF<<" "<<Rn<<" "<<b1<<" "<<Pn1<<" "<<PnF<<" "<<R<<" "<<R+NF+N<<endl;
		Time_FirstSpike=Time_FirstSpike+ISI;
		}

		//if (t>Mins)

		
		if (t>Mins)
		{
			
	    ofile3<<t<<" "<<Ica<<" "<<Ireplenish<<" "<<Scale_Replenish*(If)*(1.0-y_temp[3]/100.0)<<" "<<Scale_Replenish*(Is)*(1.0-y_temp[3]/100.0)<<" "<<y[0]<<" "<<y[2]<<" "<<y[3]<<" "<<Gpkc*m_activation*h_inactivation<<endl;
		//ofile4<<t/1000.0<<" "<<Cai<<" "<<N<<" "<<NF<<" "<<Rn<<" "<<b1<<" "<<Pn1<<" "<<PnF<<" "<<R<<" "<<R+NF+N<<endl;
		Mins=Mins+0.01;
		}
        
		if (t>Mins2 )
		{
			if (t<1.0)
			{
				Mins2=Mins2+t2;
			}
			else
			{
			
				ofile2<<t<<" "<<Ica<<" "<<Ireplenish_temp/t2<<" "<<If_temp/t2<<" "<<Is_temp/t2<<" "<<Cc_temp/t2<<" "<<Cf_temp/t2<<" "<<Cr_temp/t2<<" "<<Gpkc*m_activation*h_inactivation<<endl;
				//ofile4<<t/1000.0<<" "<<Cai<<" "<<N<<" "<<NF<<" "<<Rn<<" "<<b1<<" "<<Pn1<<" "<<PnF<<" "<<R<<" "<<R+NF+N<<endl;
				Mins2=Mins2+t2;
				Ireplenish_temp=0;
				If_temp=0;
				Is_temp=0;
				Cc_temp=0;
				Cf_temp=0;
				Cr_temp=0;
			}
		
		}
        
        
		
		

        if (t>Mins1 && t<Mins1+0.3)
		{
			
	   // ofile4<<t<<" "<<Scale_Replenish*(If)*(1.0-y_temp[3]/100.0)<<" "<<y_temp[2] * Kf / Cf1<<" "<<(1.0 / (1.0 + Mf / pow(y_temp[0], Nf)))<<" "<<(1.0-y_temp[3]/100.0)<<endl;
		ofile4<<t<<" "<<Ica<<" "<<Ireplenish<<" "<<If<<" "<<Is<<endl;
		//Mins1=Mins1+0.01;
		}

			
		for (int i=0;i<nu;i++)
		{
	
			if ((t-ti[i])>=0 && (t-ti[i])<=spikeWidth)
			{
				Release[i]=Release[i]+(y[3] * Vr) * (Ica * Kr)*t3;
			}
			if ((t-ti[i])>=0 && (t-ti[i])<=ISI)
			{
				Replenishment[i]=Replenishment[i]+Ireplenish*t3;
			}
			
		}

		if(Sucrose==true && t-150.0<-0.0001 && t-150.0> -0.0002)
		  cout<<"# of RRP before  1st surcose application"<<"     "<<"t="<<t<<"    "<<"y[3]="<<y[3]<<endl;

		if(Sucrose==true && t-300.0<0.0 && t-300.0> -0.0001 )
		  cout<<"# of RRP after  1st surcose application"<<"     "<<"t="<<t<<"    "<<"y[3]="<<y[3]<<endl;

		if(Sucrose==true && t-450.0<=0.0001 && t-450.0> 0 )
		  cout<<"# of RRP before  2nd surcose application"<<"    "<<"t="<<t<<"    "<<"y[3]="<<y[3]<<endl;



	

	
	}
	
	for (int i=0;i<nu;i++)
	{
	            
		        if (Release[0]!=0)
				{
		        

					for (int j=0; j<=i;j++)
					{
						Accumulated_release=Accumulated_release+Release[j];
						Accumulated_replenishment=Accumulated_replenishment+Replenishment[j];
					}
					cout<<"Release["<<i<<"]="<<Release[i]/Release[0]*100.0<<endl;
					cout<<"Accumulated_release="<<Accumulated_release<<endl;
					cout<<"Replenishment["<<i<<"]="<<Replenishment[i]<<endl;
					cout<<"Accumulated_replenishment="<<Accumulated_replenishment<<endl;
					ofile<<i+1<<" "<<Release[i]/Release[0]<<" "<<Accumulated_release<<" "<<Accumulated_replenishment<<" "<<Release[i]/Release[1]<<endl;
					Accumulated_release=0;
					Accumulated_replenishment=0;
				}
			
	}
	
     

	  cout<<endl;
	  cout<<endl;

	  cout<<"Release["<<"0"<<"]="<<Release[0]<<endl;
	  cout<<"Release[0]/43.0029="<<Release[0]/43.0029<<endl;
	  cout<<"Release[1]/Release[0]="<<Release[1]/Release[0]<<endl;
	  cout<<"Replenished Vesicle between 1st and 2nd stimulus"<<endl;
	  cout<<"Replenished="<<Replenish<<endl;
	  cout<<"Replenished/43.0029="<<Replenish/43.0029<<endl;
	 // string s;
     // cin >> s;
	//  cout<< s;

	  system("pause");



	double temp_If=0;
    for (double x=0;x<100;x=x+0.1)
	{   
        for (double y=0;y<100;y=y+0.1)
		{
			
			//ofile1<<x<<" "<<1.0 / (1.0 + 0.008 / pow(1.2, x))-(1.0 / (1.0 + 0.008 / pow(0.8, x)))<<endl;
		}
		
	}




	//ofile1<<endl;
    



}




double Cai_s(double t)
{

	double temp=0;

	

	for (int i=0;i<nu;i++)
	{	

		double x=t-ti[i];
		if (x>=0 && x<=spikeWidth)
		{
			Ai = 1.0 - exp(-x / Ta);
			Bi = Ci * exp(-x / Ti);
			break;
		}
		else
		{
			Ai = 0.0;
			//Ci = 1.0 - (1.0 - Bi) * exp(-(x - spikeWidth) / Tr);
		}
	}
		
	temp=Ai * Bi * Kc;

	return temp;
	}