/*
  生体システム工学の基礎
  図4.16
  gcc -o c3 c3.c -lm
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 50
#define T 40000
/* HH */
#define GNa 120.0
#define GK 36.0
#define ENa 35.0
#define EK -92.0
#define ECl -84.06686
#define Cm 1.0

/* HH */

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, j;
  double dt, dx, lambda2;
  double *u1, *u2, *u3, *m, *n, *h;
  double v, gcl, vrest, v1, im;

  u1 = calloc(N+1, sizeof(double)); 
  if(u1 == NULL) {
    fprintf(stderr, "u1: memory not enough\n");
    exit(1);
  }
  u2 = calloc(N+1, sizeof(double)); 
  if(u2 == NULL) {
    fprintf(stderr, "u2: memory not enough\n");
    exit(1);
  }
  u3 = calloc(N+1, sizeof(double)); 
  if(u3 == NULL) {
    fprintf(stderr, "u3: memory not enough\n");
    exit(1);
  }
  m = calloc(N+1, sizeof(double)); 
  if(m == NULL) {
    fprintf(stderr, "u3: memory not enough\n");
    exit(1);
  }
  n = calloc(N+1, sizeof(double)); 
  if(n == NULL) {
    fprintf(stderr, "u3: memory not enough\n");
    exit(1);
  }
  h = calloc(N+1, sizeof(double)); 
  if(h == NULL) {
    fprintf(stderr, "u3: memory not enough\n");
    exit(1);
  }

  dx = 0.5 / N;
  dt = 0.1 * dx;

  /* 初期条件 */
  for(i=0; i<=N; i++) {
    u1[i] = 0.0;
    u2[i] = 0.0;
  }

  v = 0.0;
  gcl = 0.3;
  vrest = -80.0;
  for(i=0; i<=N; i++){
    m[i] = am(v) / (am(v) + bm(v));
    h[i] = ah(v) / (ah(v) + bh(v));
    n[i] = an(v) / (an(v) + bn(v));
  }

  /* 初期値表示 */
  for(i=0; i<=N; i++) {
    printf("%lf ", u2[i]);
  }
  printf("\n");
  /* 時間発展 */
  im = 10.0;
  for(j=1; j<=T; j++) {
    for(i=1; i<N; i++) {
      m[i] = m[i] + dt * (am(u2[i])*(1-m[i]) - bm(u2[i])*m[i]);
      h[i] = h[i] + dt * (ah(u2[i])*(1-h[i]) - bh(u2[i])*h[i]);
      n[i] = n[i] + dt * (an(u2[i])*(1-n[i]) - bn(u2[i])*n[i]);
      /* 計算できるように0.0001にしたので空間方向は適切な長さの単位ではない */
      u3[i] = u2[i] + 0.0001 * dt / (dx*dx) *(u2[i+1] - 2*u2[i] + u2[i-1]) 
	- dt * GK * n[i] * n[i] * n[i] * n[i] * (u2[i] + vrest - EK)
	- dt * GNa * m[i] * m[i] * m[i] * h[i] * (u2[i] + vrest - ENa)
	- dt * gcl * (u2[i] + vrest - ECl);
      if(i==1) { /* 1ユニット目には注入電流による電位を加算する */
	u3[i] = u3[i] + dt * im;
      }
    } 
    /* 自由端 */ 
    u3[0] = u3[1];
    u3[N] = u3[N-1];

    /* コピー */
    for(i=0; i<=N; i++) {
      u1[i] = u2[i];
      u2[i] = u3[i];
    }
    /* nステップ目表示 */
    if(j%10 == 1) { /*全ステップ表示すると多いので1/10に間引く*/
      for(i=0; i<=N; i++) {
	printf("%lf ", u2[i]);
      }
      printf("\n");
    }
  }
}