/* 

************************************************************************

*  soliton.c: Solves the Kortewg-deVries Equation using a finite       *

*             difference method                                        *

*                                                                      *

*  taken from: "Projects in Computational Physics" by Landau and Paez  * 

*              copyrighted by John Wiley and Sons, New York            *      

*                                                                      *

*  written by: students in PH465/565, Computational Physics,           *

*              at Oregon State University                              *

*              code copyrighted by RH Landau                           *

*  supported by: US National Science Foundation, Northwest Alliance    *

*                for Computational Science and Engineering (NACSE),    *

*                US Department of Energy                               *

*                                                                      *

*  UNIX (DEC OSF, IBM AIX): cc soliton.c -lm                           *

*                                                                      *

*  comment: Output data is saved in 3D grid format used by gnuplot     *

************************************************************************

*/

#include <stdio.h>

#include <math.h>



#define ds 0.4				/* delta x */

#define dt 0.1				/* delta t */

#define max 2000			/* time steps */

#define mu 0.1				/* mu from KdeV equation */

#define eps 0.2				/* epsilon from KdeV eq. */ 



main()

{

   int i, j, k;

   double a1,a2,a3, fac, time, u[131][3];

   

   FILE *output;			/* save data in soliton.dat */

   output = fopen("soliton.dat","w");

   

   for(i=0; i<131; i++)			/* initial wave form */

   {

      u[i][0] = 0.5*(1.0 - tanh(0.2*ds*i - 5.0));

   }

   u[0][1]   =1.;                       /* end points */

   u[0][2]   =1.;   

   u[130][1] =0.;  

   u[130][2] =0.;

   

   fac  = mu*dt/(ds*ds*ds);             

   time = dt;

   

   for(i=1; i<130; i++)			/* first time step */

   {

      a1 = eps*dt*(u[i+1][0] + u[i][0] + u[i-1][0]) / (ds*6.0);    

           

      if((i>1) && (i<129))

      { 

           a2 = u[i+2][0] + 2.0*u[i-1][0] - 2.0*u[i+1][0] - u[i-2][0]; 

      }

      else a2 = u[i-1][0] - u[i+1][0];

         

           a3 = u[i+1][0]-u[i-1][0];

      u[i][1] = u[i][0] - a1*a3 - fac*a2/3.;       

   }   

                                 

   

   for(j=1; j<max; j++)			/* all other time steps */

   {   

      time+=dt;

      for(i=1; i<130; i++)

      {

         a1 = eps*dt*(u[i+1][1] + u[i][1] + u[i-1][1]) / (3.0*ds);

              

         if((i>1) && (i<129))

         {

            a2 = u[i+2][1] + 2.0*u[i-1][1] - 2.0*u[i+1][1] - u[i-2][1];           

         }

         else a2 = u[i-1][1] - u[i+1][1]; 

         a3      = u[i+1][1] - u[i-1][1];

         u[i][2] = u[i][0] - a1*a3 - 2.*fac*a2/3.;     

      }



      for(k=0; k<131; k++)		/* move one step ahead */

      {

         u[k][0] = u[k][1];

         u[k][1] = u[k][2];

      }

      

      if((j%200)==0)			/* plot every 200th step */

      {

         for(k=0; k<131; k+=2)

         {

            fprintf(output, "%f\n", u[k][2]);   /* gnuplot 3D format */

         }

         fprintf(output, "\n");         /* empty line for gnuplot */

      }

   } 

   printf("data stored in soliton.dat\n");

   fclose(output);

}