/*
  生体システム工学の基礎
  図4.30
  gcc -o kitem3 kitei3.c -l

 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 1000
#define T 2000
#define Lambda 0.5
#define alpha 500.0 /* 計算できるように定めた値 */

/* インピーダンスはxに対して線形に変化する */
double zzz(double x)
{
  double y;
  y = (0.01 - 0.5)/(1.0 - 0.0) * x + 0.5;
  return y;
}

/* 前庭窓に加わる力 */
double force(double t)
{
  double y, frq;
  frq = 10; /* 周波数 */
  if ((t >= 0.0) && (t <= 1 * 1.0 /frq) ) 
    y = 10 * sin(frq*2*M_PI*t);
  else
    y = 0.0;
  return y;
}

int main(int argc, char *argv[])
{
  int i, j;
  double dt, dx, lambda2;
  double *u1, *u2, *u3, *x;

  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);
  }
  x = calloc(N+1, sizeof(double)); 
  if(x == NULL) {
    fprintf(stderr, "u3: memory not enough\n");
    exit(1);
  }

  dx = 1.0 / N;
  dt = Lambda * dx;
  lambda2 = Lambda * Lambda;

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

  for(i=0; i<=N; i++) {
    x[i] = x[i] + dt * u2[i];
  }

  /* 初期値表示 */
  for(i=0; i<=N; i++) {
    printf("%lf ", x[i]);
  }
  printf("\n");

  /* 時間発展 */
  for(j=1; j<=T; j++) {
    for(i=1; i<N; i++) {
      u3[i] = 2 * (1.0 - lambda2 - alpha*dt*dt/zzz(i*dx)) * u2[i] - u1[i] + lambda2 * (u2[i+1] + u2[i-1]) - dt*(1.0 + 1.0/zzz(i*dx)) * (u2[i] - u1[i]);
    } 

    u3[0] = force(j*dt);
    u3[N] = 0.0;

    /* コピー */
    for(i=0; i<=N; i++) {
      u1[i] = u2[i];
      u2[i] = u3[i];
    }
    for(i=0; i<=N; i++) {
      x[i] = x[i] + dt * u2[i];
    }
    /* nステップ目表示 */
    for(i=0; i<=N; i++) {
      printf("%lf ", x[i]);
    }
    printf("\n");
  }
}