/*
  生体システム工学の基礎
  図4.12
  gcc -o hh hh.c -lm
 */
#include <stdio.h>
#include <math.h>

#define GNa 120.0
#define GK 36.0
#define Dt 0.025
#define ENa 35.0
#define EK -92.0
#define ECl -84.06686
#define Cm 1.0

double am(double v)
{
  double alpha;
  alpha = 0.1 * (25.0 - v)/(exp((25.0 - v)/10.0) - 1.0);
  return alpha;
}

double bm(double v)
{
  double beta;
  beta = 4.0 * exp(-v/18.0);
  return beta;
}

double ah(double v)
{
  double alpha;
  alpha = 0.07 * exp(-v/20.0);
  return alpha;
}

double bh(double v)
{
  double beta;
  beta = 1.0 / (exp((30.0 - v)/10.0) + 1.0);
  return  beta;
}

double an(double v)
{
  double alpha;
  alpha = 0.01 * (10.0 - v) / (exp(10.0 - v)/10.0 - 1.0);
  return alpha;
}

double bn(double v)
{
  double beta;
  beta = 0.125 * exp(-v/80.0);
  return beta;
}

int main(int argc, char *argv[])
{
  int i;
  double m, n, h;
  double v, vrest;
  double a, b;
  double gna, gk, gcl;
  double im, ina, ik;

  v = 0.0;
  gcl = 0.3;
  vrest = -80.0;

  m = am(v) / (am(v) + bm(v));
  h = ah(v) / (ah(v) + bh(v));
  n = an(v) / (an(v) + bn(v));

  for(i=0; i<4000; i++) {
    if(i > 200) /* 5 ms後に電流注入*/
      im = 10.0;
    else
      im = 0.0;
    m = m + (am(v)*(1-m) - bm(v)*m) * Dt;
    h = h + (ah(v)*(1-h) - bh(v)*h) * Dt;
    n = n + (an(v)*(1-n) - bn(v)*n) * Dt;
    gna = GNa * m * m * m * h;
    gk = GK * n * n * n * n;
    ina = gna*(v + vrest - ENa);
    ik = gk*(v + vrest - EK);
    v = v + (im - gna*(v + vrest - ENa) - gk*(v + vrest - EK) - gcl*(v + vrest - ECl)) / Cm * Dt;
    /* 時間，膜電流，m, h, n, ナトリウムコンダクタンス，カリウムコンダクタンス，ナトリウム電流，カリウム電流，膜電位*/
    printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", i*Dt, im, m, h, n, gna, gk, ina, ik, v + vrest);
  }
}