/* 
Uterine Single Cell Model
Author: Craig Testrow
Group: Biological Physics, The University of Manchester
Supervisor: Prof. Henggui Zhang
161214

To compile:
g++ -lm -Wall uscm.cpp -o uscm

To run the compiled program:
./uscm

*/

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <fstream>

using namespace std;
  
//Global constants and associated variables
const double R = 8314; // Molar gas constant
const double Temperature = 308.; // Temp in Kelvin (35C)
const double Faraday = 96485.; // Faraday's constant mC mM-1 // PRETEST
const double Cell_Volume = 2.65e-8; // cm3
const double Cell_Area = 7.422e-5; // cm2
const double Membrane_Capacitance = 0.0014; // mF cm-2 to give a total capacitance of about 100 pF
const double zCa = 2.0; // Calcium valence
const double zCl = -1.0; // Chlorine valence
const double zNa = 1.0; // Sodium valence
const double zK = 1.0; // Potassium valence
const double dt = 0.01; // time step (10microsecs)
const int num_cells = 1; // total number of cells
int cur_cell = 0;
const double PK(1.), PNa(0.35); // Permeabilities (rates/ease of flow) of K and Na ions

//Extracellular ion concentrations
double Ko(5.4);
double Nao(140.);
double Cao(2.5);
double Clo(146.4);
double Mgo(0.5);

//Intracellular ion concentrations
double Cai[num_cells], Nai[num_cells], Ki[num_cells], Cli[num_cells];

//SR ion concentrations
double Ca_u[num_cells], Ca_r[num_cells];

//Membrane voltage
double V[num_cells];

//Gating variables
double Ih_y[num_cells], ICaL_d[num_cells], ICaL_f1[num_cells], ICaL_f2[num_cells], ICaT_b[num_cells], ICaT_g[num_cells], INa_m[num_cells], INa_h[num_cells];
double IK1_q[num_cells], IK1_r1[num_cells], IK1_r2[num_cells], IK2_p[num_cells], IK2_k1[num_cells], IK2_k2[num_cells], IKA_s[num_cells], IKA_x[num_cells];
double V_KCa_half[num_cells], IKCa_pf[num_cells], IKCa_ps[num_cells];
double ICl_c1[num_cells], ICl_c2[num_cells], ICl_n[num_cells];

//Ryanodine receptors
double R00[num_cells],R10[num_cells],R01[num_cells],R11[num_cells];

//Current totals
double Itot;
double ICa_tot, INa_tot, IK_tot, ICl_tot;

//Other parameters
double Isus = 0.6; //sustained background current
double Istim = 0, Iinitial = 0; //reinitialised in the main program if used
double InsCa, InsNa, InsK, Insleak; // Ins currents that need to be accessed outside of function
double IsocCa, IsocNa; // Isoc currents that need to be accessed outside of function
double Iup, Itr, Irel; // SR currents that need to be accessed outside of function

//Steady states (declared here to be accessable to class_vclamp and output)
double Ih_y_steady_state;
double ICaL_d_steady_state, ICaL_f_steady_state;
double ICaT_b_steady_state, ICaT_g_steady_state;
double INa_m_steady_state, INa_h_steady_state;
double IK1_q_steady_state, IK1_r_steady_state;
double IK2_p_steady_state, IK2_k_steady_state;
double IKA_s_steady_state, IKA_x_steady_state;
double IKCa_p_steady_state;
double ICl_c_steady_state,ICl_n_steady_state;
double Inscc_fMg_steady_state;

//Time constants (declared here to be accessable to class_vclamp and output)
double Ih_y_Tau;
double ICaL_d_Tau, ICaL_f1_Tau, ICaL_f2_Tau;
double ICaT_b_Tau, ICaT_g_Tau;
double INa_m_Tau, INa_h_Tau;
double IK1_q_Tau, IK1_r1_Tau, IK1_r2_Tau;
double IK2_p_Tau, IK2_k1_Tau, IK2_k2_Tau;
double IKA_s_Tau, IKA_x_Tau;
double IKCa_pf_Tau, IKCa_ps_Tau;
double ICl_c1_Tau, ICl_c2_Tau,ICl_n_Tau;

//Kinetic cross-bridge cycle
double M, Mp, AMp, AM;
double Frac_phos, Frac_stress;

//Force
double Fp, Fx, Fa, Fs, Ft;
double lc, ls, la, lx; //lengths of cell and force elements

/* Current magnitudes */
	double Ih_mag(0.75); 
	double ICaL_mag(0.530*1.0);
	double ICaT_mag(1.333*1.0); 
	double INa_mag(0.96*1);//0.75);
	double Isoc_mag(1);
	double IsocCa_mag(0.1); 
	double IsocNa_mag(1);
	double IK1_mag(0.9*1.0);
	double IK2_mag(0.4*1.0);
	double IKA_mag(0.0075*1); 
	double IKCa_mag(10*1.0);
	double IKleak_mag(0.08);
	double Inscc_mag(1*1.0);
	double InsCa_mag(1);
	double InsNa_mag(1*1.0);
	double InsK_mag(1);
	double Insleak_mag(0.8*1);
	double ICl_mag(1); 
	double ICapump_mag(5e0); 
	double INaK_mag; // see below in time loop
	double INaCa_mag(5e0);
	double INaKCl_mag(0.02);
	double Iup_mag(1);
	double Irel_mag(1);
	double Itr_mag(1);
	

