/* 
************************************************************************
*  wave.c:   Solution of wave equation using with finite differences   *
*     du/dt+c*du/dx (c>0)
*     1: FTCS
*     2: Lax
*     3: Lax-Wendroff
*     4: 1st-order up-wind (kazakami)
************************************************************************
*/

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

#define nx 51			/* grid size */

main()
{ 
   int imode, i, j, nstep;
   double c, nu, x, dx, t, dt, u[nx], du[nx];     
   
   FILE *output;			/* save data in file wave.dat */
   output = fopen("wave.dat","w");

   printf("wave equation\n");
   printf("du/dt+c*du/dx (c>0)\n");
   printf("1: FTCS\n");
   printf("2: Lax\n");
   printf("3: Lax-Wendroff\n");
   printf("4: 1st-order up-wind (kazakami)\n");
   printf("Input: imode ?\n");
   scanf("%d", &imode);
   printf("Input: nu ?\n");
   scanf("%lf", &nu);
   printf("Input: step number ?\n");
   scanf("%d", &nstep);

   c=1.0;
   dx=1.0/nx;
   dt=nu*dx/c;

   t=0.0;

   for(j=0; j<nx; j++) {
        x=dx*j;
        if(x < 0.5) 
	  {
            u[j]=1.0;
	  }
        else
	  {
	    u[j]=0.0;
	  }
      }

   for(j=0 ; j<nx; j++)  
     {
       fprintf(output, "%f\n", u[j]);
     }
   fprintf(output, "\n");			/* empty line for gnuplot */
   
   for(i=1; i<=nstep; i++)        	/* loop over max timesteps */
   {
      for(j=1; j<(nx-1); j++)         /* loop over space */
      {
	du[j]=-nu/2.0*(u[j+1]-u[j-1]);
	if(imode == 1) 
	  du[j]=du[j];
	else if(imode == 2) 
	  du[j]=du[j]+1.0/2.0*(u[j+1]-2*u[j]+u[j-1]);
	else if(imode == 3) 
	  du[j]=du[j]+nu*nu/2.0*(u[j+1]-2*u[j]+u[j-1]);
	else if(imode == 4) 
	  du[j]=du[j]+nu/2.0*(u[j+1]-2*u[j]+u[j-1]);
      }
      for(j=1; j<(nx-1); j++)         /* loop over space */
      {
	u[j] = u[j] + du[j];
      }
      
      for(j=0 ; j<nx; j++)  
	{
	  fprintf(output, "%f\n", u[j]);
	}
      fprintf(output, "\n");			/* empty line for gnuplot */
 
   }

   printf("data stored in heat.dat\n");
   fclose(output);

}