/*********************************************************/
/* One-dimensional time dependent Schrodinger Equation.  */
/* The scattering of a wavepack on a square barrier.     */
/* Use Crank-Nicholson/Cayley algorithm...               */
/* Stable, Norm Conserving.     Li Ju. March.28,1995     */
/*********************************************************/

#include <iostream.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "complex.h"
/* The complex operation packet */

#define MESH 1000    /* # of space intevals from 0-1 */
#define TAU  2E-6    /* Time inteval                 */
#define PI   3.1415925

Complex w[MESH+1];
/* w[] is where the wave function is kept. */

float Width,Height,K;
/* Width is the width of the barrier */
/* Height is the height              */
/* K is the wave vector              */

Complex HEAD,TAIL;
/* Dirichlet boundary condition -- might be time dependent */

/*
** Potential function
*/
float V(int i)
{
  if (abs(i-500)<Width/2) return (Height);
  else return(0.);
}

/*
** Initialize the wave pack.
*/
void InitWave()
{
  int i;
  double module;
  for (i=100;i<350;i++)
    {
      module = 1.;
      w[i].real = module*cos(K*i/MESH);
      w[i].img  = module*sin(K*i/MESH);
    }

/*  for (i=1;i<(MESH-Width)/2-1;i++)
    {
      module = exp(-0.5*((i-250.)/40)*((i-250.)/40.));
      w[i].real = module*cos(K*i/MESH);
      w[i].img  = module*sin(K*i/MESH);
    }
*/
  /* The wavepack is Guassian, character width = 40/1000 */
  /* It's also moving,  wave vector = 50*PI              */
}

/* 
** Cayley propagator.
** Pass in w, gives out w.
*/
void Cayley()
{
 static Complex alpha[MESH+1],beta[MESH+1],gamma[MESH+1];
 Complex AP(0,-0.5*MESH*MESH*TAU);
 Complex A0,A00(1,MESH*MESH*TAU),A0V(0,0.5*TAU);
 Complex bi;

 int i,N=MESH;

 alpha[N-1] = 0;
 beta[N-1]  = TAIL;
 
 for (i=N-1;i>=1;i--)
   {
     A0 = A00+A0V*V(i);
     bi = (2-A0)*w[i]-AP*(w[i-1]+w[i+1]);
     gamma[i] = -1/(A0+AP*alpha[i]);
     alpha[i-1] = gamma[i]*AP;
     beta[i-1]  = gamma[i]*(AP*beta[i]-bi);
   }

 w[0] = HEAD;
 for (i=0;i<=N-1;i++)
   w[i+1] = alpha[i]*w[i]+beta[i];
}


int main()
{
 FILE * OUTPUT;
 long Steps,Shot,i,j;

/* Total time step,Barrier width,height */
 
 scanf ("%ld %ld %f %f %f",&Steps,&Shot,&Width,&Height,&K);
 
 OUTPUT = fopen("a.out","w");
 InitWave();  

 for (i=0;i<=Steps;i++)
   {
     Cayley();
     if (i%Shot == 0)
      {
	for (j=0;j<=MESH-1;j++)
	  fprintf (OUTPUT,"%lf\n",w[j].Abs());
	printf ("Running to Step %ld\n",i);
      }
   }
 return(1);
}