#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<string.h>
#define max 89

/* global variables */
   double v[91][91],p2[91][91]; 
   FILE *pf,*outparms;

void potential()
{ //potential slit
  int i,j;
  for(j=0;j<=max+1;j++){
    for(i=0;i<=max+1;i++){
      if((j==36)&&(i<=40)||(j==36)&&(i>=51))v[i][j]=0.5;
      else v[i][j]=0.0;
    }//i
  }//for j  
}

void output(int n)
{
  //output name for different times
  int i,j;
  char stri[8],strin1[3],*slitstr; 
  slitstr="slit";
  strcpy(stri,slitstr);
  if(n<10){
    strin1[0]=48;
    strin1[1]=48;
    strin1[2]=n+48;  
  }
  if(n>9 &&n<100){
    strin1[0]=48;
    strin1[1]=n/10+48;
    strin1[2]=n%10+48;
  }
  if(n>99 && n<1000){
    strin1[0]=n/100+48;
    strin1[1]=(n%100)/10+48;
    strin1[2]=n%10+48;
  }
  strcat(stri,strin1);
  pf=fopen(stri, "w");
  //write probabilities plus potential
  for(j=0;j<=max;j=j+3){
    for(i=0;i<=max;i=i+2){
      fprintf(pf,"%lf\n",p2[i][j]+v[i][j]);
    }//i
    fprintf(pf,"\n"); 
  }//j  
  fclose(pf);
}

main(){
  double psr[91][91][2],psi[91][91][2];
  double a1,a2,dt,dx,k0x,k0y,x0,y0,x,y,rex,iex,x2,y2;
  int i,j,n,tiempo;
  time_t timein, timefin;
  timein=time(NULL);
  printf("Input positive integer between 1 and 800, \n");
  printf("it is the number of total time steps\n");
  scanf("%d",&tiempo);
  printf("time selected %d\n",tiempo);
  outparms=fopen("slitparm","w");

  dx=0.2;
  dt=0.0025/(dx*dx);
  k0x=0.0;
  k0y=2.5;
  x0=0.0;
  y0=-7.0;
  y=-9.0;
  potential();
  //initial wave packet
  for(j=0;j<=max+1;j++){ 
    x=-9.0;
    for(i=0;i<=max+1;i++){
      rex=cos(k0x*x+k0y*y);
      iex=sin(k0x*x+k0y*y);
      x2=(x-x0)*(x-x0);
      y2=(y-y0)*(y-y0);
      a1=exp(-0.5*(x2+y2));
      psr[i][j][0]=a1*rex;
      psi[i][j][0]=a1*iex;
      x=x+dx;
    }//i
    y=y+dx;
  }//for j
  //propagate in time
  for(n=1;n<=tiempo;n++){
    for(j=1;j<=max;j++){
      for(i=1;i<=max;i++){
        a2=v[i][j]*psi[i][j][0]+2.0*dt*psi[i][j][0];
        a1=psi[i+1][j][0]+psi[i-1][j][0]+psi[i][j+1][0]+psi[i][j-1][0];
        psr[i][j][1]=psr[i][j][0]-dt*a1+2.0*a2;
        if((n%20==0)||(n==1))p2[i][j]=psr[i][j][0]*psr[i][j][0]
                 +psi[i][j][0]*psi[i][j][0];
      }//for i
      psr[0][j][1]=psr[1][j][1];
      psr[max+1][j][1]=psr[max][j][1];
    }//for j
    //imaginary part w.f.
    for(j=1;j<=max;j++){
      for(i=1;i<=max;i++){
        a2=v[i][j]*psi[i][j][1]+2.0*dt*psr[i][j][1];
        a1=psr[i+1][j][1]+psr[i-1][j][1]+psr[i][j-1][1]+psr[i][j+1][1];
        psi[i][j][1]=psi[i][j][0]+dt*a1-2.0*a2;
      }//for i
      //derivative 0 at edges
      psi[0][j][1]=psi[1][j][1];
      psi[max+1][j][1]=psi[max][j][1];
    }//j
    //new iterations are now the old ones, recycling
      for(j=0;j<=max+1;j++){
        for(i=0;i<=max+1;i++){
          psi[i][j][0]=psi[i][j][1];
          psr[i][j][0]=psr[i][j][1];
	}// for i
      }//for j  

     if((n%20==0)||(n==1))output(n);
  }//for n
 
    timefin=time(NULL);
    fprintf(outparms,"\n Elapsed time in seconds %ld\n",timefin-timein);
    fclose(outparms);
}//program