\File{oppack.c},{10:58},{Jan  5 1998}
\L{\LB{\C{}\/* \$Source: \/mit\/bert\/Cvsroot\/Thesis\/src\/exact\/oppack.c,v \$}}
\L{\LB{ *}}
\L{\LB{ * generation of packets and operators}}
\L{\LB{ *\/\CE{}}}
\L{\LB{}}
\L{\LB{\K{\#include} \<unistd.h\>}}
\L{\LB{}}
\L{\LB{\K{\#include} \S{}\"complex.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"oppack.h\"\SE{}}}
\L{\LB{\K{\#include} \S{}\"param.h\"\SE{}}}
\L{\LB{}}
\L{\LB{\K{static} \K{char} rcsid[] = \S{}\"\$Id: oppack.c,v 1.17 1998\/01\/05 15:58:56 bert Exp \$\"\SE{};}}
\L{\LB{}}
\L{\LB{\C{}\/* NOTE: The code in this file is executed only a few times during program}}
\L{\LB{ *   execution.  As a result, sanity and readability has often been}}
\L{\LB{ *   stressed over efficiency.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{}}
\L{\LB{\C{}\/** syntactic sugar **\/\CE{}}}
\L{\LB{}}
\L{\LB{\K{\#define} square(a)   ((a)*(a))}}
\L{\LB{\K{\#define} cube(a)     ((a)*(a)*(a))}}
\L{\LB{}}
\L{\LB{\C{}\/** discretization stuff **\/\CE{}}}
\L{\LB{}}
\L{\LB{\C{}\/* Turn a discretized energy\-operator (V or K) matrix M into a propagator}}
\L{\LB{ *   coefficient matrix given by U[i,j] = exp(\-i M[i,j] dt\/\!hbar).}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   M \-\-     complex array containing (real) values for the operator.}}
\L{\LB{ *              The value will be modified in place.}}
\L{\LB{ *   n \-\-     dimension of matrix.  (must be a power of 2 because of FFT.)}}
\L{\LB{ *   dt \-\-    discrete time step.}}
\L{\LB{ * Last one isn\'t type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{void} make\_propagator (\K{complex} *M, \K{long} n, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  \K{long} i;}}
\L{\LB{}}
\L{\LB{  \C{}\/* we do this for all rows and all columns, so we don\'t need two loops *\/\CE{}}}
\L{\LB{  \K{for} (i=0; i\<(n*n); i++) \{}}
\L{\LB{    M[i] = Cexp(\-I * ((\K{real}) dt \/ hbar) * M[i]);}}
\L{\LB{  \}}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/* Generate a matrix of values of a function at regularly spaced points}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   M \-\-     complex array for return matrix.}}
\L{\LB{ *              (must be already allocated and be at least n\^2 big.)}}
\L{\LB{ *              The return will only contain real values.}}
\L{\LB{ *   n \-\-     dimension of matrix.  (must be a power of 2 because of FFT.)}}
\L{\LB{ *   delta \-\- discrete space step.}}
\L{\LB{ *   tran \-\-  coordinate translation (or NULL for no translation)}}
\L{\LB{ *   F \-\-     potential function (takes X and Y args, returns V)}}
\L{\LB{ * Last two aren\'t of type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{void} sampled\_matrix (\K{complex} *M, \K{long} n, \K{double} delta, }}
\L{\LB{}\Tab{16}{     translation *tran, xy\_function F)}}
\L{\LB{\{}}
\L{\LB{  \K{long} i,j;}}
\L{\LB{  \K{static} translation no\_translation = \{ 0., 0. \};}}
\L{\LB{}}
\L{\LB{  \K{if} (tran == NULL)}}
\L{\LB{    tran = \&no\_translation;}}
\L{\LB{}}
\L{\LB{  \K{for} (i = \-(n\/2) ; i \< (n\/2); i++) \{}}
\L{\LB{    \K{for} (j = \-(n\/2) ; j \< (n\/2); j++) \{}}
\L{\LB{      \C{}\/* adding n and taking modulo n leaves positive indices alone,}}
\L{\LB{       * but maps negative ones to the correct array index.}}
\L{\LB{       *\/\CE{}}}
\L{\LB{      M[n * ((i+n)\% n) + ((j+n)\% n)] = F(i*delta + tran\-\>x, j*delta + tran\-\>y);}}
\L{\LB{    \}}}
\L{\LB{  \}}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/** potential energy functions **\/\CE{}}}
\L{\LB{}}
\L{\LB{\C{}\/* Return the value of Delfi potential at specified coordinates.}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   x, y \-\-  coordinates}}
\L{\LB{ * They aren\'t of type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{real} V\_delfi (\K{double} x, \K{double} y)}}
\L{\LB{\{}}
\L{\LB{  \K{return}  square(square(x)+square(y))\/4. + mu*cube(x)*(y\-1)}}
\L{\LB{          \- y*(square(x)+square(y)) + (square(x)+3*square(y))\/2.;}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/* Return the value of Nelson potential at specified coordinates.}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   x, y \-\-  coordinates}}
\L{\LB{ * They aren\'t of type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{real} V\_nelson (\K{double} x, \K{double} y)}}
\L{\LB{\{}}
\L{\LB{  \K{return}  vcoef * ( mu * square(x)\/2 + square(y \- square(x)\/2) );}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/* Return the value of harmonic potential at specified coordinates.}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   x, y \-\-  coordinates}}
\L{\LB{ * They aren\'t of type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{real} V\_harmonic (\K{double} x, \K{double} y)}}
\L{\LB{\{}}
\L{\LB{  \K{return}  vcoef * ( square(x) + square(y \/ mu) );}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/* Return the value of a constant potential at specified coordinates.}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   x, y \-\-  coordinates}}
\L{\LB{ * They aren\'t of type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{real} V\_constant (\K{double} x, \K{double} y)}}
\L{\LB{\{}}
\L{\LB{  \K{return}  vcoef;}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/** kinetic energy function **\/\CE{}}}
\L{\LB{}}
\L{\LB{\C{}\/* Return the value of kinetic energy for specified momenta}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   px, py \-\- momenta (NOT wave numbers!)}}
\L{\LB{ * They aren\'t of type \"real\" because of compiler bugs wrt long double.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{real} K\_p2 (\K{double} px, \K{double} py)}}
\L{\LB{\{}}
\L{\LB{  \K{return}  ( square(px) + square(py) ) \/ (2 * mass);}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/* Return the value of kinetic energy at specified points.}}
\L{\LB{ *}}
\L{\LB{ * This expression was derived taking discretization into account, so the}}
\L{\LB{ * simplest units to use are the ones that make values of momentum coordi\-}}
\L{\LB{ * nates cover the interval [\-pi\/2,pi\/2], rather than any \"sane\" units.}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   ax, ay \-\- momentum coordinates scaled to cover [\-pi\/2,pi\/2]}}
\L{\LB{ * They aren\'t of type \"real\" because of compiler bugs wrt long double.  *\/\CE{}}}
\L{\LB{\K{static} \K{real} K\_sin2 (\K{double} ax, \K{double} ay)}}
\L{\LB{\{}}
\L{\LB{  \K{return} ( square(sin(ax)) + square(sin(ay)) );}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/** convenience wrappers **\/\CE{}}}
\L{\LB{}}
\L{\LB{\K{void} generate\_Vprop\_delfi (\K{complex} *V, \K{long} n, \K{double} delta, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  sampled\_matrix(V, n, delta, \&translate, V\_delfi);}}
\L{\LB{  make\_propagator(V, n, dt);}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\K{void} generate\_Vprop\_nelson (\K{complex} *V, \K{long} n, \K{double} delta, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  sampled\_matrix(V, n, delta, \&translate, V\_nelson);}}
\L{\LB{  make\_propagator(V, n, dt);}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\K{void} generate\_Vprop\_harmonic (\K{complex} *V, \K{long} n, \K{double} delta, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  sampled\_matrix(V, n, delta, \&translate, V\_harmonic);}}
\L{\LB{  make\_propagator(V, n, dt);}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\K{void} generate\_Vprop\_constant (\K{complex} *V, \K{long} n, \K{double} delta, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  sampled\_matrix(V, n, delta, \&translate, V\_constant);}}
\L{\LB{  make\_propagator(V, n, dt);}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\K{void} generate\_Kprop\_p2 (\K{complex} *K, \K{long} n, \K{double} delta, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  \C{}\/* if spacing of x;y is delta, spacing of kx;ky is (2 pi)\/(N delta),}}
\L{\LB{   * and spacing of px;py is hbar times that.}}
\L{\LB{   *\/\CE{}}}
\L{\LB{  sampled\_matrix(K, n, hbar*(2*M\_PI)\/(n*delta), NULL, K\_p2);}}
\L{\LB{  make\_propagator(K, n, dt);}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/* This one is a bit more messy, since it explicitly takes discretization}}
\L{\LB{ * into account rather than just dealing with a scaled coordinate system.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{void} generate\_Kprop\_sin2 (\K{complex} *K, \K{long} n, \K{double} delta, \K{double} dt)}}
\L{\LB{\{}}
\L{\LB{  \C{}\/* units used map the momenta onto the interval [\-pi\/2,pi\/2] ... *\/\CE{}}}
\L{\LB{  sampled\_matrix(K, n, M\_PI\/n, NULL, K\_sin2);}}
\L{\LB{  \C{}\/* ... and we also need to scale things in different ways *\/\CE{}}}
\L{\LB{  make\_propagator(K, n, dt * 2 * square(hbar) \/ square(delta));}}
\L{\LB{\}}}
\L{\LB{}}
\L{\LB{\C{}\/** wave packet **\/\CE{}}}
\L{\LB{}}
\L{\LB{\C{}\/* generate a Gaussian wave packet in coordinate form.}}
\L{\LB{ * We could use sampled\_matrix() to do this, but it takes a zillion parameters}}
\L{\LB{ *}}
\L{\LB{ * arguments:}}
\L{\LB{ *   pack \-\-    complex array for return}}
\L{\LB{ *                (must be already allocated and be at least n\^2 big.)}}
\L{\LB{ *   n \-\-       dimension of matrix.  (must be a power of 2 because of FFT.)}}
\L{\LB{ *   rx,ry \-\-   components of initial avg position.}}
\L{\LB{ *   kx,ky \-\-   components of initial avg wave vector.}}
\L{\LB{ *   delta \-\-   discrete space step.}}
\L{\LB{ *   eps \-\-     half\-width of Gaussian (in coordinate form).}}
\L{\LB{ *   norm \-\-    normalization constant.}}
\L{\LB{ * Last six aren\'t long double because of compiler bugs.}}
\L{\LB{ *\/\CE{}}}
\L{\LB{\K{void} generate\_packet (\K{complex} *pack, \K{long} n,}}
\L{\LB{}\Tab{16}{      \K{double} rx, \K{double} ry,  \K{double} kx, \K{double} ky,}}
\L{\LB{}\Tab{16}{      \K{double} delta, \K{double} eps, \K{double} norm)}}
\L{\LB{\{}}
\L{\LB{  \K{long} i,j;}}
\L{\LB{  \K{complex} ikr, r2;}}
\L{\LB{  \K{real} ncoef;}}
\L{\LB{}}
\L{\LB{  \C{}\/* this makes the total normalization of the packet be == norm *\/\CE{}}}
\L{\LB{  ncoef = norm * sqrt(0.5 \/ M\_PI) * delta \/ eps;}}
\L{\LB{}}
\L{\LB{  \K{for} (i = \-(n\/2) ; i \< (n\/2); i++) \{}}
\L{\LB{    \K{for} (j = \-(n\/2) ; j \< (n\/2); j++) \{}}
\L{\LB{      \C{}\/* \!psi = ncoef * exp(i*k0 \!cdot r) * exp(\-1\/4*(r\-r0)\^2\/\!epsilon\^2) *\/\CE{}}}
\L{\LB{      ikr = delta*(i*kx + j*ky)*I;}}
\L{\LB{      r2 = \-.25*(square(delta*i\-rx) + square(delta*j\-ry))\/square(eps);}}
\L{\LB{      \C{}\/* adding n and taking modulo n leaves positive indices alone,}}
\L{\LB{       * but maps negative ones to the correct array index.}}
\L{\LB{       *\/\CE{}}}
\L{\LB{      pack[n * ((i+n) \% n) + ((j+n) \% n)] = ncoef * Cexp(ikr + r2);}}
\L{\LB{    \}}}
\L{\LB{  \}}}
\L{\LB{\}}}
