#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <fstream>
                                                                                          
// Parameters for changing Eleak
double swtch=1,swtch2=1,x1=0,x2=100000,Y1=0,Y2=6,aa1=0,aa2=0,hh1=0,hh2=0;

// general parameters
double C=20;
double gL=2.8,gsynE=.1,gsynI=.60;
double EsynE=-10,EsynI=-75;
double EK=-85,ENa=50;
double Eleak5=-50;
//double Eleak1=-50, Eleak2=-57, Eleak3=-62, Eleak4=-67, Eleak5=-65;
// double beta=5,E0=-57; -55 + 2 is also good
double beta=4.5,E0=-54.5;
double Eleak1=-24;
double Eleak2=E0;
double Eleak3=E0-2;
//double Eleak3=E0-1*beta;
double Eleak4=E0-2*beta;
// parameters of NaP and KA currents
//double v12n=-30,v12hp=-55,v12mp=-40,v12hk=-92,v12mk=-22;
//double kn=-4,khp=10,kmp=-6,khk=10,kmk=-6;
double v12n=-30,v12hp=-59.,v12mp=-40,v12hk=-92,v12mk=-22;
double kn=-4,khp=10,kmp=-6,khk=10,kmk=-6;
double tauhp=5000,gNaP=5,tauhk=300,gka=0;//140;
double g1=5,g2=5;

// parameters of output function
double lowerlimE=-50,upperlimE=0;
double lowerlimI=-50,upperlimI=0;

inline double heav(double x) { return x>0?1:0; }
inline double foute(double v) { return heav(v-lowerlimE)*(v-lowerlimE)/(upperlimE-lowerlimE)+heav(v-upperlimE)*(1-(v-lowerlimE)/(upperlimE-lowerlimE)); }
inline double fouti(double v) { return heav(v-lowerlimI)*(v-lowerlimI)/(upperlimI-lowerlimI)+heav(v-upperlimI)*(1-(v-lowerlimI)/(upperlimI-lowerlimI)); }

double a=0.8;
double a1=a;
double a2=a;
double a3=a;
double a4=a;
double a5=2;

double a12=a1;
double a13=a1;
double a14=a1;

double a21=a2;
double a23=a2;
double a24=a2;

double a31=a3;
double a32=a3;
double a34=a3;

double a41=a4;
double a42=a4;
double a43=a4;

double brek=1;
// Initial Conditions
double v1=-60,v2=-60,v3=-60,v4=-60,v5=-60;
double h1=0.35,h2=0.35,h3=0.35,h4=0.35,h5=0.35,count=0,hk1=0.05,hk2=0.05,hk3=0,hk4=0;

// HH gating functions
inline double ninf(double v) { return 1./(1+exp((v-v12n)/kn)); }
inline double mpinf(double v) { return 1./(1+exp((v-v12mp)/kmp)); }
inline double hpinf(double v) { return 1./(1+exp((v-v12hp)/khp)); }
inline double tauinf(double v) { return tauhp/cosh((v-v12hp)/(khp*2)); }
inline double mkinf(double v) { return 1./(1+exp((v-v12mk)/kmk)); }
inline double hkinf(double v) { return 1./(1+exp((v-v12hk)/khk)); }

// Leak
double coeff1=0.05;
inline double block(double count) {return (1-.15*(heav(count-hh1*1000)-heav(count-hh2*1000))) ;}
inline double leak1(double count) { return Eleak1*(count*(Y2-Y1)/x2+Y1); }
inline double leak2(double count) { return Eleak2*(count*(Y2-Y1)/x2+Y1); }
inline double leak3(double count) { return Eleak3*(count*(Y2-Y1)/x2+Y1); }
inline double leak4(double count) { return block(count)*Eleak4*(count*(Y2-Y1)/x2+Y1); }


// Persisten sodium (NaP), Potassium (Kdr), Adaptive (AD), and leak (L) current descriptions
inline double gN(double count) {return count*(g2-g1)/x2+g1; }
inline double Inap1(double v,double h1) { return gN(count)*mpinf(v)*h1*(v-ENa); }
inline double Inap2(double v,double h2) { return gN(count)*mpinf(v)*h2*(v-ENa); }
inline double Il1(double v) { return gL*(v-leak1(count)); }
inline double Il2(double v) { return gL*(v-leak2(count)); }
inline double Il3(double v) { return gL*(v-leak3(count)); }
inline double Il4(double v) { return gL*(v-leak4(count)); }
inline double Il5(double v) { return gL*(v-Eleak5); }

inline double aa(double count) { return count*(aa2-aa1)/x2+aa1; }


inline double Ika1(double v,double hk1) { return gka*mkinf(v)*hk1*(v-EK); }
inline double Ika2(double v,double hk2) { return gka*mkinf(v)*hk2*(v-EK); }

double dt=1,total=100000; 


