/*
 * wave1-initial.c
 */

#include <stdio.h>
#include <math.h>

double pi = M_PI;

void prompt()
{
  /* 初期値を選ぶ */
  printf("1: 定常波(1) 2: 定常波(3) 3: 2つのモードの和\n");
  printf("4: 割れる三角山 5: 動く三角山\n");
  printf("6: 割れる滑らか山 7: 動く滑らか山  8: 二つの動く滑らか山\n");
  printf("9: ホイヘンスの原理反例\n");
  printf("10: 割れる滑らか山(2) 11: 動く滑らか山(2)  12:  二つの動く山(2)\n");
  printf("13: 割れる滑らか山(3) 14: 動く滑らか山(3)  15:  二つの動く山(3)\n");
  printf("nfunc(初期値の種類 1..15)=");
}


/* 初期値 (nfunc=3,4) を作るための道具 */
double sqr(double x)
{
  return x * x;
}

double mountain(double x)
{
  double a = 0.4, b = 0.6;
  if (x <= a || x >= b)
    return 0;
  else {
    double c = sqr(sqr(2 / (b - a)));
    return c * sqr((x - a) * (x - b));
  }
}

double d_mountain(double x)
{
  double a = 0.4, b = 0.6;
  if (x <= a || x >= b)
    return 0;
  else {
    double c = sqr(sqr(2 / (b - a)));
    return 2 * c * (x - a) * (x - b) * (2 * x - (a + b));
  }
}

/* 初期値 (nfunc=5) を作るための道具 */

/* これは初期値ではない！ */
double f(double x)
{
  if (x <= 0)
    return 0;
  else
    return exp(- 1.0 / x);
}

double d_f(double x)
{
  if (x <= 0)
    return 0;
  else
    return exp(- 1.0 / x) / sqr(x);
}

/* これは初期値ではない！ */
double g(double x)
{
  return f(x) / (f(x) + f(1 - x));
}

double d_g(double x)
{
  return (f(1-x)*d_f(x)+d_f(1-x)*f(x))/sqr(f(x)+f(1-x));
}

/*
 *             1   (|x|<a)
 *  h_ab(x) ={ 
 *             0   (|x|>b)
 *
 */
double h_ab(double x, double a, double b)
{
  return g((x + a) / (a - b)) * g((- x + a) / (a - b));
}

double d_h_ab(double x, double a, double b)
{
  return (d_g((x + a) / (a - b)) * g((- x + a) / (a - b))
	  - g((x + a) / (a - b)) * d_g((- x + a) / (a - b))) / (a - b);
}

/* 初期値用の道具 */
double sankaku(double x)
{
    if (x > 0.375 && x <= 0.5)
      return 8.0 * x - 3.0;
    else if (x > 0.5 && x < 0.625)
      return -8.0 * x + 5.0;
    else
      return 0.0;
}

double d_sankaku(double x)
{
    if (x > 0.375 && x < 0.5)
      return 8.0;
    else if (x > 0.5 && x < 0.625)
      return - 8.0;
    else
      return 0.0;
}

/* */

double C = 0.01;

double cos_mountain(double x)
{
  if (fabs(x) < pi * C)
    return 0.5 * (cos(x / C) + 1);
  else
    return 0;
}

double d_cos_mountain(double x)
{
  if (fabs(x) < pi * C)
    return - 0.5 * sin(x / C) / C;
  else
    return 0;
}

/* */
double A = 0.05, B = 0.01;

/* 初期値 φ=u(・,0) */
double phi(double x, int nfunc)
{
  switch (nfunc) {
  case 1:
    return sin(pi * x);
  case 2:
    return sin(3 * pi * x);
  case 3:
    return 0.5 * (sin(pi * x) + sin(3 * pi * x));
  case 4:
    return sankaku(x);
  case 5:
    return sankaku(x);
  case 6:
    {
      return cos_mountain(x - 0.5);
    }
  case 7:
    {
      return cos_mountain(x - 0.5);
    }
  case 8:
    {
      return 0.4 * cos_mountain(x - 0.2) + 0.2 * cos_mountain(x - 0.8);;
    }
  case 9:  /* ホイヘンス反例 */
    return mountain(x);
  case 10:
    {
      return h_ab(x - 0.5, A, B);
    }
  case 11:
    {
      return h_ab(x - 0.5, A, B);
    }
  case 12:
    {
      return 0.4 * h_ab(x - 0.2, A, B) + 0.2 * h_ab(x - 0.8, A, B);
    }
  case 13:
    return mountain(x);
  case 14:
    return mountain(x);
  case 15:
    return 0.4 * mountain(x + 0.3) + 0.2 * mountain(x - 0.3);
  default:
    return 0;
  }
}

/* 初期値 ψ=ut(・,0) */
double psi(double x, int nfunc)
{
  switch (nfunc) {
  case 1:
    return 0.0;
  case 2:
    return 0.0;
  case 3:
    return 0.0;
  case 4:
    return 0.0;
  case 5:
    return - d_sankaku(x);
  case 6:
    return 0.0;
  case 7:
    return - d_cos_mountain(x - 0.5);
  case 8:
    {
      return - 0.4 * d_cos_mountain(x - 0.2) + 0.2 * d_cos_mountain(x - 0.8);;
    }
  case 9: /* ホイヘンス反例 */
    return 5 * mountain(x);
  case 10:
    return 0.0;
  case 11:
    {
      return - d_h_ab(x - 0.5, A, B);
    }
  case 12:
    {
      return - 0.4 * d_h_ab(x - 0.2, A, B) + 0.2 * d_h_ab(x - 0.8, A, B);
    }
  case 13:
    return 0.0;
  case 14:
    return - d_mountain(x);
  case 15:
    return - 0.4 * d_mountain(x + 0.3) + 0.2 * d_mountain(x - 0.3);
  default:
    return 0.0;
  }
}

/*
 * ベクトル b をベクトル a にコピー
 */
void copy_vector(int N, double *a, double *b)
{
  int i;
  for (i = 0; i < N; i++)
    a[i] = b[i];
  return;
}