/*********************************************************************/
/* Numerical hydrodynamics in one-dimensional space                  */
/*         sample program using Lax-Wendroff method                  */
/*                                                                   */
/*                                      11/18/93  by R.Matsumoto     */
/*                         modified     09/25/97  by R.M.            */
/*********************************************************************/

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

#define NX 1002
#define SAFETY 0.4

/* global variables */

float   dt,tm;
float   qav,gm;
int     ns,nstop;
int	ix;
float 	x[NX],dx[NX],dxm[NX],w0[NX],w1[NX],dx0;
float   ro[NX],rvx[NX],ee[NX],pr[NX],g;
FILE    *fp1,*fp2,*fp3;

main()
{
	int i,nstep;

        tm = 0.0;
        ns   = 0;

	read_parameters();
        initialize_fp();
        generate_grid();
        shock_tube_model();
        minimum_dt();
	printf("dt = %f\n",dt);
        for(nstep=1; nstep<=nstop; nstep++) {
                ns = ns+1;
                tm = tm+dt;
                Lax_Wendroff_1d();
                apply_boundary_condition();
                minimum_dt();
        }
        printf("ns = %d  time = %f\n",ns,tm);
        graph_output();
        fclose(fp1);
        fclose(fp2);
        fclose(fp3); 
}

read_parameters()
{
       scanf("%f %f %d %d",&qav,&gm,&ix,&nstop);
       dx0   = 100.0/(ix-1);
}

initialize_fp()
{
        if ((fp1=fopen("output.density","w+"))==NULL) {
            fprintf(stderr,"cannot open output.density\n");
            exit(-2);
        };
        if ((fp2=fopen("output.velocity","w+"))==NULL) {
            fprintf(stderr,"cannot open output.velocity\n");
            exit(-2);
        };
        if ((fp3=fopen("output.pressure","w+"))==NULL) {
            fprintf(stderr,"cannot open output.pressure\n");
            exit(-2);
        };
}

generate_grid() 
{
	int i;

        dx[0] = dx0;
        x[0]  = 0.0;

        for(i=0; i<ix; i++) {
                x[i+1] =x[i]+dx[i];
                dx[i+1]=dx[i];
        }

        for(i=1; i<ix; i++) {
                dxm[i]=0.5*(dx[i-1]+dx[i]);
                w1[i]=0.5*dx[i-1]/dxm[i];
                w0[i]=0.5*dx[i]/dxm[i];
        }
}

shock_tube_model()
{
        int i;

	g=0.0;

        for(i=0; i<=ix; i++) {
               if (x[i] <= 20.0) {
                        ro[i] = 50.0;
                        rvx[i]=100.0;
                        pr[i] =100.0;
               }
               else {
                        ro[i] = 10.0;
                        rvx[i]= 0.0;
                        pr[i] = 0.0;
               }
        }
}

graph_output()
{
        int i;

        for(i=0; i<ix; i++) fprintf(fp1,"%9.4f %9.4f\n",x[i],ro[i]);
        for(i=0; i<ix; i++) fprintf(fp2,"%9.4f %9.4f\n",x[i],rvx[i]/ro[i]);
        for(i=0; i<ix; i++) fprintf(fp3,"%9.4f %9.4f\n",x[i],pr[i]);       
}

minimum_dt()
{
        int i;
        float dtq[NX],vx,dtmin;

        for(i=1; i<ix; i++) {
                vx = rvx[i]/ro[i];
                dtq[i]=SAFETY*dx[i]/sqrt(vx*vx+gm*pr[i]/ro[i]+1.0e-10);
        }

	dt = 1.0e10;
        for(i=1; i<ix; i++) 
                if (dtq[i] < dt) dt=dtq[i];
}

