#include <stdio.h>
#include <math.h>
#include <graphics.h>
#include <conio.h>
 
Void DrawFunction (int Npoint, double data[], double H) ;/* Display results */
 
Void MakePotential (int Npoint, double x0, double dx, double V0, double v[], double vd[]) ; /*potential*/
 
Void MakePacket (int Npoint, double x0, double dx, double k0, double *norm1, double yr[], double yi[]) ; /*make wavepacket*/
 
Void TimeIntegral (int Npoint, double dx2, double dt, double *norm2, double yr[], double yi[], double yd[], double v[], double a[], double b[]);
 
        double a[4] = {0.0, 0.5, 0.5, 1.0};
        double b[4] = {1/6.0, 1/3.0, 1/3.0, 1/6.0};
        double yr[401],yi[401],yd[401],v[401],vd[401];
        double kr[401],ki[401],dyr[401],dyi[401],ykr[401],yki[401];
void main(void)
{
        int      t, I, itest, Nmax;
        double  V0, x0, x, dx, dx2, dt, k0, E, norm1, norm2;
        int      graphriver = DETECT, graphmode;
 
        initgraph( &graphriver, &graphmode , “ “);
        setviewport(200, 50, 600, 350, 1);
        rectangle(0,0,400,300);
        dx  = 0.02;     
Nmax = 400;
        dx2 = dx*dx;
        dt  = dx2/2.0;
        x0  = 0.0; 
        V0  = 400.0;
 
         MakePotential (Nmax, x0, dx, V0, v, vd);
         DrawFunction (Nmax, vd, 0);
        Gotoxy(1,1);
        printf(“Input wave number (k0<25) ! “);
        scanf(“%1f”, &k0);   E=k0*ko;
         printf(“ Time  Norm”);
         norm 1=0;
         MakePacket(Nmax, x0, dx, k0, &norm1, yr, yi);
         norm2=norm1;
    for( t=0; t<=600; t++)
       {   gotoxy(1,3); printf(“%5d %8.41f”,t, norm2/norm1); printf(“%d \n”,t);
 
    TimeIntegral (Nmax, dx2,dt, %norm2, yr,yi,yd,v,a,b);
 
    itest=fmod(t,10);
    if (itest==0) { clearviewport();
          rectangle(0,0,400,300);
          setcolor(4); DrawFunction (Nmax,yr,E/V0); 
          setcolor(1); DrawFunction (Nmax,yr,E/V0); 
          setcolor(7); DrawFunction (Nmax,yr,E/V0);
          DrawFunction (Nmax,vd,0);    }
}
     printf(“ Calculation ends,if ready, please return key!”);
     getche();
     closegraph();
}      /* -- End of main function -- */
 
void DrawFunction(int Npoint,double data[], double H)
       {       int n, iy +
               moveto(0,250-(data[0]+H)*100);
               for( n=0; n<=Npoint-1; n++)
               { iy=250-(data[n]+H*100;     lineto(n, iy);
               }
       }
 
void MakePotential (int Npoint,double x0, double dx, double V0,double v[], double vd[])
       {  int    n;
          double x;
             x=x0 ;
             for ( n=0; n<=Npoint-1; n++ )
              { if (n<220) vd[n]=0 ;else if (n>230) vd[n]=0;
               else vd[n]=1.0;
               v[n]=V0*vd[n];    x+=dx;  }
       }
 
void MakePacket (int Npoint,double x0, double dx, double k0, double *norm1, double yr[], double yi[])
       {  int    n;
          double x;
             x=x0 ;
             for ( n=0; n<=Npoint-1; n++ )
              { yr[n]= cos(k0*x}*exp(-(x-2.5)*(x-2.5)*3.0) ;
               yi[n]= sin(k0*x}*exp(-(x-2.5)*(x-2.5)*3.0) ;
               x+=dx;
               *norm1=norm1+yr[n]*yr[n]+yi[n]*yi[n]; }
       }
 
void TimeIntegralt (int Npoint,double dx2, double dt, double *norm2, double yr[], double yi[], double yd[], double v[], 
                   double a[], double b[])
       {  int    nx,ij ;
             for ( nx=0; nx<=Npoint; nx++ )
                 { dyr[nx] = 0.0; dyi[nx] = 0.0 ;
                  kr[nx]  = 0.0;  ki[nx] = 0.0; }
 
             for ( ij=0 ;ij<=3 ;ij++)    /* I loop */
                 {
                  for ( nx=0; nx<=Npoint; nx++ )
                     { ykr[nx]= yr[nx] + a[ij]*kr[nx] ;
                       yki[nx]= yi[nx] + a[ij]*ki[nx] ; }
 
                  for ( nx=1; nx<=Npoint-1; nx++ )   
                     { kr[nx]=-((yki[nx+1]-2.0*yki[nx]+yki[nx-1])/dx2-v[nx]*yki[nx])*dt;
                      ki[nx]=-((ykr[nx+1]-2.0*ykr[nx]+ykr[nx-1])/dx2-v[nx]*yki[nx])*dt;
                     }
 
                  for ( nx=0; nx<=Npoint; nx++ )
                     { dyr[nx] = dyr[nx] + b[ij]*kr[nx] ;
dyi[nx] = dyi[nx] + b[ij]*ki[nx] ; }
                 }                         /* i-loop end */
 
                *norm2=0.0;
                for ( nx=0; nx<=Npoint; nx++ )
                   { yr[nx] = yr[nx] + dyr[nx];
                    yi[nx] = yi[nx] + dyi[nx];
                    yd[nx] = yr[nx]*yr[nx]+yi[nx]*yi[nx];
                    *norm2=*norm2+yd[nx] ;
                   }
}