#include <stdio.h>
#include "patch.h"
#include <math.h>

//#define DEBUG_OUTPUT

/**************************************************************
 *
 *											  Constants 
 *
 **************************************************************/

const double R = 8.3144;
const double F = 96485;

/**************************************************************
 *
 *                   Diagnostic Functions
 *
 **************************************************************/

#ifdef DEBUG_OUTPUT

void debugDump(TParameter p,TStimulus s,TModel m) {
	FILE *file;
	if((file = fopen("debug.txt","wt")) != NULL) {
		fprintf(file,"PATAMETERS:\n");
		fprintf(file,"INTEGRATION: TSPAN = [%.1fms %.1fms], h = %.1fus, ES = %d\n",p.t0*1e3,p.t1*1e3,p.h*1e6,p.EarlyStop);
		fprintf(file,"SAVE       : SAVE = %d, Fs = %d\n\n",p.save,p.Fs);

		fprintf(file,"MODEL:\n");
		fprintf(file,"ELECTRICAL: [Cn = %.2fpF Ci = %.2fpF Cm = %.2fpF Ril = %.1fMohm eNa = %.2fmV eK = %.2fmV\n",
			m.Cn*1e12,m.Ci*1e12,m.Cm*1e12,m.Ril*1e-6,m.eNa*1e3,m.eK*1e3);
		fprintf(file,"NODAL CONDUCTANCE     : gNaf = %.2fnS, gNap = %.2fnS, gKf = %.2fnS, gKs = %.2fnS\n",
			m.gNode.gNaf*1e9,m.gNode.gNap*1e9,m.gNode.gKf*1e9,m.gNode.gKs*1e9);
		fprintf(file,"INTER-NODAL CONDUCTACE: gNaf = %.2fnS gKs = %.2fnS gL = %.2fnS eL = %.2fmV\n",
			m.gInter.gNaf*1e9,m.gInter.gKs*1e9,m.gInter.gL*1e9,m.gInter.eL*1e3);
		
		fprintf(file,"RATE A\n");
		fprintf(file,"am : %10.2f  %10.2f %10.2f\n",m.kin.am.A,m.kin.am.B,m.kin.am.C);
		fprintf(file,"ah : %10.2f  %10.2f %10.2f\n",m.kin.ah.A,m.kin.ah.B,m.kin.ah.C);
		fprintf(file,"ap : %10.2f  %10.2f %10.2f\n",m.kin.ap.A,m.kin.ap.B,m.kin.ap.C);
		fprintf(file,"an : %10.2f  %10.2f %10.2f\n",m.kin.an.A,m.kin.an.B,m.kin.an.C);
		fprintf(file,"as : %10.2f  %10.2f %10.2f\n",m.kin.as.A,m.kin.as.B,m.kin.as.C);
		
		fprintf(file,"RATE B\n");
		fprintf(file,"bm : %10.2f  %10.2f %10.2f\n",m.kin.bm.A,m.kin.bm.B,m.kin.bm.C);
		fprintf(file,"bh : %10.2f  %10.2f %10.2f\n",m.kin.bh.A,m.kin.bh.B,m.kin.bh.C);
		fprintf(file,"bp : %10.2f  %10.2f %10.2f\n",m.kin.bp.A,m.kin.bp.B,m.kin.bp.C);
		fprintf(file,"bn : %10.2f  %10.2f %10.2f\n",m.kin.bn.A,m.kin.bn.B,m.kin.bn.C);
		fprintf(file,"bs : %10.2f  %10.2f %10.2f\n",m.kin.bs.A,m.kin.bs.B,m.kin.bs.C);
		
		fprintf(file,"INITIAL CONDITIONS:\n");
		fprintf(file,"POTENTIALS  : En = %.2fmV Ei = %.2fmV\n",m.E0*1e3,m.Ei0*1e3);
		fprintf(file,"NODAL       : m = %.4f h = %.4f p = %.4f n = %.4f s = %.4f\n",
			m.Knode0.m,m.Knode0.h,m.Knode0.p,m.Knode0.n,m.Knode0.s);
		fprintf(file,"INTER-NODAL : m = %.4f h = %.4f s = %.4f\n",m.Kinter0.m,m.Kinter0.h,m.Kinter0.s);

		fprintf(file,"STIMULUS: \n");
		fprintf(file,"RAW DATA : %d %e %e %e %e %e\n",s.type,s.Is,s.Ts,s.Ip,s.Ip,s.Tisi);
		fprintf(file,"STIMULUS : ");
		switch(s.type) {
		case RECT   : fprintf(file,"PULSE Is = %.2fnA Ts = %.2fms\n",s.Is*1e9,s.Ts*1e3); break;
		case RAMP   : fprintf(file,"RAMP  Is = %.2fnA Ts = %.2fms\n",s.Is*1e9,s.Ts*1e3); break;
		case CPULSE : fprintf(file,"C-PULSE Ic = %.2fnA Tc = %.2fms\n Tisi = %.2fms Is = %.2fnA Ts = %.2fms\n",
										s.Ip*1e9,s.Tp*1e3,s.Tisi*1e3,s.Is*1e9,s.Ts*1e3); break;
		case PRECT  : fprintf(file,"R-PREPULSE Ip = %.2fnA Tp = %.2fms\n Tisi = %.2fms Is = %.2fnA Ts = %.2fms\n",
										s.Ip*1e9,s.Tp*1e3,s.Tisi*1e3,s.Is*1e9,s.Ts*1e3); break;
		}

		fclose(file);
	} 
}

