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

#define LMAX 256
#define PI 3.14159265

void calc_F(double F[], double f[], double V[], double dx);
void RK4(double y_r[], double y_i[], double V[], double dx, double dt);

int main(int argc, char* argv[])
{
	double y_r[LMAX], y_i[LMAX], V[LMAX], Vo;
	double dx, dt, k;
	double center, delta, gauss; /*gaussian wave packet parameters */
	int ix, time, time_max, w_time;

	Vo = 0.5;
	/* set potential */
	for(ix=0;ix LMAX/2 && ix < LMAX/2 + 10 ) V[ix]=Vo;
		else V[ix] = 0.0;
	}
    
	dx = 1.; k = 2.*PI / (dx*10);
	center = 50.; delta = 10.;
	/* set initial wave function */
        for(ix=0;ix < LMAX;ix++){
       gauss = exp(- (ix*dx - center)*(ix*dx - center)/delta/delta);
	   y_r[ix] = gauss * cos (k*ix*dx);
	   y_i[ix] = gauss * sin (k*ix*dx);
	}
		
	    printf("0\n");
		for(ix=0;ix<LMAX;ix++)printf("%f\t",V[ix]);
		printf("\n");          
		for(ix=0;ix<LMAX;ix++)printf("%f\t",y_r[ix]);
		printf("\n");
		for(ix=0;ix<LMAX;ix++)printf("%f\t",y_i[ix]);
		printf("\n");
		printf("\n");		   

	w_time = 100; dt=0.1;time_max=2000;
	for(time=1;time<=time_max;time++){
       RK4(y_r, y_i, V, dx, dt);
	   if(time%w_time==0){
		  printf("%d\n", time);
		  for(ix=0; ix<LMAX; ix++)printf("%f\t",y_r[ix]);
		  printf("\n");
		  for(ix=0; ix<LMAX; ix++)printf("%f\t",y_i[ix]);
		  printf("\n");
		  printf("\n");
	   }
	}
	return 0;
}

void calc_F(double F[], double f[], double V[], double dx){
	int ix;
	for(ix=1; ix<LMAX-1; ix++){
		F[ix]= (f[ix+1] - 2.* f[ix] + f[ix-1])/dx/dx - V[ix]*f[ix];
	}

	return;
}

void RK4(double y_r[], double y_i[], double V[], double dx, double dt){
    static double  F[LMAX];
	static double k1_r[LMAX], k1_i[LMAX], k2_r[LMAX], k2_i[LMAX],
		          k3_r[LMAX], k3_i[LMAX], k4_r[LMAX], k4_i[LMAX];
	static double y_k_i[LMAX], y_k_r[LMAX];
	static int ix;
	
	/* k1 calculation */
	calc_F(F, y_r, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k1_i[ix] = F[ix]*dt;
	}
	calc_F(F, y_i, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k1_r[ix] = -F[ix]*dt;
	}

	/* k2 calculation */
	for(ix=0; ix<LMAX; ix++){
		y_k_i[ix] = y_i[ix] + 0.5*k1_i[ix]; /* y + 0.5*k1 */
		y_k_r[ix] = y_r[ix] + 0.5*k1_r[ix];
	}
	calc_F(F, y_k_r, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k2_i[ix] = F[ix]*dt;
	}
	calc_F(F, y_k_i, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k2_r[ix] = -F[ix]*dt;
	}

		/* k3 calculation */
	for(ix=0; ix<LMAX; ix++){
		y_k_i[ix] = y_i[ix] + 0.5*k2_i[ix]; /* y + 0.5*k2 */
		y_k_r[ix] = y_r[ix] + 0.5*k2_r[ix];
	}
	calc_F(F, y_k_r, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k3_i[ix] = F[ix]*dt;
	}
	calc_F(F, y_k_i, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k3_r[ix] = -F[ix]*dt;
	}

		/* k4 calculation */
	for(ix=0; ix<LMAX; ix++){
		y_k_i[ix] = y_i[ix] + k3_i[ix]; /* y + k3 */
		y_k_r[ix] = y_r[ix] + k3_r[ix];
	}
	calc_F(F, y_k_r, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k4_i[ix] = F[ix]*dt;
	}
	calc_F(F, y_k_i, V, dx);
    for(ix=1; ix<LMAX-1; ix++){
		k4_r[ix] = -F[ix]*dt;
	}

	for(ix=0; ix<LMAX; ix++){
		y_i[ix] += (k1_i[ix] + 2.*(k2_i[ix] + k3_i[ix]) + k4_i[ix])/6.; /* y + k3 */
		y_r[ix] += (k1_r[ix] + 2.*(k2_r[ix] + k3_r[ix]) + k4_r[ix])/6.;
	}

	return;
}