// Differential equations 1 fast (V) and slow (h or m) for each cell
void step(double t)
{
	double v1p=(-Inap1(v1,h1)-Ika1(v1,hk1)-Il1(v1)-gsynE*brek*(aa(count)*foute(v2)+aa(count)*foute(v3)+aa(count)*foute(v4))*(v1-EsynE))/C;
	double v2p=(-Inap2(v2,h2)-Ika2(v2,hk2)-Il2(v2)-gsynE*brek*(aa(count)*0*foute(v1)+aa(count)*foute(v3)+aa(count)*foute(v4))*(v2-EsynE))/C;
	double h1p=(hpinf(v1)-h1)/tauinf(v1);
	double h2p=(hpinf(v2)-h2)/tauinf(v2);
	double hk1p=(hkinf(v1)-hk1)/tauhk;
	double hk2p=(hkinf(v2)-hk2)/tauhk;

	double v3p=(-Inap1(v3,h3)-Ika1(v3,hk3)-Il3(v3)-gsynE*brek*(aa(count)*0*foute(v1)+aa(count)*foute(v2)+aa(count)*foute(v4))*(v3-EsynE))/C;
	double v4p=(-Inap2(v4,h4)-Ika2(v4,hk4)-Il4(v4)-gsynE*brek*(aa(count)*0*foute(v1)+aa(count)*foute(v2)+aa(count)*foute(v3))*(v4-EsynE))/C;
	double h3p=(hpinf(v3)-h3)/tauinf(v3);
	double h4p=(hpinf(v4)-h4)/tauinf(v4);
	double hk3p=(hkinf(v3)-hk3)/tauhk;
	double hk4p=(hkinf(v4)-hk4)/tauhk;

	double v5p=(-Inap1(v5,h5)-Ika1(v1,hk1)-Il5(v5)-gsynE*brek*(a5*foute(v1)+a5*foute(v2)+a5*foute(v3)+a5*foute(v4))*(v5-EsynE))/C;
	double h5p=(hpinf(v5)-h5)/tauinf(v5);

	v1+=dt*v1p; v2+=dt*v2p; h1+=dt*h1p; h2+=dt*h2p; hk1+=dt*hk1p; hk2+=dt*hk2p;
	v3+=dt*v3p; v4+=dt*v4p; h3+=dt*h3p; h4+=dt*h4p; hk3+=dt*hk3p; hk4+=dt*hk4p;
	v5+=dt*v5p; h5+=dt*h5p;
	count=t;
}

using namespace std;

int main(int argc,char **argv)
{
	int freq=10,poi=0,ini=0,vF=0,vE=0,du=0,dum=0,vZ=0;
	double th=-43,ed=0,ef=0,per=0;
	for(int i=1;i<argc;i++)
	{
		if(strcmp(argv[i],"-dt")==0) dt=atof(argv[++i]);
		else if(strcmp(argv[i],"-f")==0) freq=atoi(argv[++i]);
		else if(strcmp(argv[i],"-p")==0) poi=1;
		else if(strcmp(argv[i],"-i")==0) ini=1;
		else if(strcmp(argv[i],"-v")==0) vF=1;
		else if(strcmp(argv[i],"-w")==0) vE=1;
		else if(strcmp(argv[i],"-z")==0) vZ=1;
		else if(strcmp(argv[i],"-z")==0) du=1;
		else if(strcmp(argv[i],"-l")==0) total=x2=atof(argv[++i])*1000;
		else if(strcmp(argv[i],"-a")==0) { aa1=atof(argv[++i]); aa2=atof(argv[++i]); }
		else if(strcmp(argv[i],"-g")==0) { g1=atof(argv[++i]); g2=atof(argv[++i]); }
		else if(strcmp(argv[i],"-y")==0) { Y1=atof(argv[++i]); Y2=atof(argv[++i]); }
		else if(strcmp(argv[i],"-h")==0) { hh1=atof(argv[++i]); hh2=atof(argv[++i]); }
		else { cerr<<"unknown option"<<endl; return 1; }
	}
	double tpre=-1,t2pre=-1;
	if(ini) ifstream("bart.ini")>>v1>>v2>>v3>>v4>>h1>>h2>>h3>>h4;
	double v1m=v1,v2m=v2,v3m=v3,v4m=v4;
	for(int i=0;i*dt<total;i++)
	{
		double vpre=v2,vpre4=v4,vpre3=v3,t=i*dt;
		step(t);
		v1m=max(v1m,v1);v2m=max(v2m,v2);v3m=max(v3m,v3);v4m=max(v4m,v4);
		if(!poi && i%freq==0) cout<<v1<<'\t'<<v2<<'\t'<<v3<<'\t'<<v4<<'\t'<<t<<'\t'<<h1<<'\t'<<h2<<'\t'<<h3<<'\t'<<h4<<'\t'<<t<<endl;
		if(vF && poi && vpre<=th && v2>th) 
		{
			if(t2pre>0)
			tpre=t;
		}
		if(vF && poi && vpre>=th && v2<th) 
		{
			if(tpre>0) cout<<(t-t2pre)<<'\t'<<((t-tpre)/(t-t2pre))<<'\t'<<((t-t2pre)/1000)<<'\t'<<v3m<<'\t'<<t<<'\t'<<v1m+v2m+v3m+v4m<<endl;
			t2pre=t; v2m=v2;
		}
		if(vE && poi && vpre4<=th && v4>th) 
		{
			if(t2pre>0)
			tpre=t;
		}
		if(vE && poi && vpre4>=th && v4<th) 
		{
			if(tpre>0) cout<<(t-t2pre)<<'\t'<<((t-tpre)/(t-t2pre))<<'\t'<<((t-t2pre)/1000)<<'\t'<<v3m<<'\t'<<t<<'\t'<<v1m+v2m+v3m+v4m<<endl;
			t2pre=t; v4m=v4;v3m=v3;v2m=v2;v1m=v1;
		}
		if(vZ && poi && vpre3<=th && v3>th) 
		{
			if(t2pre>0)
			tpre=t;
		}
		if(vZ && poi && vpre3>=th && v3<th) 
		{
			if(tpre>0) cout<<(t-t2pre)<<'\t'<<((t-tpre)/(t-t2pre))<<'\t'<<((t-t2pre)/1000)<<'\t'<<v3m<<'\t'<<t<<endl;
			t2pre=t; v3m=v3;
		}

	}
	ofstream("bart.ini")<<v1<<'\t'<<v2<<'\t'<<v3<<'\t'<<v4<<'\t'<<v5<<'\t'<<h1<<'\t'<<h2<<'\t'<<h3<<'\t'<<h4<<'\t'<<h5<<endl;
	return 0;
}