#endif

/**************************************************************
 *
 *                     Stimuli Functions
 *
 **************************************************************/

double stimulus(double t,TStimulus s) {
	switch(s.type) {
		case RECT      : 
			if(t < 0) 					return 0 + s.Idc;
			else if(t <= s.Ts)	return s.Is + s.Idc;
			else								return 0 + s.Idc;
		case RAMP      :
			if(t < 0)						return 0 + s.Idc;
			else if(t <= s.Ts)	return s.Is*(t/s.Ts) + s.Idc;
			else								return 0 + s.Idc;
		case CPULSE    : {
			double ic;
			double is;
			if(t < 0) 								ic = 0;
			else if(t <= s.Tp)				ic = s.Ip;
			else											ic = 0;
			if(t-s.Tisi < 0) 					is = 0;
			else if(t-s.Tisi <= s.Ts)	is = s.Is;
			else											is = 0;
			return is+ic+s.Idc;
										 } break;
		case PRECT     : {
			if(t < 0)								return 0 + s.Idc;
			else if(t <= s.Tp)			return s.Ip + s.Idc;
			else if(t-s.Tp <= s.Ts)	return s.Is + s.Idc;
			else										return 0 + s.Idc;
										 } break;
		case EXPR			 : {
			if(t < 0)								return 0 + s.Idc;
			else										return s.Is*(1-exp(-(t/s.Ts))) + s.Idc;
										 } break;
		case FPRECT1    : {
			if(t < 0)               return 0;
			else if(t <= s.Ts*s.Tp)      return s.Ip*s.Is;
			else if((s.Ts*s.Tp < t) && (t <= s.Ts)) return s.Is;
			else                    return 0;
										 } break;
		case FPRECT2    : {
			if(t < 0)               return 0;
			else if(t <= s.Tp)      return s.Ip*s.Is;
			else if((s.Tp < t) && (t <= s.Ts+s.Tp)) return s.Is;
			else                    return 0;
										 } break;
		default				 : return 0;
	}
}

/**************************************************************
 *
 *                      Core Functions
 *
 **************************************************************/

double rate1(double E,TRate r) { return 1e3*r.A*(E-r.B)/(1-exp((r.B-E)/r.C)); }
double rate2(double E,TRate r) { return 1e3*r.A*(r.B-E)/(1-exp((E-r.B)/r.C)); }
double rate3(double E,TRate r) { return 1e3*r.A/(1+exp((r.B-E)/r.C)); }

double dmdt(double m,double E,TKinetics k) { return rate1(E,k.am)*(1-m)-rate2(E,k.bm)*m; }
double dhdt(double h,double E,TKinetics k) { return rate2(E,k.ah)*(1-h)-rate3(E,k.bh)*h; }
double dpdt(double p,double E,TKinetics k) { return rate1(E,k.ap)*(1-p)-rate2(E,k.bp)*p; }
double dndt(double n,double E,TKinetics k) { return rate1(E,k.an)*(1-n)-rate2(E,k.bn)*n; }
double dsdt(double s,double E,TKinetics k) { return rate1(E,k.as)*(1-s)-rate2(E,k.bs)*s; }

TIon_i Iion_i(double Ei,TKinetics_i k,TModel model) {
	TIon_i i;
	i.iNaf  = model.gInter.gNaf*k.m*k.m*k.m*k.h*(Ei-model.eNa);
	i.iKs   = model.gInter.gKs*k.s*(Ei-model.eK);
	i.iL    = model.gInter.gL*(Ei-model.gInter.eL);
	return i;
}

