#include <stdio.h>
#include <stdlib.h>

/* 解説: strcmp の関数定義が含まれている */
#include <string.h>

/*
つぎの1階常微分方程式を、初期値と係数を与えて数値的近似的に解く。
dx/dt=a0+a1*x+...+an*x^n

これは、非線形になりうるので、一般に解析的にとくのは難しい。

オイラー法 (Euler)
   Δxx(t+Δt)-x(t)
   Δx≒(dx/dt)Δt
∴ x(t+Δt)≒x(t)+(dx/dt)Δt

*/

#define X0 0.0
#define T0 0.0
#define DeltaT 1e-2
#define Rept 200
#define ReptMax 1000

/* 「冗舌モードを示す」グローバル変数 */
int verbose;
/* 「デバッグモードを示す」グローバル変数 */
int debug;

/* 入力関数(サブルーチン) */
int
input( coef, maxorder, file ) 
     double *coef;
     int maxorder;
     char *file;
{
  FILE *fd;
  int i, r, n;

  fd=fopen( file, "r" );
  if( fd==NULL ){
    fprintf(stderr,"ファイル %s をオープン出来ませんでした。\n", file);
    return -1;
  }

  /* 次数を読む */
  r=fscanf( fd, "%d", &n );
  if( r!=1 ){
    fprintf(stderr,"次数を読め出来ませんでした。\n");
    fclose(fd);
    return -1;
  }
  if( n<0 || n>maxorder ){
    fprintf(stderr,"次数が正しくありあません。(最大%d)\n",maxorder);
    fclose(fd);
    return -1;
  }
  if( debug ){/* デバッグモードならファイルから読んだ次数を表示 */
    fprintf(stderr, "debug:n=%d\n",n);
  }
  for( i=0 ; i<=n ; i++ ){
    r=fscanf( fd, "%lf", &coef[n-i] );/* 高い次数から読む */
    if( r!=1 ){
      fprintf( stderr, "%d次の係数が読めません。\n",n-i);
      fclose(fd);
      return -1;
    }
    if( debug ){/* デバッグモードならファイルから読んだ係数を表示 */
      fprintf(stderr, "debug:%g(%d)\n",coef[n-i],n-i);
    }
  }
  fclose( fd );
  return n;/* 次数を返す */
}

/* 多項式を計算する */
double
poly( x, coef, order )
double x;
double coef[];
int order;
{
  double y;
  int i;

  i=order;/* iを最高次数に初期化する */

  y=coef[i];
  while( --i>=0 ){
    /*解説 : --i は、変数iを1だけ減算し、減算後の値をもつ 
            これに対して i--、変数iを1だけ減算し、減算前の値をもつ */
    /* 次数を減らし次数が0以上であるあいだ繰り返す */
    y=y*x+coef[i];
  }
  return y;
}

/* 多項式を計算する */
polyout( coef, order, fd )
double coef[];
int order;
FILE *fd;
{
  int i;
  /* 係数を多項式の形で出力する */
  for( i=0 ; i<=order ; i++ ){
    fprintf(fd, " %+g*x^%d", coef[order-i], order-i );
  }
}

/* 最も単純な「オイラー法」(Euler法)による近似解を求める */
/* 1段階だけ */
double
solve( coef, order, dt, x )
double *coef, x, dt;
int order;
{
  double y;
  x=x+dt*poly( x, coef, order );
  return x;
}

/* 解説:
   argc は、コマンドラインの項目数(コマンド名＝プログラム名を含む) 
   argv[i] はコマンドラインのi番目の項目を示す文字列
           (つまり先頭文字へのポインタ)
   argv[0] はコマンド名 */

#define N 10

