/*--------------------------------------------------------------------------
   Author: Thomas Nowotny
  
   Institute: Institute for Nonlinear Dynamics
              University of California San Diego
              La Jolla, CA 92093-0402
  
   email to:  tnowotny@ucsd.edu
  
   initial version: 2005-08-19
  
--------------------------------------------------------------------------*/

#ifndef CN_LPANEURON_CC
#define CN_LPANEURON_CC

#include "CN_neuron.cc"

LPAneuron::LPAneuron(int inlabel, double *the_p= LPA_p):
  neuron(inlabel, LPA_IVARNO, LPANEURON, the_p, LPA_PNO)
{
}

LPAneuron::LPAneuron(int inlabel, tnvector<int> inpos, double *the_p= LPA_p):
  neuron(inlabel, LPA_IVARNO, LPANEURON, inpos, the_p, LPA_PNO)
{
}

inline double LPAneuron::E(double *x)
{
  return x[idx];
}

#define efunc(X,Y,Z) (1.0/(1.0+expl(((X)+(Y))/(Z))))
#define E x[idx]

#define Capacit p[0]
#define tau_Ca p[1]
#define f p[2]
#define C_Ca_0 p[3]
#define E_Na p[4]
#define E_K p[5]
#define E_H p[6]
#define E_leak p[7]
#define g_H p[8]
#define g_leak p[9]
#define g_Na p[10]
#define g_CaT p[11]
#define g_CaS p[12]
#define g_A p[13]
#define g_KCa p[14]
#define g_Kd p[15]
#define Area p[16]

void LPAneuron::currents(ostream &os, double *x)
{
  _E_Ca=(12.5*log(3000.0/x[idx+12]));
  
  I_Na=  g_Na*ipower(x[idx+1],3)*x[idx+2]*(E-E_Na)*Area;
  I_CaT= g_CaT*ipower(x[idx+3],3)*x[idx+4]*(E-_E_Ca)*Area;
  I_CaS= g_CaS*ipower(x[idx+5],3)*x[idx+6]*(E-_E_Ca)*Area;
  I_A=   g_A*ipower(x[idx+7],3)*x[idx+8]*(E-E_K)*Area;
  I_KCa= g_KCa*ipower(x[idx+9],4)*(E-E_K)*Area;
  I_Kd=  g_Kd*ipower(x[idx+10],4)*(E-E_K)*Area;
  I_H=   g_H*x[idx+11]*(E-E_H)*Area;
  I_leak= g_leak*(E-E_leak)*Area;
  
  os << I_Na << " ";
  os << I_CaT << " ";
  os << I_CaS << " ";
  os << I_A << " ";
  os << I_KCa << " ";
  os << I_Kd << " ";
  os << I_H << " ";
  os << I_leak << " ";
}