TIon_n Iion_n(double En,TKinetics_n k,TModel model) {
	TIon_n i;
	i.iNaf  = model.gNode.gNaf*k.m*k.m*k.m*k.h*(En-model.eNa);
	i.iNap  = model.gNode.gNap*k.p*k.p*k.p*(En-model.eNa);
	i.iKf   = model.gNode.gKf*k.n*k.n*k.n*k.n*(En-model.eK);
	i.iKs   = model.gNode.gKs*k.s*(En-model.eK);
	return i;
}

void prnResult(FILE *file,
							 double t,double Istim,double En,double Ei,
							 TKinetics_n kNode,TKinetics_i kInter,
							 TIon_n In,TIon_i Ii) {
	fprintf(file,"%.3e  %.3e %.2e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n",
						   t*1e3,En*1e3,Ei*1e3,Istim*1e9,
							 kNode.m,kNode.h,kNode.p,kNode.n,kNode.s,		
							 kInter.m,kInter.h,kInter.s,
							 In.iNaf*1e9,In.iNap*1e9,In.iKf*1e9,In.iKs*1e9,
							 Ii.iNaf*1e9,Ii.iKs*1e9,Ii.iL*1e9);
}

void prnX0(double En,double Ei,TKinetics_n kNode,TKinetics_i kInter) {
	FILE *file;
	if((file = fopen("x0.out","wt")) != NULL) {
		fprintf(file,"%e %e %e %e %e %e %e %e %e %e",
			En,Ei,kNode.m,kNode.h,kNode.p,kNode.n,kNode.s,kInter.m,kInter.h,kInter.s);
		fclose(file);
	}
}
int simPatch(TModel m,TStimulus s,TParameter p) {
	TIon_i Ii;
	TIon_n In;
	double Istim;
	TKinetics_i kInter = m.Kinter0;
	TKinetics_n kNode  = m.Knode0; 
	double En0;
	double Ei0;
	double En1 = m.E0;
	double Ei1 = m.Ei0;
	double dEn;

	double t = p.t0;
	int		 step = p.Fs;
	int		 stop = 0;
	int		 ap   = 0;
	int    state = IDLE;
	double tAP  = -1;

	FILE   *file;

#ifdef DEBUG_OUTPUT
	debugDump(p,s,m);
#endif
	if(p.save) file = fopen("data.out","wt");
	else file = NULL;

	while((t <= p.t1) && !stop) {
		Istim = stimulus(t,s);
		Ii = Iion_i(Ei1,kInter,m); 
		In = Iion_n(En1,kNode,m); 
		dEn = -(In.iKf+In.iKs+In.iNaf+In.iNap+Istim - (Ei1 - En1)/m.Ril)/(m.Cn + m.Cm);
		En0 = En1+p.h*dEn;
		Ei0 = Ei1-p.h*(Ii.iKs+Ii.iL+Ii.iNaf + (Ei1 - En1)/m.Ril - m.Cm*dEn)/m.Ci;
		kInter.m = kInter.m + p.h*dmdt(kInter.m,Ei1*1e3,m.kin);
		kInter.h = kInter.h + p.h*dhdt(kInter.h,Ei1*1e3,m.kin);
		kInter.s = kInter.s + p.h*dsdt(kInter.s,Ei1*1e3,m.kin);
		kNode.m = kNode.m + p.h*dmdt(kNode.m,En1*1e3,m.kin);
		kNode.h = kNode.h + p.h*dhdt(kNode.h,En1*1e3,m.kin);
		kNode.p = kNode.p + p.h*dpdt(kNode.p,En1*1e3,m.kin);
		kNode.n = kNode.n + p.h*dndt(kNode.n,En1*1e3,m.kin);
		kNode.s = kNode.s + p.h*dsdt(kNode.s,En1*1e3,m.kin);

		if((step == p.Fs) && (file != NULL) && p.save) {
			prnResult(file,t,Istim,En1,Ei1,kNode,kInter,In,Ii);
			step = 1;		
		}

		//Update before next step
		En1 = En0; Ei1 = Ei0; step++; t = t + p.h;
		
		//Detect action potential
		if(state == IDLE) {
			if(En1 > -30e-3) {
				state = DETECTED;
				ap++;
				if(p.EarlyStop == ap) stop = 1;
			}		
		} else {
			if(En1 < -30e-3) {
				state = IDLE;
			}
		}		
	}

	if((file != NULL) && p.save) {
		prnResult(file,t,Istim,En1,Ei1,kNode,kInter,In,Ii);
	}
	if(p.save)
		prnX0(En1,Ei1,kNode,kInter);

	if(file != NULL) fclose(file);
	return ap;
}