//RodentPreBotCNeuron_dy.c
//Jessica Parker, October 5, 2024
//
//This is the equation file. The set of differential equations for this model are defined here.

#include "RodentPreBotCNeuron.h"

int func(double t, const double y[], double f[],void *params) //params is the pointer to the array of parameter values
{

  //dV1
  /*  f[0]=-(gNaF*1.0/(1.0+exp((Vclamp-VmNaF)/kmNaF))*1.0/(1.0+exp((Vclamp-VmNaF)/kmNaF))*1.0/(1.0+exp((Vclamp-VmNaF)/kmNaF))*y[2]*(Vclamp-rtf*log(Nae/y[3]))   //Fast Sodium from Rybak et al. 2003 experimental
	 +gNaP*y[4]*(Vclamp-rtf*log(Nae/y[3]))   //Slowly-inactivating component of persistent sodium, mNaP and hNaP from two different sources
	 +gK*y[1]*y[1]*y[1]*y[1]*(Vclamp-EK)	 //Potassium from Barrett-Jolley et al. 2000 except maybe the time constant might be from Itri et al. 2005
	 +gL*(EL-EK)/(ENa-EK)*(Vclamp-rtf*log(Nae/y[3]))      //Sodium leak
	 +gL*(EL-ENa)/(EK-ENa)*(Vclamp-EK)      //Potassium leak
	 +gModWE*(Vclamp-0.0)+gModWI*(Vclamp+75.0)-Iinj)/capac;	 */

  //mK
  f[0]=(1.0/(1.0+exp((Vclamp-VmK)/kmK))-y[0])/taumK;

  //hNaF
  f[1]=(1.0/(1.0+exp((Vclamp-VhNaF)/khNaF))-y[1])/(tauhNaF/cosh((Vclamp-VhNaF)/12.8));

  //[Na+]i
  f[2] = (-1.0/(vlm*frdy))*(gL*(EL-EK)/(ENa-EK)*(Vclamp-rtf*log(Nae/y[2]))   //ILNa
			    +gNaF*1.0/(1.0+exp((Vclamp-VmNaF)/kmNaF))*1.0/(1.0+exp((Vclamp-VmNaF)/kmNaF))*1.0/(1.0+exp((Vclamp-VmNaF)/kmNaF))*y[1]*(Vclamp-rtf*log(Nae/y[2])) //INaF
			    +gNaP*y[3]*(Vclamp-rtf*log(Nae/y[2]))   //Persistent Sodium
			    +3.0*Ipmpmx/(1.0+pow(Naih*exp(alpha*Vclamp/rtf)/y[2],nhillNa)));   //3*IPump  

  //mNaP
  f[3] = (1.0/(1.0+exp((Vclamp-VmNaP)/kmNaP))-y[3])/taumNaP; //From Rybak et al. 2003 experimental

  return GSL_SUCCESS;
}