main( argc, argv )
     int argc;
     char *argv[];
{
  int i;
  char *file;

  int jisuu;
  double keisuu[N+1];
  double kekka[ReptMax+1];
  double t, x;

  double t0, x0;/* 初期値 */ 
  double delta_t; /* 時間の刻幅 */
  int rept;/*繰り返し数*/
  int r;

/* 標準値の代入 */
  x0=X0;
  t0=T0;
  delta_t=DeltaT;
  rept=Rept;

  verbose=0;   /* verbose モード: 標準値は「verbose でない」*/
  file=NULL;   /* NULL はファイル名が指定されてい
                  ない状態を示すことにする */

  /* コマンドラインの解析 */
  for( i=1 ; i<argc ; i++ ){
    /* 解説: strcmp( a, b ) は文字列 a と b を比較する。
       文字列が等しい時 0 を返す。
      (文字列としての大小もわかるがここでは使わない) */
    if( strcmp( argv[i], "-v")==0 ){ /*バーバスモードか*/
      verbose=1;
    }else
    if( strcmp( argv[i], "-d")==0 ){ /*デバッグモードか*/
      debug=1;
    }else
    if( strcmp( argv[i], "-i")==0 && i+1<argc ){/*初期値の指定か*/
      i++;
      r=sscanf( argv[i], "%lf,%lf", &x0, &t0 );
      if( r<1 ){
	fprintf(stderr,"'%s' 不正な初期値を無視します\n", argv[i]);
	x0=X0;t0=T0;/* 標準値に戻す */
	continue;
      }
      if( r==1 ){
	t0=0.0;/* ,t0 が省略されたときは 0 とする */
	fprintf(stderr,",t0 が省略されたので、t0=%g とします。\n",t0);
      }
    }else
    if( strcmp( argv[i], "-dt" )==0 && i+1<argc ){/* 時間きざみの指定か */
      i++;
      r=sscanf( argv[i], "%lf", &delta_t );
      if( r<1 || delta_t<=0.0 ){
	fprintf(stderr,"'%s' 不正な時間幅を無視します\n", argv[i]);
	delta_t=DeltaT;/* 標準値に戻す */
	continue;
      }
    }else
    if( strcmp( argv[i], "-r" )==0 && i+1<argc ){
      i++;
      r=sscanf( argv[i], "%d", &rept );
      if( r<1 || rept<=0 || rept>ReptMax ){
	fprintf(stderr,"'%s' 不正な繰り返し数を無視します(最大%d)\n",
		argv[i],ReptMax );
	rept=Rept;/* 標準値に戻す */
	continue;
      }
    }else{
      if( file==NULL ){
	/* ファイル名がこれまで指定されていないとき */
	file=argv[i];
      }else{
	/* ファイル名がすでに指定されている時 */
	fprintf(stderr,
"ファイルはすでに'%s'に指定されています。\nパラメタ %d '%s' は無視します。\n",
		file, i,argv[i]);
      }
    }
  }
  if( file==NULL ){  /* ファイルが指定されていない時は使い方を表示する */
    fprintf(stderr,"使い方: sample12 [-v][-i x0,t0][-r repeat][-d delta_t] file\n");
    fprintf(stderr,"           fileには、次数と係数を高次から順に入れる。\n");
    exit(1);
  }

  /* 係数を入力する */
  jisuu=input( keisuu, N, file );
  if( jisuu==-1 )exit(1);

  if( verbose ){/* バーバスモードなら */
    /* 係数を微分方程式の形で出力する */
    fprintf(stderr,"微分方程式:dx/dt=");
    polyout( keisuu, jisuu, stderr );
    /* 条件を出力する */
    fprintf( stderr,"\n初期値 x=%g(t=%g)、時間きざみ %g、繰り返し数 %d\n",
	    x0, t0, delta_t, rept );
  }

  /* 解く */
  /* (ここでは結果を kekka に一旦格納してから出力しているが
      格納せずに毎回出力する方法もある。                     */

  /* 初期値を格納する */
  i=0;
  x=x0;
  kekka[i]=x;

  for( i=1 ; i<=rept ; i++ ){    /* rept 回つぎの値の計算を繰り返す */
    x=solve( keisuu, jisuu, delta_t, x );
    kekka[i]=x;
  }

  /* 出力する */
  /* ファイル 初期値などのパラメタを書く */
  /* ファイルにしまうとき、あとでどういうデータかわかる */
  /* ( 例えば、sample12 sample12.dat > sample12.x  ) */
  /* 行の先頭に # をつける事により gnuplot で描くときにこれらの
     行は無視される。( gnuplot で plot "sample12.x" とできる ) */
  printf("# file=%s\n", file);

  /* 係数を微分方程式の形で出力する */
  printf("# dx/dt=");
  polyout( keisuu, jisuu, stdout );
  putchar('\n');

  /* 初期値 時間キザミ 繰り返し数を出力する */
  printf("# t0=%g x0=%g deltaT=%g rept=%d\n", t0, x0, delta_t, rept );

  /* データを出力する */  
  t=t0;                                /* 時間を出力するための初期化 */
  for( i=0 ; i<=rept ; i++ ){
    printf( "%g %g\n", t, kekka[i] );
    t+=delta_t;                        /* 時間を出力するための加算 */
  }
  exit(0);
}