Lax_Wendroff_1d()
{
        int i;

        float roh[NX],rvxh[NX],eeh[NX];
        float dro[NX],drvx[NX],dee[NX];
        float f[NX],s[NX],vx[NX];

/* initialize dro etc */

        for(i=0; i<=ix; i++) {
                dro[i] = 0.0;
                drvx[i]= 0.0;
                dee[i] = 0.0;
              }

/* half step */

	for(i=0; i<=ix; i++) { 
		vx[i]=rvx[i]/ro[i];
		ee[i]=pr[i]/(gm-1.0)+0.5*rvx[i]*vx[i];
	}

        for(i=0; i<=ix; i++) {
                f[i] = rvx[i];
                s[i] = 0.0;
        }
        half_step(ro,roh,f,s);
        
        for(i=0; i<=ix; i++) {
                f[i] = rvx[i]*vx[i]+pr[i];
                s[i] = ro[i]*g;
        }
        half_step(rvx,rvxh,f,s);

        for(i=0; i<=ix; i++) {
                f[i] = rvx[i]*(gm*pr[i]/(gm-1.0)/ro[i]+0.5*vx[i]*vx[i]);
                s[i] = rvx[i]*g;
        }
        half_step(ee,eeh,f,s);

/* full step */

	for(i=0; i<ix; i++) {
		vx[i]=rvxh[i]/roh[i];
                pr[i]=(gm-1.0)*(eeh[i]-0.5*rvxh[i]*vx[i]);
	}

        for(i=0; i<ix; i++) {
                f[i] = rvxh[i];
                s[i] = 0.0;
        }
        full_step(dro,f,s);

        for(i=0; i<ix; i++) {
                f[i] = rvxh[i]*vx[i]+pr[i];
                s[i] = roh[i]*g;
        }
        full_step(drvx,f,s);

        for(i=0; i<ix; i++) {
                f[i] = rvxh[i]*(gm*pr[i]/(gm-1.0)/roh[i]+0.5*vx[i]*vx[i]);
                s[i] = rvxh[i]*g;
        }
        full_step(dee,f,s);

/* artificial viscosity */

        for(i=0; i<ix; i++) {
                f[i] = qav*dx[i]*fabs(rvx[i+1]/ro[i+1]-rvx[i]/ro[i]);
	}

        artificial_viscosity(ro,dro,f);
        artificial_viscosity(rvx,drvx,f);
        artificial_viscosity(ee,dee,f);

/* update internal points */

        for(i=1; i<ix; i++) {
                ro[i] = ro[i] +dro[i];
                rvx[i]= rvx[i]+drvx[i];
                ee[i] = ee[i] +dee[i];
        }

/* convert from total energy to pressure */

        for(i=1; i<ix; i++) 
                pr[i] = (gm-1.0)*(ee[i]-0.5*rvx[i]*rvx[i]/ro[i]);
                      
}


apply_boundary_condition()
{

/* free boundary */

        ro[0] = ro[1];
        rvx[0]= rvx[1];
        pr[0] = pr[1];

        ro[ix] = ro[ix-1];
        rvx[ix]= rvx[ix-1];
        pr[ix] = pr[ix-1];
}


half_step(u,uh,f,s)
float *u,*uh,*f,*s;
{
        int i;
	float dfdx,sh;

        for(i=0; i<ix; i++) {
		dfdx = (f[i+1]-f[i])/dx[i];
		sh   = 0.5*(s[i+1]+s[i]);
                uh[i]= 0.5*(u[i+1]+u[i])-0.5*dt*(dfdx+sh);
	}
}

full_step(du,f,s)
float *du,*f,*s;
{
        int i;
	float dfdx,sh;

        for(i=1; i<ix; i++) {
                dfdx  = (f[i]-f[i-1])/dxm[i];
                sh    = w0[i]*s[i-1]+w1[i]*s[i];
                du[i] = du[i]-dt*(dfdx+sh);
        }
}

artificial_viscosity(u, du, qx)
float *u,*du,*qx;
{
        int i;

        for(i=1; i<ix; i++) {
                du[i] = du[i]+0.5*dt*(qx[i]*(u[i+1]-u[i])/dx[i]
                                 -qx[i-1]*(u[i]-u[i-1])/dx[i-1])/dxm[i];
              }
}