void LPAneuron::derivative(double *x, double *dx)
{
  Isyn= 0.0;
  forall (den_it) {
    Isyn+= den_it->c_value()->Isyn(x);
  }

  _E_Ca=(12.5*log(3000.0/x[idx+12]));
  
  I_Na=  g_Na*ipower(x[idx+1],3)*x[idx+2]*(E-E_Na);
  I_CaT= g_CaT*ipower(x[idx+3],3)*x[idx+4]*(E-_E_Ca);
  I_CaS= g_CaS*ipower(x[idx+5],3)*x[idx+6]*(E-_E_Ca);
  I_A=   g_A*ipower(x[idx+7],3)*x[idx+8]*(E-E_K);
  I_KCa= g_KCa*ipower(x[idx+9],4)*(E-E_K);
  I_Kd=  g_Kd*ipower(x[idx+10],4)*(E-E_K);
  I_H=   g_H*x[idx+11]*(E-E_H);
  I_leak= g_leak*(E-E_leak);
  
  dx[idx]= Isyn/Capacit -(I_Na + I_Kd + I_CaT + I_CaS + I_A + I_KCa +
			       I_H + I_leak)*Area/Capacit;
  
  // differential eqn for mNa
  _minf=efunc(E,25.5,-5.29);
  _taum=2.64-2.52*efunc(E,120.0,-25.0);
  dx[idx+1]=(_minf-x[idx+1])/_taum;
  // differential eqn for hNa
  _hinf= efunc(E,48.9,5.18);
  _tauh= (1.34*efunc(E,62.9,-10.0))*(1.5+efunc(E,34.9,3.6));
  if (_tauh < 1e-2) _tauh= 1e-2; // this is a horrible numerical instability
  // if the neuron is hyperpolarized, _tauh goes to 0 amplifying the already
  // clear disadvantageous rounding errors in (order of 1 - order of 1)
  // in mNa *a lot*. This leads to an explosion and eventually nan ...
  //  cerr << _hinf << " " << _tauh << " ";
  dx[idx+2]=(_hinf-x[idx+2])/_tauh;

  // differential eqn for mCa1
  _minf= efunc(E,27.1,-7.2);
  _taum= 43.4-42.6*efunc(E,68.1,-20.5);
  dx[idx+3]=(_minf-x[idx+3])/_taum;
  // differential eqn for hCa1
  _hinf= efunc(E,32.1,5.5);
  _tauh= 210.0-179.6*efunc(E,55.0,-16.9);
  dx[idx+4]=(_hinf-x[idx+4])/_tauh;

  // differential eqn for mCa2
  _minf= efunc(E,33.0,-8.1);
  _taum= 2.8+(14.0/(exp((E+27.0)/10.0)+exp((E+70.0)/-13.0)));
  dx[idx+5]=(_minf-x[idx+5])/_taum;
  // differential eqn for hCa2
  _hinf= efunc(E,60,6.2);
  _tauh= 120+(300.0/(exp((E+55.0)/9.0)+exp((E+65.0)/-16)));
  dx[idx+6]=(_hinf-x[idx+6])/_tauh;

  // differential eqn for mA
  _minf= efunc(E,27.2,-8.7);
  _taum= 23.2-20.8*efunc(E,32.9,-15.2);
  dx[idx+7]=(_minf-x[idx+7])/_taum;
  // differential eqn for hA
  _hinf= efunc(E,56.9,4.9);
  _tauh= 77.2-58.4*efunc(E,38.9,-26.5);
  dx[idx+8]=(_hinf-x[idx+8])/_tauh;

  // differential eqn for mKCa
  _minf= (x[idx+12]/(x[idx+12]+3.0))*efunc(E,28.3,-12.6);
  _taum= 180.6-150.2*efunc(E,46.0,-22.7);
  dx[idx+9]=(_minf-x[idx+9])/_taum;
  
  // differential eqn for mKd
  _minf= efunc(E,12.3,-11.8);
  _taum= 14.4-12.8*efunc(E,28.3,-19.2);
  dx[idx+10]=(_minf-x[idx+10])/_taum;
  
  // differential eqn for mh
  _minf= efunc(E,75.0,5.5);
  _taum= 2.0/(exp((E+169.7)/-11.6)+exp((E-26.7)/14.3));
  if (_taum < 1e-2) _taum= 1e-2; // same instability as for _tauh_Na ...
  dx[idx+11]=(_minf-x[idx+11])/_taum;

  // differential eqn for o, prob. of Ca mediated K channel activation 
  dx[idx+12]=(-f*(I_CaT+I_CaS)-x[idx+12]+C_Ca_0)/tau_Ca;

  // the Isyn is calculated elsewhere
  dx[idx+13]= 0.0; 
}

#undef efunc
#undef E
#undef Capacit
#undef tau_Ca
#undef f
#undef C_Ca_0
#undef E_Na
#undef E_K
#undef E_H
#undef E_leak
#undef g_H
#undef g_leak
#undef g_Na
#undef g_CaT
#undef g_CaS
#undef g_A
#undef g_KCa
#undef g_Kd
#undef Area

#endif