// Cell class
	class cell
{
	public:
	// Constructor to initialise all current variables, etc
	cell(void)
	{
	//Initial conditions for parameters
		
		
		for (int i=0;i<num_cells;i++)
		{
			//Initialise gating variables
			Ih_y[i]=0.00223178;
			ICaL_d[i]=0.0125799;
			ICaL_f1[i]=0.942219;
			ICaL_f2[i]=0.942218;
			ICaT_b[i]=0.169361;
			ICaT_g[i]=0.0318277;
			INa_m[i]=0.124663;
			INa_h[i]=0.371776;
			IK1_q[i]=0.07298;
			IK1_r1[i]=0.159599;
			IK1_r2[i]=0.129642;
			IK2_p[i]=0.0782046;
			IK2_k1[i]=0.995883;
			IK2_k2[i]=0.95664;
			IKA_s[i]=0.0366284;
			IKA_x[i]=0.653819;
			V_KCa_half[i]=-41.7 * log10(Cai[i]) - 128.2;
			IKCa_pf[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
			IKCa_ps[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
			ICl_c1[i]=0.000864388;
			ICl_c2[i]=0.000864392;
			ICl_n[i]=0.445644;
			
			//Initialise each element of concentration arrays
			Cai[i]=8.97e-5;
			Nai[i]=4.;
			Ki[i]=140;
			Cli[i]=130;
			Ca_u[i]=0.0273854;
			Ca_r[i]=0.0232185;
			
			//Initialise each element of membrane voltage array
			V[i]=-52.541;
			
			//Initialise ryanodine receptors
			R00[i]=6.99E-09;//0.25;
			R01[i]=7.85E-12;//.25;
			R10[i]=1.88E-11;//.25;
			R11[i]=2.11E-14;//.25;

		/*	//Initialise gating variables
			Ih_y[i]=0.00212991;
			ICaL_d[i]=0.0133195;
			ICaL_f1[i]=0.938986;
			ICaL_f2[i]=0.938985;
			ICaT_b[i]=0.179073;
			ICaT_g[i]=0.0308496;
			INa_m[i]=0.129659;
			INa_h[i]=0.360711;
			IK1_q[i]=0.0741461;
			IK1_r1[i]=0.151067;
			IK1_r2[i]=0.108987;
			IK2_p[i]=0.0794842;
			IK2_k1[i]=0.995353;
			IK2_k2[i]=0.954862;
			IKA_s[i]=0.0385646;
			IKA_x[i]=0.644685;
			V_KCa_half[i]=-41.7 * log10(Cai[i]) - 128.2;
			IKCa_pf[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
			IKCa_ps[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
			ICl_c1[i]=0.00089031;
			ICl_c2[i]=0.000890317;
			ICl_n[i]=0.440717;
			
			//Initialise each element of concentration arrays
			Cai[i]=9.43e-5;
			Nai[i]=4.;
			Ki[i]=140;
			Cli[i]=130;
			Ca_u[i]=2.64648;
			Ca_r[i]=2.24248;
			
			//Initialise each element of membrane voltage array
			V[i]=-52.1358;
			
			//Initialise ryanodine receptors
			R00[i]=1;//0.25;
			R01[i]=0;//.25;
			R10[i]=0;//.25;
			R11[i]=0;//.25;
		*/	
			//Initialise force elements
			//These do not need to be arrays as they're never used beyond single cell
			M = 0.97;
			Mp = 0.01;
			AMp = 0.01;
			AM = 1-(M+Mp+AMp);
			ls = 30.08;
			la = 25.3;
		  
		}
		
	}
	
//Current functions

/* ##############################################
 * #                                            #
 * #  Hyperpolarisation-activated Current - Ih  #
 * #                                            #
 * ############################################## */

double Ih ()
{	
	double gh = 0.0542;
	double Eh = (R * Temperature / Faraday) * log((Ko + (PNa / PK) * Nao) / (Ki[cur_cell] + (PNa / PK) * Nai[cur_cell]));
		
	// Activation
	Ih_y_steady_state = 1/(1+exp((V[cur_cell]+105.39)/8.66)); 
	if (V[cur_cell]>160)
	{Ih_y_Tau = 1/(3.5e-6*(1-0.0487*V[cur_cell]+pow(-0.0487*V[cur_cell],2)/2
		 +pow(-0.0487*V[cur_cell],3)/6)+0.04*(1+0.0521*V[cur_cell]+pow(0.0521*V[cur_cell],2)/2+pow(0.0521*V[cur_cell],3)/6));} 
	else
	{Ih_y_Tau = 1/(3.5e-6*exp(-0.0497*V[cur_cell])+0.04*exp(0.0521*V[cur_cell]));}
	Ih_y[cur_cell] = Ih_y[cur_cell] + dt * ((Ih_y_steady_state - Ih_y[cur_cell]) / Ih_y_Tau);

	return (gh * pow(Ih_y[cur_cell],1) * (V[cur_cell] - Eh));
}

/* ###################################
 * #                                 #
 * #  L-Type Calcium Current - ICaL  #
 * #                                 #
 * ################################### */
	
double ICaL (double nif)
{	
	double gCaL(0.6/* *0.4*/);//reduced 60% for estradiol
	double ECaL = 45;
	const double KmCaL(0.0006);
	
	// Activation
	ICaL_d_steady_state = 1/(1+exp(-(V[cur_cell]+22)/7));
	ICaL_d_Tau = 2.29+5.7/(1+pow((V[cur_cell]+29.97)/9,2));
	ICaL_d[cur_cell] = ICaL_d[cur_cell] + dt * ((ICaL_d_steady_state - ICaL_d[cur_cell]) / ICaL_d_Tau);

	// Inactivation
	ICaL_f_steady_state = 1/(1+exp((V[cur_cell]+33)/7));
	ICaL_f1_Tau = 12;
	ICaL_f2_Tau = 90.97-(90.97/((1+exp((V[cur_cell]+13.96)/45.38))*(1+exp(-(V[cur_cell]+9.5)/3.39)))); 
	ICaL_f1[cur_cell] = ICaL_f1[cur_cell] + dt * ((ICaL_f_steady_state - ICaL_f1[cur_cell]) / ICaL_f1_Tau);
	ICaL_f2[cur_cell] = ICaL_f2[cur_cell] + dt * ((ICaL_f_steady_state - ICaL_f2[cur_cell]) / ICaL_f2_Tau);

	const double IC50_nif =1e-6;
	double Hill=1/(1+pow((nif/IC50_nif),1));
	gCaL = gCaL*Hill;
	
	return (gCaL * pow(ICaL_d[cur_cell],2.0) * (1 / (1 + pow((Cai[cur_cell]/KmCaL),1))) * pow((0.8 * ICaL_f1[cur_cell] + 0.2 * ICaL_f2[cur_cell]),1))*(V[cur_cell] - ECaL);
}


/* ###################################
 * #                                 #
 * #  T-Type Calcium Current - ICaT  #
   #                                 #
   ################################### */

double ICaT ()
{
	double gCaT(0.0174/* *0.4*/);//reduced 60% for estradiol
	double ECaT = 42;
	
	// Activation	
	ICaT_b_steady_state = 1/(1+exp(-(V[cur_cell]+43.)/6.)); //best fit to data 14xxxx. This provides closest IV curve, despite being odd fit to activation
	ICaT_b_Tau = (0.45+ 3.9/(1+((V[cur_cell]+66)/26)*((V[cur_cell]+66)/26)));
	ICaT_b[cur_cell] = ICaT_b[cur_cell] + dt * ((ICaT_b_steady_state - ICaT_b[cur_cell]) / ICaT_b_Tau);

	// Inactivation
	ICaT_g_steady_state = 0.02+0.98/(1+exp((V[cur_cell]+72.98)/4.64));
	ICaT_g_Tau = 150 - 150/((1+exp((V[cur_cell]-417.43)/203.18))*(1+exp(-(V[cur_cell]+61.11)/8.07)));
	ICaT_g[cur_cell] = ICaT_g[cur_cell] + dt * ((ICaT_g_steady_state - ICaT_g[cur_cell]) / ICaT_g_Tau);

	return (gCaT * pow(ICaT_b[cur_cell],2) * (pow(ICaT_g[cur_cell],1)) * (V[cur_cell] - ECaT));
}


/* ###############################
 * #                             # 
 * #  Fast Sodium Current - INa  #
   #                             #
   ############################### */

double INa ()
{
	double gNa(0.12);
	double ENa = (R * Temperature / Faraday) * log(Nao / Nai[cur_cell]);
    
	// Activation
	INa_m_steady_state = 1/(1+exp(-(V[cur_cell]+35.)/9)); 
	INa_m_Tau = 0.25+7/(1+exp((V[cur_cell]+38)/10));
	INa_m[cur_cell] = INa_m[cur_cell] + dt * ((INa_m_steady_state - INa_m[cur_cell]) / INa_m_Tau);

	// Inactivation
	INa_h_steady_state = 1/(1+exp((V[cur_cell]+57)/8.5));
	INa_h_Tau = 0.9+1002.85/(1+((V[cur_cell]+47.5)/1.5)*((V[cur_cell]+47.5)/1.5));
	INa_h[cur_cell] = INa_h[cur_cell] + dt * ((INa_h_steady_state - INa_h[cur_cell]) / INa_h_Tau);		

	return (gNa * pow(INa_m[cur_cell],3) * pow(INa_h[cur_cell],1) * (V[cur_cell] - ENa));
}


/* ########################################################
 * #                                                      #
   #  Store-Operated Non-Selective Cation Current - Isoc  #
   #                                                      #
   ########################################################*/ 

double Isoc ()
{
	double ECaL = 45;
	double ECaT = 42;
	double ECa = ECaL + ECaT;
	double ENa = (R * Temperature / Faraday) * log(Nao / Nai[cur_cell]);
	double gsocCa(8.3e-5);
	double gsocNa(5.75e-4);
	const double Kmsoc(0.0001);
	
	double Psoc = 1 / (1 + Cai[cur_cell] / Kmsoc);
	IsocNa = gsocNa * Psoc * (V[cur_cell] - ENa);
	IsocCa = gsocCa * Psoc *  (V[cur_cell] - ECa);

	return (IsocCa + IsocNa);
}

/* ################################################
   #                                              #
   #  Delayed Rectifying Potassium Current - IK1  #
   #                                              #
   ################################################ */

double IK1 ()
{
	double gK1(0.2275/* *0.7*/);//reduced 30% for estradiol
	double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);

	// Activation
	IK1_q_steady_state = 1/(1+exp(-(V[cur_cell]-7.7)/23.7));
	IK1_q_Tau = (500/(1+((V[cur_cell]+60.71)/15.79)*((V[cur_cell]+60.71)/15.79)));
	IK1_q[cur_cell] = IK1_q[cur_cell] + dt * ((IK1_q_steady_state - IK1_q[cur_cell]) / IK1_q_Tau);

	// Inactivation
	IK1_r_steady_state = 1/(1+exp((V[cur_cell]+63)/6.3));
	IK1_r1_Tau = (5e3/(1+((V[cur_cell]+62.71)/35.86)*((V[cur_cell]+62.71)/35.86))); 
	IK1_r2_Tau = (3e4+2.2e5/(1+exp((V[cur_cell]+22)/4)));
	IK1_r1[cur_cell] = IK1_r1[cur_cell] + dt * ((IK1_r_steady_state - IK1_r1[cur_cell]) / IK1_r1_Tau);
	IK1_r2[cur_cell] = IK1_r2[cur_cell] + dt * ((IK1_r_steady_state - IK1_r2[cur_cell]) / IK1_r2_Tau);

	return (gK1 * pow(IK1_q[cur_cell],1) * pow((0.38 * IK1_r1[cur_cell] + 0.62 * IK1_r2[cur_cell]),1) * (V[cur_cell] - EK));
}


/* ################################################
 * #                                              #
 * #  Delayed Rectifying Potassium Current - IK2  #
   #                                              #
   ################################################ */

double IK2 ()
{
	double gK2(0.07875/* *0.5*/);//reduced 50% for estradiol
	double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
	
	// Activation
	IK2_p_steady_state = 1/(1+exp(-(V[cur_cell]-4.2)/23));
	IK2_p_Tau = 100/(1+((V[cur_cell]+64.1)/28.67)*((V[cur_cell]+64.1)/28.67));
	IK2_p[cur_cell] = IK2_p[cur_cell] + dt * ((IK2_p_steady_state - IK2_p[cur_cell]) / IK2_p_Tau);

	// Inactivation
	IK2_k_steady_state = 1/(1+exp((V[cur_cell]+21.2)/5.7));
	IK2_k1_Tau = 1e6-1e6/((1+exp((V[cur_cell]-315)/50))*(1+exp(-(V[cur_cell]+74.9)/8)));
	IK2_k2_Tau = 2.5e6-2.5e6/((1+exp((V[cur_cell]-132.87)/25.4))*(1+exp(-(V[cur_cell]+24.92)/2.68)));
	IK2_k1[cur_cell] = IK2_k1[cur_cell] + dt * ((IK2_k_steady_state - IK2_k1[cur_cell]) / IK2_k1_Tau);
	IK2_k2[cur_cell] = IK2_k2[cur_cell] + dt * ((IK2_k_steady_state - IK2_k2[cur_cell]) / IK2_k2_Tau);

	return (gK2 * IK2_p[cur_cell] * (0.75 * IK2_k1[cur_cell] + 0.25 * IK2_k2[cur_cell]) * (V[cur_cell] - EK));
}


/* ##############################################
 * #                                            #
 * #  A-Type Transient Potassium Current - IKA  #
   #                                            #
   ############################################## */

double IKA ()
{
	double gKA(0.189/* *0.4*/);//reduced 60% for estradiol
	double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
	
	// Activation
	IKA_s_steady_state = 1/(1+exp(-(V[cur_cell]+27.79)/7.57));
	IKA_s_Tau = 17/(1+((V[cur_cell]+20.52)/35)*((V[cur_cell]+20.52)/35));
	IKA_s[cur_cell] = IKA_s[cur_cell] + dt * ((IKA_s_steady_state - IKA_s[cur_cell]) / IKA_s_Tau);

	// Inactivation
	IKA_x_steady_state = 0.3293+0.71165/(1+exp((V[cur_cell]+53.919)/7.8111));
	IKA_x_Tau = 7.5+10/(1+pow((V[cur_cell]+34.18)/120,2));
	IKA_x[cur_cell] = IKA_x[cur_cell] + dt * ((IKA_x_steady_state - IKA_x[cur_cell]) / IKA_x_Tau);

	return (gKA * pow(IKA_s[cur_cell],1) * exp(pow(IKA_x[cur_cell],1)) * (V[cur_cell] - EK));
}


/* ################################################
 * #                                              #
 * #  Calcium-Activated Potassium Current - IKCa  #
 * #                                              #
 * ################################################ */

double IKCa ()
{
	double PBKCa = 3.9e-13;
	double NBKCa = 5.94e6;
	double iKCa;
	
	// Activation
	V_KCa_half[cur_cell]=-41.7 * log10(Cai[cur_cell]) - 128.2;
	IKCa_p_steady_state = 1 / (1 + exp(-(V[cur_cell] - V_KCa_half[cur_cell])/18.25));
	IKCa_pf_Tau = 0.84;
	IKCa_ps_Tau = 35.9;
	
	// Need an approximation to avoid /zero at higher voltages. Based on taylor expansion of exponential
	if(V[cur_cell]>160)
	{
		iKCa = 1e4 * PBKCa * V[cur_cell] * (Faraday * Faraday / (R * Temperature)) *
			((Ko - Ki[cur_cell] * (1+(V[cur_cell] * Faraday / (R * Temperature))+pow(V[cur_cell] * Faraday / (R * Temperature),2)/2+pow(V[cur_cell] * Faraday / (R * Temperature),3)/6)) /
			(1-(1+(V[cur_cell] * Faraday / (R * Temperature))+pow(V[cur_cell] * Faraday / (R * Temperature),2)/2+pow(V[cur_cell] * Faraday / (R * Temperature),3)/6)));
	}
	else
	{
		iKCa = 1e4 * PBKCa * V[cur_cell] * (Faraday * Faraday / (R * Temperature)) * ((Ko - Ki[cur_cell] * exp(V[cur_cell] * Faraday / (R * Temperature))) /
			          (1 - exp(V[cur_cell] * Faraday / (R * Temperature))));
	}
	
	IKCa_pf[cur_cell] = IKCa_pf[cur_cell] + dt * ((IKCa_p_steady_state - IKCa_pf[cur_cell]) / IKCa_pf_Tau);
	IKCa_ps[cur_cell] = IKCa_ps[cur_cell] + dt * ((IKCa_p_steady_state - IKCa_ps[cur_cell]) / IKCa_ps_Tau);
	double PKCa = 0.17 * IKCa_pf[cur_cell] + 0.83 * IKCa_ps[cur_cell];
	
	return (Cell_Area * NBKCa * PKCa * iKCa / 100);
}


/* ###################################################
 * #                                                 #
 * #  Background Potassium Leakage Current - IKleak  #
 * #                                                 #
 * ###################################################*/ 

double IKleak ()
{
	double gKleak(0.005);
	double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
	return (gKleak * (V[cur_cell] - EK));
}


/* ##########################################
 * #                                        #
 * #  Non-Selective Cation Current - Inscc  #
 * #                                        #
 * ########################################## */

double Inscc ()
{
	const double gns(0.0123), g_leak(0.009685);
	const double KdMg(0.28);
	double PNa_PCs(0.9), PK_PCs(1.3), PCa_PCs(0.89); // Permeability ratios for Inscc

	double Ens = (R * Temperature / Faraday) *
	      log((PNa_PCs * Nao + PK_PCs * Ko + 4 * PCa_PCs * Cao * (Cao / (1 + exp((V[cur_cell] * Faraday) / (R * Temperature))))) / //should there be a permeability of calcium
	          (PNa_PCs * Nai[cur_cell] + PK_PCs * Ki[cur_cell] + 4 * PCa_PCs * Cai[cur_cell] * (Cai[cur_cell] / (1 + exp((V[cur_cell] * Faraday) / (R * Temperature)))))); //factor in here near PCa_PCs? if so, what is it?

	Inscc_fMg_steady_state = 0.1 + 0.9 / ((1 + pow((Mgo / KdMg),1.3)));

	double gnscc_Ca = (0.03 / (1 + pow((150 / (Cao + 1e-8)),2)))/0.000525;
	double gnscc_Na = (0.03 / (1 + pow((150 / (Nao + 1e-8)),2)))/0.0123;
	double gnscc_K = (0.03 / (1 + pow((150 / (Ko + 1e-8)),2)))/0.0123;

	InsCa = InsCa_mag*0.15*Inscc_fMg_steady_state * gnscc_Ca * gns * (V[cur_cell]-Ens);
	InsNa = InsNa_mag*Inscc_fMg_steady_state * gnscc_Na * gns * (V[cur_cell]-Ens);
	InsK = InsK_mag*1.19*Inscc_fMg_steady_state * gnscc_K * gns * (V[cur_cell]-Ens);
	Insleak = Insleak_mag*g_leak * (V[cur_cell] - Ens);

	return (InsCa + InsNa + InsK + Insleak);
}


/* ##############################################
 * #                                            #
 * #  Calcium-Activated Chloride Current - ICl  #
 * #                                            #
 * ############################################## */

double ICl ()
{
	double gCl(0.1875);
	double ECl_pos = -30, ECl_neg = -5;
	double KdCl(256.98e-3);
	double gCl_value;
	
	// Activation
	ICl_c_steady_state = 1/(1+exp(-(V[cur_cell]-44.08)/13.7));
	ICl_c1_Tau = 24.-16.88/(1+exp(-(V[cur_cell]-48.14)/13.));
	ICl_c2_Tau = 32.+81.54/(1+exp(-(V[cur_cell]+43.51)/39.5));
	ICl_c1[cur_cell] = ICl_c1[cur_cell] + dt * ((ICl_c_steady_state - ICl_c1[cur_cell]) / ICl_c1_Tau);
	ICl_c2[cur_cell] = ICl_c2[cur_cell] + dt * ((ICl_c_steady_state - ICl_c2[cur_cell]) / ICl_c2_Tau);

	// Inactivation
	ICl_n_steady_state = 1 / (1 + exp((V[cur_cell] + 56.97) / 20.29));
	ICl_n_Tau = 0.69;
	ICl_n[cur_cell] = ICl_n[cur_cell] + dt * ((ICl_n_steady_state - ICl_n[cur_cell]) / ICl_n_Tau);

	if (V[cur_cell] <= -20.0)
	{	gCl_value = gCl * (19.26 + 365.8 * exp(0.0564 * V[cur_cell]));}
	else 
	{	gCl_value = gCl * (0.752 + 0.0349 * exp(-0.414 * V[cur_cell]));}

	if (V[cur_cell] > 0)
	{	
		return gCl_value * (Cai[cur_cell] / (Cai[cur_cell] + KdCl)) * pow(ICl_c1[cur_cell],0.5) * (V[cur_cell] - ECl_pos);
	}
	else
	{	
		return gCl_value * (Cai[cur_cell] / (Cai[cur_cell] + KdCl)) * pow(ICl_c2[cur_cell],1) * pow(ICl_n[cur_cell],1) * (V[cur_cell] - ECl_neg);
	}

}


/* ####################################
 * #                                  #
 * #  Calcium-Pump Current - ICapump  #
 * #                                  #
 * #################################### */

double ICapump ()
{
	const double JpCa(7e-7); //Max flux in mM/ms
	const double KpCa(0.0005); //Half saturation of cai in mM
	const double hpCa(2); //Hill coefficient
	return ((zCa*Faraday*Cell_Volume)/(Membrane_Capacitance*Cell_Area))*(JpCa/(1+pow((KpCa/Cai[cur_cell]),hpCa)));
}


/* ###############################################
 * #                                             #
 * #  Sodium/Potassium Exchanger Current - INaK  #
 * #                                             #
 * ############################################### */

double INaK ()
{
	const double gNaK(1.7);
	const double KdK(2), KdNa(22);
	const double nK(1.5), nNa(2); //Hill coefficients
	double kNaK = 1 / (1 + pow((KdK / Ko),nK));
	double nNaK = 1 / (1 + pow((KdNa / Nai[cur_cell]),nNa));
	double fNaK = 1 / (1 + 0.125 * exp(-(0.1 * V[cur_cell] * Faraday) / (R * Temperature)) +
		      0.00219 * exp(Nao / 49.71) * exp(-(1.9 * V[cur_cell] * Faraday) / (R * Temperature)));
	
	return (gNaK * fNaK * kNaK * nNaK);
}


/* ##############################################
 * #                                            #
 * #  Sodium/Calcium Exchanger Current - INaCa  #
 * #                                            #
 * ############################################## */

double INaCa ()
{
	double JNaCa(3.5e-6); // Maximum flux in mM/ms
	const double KmAllo(0.0003), KmNai(30), KmCai(0.007), KmNao(87.5), KmCao(1.3);
	const double nallo(4);
	const double NaCa_energy_barrier(0.35); // Gamma, energy barrier of Na/Ca exchanger in the electric field
	const double ksat(0.27);  // Saturation if INaCa at negative potentials
	
	double VF_RT = (V[cur_cell]*Faraday)/(R*Temperature);
	double barrier = exp(NaCa_energy_barrier * VF_RT);
	double barrier_1 = exp((NaCa_energy_barrier-1) * VF_RT);

	double Allo = 1 / (1 + pow((KmAllo / Cai[cur_cell]),nallo));

	double E1 = JNaCa * pow(Nai[cur_cell],3) * Cao * barrier - JNaCa * pow(Nao,3) * Cai[cur_cell] * barrier_1;
    double E2 = 1 + ksat * barrier_1;
	double E3 = KmCao * pow(Nai[cur_cell],3) + pow(KmNao,3) * Cai[cur_cell] + pow(KmNai,3) * Cao * (1 + (Cai[cur_cell] / KmCai));
	double E4 = Cao * pow(Nai[cur_cell],3) + pow(Nao,3) * Cai[cur_cell] + pow(Nao,3) * KmCai * (1 + (pow(Nai[cur_cell],3)/KmNai));
	double delta_E = E1 / (E2 * (E3 + E4));

	double INaCa_JNaCa = Allo * delta_E;

	return ((INaCa_JNaCa * zCa * Cell_Volume * Faraday) / (Cell_Area * Membrane_Capacitance));
}


/* #############################################################
 * #                                                           #
 * #  Sodium-Potassium-Chloride Co-Transport Current - INaKCl  #
 * #                                                           #
 * ############################################################# */

double INaKCl ()
{
	double LNaKCl(1.79e-8);
	double RNaKCl(1); // Co-Transport regulation

	return (-RNaKCl * zCl * Cell_Area * LNaKCl * R * Faraday * Temperature *
			log(((Nao * Ko) / (Nai[cur_cell] * Ki[cur_cell])) *pow( (Clo / Cli[cur_cell]),2)));
}

/* ##############################
   #                            #
   # Kinetic Cross-Bridge Cycle #
   #                            #
   ############################## */

void KCBC ()
{
	double Wortmannin = 1; //1=off, 0=on
	const double MLCK_max = 0.84;
	//KCa_MLCK is a key component here due to it being relative to the Cai and then raised to a power
	//Reducing it causes KCBC to be more sensistive to calcium levels
	const double K_MLCK = 721.746e-6,KCa_MLCK = 6*180e-6;
	const double pM=1, nM=8.7613;
	double MLCK = Wortmannin * MLCK_max/(1+pow((Cai[cur_cell]/K_MLCK),pM));
	double K1 = MLCK/(1+pow((KCa_MLCK/Cai[cur_cell]),nM));
	double K6 = K1;
	const double K2 = 0.4, K3 = 1.8, K4 = 0.1, K5 = 0.4, K7 = 0.045; //Yang 2003
	
	M = M+dt*(-K1*M+K2*Mp+K7*AM);
	Mp = Mp+dt*(K4*AMp+K1*M-(K2+K3)*Mp);
	AMp = AMp+dt*(K3*Mp+K6*AM-(K4+K5)*AMp);
	AM = AM+dt*(K5*AMp-(K7+K6)*AM);
	
	Frac_phos = AMp + Mp;
	Frac_stress = AMp + AM;
}

/* ####################
 * #                  #
 * # Force Production #
 * #                  #
 * #################### */

void force ()
{
	lc=120; // cell length 
	double kp=0.1, kx1=12.5, kx2=8.8, ks=0.2, mu_s=0.01, lc0=40, ls0=30;
	double lopt=100, fAMp=1.3, fAM=85.5, vx=5, beta=7.5, alpha_p=0.1, alpha_s=4.5;
	
	KCBC();
	
	Fp = kp*(exp(alpha_p*(lc-lc0)/lc0)-1); 
	Fx = (kx1*AMp+kx2*AM)*lx*exp(-beta*(pow(((la-lopt)/lopt),2)));
	Fa = (fAMp*AMp*(vx+la)+fAM*AM*la)*exp(-beta*(pow(((-la-lopt)/lopt),2))); 
	Fs = Fa; 
	Ft = Fa+Fp;

	// Following lengths are initialised in main program
	ls = ls+dt*((1/mu_s)*(Fs-ks*(exp(alpha_s*(ls-ls0)/ls0)-1))); 
	la = la+dt*((Fa/exp(-beta*(pow(((la-lopt)/lopt),2)))-(fAMp*AMp*vx))/(fAM*AM+fAMp*AMp));  
	lx=(lc-la-ls);
}


/* ######################
 * #                    #
 * # Ryanodine receptor #
 * #                    #
 * ###################### */

void RR ()
{
	const double Kr1=2500,Kr2=1.05,K_r1=0.0076,K_r2=0.084; //according to Yang 2003
	
	R00[cur_cell]+=dt*(-(Kr1*Cai[cur_cell]+Kr2)*Cai[cur_cell])*R00[cur_cell];
   	R10[cur_cell]+=dt*(Kr1*Cai[cur_cell]*Cai[cur_cell]*R00[cur_cell]-(K_r1+Kr2*Cai[cur_cell])*R10[cur_cell]+K_r2*R11[cur_cell]);
   	R01[cur_cell]+=dt*(Kr2*Cai[cur_cell]*R00[cur_cell]+K_r1*R11[cur_cell]-(K_r2+Kr1*Cai[cur_cell]*Cai[cur_cell])*R01[cur_cell]);
   	R11[cur_cell]+=dt*(Kr2*Cai[cur_cell]*R10[cur_cell]-(K_r1+K_r2)*R11[cur_cell]+Kr1*Cai[cur_cell]*Cai[cur_cell]*R01[cur_cell]);
	// Initialised above
}


/* ##################################################
 * #                                                #
 * # Sarcoplasmic Reticulum control and SR currents #
 * #                                                #
 * ################################################## */

void ISR ()
{
	double Iup_max=6.68*0.1, Kmup=1e-3*10, vol_u=1.33e-9, vol_r=1.33e-11, tau_tr=180,Frel=0.2,Rleak=75.07e-4,tau_rel=0.015,K_CSQN=0.8,CSQN=15,conc_red=1e6; 
	RR();

	Iup = Iup_max * (Cai[cur_cell]/(Cai[cur_cell]+Kmup));
	Itr = ((Ca_u[cur_cell]-Ca_r[cur_cell])*zCa*Faraday*vol_u*conc_red)/tau_tr; 
	Irel = Frel*((R10[cur_cell]*R10[cur_cell]+Rleak)*(Ca_r[cur_cell]-Cai[cur_cell])*zCa*Faraday*vol_r*conc_red)/tau_rel;

	Ca_u[cur_cell] += dt*((Iup-Itr)/(zCa*Faraday*vol_u*conc_red));
	Ca_r[cur_cell] += dt*((Itr-Irel)/(zCa*Faraday*vol_r*conc_red))*pow(1+(CSQN*K_CSQN)/pow(K_CSQN+Ca_r[cur_cell],2),-1);
}
	
	
/* ##################
 * #                #
 * # Calcium buffer # 
 * #                #
 * ################## */


double CaBuffer ()
{
	double FBCa = 2.0e1, Scm = 0.3, Bf = 0.2, Kd = 2.6e-4, Kdb = 5.39810e-4;
	
	double frac1 = (Scm*Kd)/pow(Kd+Cai[cur_cell],2);
	double frac2 = (Bf*Kdb)/pow(Kdb+Cai[cur_cell],2);

	return FBCa*pow(1+frac1+frac2,-1);
}

/* ########################
 * #                      #
 * # Total current - Itot # 
 * #                      #
 * ######################## */

double Itotal(double nif)
{
	return Ih_mag*Ih()+ICaL_mag*ICaL(nif)+ICaT_mag*ICaT()+INa_mag*INa()
	+IK1_mag*IK1()+IK2_mag*IK2()+IKA_mag*IKA()+IKCa_mag*IKCa()
	+Inscc_mag*Inscc()+ICl_mag*ICl()+ICapump_mag*ICapump()+INaK_mag*INaK()
	+INaCa_mag*INaCa()+Isus+Isoc_mag*Isoc()+IKleak_mag*IKleak()
	+INaKCl_mag*INaKCl();
}

/* #############################################
 * #                                           #
 * # Intracellular calcium concentration - Cai # 
 * #                                           #
 * ############################################# */

double ConCai(double nif)
{
	ISR();
	
	ICa_tot = ICaL_mag*ICaL(nif) + ICaT_mag*ICaT() - 2*INaCa_mag*INaCa()
	+ InsCa_mag*InsCa + IsocCa_mag*IsocCa 
	+ ICapump_mag*ICapump() + Iup_mag*Iup - Irel_mag*Irel;
	
	Cai[cur_cell] +=dt * (-(((Cell_Area * Membrane_Capacitance)
	/(zCa*Cell_Volume*Faraday))
	*(ICa_tot)*CaBuffer()));
	
	return Cai[cur_cell];
}

/* ########################################
 * #                                      #
 * # Intracellular potassium current - Ki # 
 * #                                      #
 * ######################################## */

void IK()
{
	ISR();
	
	IK_tot = IK1_mag*IK1() + IK2_mag*IK2() + IKA_mag*IKA()
	+ IKCa_mag*IKCa() - 2*InsK_mag*InsK 
	+ IKleak_mag*IKleak() + INaKCl_mag*INaKCl() - 2*INaK_mag*INaK();
}

}; // end of cell class


int main () 
{	
	// Generate output files
	ofstream output_file1, output_file2, output_file3;
	// Open the files
	output_file1.open ("Output_AP.dat");
	output_file2.open ("Output_sstau.dat");
	output_file3.open ("Output_force.dat");
	ofstream IC;
	IC.open("IC.txt");

	//Column titles in output files
	// Currents		
	output_file1 << 				
	"time" <<"\t"<< "v" <<"\t"<< "Cai" <<"\t"<< "Itot" <<"\t"<< "ICa_tot" <<"\t"<< "IK_tot" <<"\t"<<
	"Ih" <<"\t"<< "ICaL" <<"\t"<< "ICaT" <<"\t"<< "INa" <<"\t"<< "IK1" <<"\t"<< "IK2" <<"\t"<< "IKA" <<"\t"<< 
	"IKCa" <<"\t"<< "Inscc" <<"\t"<< "ICl" <<"\t"<< "ICapump" <<"\t"<< "INaK" <<"\t"<< "INaCa"  <<"\t"<< 
	"Isus" <<"\t"<< "Isoc" <<"\t"<< "INaKCl" <<"\t"<< "IKleak" << "\t" << 
	"InsCa" << "\t" << "InsNa" << "\t" << "InsK" <<"\t"<< "Insleak" << "\t"
	"Iup" <<"\t"<< "Irel" <<"\t"<< "Ca_u" <<"\t"<< "Ca_r" <<"\t"<<"Itr"<< "\t"<<
	"cur_cell" <<"\t"<< "Istim" <<"\t"<< "Ft" <<"\n";
	
	/*SS and time const*/   
	output_file2 <<
	"time" << "\t" << "v" << "\t" << "cai" << "\t" << "Ih_y_steady_state" << "\t" << "Ih_y_Tau" << "\t" <<
	"ICaL_d_steady_state"  << "\t" << "ICaL_d_Tau" << "\t" << 
	"ICaL_f1_steady_state" << "\t" << "ICaL_f1_Tau" << "\t" << "ICaL_f2_Tau" << "\t" <<
	"ICaT_b_steady_state" << "\t" << "ICaT_b_Tau" << "\t" << "ICaT_g_steady_state" << "\t" << "ICaT_g_Tau" << "\t" <<
	"INa_m_steady_state" << "\t" << "INa_m_Tau" << "\t" << "INa_h_steady_state" << "\t" << "INa_h_Tau" << "\t" << 
	"IK1_q_steady_state" << "\t" << "IK1_q_Tau" << "\t" << "IK1_r_steady_state" << "\t" << "IK1_r1_Tau" << "\t" << "IK1_r2_Tau" << "\t" << 
	"IK2_p_steady_state" << "\t" << "IK2_p_Tau" << "\t" << "IK2_k_steady_state" << "\t" << "IK2_k1_Tau" << "\t" << "IK2_k2_Tau" << "\t" <<
	"IKA_s_steady_state" << "\t" << "IKA_s_Tau" << "\t" << "IKA_x_steady_state" << "\t" << "IKA_x_Tau" << "\t" <<
	"IKCa_p_steady_state" << "\t" << "IKCa_pf_Tau" << "\t" << "IKCa_ps_Tau" << "\t" <<
	"Inscc_fMg_steady_state" << "\t" <<  
	"ICl_c_steady_state" << "\t" << "ICl_c1_Tau" << "\t" << "ICl_c2_Tau" << "\t" << "ICl_n_steady_state" << "\t" << "ICl_n_Tau" << "\n";
	
	/*Force*/		
	output_file3 << 
	"time" << "\t" << "v" << "\t" << "Cai" << "\t" <<
	"M" << "\t" << "AM" << "\t" << "Mp" << "\t" << "AMp" << "\t" <<
	"Frac_phos" << "\t" << "Frac_stress" << "\t" <<
	"Ft" << "\t" << "Fs" << "\t" << "Fa" << "\t" << "Fx" << "\t" << "Fp" << "\t" <<
	"Istim" << "\t" << "ICa_tot" << "\t" <<
	"lc" << "\t" << "ls" << "\t" << "la" << "\t" << "lx" <<  "\n";	
	
	// Generate cells
	cell myocyte[num_cells];
	// Intialise time counter for outputs
	int counter = 0;
	
	// Set nifedipine conc
	double nif = 1e-12;

	// Set Stim condutions
//	int pulse_time=5e3;
//	int pulse_duration=100;
//	int pulse_counter=0;
	
	//Time loop (remember time is in ms, dt is 0.01 ms. So after 100 time loops Time == 1. E.g. Time == 1e3 == 1 second)
	for (double Time=0; Time<=25e3; (Time += dt) && (counter++))
	{	
		//if (counter==10e5) ICl_mag=ICl_mag*2;

		// Check voltage and set INaK_mag
		if (V[cur_cell]>0) {INaK_mag = 2;}
		else {INaK_mag = 0;}
		
		
		//----------------Stimuli---------------------

		//Default stim to zero (or reset for loops)
		Istim=0;
		//Stimulus current (simple single pulse)
		if ((Time > 5e3 && Time <= 15e3)) 
		{
			//if (Time>9850 && Time<=9860) Istim = -0.400; //for short pulses during AP change Istim to higher value
			//else Istim = -0.400;
			Istim = -0.400;
		}//0.05 //Iclamp;
		else {Istim = 0.;}

		//Repetative stimulus
/*		if (Time > 5e3 && pulse_counter<10 && (counter*dt)>pulse_time && (counter*dt)<=pulse_time+pulse_duration)
		{
			Istim = -0.400;
		}//0.05 //Iclamp;
		else {Istim = 0.;}
		//if (counter%10000==0) cout << counter*dt << "\t" << pulse_time << "\t" << pulse_duration << endl;
		if ((counter*dt)==pulse_time+pulse_duration)
		{
			pulse_time+=400;
			pulse_counter++;
		}
*/
		//--------------------------------------------
		
		//Calculate intracellular calcium
		myocyte[cur_cell].ConCai(nif);
		//Calculate net current through membrane
		Itot = myocyte[cur_cell].Itotal(nif);
		//Calculate forces
		myocyte[cur_cell].force();
		//Calculate total K current
		myocyte[cur_cell].IK();

		//------------Voltage protocols---------------
		
		//Voltage protocol (forward Euler)
		V[cur_cell] += (-(Itot + Istim) * dt)/Membrane_Capacitance;
		
		//Voltage clamp pulses
	/*	if (Time > 5e3 && pulse_counter<10 && (counter*dt)>pulse_time && (counter*dt)<=pulse_time+pulse_duration)
		{
			V[cur_cell] = -30;
		}
		else V[cur_cell] = -60.;
		if ((counter*dt)==pulse_time+pulse_duration)
		{
			pulse_time+=400;
			pulse_counter++;
		}
	*/
		//--------------------------------------------
		

		//Output to terminal
		if (counter%10000==0) // 100 for every ms, 1000 for every 10 ms
		{
			cout << Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << endl;
		}
		
		//Output data to files
		if (counter%100==0 /*&& Time>11e3*/) // 100 for every ms, 1000 for every 10 ms
		{
			//Currents
			output_file1 <<
			Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << "\t" << 
			Itot << "\t" <<
			ICa_tot << "\t" <<
			IK_tot << "\t" <<
			Ih_mag*myocyte[cur_cell].Ih() << "\t" <<  
			ICaL_mag*myocyte[cur_cell].ICaL(nif) << "\t" <<  
			ICaT_mag*myocyte[cur_cell].ICaT() << "\t" <<  
			INa_mag*myocyte[cur_cell].INa() << "\t" <<  
			IK1_mag*myocyte[cur_cell].IK1() << "\t" << 
			IK2_mag*myocyte[cur_cell].IK2() << "\t" << 
			IKA_mag*myocyte[cur_cell].IKA() << "\t" << 
			IKCa_mag*myocyte[cur_cell].IKCa() << "\t" << 
			Inscc_mag*myocyte[cur_cell].Inscc() << "\t" << 
			ICl_mag*myocyte[cur_cell].ICl() << "\t" << 
			ICapump_mag*myocyte[cur_cell].ICapump() << "\t" <<
			INaK_mag*myocyte[cur_cell].INaK() << "\t" <<
			INaCa_mag*myocyte[cur_cell].INaCa() << "\t" << 
			Isus << "\t" <<
			Isoc_mag*myocyte[cur_cell].Isoc() << "\t" <<
			INaKCl_mag*myocyte[cur_cell].INaKCl() << "\t" <<
			IKleak_mag*myocyte[cur_cell].IKleak() << "\t" <<
			InsCa_mag*InsCa << "\t" <<
			InsNa_mag*InsNa << "\t" <<
			InsK_mag*InsK << "\t" <<
			Insleak_mag*Insleak << "\t" <<
			Iup_mag*Iup << "\t" <<
			Irel_mag*Irel << "\t" <<
			Ca_u[cur_cell] << "\t" <<
			Ca_r[cur_cell] << "\t" <<
			Itr_mag*Itr << "\t" <<
			cur_cell << "\t" << 
			Istim << "\t" <<
			Ft << "\n";
		
			/*SS and time const*/   	
			output_file2 << Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << "\t" << Ih_y_steady_state << "\t" << Ih_y_Tau << "\t" <<
			ICaL_d_steady_state  << "\t" << ICaL_d_Tau << "\t" << 
			ICaL_f_steady_state << "\t" << ICaL_f1_Tau << "\t" << ICaL_f2_Tau << "\t" << 
			ICaT_b_steady_state << "\t" << ICaT_b_Tau << "\t" << ICaT_g_steady_state << "\t" << ICaT_g_Tau << "\t" <<
			INa_m_steady_state << "\t" << INa_m_Tau << "\t" << INa_h_steady_state << "\t" << INa_h_Tau << "\t" << 
			IK1_q_steady_state << "\t" << IK1_q_Tau << "\t" << IK1_r_steady_state << "\t" << IK1_r1_Tau << "\t" << IK1_r2_Tau << "\t" << 
			IK2_p_steady_state << "\t" << IK2_p_Tau << "\t" << IK2_k_steady_state << "\t" << IK2_k1_Tau << "\t" << IK2_k2_Tau << "\t" <<
			IKA_s_steady_state << "\t" << IKA_s_Tau << "\t" << IKA_x_steady_state << "\t" << IKA_x_Tau << "\t" <<
			IKCa_p_steady_state << "\t" << IKCa_pf_Tau << "\t" << IKCa_ps_Tau << "\t" <<
			Inscc_fMg_steady_state << "\t" <<  
			ICl_c_steady_state << "\t" << ICl_c1_Tau << "\t" << ICl_c2_Tau << "\t" << ICl_n_steady_state << "\t" << ICl_n_Tau << "\n";
						
			/*Force*/			
			output_file3 << Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << "\t" <<
			M << "\t" << AM << "\t" << Mp << "\t" << AMp << "\t" <<
			Frac_phos << "\t" << Frac_stress << "\t" <<
			Ft << "\t" << Fs << "\t" << Fa << "\t" << Fx << "\t" << Fp<< "\t" <<
			Istim << "\t" << ICa_tot << "\t" <<
			lc << "\t" << ls << "\t" << la << "\t" << lx <<  "\n";
		}

		if (counter>1100e4 && counter<=1101e4)//Time>(Time-1e2))
		{
			IC << "Time" << "\t" << "Ih_y[0]" << "\t" << "ICaL_d[0]" << "\t" << "ICaL_f1[0]" << "\t" << "ICaL_f2[0]" << "\t" <<"ICaT_b[0]" 
			<< "\t" << "ICaT_g[0]" << "\t" << "INa_m[0]" << "\t" << "INa_h[0]" << "\t" << "IK1_q[0]" 
			<< "\t" << "IK1_r1[0]" << "\t" << "IK1_r2[0]" << "\t" << "IK2_p[0]" << "\t" << "IK2_k1[0]" << "\t" << "IK2_k2[0]"
			<< "\t" << "IKA_s[0]" << "\t" << "IKA_x[0]"
			<< "\t" << "ICl_c1[0]" << "\t" << "ICl_c2[0]" << "\t" << "ICl_n[0]" << "\t" << "V[0]" 
			<< "\t" << "Cai[0]" << "\t" << "Ca_u[0]" << "\t" << "Ca_r[0]" << "\t" << "R00[0]" 
			<< "\t" << "R01[0]" << "\t" << "R10[0]" << "\t" << "R11[0]" << "\t" << "V_KCa_half" 
			<< "\t" << "IKCa_pf[0]" << "\t" << "IKCa_ps[0]" << endl;
			
			IC << Time << "\t" << Ih_y[cur_cell] << "\t" << ICaL_d[cur_cell] << "\t" << ICaL_f1[cur_cell] << "\t" << ICaL_f2[cur_cell] << "\t" <<ICaT_b[cur_cell] 
			<< "\t" << ICaT_g[cur_cell] << "\t" << INa_m[cur_cell] << "\t" << INa_h[cur_cell] << "\t" << IK1_q[cur_cell] 
			<< "\t" << IK1_r1[cur_cell] << "\t" << IK1_r2[cur_cell] << "\t" << IK2_p[cur_cell] << "\t" << IK2_k1[cur_cell] << "\t" << IK2_k2[cur_cell]
			<< "\t" << IKA_s[cur_cell] << "\t" << IKA_x[cur_cell]
			<< "\t" << ICl_c1[cur_cell] << "\t" << ICl_c2[cur_cell] << "\t" << ICl_n[cur_cell] << "\t" << V[cur_cell] 
			<< "\t" << Cai[cur_cell] << "\t" << Ca_u[cur_cell] << "\t" << Ca_r[cur_cell] << "\t" << R00[cur_cell] 
			<< "\t" << R01[cur_cell] << "\t" << R10[cur_cell] << "\t" << R11[cur_cell] << "\t" << V_KCa_half[cur_cell] 
			<< "\t" << IKCa_pf[cur_cell] << "\t" << IKCa_ps[cur_cell] << endl;
		}
	} //end of time loop
	
	
	
	
	

  //cell A;
  
} // end of main

/*Notes:
 * Elements cell class must contain:
 * channel functions to work out current through membrane
 * membrane voltage
 * intracellular ion handling
 *
 * Elements outside of class:
 * diffusion between cells
 * number of cells etc
 *
 * */