/*********************************************************/ /* 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); }