\File{twodim.c},{10:59},{Jan  5 1998}
\L{\LB{\C{}\/* \$Source: \/mit\/bert\/Cvsroot\/Thesis\/src\/exact\/twodim.c,v \$}}
\L{\LB{ *}}
\L{\LB{ * 2\-dimensional solution of Schr{\accent"7Fo}dinger\'s equation for a wavefunction}}
\L{\LB{ *\/\CE{}}}
\L{\LB{}}
\L{\LB{\K{\#include} \<stdio.h\>}}
\L{\LB{\K{\#include} \<stdlib.h\>}}
\L{\LB{\K{\#include} \<string.h\>}}
\L{\LB{\K{\#include} \<errno.h\>}}
\L{\LB{\K{\#include} \<sys\/types.h\>}}
\L{\LB{\K{\#include} \<time.h\>}}
\L{\LB{\K{\#include} \S{}\"complex.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"fft.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"oppack.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"propagate.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"expect.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"util.h\"\SE{}}}
\L{\LB{}}
\L{\LB{\K{\#define} TWODIM}}
\L{\LB{\K{\#include} \S{}\"param.h\"\SE{}}}
\L{\LB{}}
\L{\LB{\K{\#ifdef} HAVE\_PROF\_H}}
\L{\LB{\K{\#include} \<prof.h\>}}
\L{\LB{\K{\#else}}}
\L{\LB{\K{\#define} MARK(a)}}
\L{\LB{\K{\#endif}}}
\L{\LB{}}
\L{\LB{\K{static} \K{char} rcsid[] = \S{}\"\$Id: twodim.c,v 1.18 1998\/01\/05 15:59:01 bert Exp \$\"\SE{};}}
\L{\LB{}}
\L{\LB{\C{}\/***** main program *****\/\CE{}}}
\L{\LB{}}
\L{\LB{\K{int} main(\K{int} argc, \K{char} **argv)}}
\L{\LB{\{}}
\L{\LB{  \K{int} completed;}}
\L{\LB{}}
\L{\LB{  \C{}\/* \"real\" and \"complex\" types are defined in complex.h *\/\CE{}}}
\L{\LB{  \K{complex} *packet,*ipacket, *propV2,*propV, *propK;}}
\L{\LB{  \K{real} norm1, norm2;}}
\L{\LB{  \K{real} curr\_max\_sq, max\_max\_sq;}}
\L{\LB{  \K{complex} over;}}
\L{\LB{  \K{char} tmpname[128];}}
\L{\LB{}}
\L{\LB{  \C{}\/* this is for speed calculations *\/\CE{}}}
\L{\LB{  \K{time\_t} t\_start;}}
\L{\LB{}}
\L{\LB{  \C{}\/* read the command line arguments *\/\CE{}}}
\L{\LB{}}
\L{\LB{  read\_arguments(argc, argv);}}
\L{\LB{}}
\L{\LB{  \C{}\/* allocate everything *\/\CE{}}}
\L{\LB{}}
\L{\LB{  packet =}\Tab{16}{(\K{complex}*) safe\_malloc(size*size*\K{sizeof}(\K{complex}),\S{}\"packet\"\SE{});}}
\L{\LB{  propK =}\Tab{16}{(\K{complex}*) safe\_malloc(size*size*\K{sizeof}(\K{complex}),\S{}\"propK\"\SE{});}}
\L{\LB{  propV2 =}\Tab{16}{(\K{complex}*) safe\_malloc(size*size*\K{sizeof}(\K{complex}),\S{}\"propV2\"\SE{});}}
\L{\LB{  propV =}\Tab{16}{(\K{complex}*) safe\_malloc(size*size*\K{sizeof}(\K{complex}),\S{}\"propV\"\SE{});}}
\L{\LB{  ipacket =}\Tab{16}{(\K{complex}*) safe\_malloc(size*size*\K{sizeof}(\K{complex}),\S{}\"ipacket\"\SE{});}}
\L{\LB{}}
\L{\LB{  \C{}\/* allocate and initialize data *\/\CE{}}}
\L{\LB{}}
\L{\LB{  fprintf(stdout, \S{}\"initializing FFT parameters...\"\SE{});  fflush(stdout);}}
\L{\LB{  fft2\_init(size);}}
\L{\LB{  fprintf(stdout, \S{}\" done!\!n\"\SE{});}}
\L{\LB{}}
\L{\LB{  generate\_packet(packet, size,  rx, ry,  kx, ky,  delta, eps, 1.);}}
\L{\LB{  copy\_data (packet, ipacket, size);}}
\L{\LB{  fprintf(stdout,\S{}\"packet initialized (\%ux\%u)\!n\"\SE{}, size,size);}}
\L{\LB{  \K{if} (intermediate \|\| (debug\>2)) \{}}
\L{\LB{    dump\_Re(packet, size, SCRATCH \S{}\"packet0\_r\"\SE{});}}
\L{\LB{    dump\_Im(packet, size, SCRATCH \S{}\"packet0\_i\"\SE{});}}
\L{\LB{  \}}}
\L{\LB{}}
\L{\LB{  generate\_Kprop (propK, size, delta, dt);}}
\L{\LB{  fprintf(stdout,\S{}\"propK initialized\!n\"\SE{});}}
\L{\LB{  \K{if} (debug\>2) \{}}
\L{\LB{    dump\_Re(propK, size, SCRATCH \S{}\"propK\_r\"\SE{});}}
\L{\LB{    dump\_Im(propK, size, SCRATCH \S{}\"propK\_i\"\SE{});}}
\L{\LB{  \}}}
\L{\LB{}}
\L{\LB{  generate\_Vprop (propV2, size, delta, dt\/2);}}
\L{\LB{  fprintf(stdout,\S{}\"propV2 initialized\!n\"\SE{});}}
\L{\LB{  \K{if} (debug\>2) \{}}
\L{\LB{    dump\_Re(propV2, size, SCRATCH \S{}\"propV2\_r\"\SE{});}}
\L{\LB{    dump\_Im(propV2, size, SCRATCH \S{}\"propV2\_i\"\SE{});}}
\L{\LB{  \}}}
\L{\LB{}}
\L{\LB{  generate\_Vprop (propV, size, delta, dt);}}
\L{\LB{  fprintf(stdout,\S{}\"propV initialized\!n\"\SE{});}}
\L{\LB{  \K{if} (debug\>2) \{}}
\L{\LB{    dump\_Re(propV, size, SCRATCH \S{}\"propV\_r\"\SE{});}}
\L{\LB{    dump\_Im(propV, size, SCRATCH \S{}\"propV\_i\"\SE{});}}
\L{\LB{  \}}}
\L{\LB{}}
\L{\LB{  norm1 = sum\_sq(packet,size);}}
\L{\LB{  printf(\S{}\"probability at t=0: \%20.16e\!n\"\SE{}, norm1);}}
\L{\LB{}}
\L{\LB{  \C{}\/* set up calculation of max amplitude, expectation values, and speed *\/\CE{}}}
\L{\LB{}}
\L{\LB{  max\_max\_sq = max\_sq(packet, size);}}
\L{\LB{  printf(\S{}\"initial max \|psi\| = \%20.16g\!n\"\SE{}, sqrt(max\_max\_sq));}}
\L{\LB{\K{\#ifdef} EXPECT\_IN\_TWODIM}}
\L{\LB{  printf(\S{}\"initial \<x\> = \%20.16g\!n\"\SE{}, expect\_X(packet, size, delta));}}
\L{\LB{  printf(\S{}\"initial \<y\> = \%20.16g\!n\"\SE{}, expect\_Y(packet, size, delta));}}
\L{\LB{\K{\#endif} \C{}\/* EXPECT\_IN\_TWODIM *\/\CE{}}}
\L{\LB{}}
\L{\LB{  t\_start = time(NULL);}}
\L{\LB{}}
\L{\LB{  \C{}\/* set propagation counter variables *\/\CE{}}}
\L{\LB{}}
\L{\LB{  \K{if} (intermediate == 0)}}
\L{\LB{    intermediate = total\_iter \- start\_iter;  \C{}\/* ie, calculate all at once *\/\CE{}}}
\L{\LB{}}
\L{\LB{  \C{}\/* do the initial propagation *\/\CE{}}}
\L{\LB{}}
\L{\LB{  completed = start\_iter;}}
\L{\LB{}}
\L{\LB{  \K{if} (start\_iter != 0) \{}}
\L{\LB{    propagate (packet, propV2, propV, propK, size, start\_iter, 0);}}
\L{\LB{}}
\L{\LB{    \C{}\/* intermediate results *\/\CE{}}}
\L{\LB{    sprintf(tmpname, SCRATCH \S{}\"packet\%d\_r\"\SE{}, completed);}}
\L{\LB{    dump\_Re(packet, size, tmpname);}}
\L{\LB{    sprintf(tmpname, SCRATCH \S{}\"packet\%d\_i\"\SE{}, completed);}}
\L{\LB{    dump\_Im(packet, size, tmpname);}}
\L{\LB{}}
\L{\LB{    \C{}\/* parameters for plots *\/\CE{}}}
\L{\LB{    curr\_max\_sq = max\_sq(packet, size);}}
\L{\LB{    \K{if} (max\_max\_sq \< curr\_max\_sq)}}
\L{\LB{      max\_max\_sq = curr\_max\_sq;}}
\L{\LB{    printf(\S{}\"max \|psi\| = \%20.16g\!n\"\SE{}, sqrt(curr\_max\_sq));}}
\L{\LB{\K{\#ifdef} EXPECT\_IN\_TWODIM}}
\L{\LB{    printf(\S{}\"\<x\> = \%20.16g\!n\"\SE{}, expect\_X(packet, size, delta));}}
\L{\LB{    printf(\S{}\"\<y\> = \%20.16g\!n\"\SE{}, expect\_Y(packet, size, delta));}}
\L{\LB{\K{\#endif} \C{}\/* EXPECT\_IN\_TWODIM *\/\CE{}}}
\L{\LB{  \}}}
\L{\LB{}}
\L{\LB{  \C{}\/* do the \"main\" propagation *\/\CE{}}}
\L{\LB{}}
\L{\LB{  \K{while} (completed \< total\_iter) \{}}
\L{\LB{    propagate (packet, propV2, propV, propK, size, intermediate, completed);}}
\L{\LB{    completed += intermediate;}}
\L{\LB{}}
\L{\LB{    \C{}\/* intermediate results *\/\CE{}}}
\L{\LB{    sprintf(tmpname, SCRATCH \S{}\"packet\%d\_r\"\SE{}, completed);}}
\L{\LB{    dump\_Re(packet, size, tmpname);}}
\L{\LB{    sprintf(tmpname, SCRATCH \S{}\"packet\%d\_i\"\SE{}, completed);}}
\L{\LB{    dump\_Im(packet, size, tmpname);}}
\L{\LB{}}
\L{\LB{    \C{}\/* parameters for plots *\/\CE{}}}
\L{\LB{    curr\_max\_sq = max\_sq(packet, size);}}
\L{\LB{    \K{if} (max\_max\_sq \< curr\_max\_sq)}}
\L{\LB{      max\_max\_sq = curr\_max\_sq;}}
\L{\LB{    printf(\S{}\"max \|psi\| = \%20.16g\!n\"\SE{}, sqrt(curr\_max\_sq));}}
\L{\LB{\K{\#ifdef} EXPECT\_IN\_TWODIM}}
\L{\LB{    printf(\S{}\"\<x\> = \%20.16g\!n\"\SE{}, expect\_X(packet, size, delta));}}
\L{\LB{    printf(\S{}\"\<y\> = \%20.16g\!n\"\SE{}, expect\_Y(packet, size, delta));}}
\L{\LB{\K{\#endif} \C{}\/* EXPECT\_IN\_TWODIM *\/\CE{}}}
\L{\LB{  \}}}
\L{\LB{}}
\L{\LB{  fprintf(stdout,\S{}\"... done.\!n\"\SE{});}}
\L{\LB{}}
\L{\LB{  \K{if} (total\_iter \> 0)}}
\L{\LB{    printf(\S{}\"speed: approximately \%f seconds per iteration\!n\"\SE{},}}
\L{\LB{}\Tab{8}{   ((\K{double})(time(NULL)\-t\_start)) \/ total\_iter);}}
\L{\LB{}}
\L{\LB{  \C{}\/* parameters for plots *\/\CE{}}}
\L{\LB{  curr\_max\_sq = max\_sq(packet, size);}}
\L{\LB{  \K{if} (max\_max\_sq \< curr\_max\_sq)}}
\L{\LB{    max\_max\_sq = curr\_max\_sq;}}
\L{\LB{  printf(\S{}\"max max \|psi\| = \%20.16g\!n\"\SE{}, sqrt(max\_max\_sq));}}
\L{\LB{\K{\#ifdef} EXPECT\_IN\_TWODIM}}
\L{\LB{  printf(\S{}\"final \<x\> = \%20.16g\!n\"\SE{}, expect\_X(packet, size, delta));}}
\L{\LB{  printf(\S{}\"final \<y\> = \%20.16g\!n\"\SE{}, expect\_Y(packet, size, delta));}}
\L{\LB{\K{\#endif} \C{}\/* EXPECT\_IN\_TWODIM *\/\CE{}}}
\L{\LB{}}
\L{\LB{  \C{}\/* propagation statistics *\/\CE{}}}
\L{\LB{  norm2 = sum\_sq(packet,size);}}
\L{\LB{  printf(\S{}\"probability at t=T: \%20.16e\!n\"\SE{}, norm2);}}
\L{\LB{  printf(\S{}\"probability conservation error: \%+21.16e\!n\"\SE{},}}
\L{\LB{}\Tab{8}{ (norm2\/norm1)\-1);}}
\L{\LB{}}
\L{\LB{  over = overlap (ipacket, packet, size);}}
\L{\LB{  over \/= norm1;}}
\L{\LB{  printf(\S{}\"overlap at t=\%f:   \%14.10f (\%10.6e + \%10.6e i)\!n\"\SE{},}}
\L{\LB{}\Tab{8}{ dt*total\_iter, Csqabs(over), Re(over),Im(over));}}
\L{\LB{}}
\L{\LB{  fft2\_end();}}
\L{\LB{}}
\L{\LB{  free(packet);}}
\L{\LB{  free(propK);}}
\L{\LB{  free(propV2);}}
\L{\LB{  free(propV);}}
\L{\LB{  \K{return} 0;}}
\L{\LB{\}}}
