//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((y[0]-VmNaF)/kmNaF))*1.0/(1.0+exp((y[0]-VmNaF)/kmNaF))*1.0/(1.0+exp((y[0]-VmNaF)/kmNaF))*y[2]*(y[0]-rtf*log(Nae/Nai))   //Fast Sodium 
	 +gNaP*y[3]*(y[0]-rtf*log(Nae/Nai))   //Persistent sodium
	 +gK*y[1]*y[1]*y[1]*y[1]*(y[0]-EK)	 //Potassium 
	 +gL*(EL-EK)/(ENa-EK)*(y[0]-rtf*log(Nae/Nai))      //Sodium leak
	 +gL*(EL-ENa)/(EK-ENa)*(y[0]-EK)      //Potassium leak
	 +Ipmpmx/(1.0+pow(Naih*exp(alpha*y[0]/rtf)/Nai,nhillNa))    //Pump Current 
	 +gModWE*(y[0]-0.0)+gModWI*(y[0]+75.0)-Iinj)/capac;	 

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

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

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

  return GSL_SUCCESS;
}