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

#define N 200           //系の範囲
#define a 0.5           //係数aの値
#define b 0.1           //係数bの値
#define dt 0.05         //時間的刻み幅
#define dx 0.3          //系の刻み幅
int main()
{
        int t,i,tmax;
        double u[3][N];

        FILE *output;
        output=fopen("output.data","w");

        tmax=1000;              //時間的最大繰り返し回数

/*初期状態の波。この関数を使うと浅瀬へ向かう穏やかな波が表現できる。*/
        for(i=0;i<N;i++) {
                u[0][i]=(1.0-tanh((i-100.0)/1.5))/2.0;
        }

        for(i=2;i<N-2;i++) {
                u[1][i]=u[0][i]-a*dt/(dx)*u[0][i]*(u[0][i+1]-u[0][i-1])-3.0*b*dt/(8.0*dx*dx*dx)*(u[0][i+2]-u[0][i-2]-2.0*u[0][i+1]+2.0*u[0][i-1]);
        }
        u[1][0]=u[1][1]=1.0;            //境界条件
        u[1][N-1]=u[1][N-2]=0.0;

        for(t=0;t<tmax;t++) {
                for(i=2;i<N-2;i++) {
                        u[2][i]=u[0][i]-a*dt/dx*u[1][i]*(u[1][i+1]-u[1][i-1])-3.0*b*dt/(4.0*dx*dx*dx)*(u[1][i+2]-u[1][i-2]-2.0*u[1][i+1]+2.0*u[1][i-1]);
                }
                u[2][0]=u[2][1]=1.0;    //境界条件
                u[2][N-1]=u[2][N-2]=0.0;

/*次の計算のために配列の入れ替え。ここで一番古い情報は失われる。*/
                for(i=0;i<N;i++) {
                        u[0][i]=u[1][i];
                        u[1][i]=u[2][i];
                }

/*一定おきに結果をファイルに書き込むことにしている。*/
                if(t%10==0) {
                        for(i=0;i<N;i++) {
                                fprintf(output,"%d %d %f\n",t/10,i,u[2][i]);
                        }
                        fprintf(output,"\n");
                }

        }

        fclose(output);

        return 0;
}