//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
{

  //V
  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/y[3]))   //Fast Sodium 
	 +gNaP*y[4]*(y[0]-rtf*log(Nae/y[3]))   //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/y[3]))      //Sodium leak
	 +gL*(EL-ENa)/(EK-ENa)*(y[0]-EK)      //Potassium leak
	 +Ipmpmx/(1.0+pow(Naih*exp(alpha*y[0]/rtf)/y[3],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));

  //[Na+]i
  f[3] = (-1.0/(vlm*frdy))*(gL*(EL-EK)/(ENa-EK)*(y[0]-rtf*log(Nae/y[3]))   //Sodium Leak
			    +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/y[3])) //INaF
			    +gNaP*y[4]*(y[0]-rtf*log(Nae/y[3]))   //INaP
			    +3.0*Ipmpmx/(1.0+pow(Naih*exp(alpha*y[0]/rtf)/y[3],nhillNa)));   //3*IPump  

  //mNaP
  f[4] = (1.0/(1.0+exp((y[0]-VmNaP)/kmNaP))-y[4])/taumNaP; 

  return GSL_SUCCESS;
}