{VERSION 6 0 "IBM INTEL NT" "6.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
267 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 
2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 20 "Multidimensional FFT" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 171 "A test o
f Maple's new discrete Fourier transform package (available as of Mapl
e9). First, we test it in one dimension and demonstrate the effect of \+
two mesh distributions." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT -1 507 "First case: a one-dimensional mesh that includes th
e origin, and is similar in both the coordinate and momentum domain. W
e restrict ourselves to meshes with size in powers of two (as in the o
ld FFT program), even though the algorithm is capable of handling othe
r meshes as well. We will use the interpretation of the Fourier transf
orm as familiar from quantum mechanics, i.e., connecting the coordinat
e and momentum representations. The alternative would be to look at ti
me signals in the frequency domain. " }}{PARA 0 "" 0 "" {TEXT -1 160 "
We do not overemphasize the QM aspect, i.e., we work in units where Pl
anck's constant (hbar) equals unity. Our domains of integration will b
e symmetric both in " }{TEXT 19 1 "x" }{TEXT -1 8 " and in " }{TEXT 
19 1 "p" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT 263 22 "1. One-dimensional FFT" }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "restart; Digits:=15;
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "with(DiscreteTransforms
);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N:=64;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "dx:=0.3;" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 24 "dp:=evalhf(2*Pi/(N*dx));" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 120 "Our choice of parameters resulted in a similar spacing
 in coordinate and momentum. We pick a function to be transformed." }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "fx:=1/(1+x^2)^2;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Z:=Array(1..N,datatype=compl
ex[8]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 327 "In the first 1d imple
mentation we introduce a coordinate mesh which includes the origin. Th
e mesh is arranged in such a way that one starts with zero, goes to th
e largest represented positive value, then continues with the largest \+
negative value (which is one unit further than the largest positive) a
nd returns to the origin. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
23 "for i from 1 to N/2 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 162 "pv[
i]:=evalhf((i-1)*dp); pv[N-i+1]:=evalhf(-i*dp); xv[i]:=evalhf((i-1)*dx
); Z[i]:=eval(fx,x=xv[i]); xv[N-i+1]:=evalhf(-i*dx);  Z[N-i+1]:=eval(f
x,x=xv[N-i+1]): od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "plot
([seq([xv[i],Z[i]] ,i=1..N)],x=-10..10,style=point);" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 75 "Note that our mesh spends a lot of points where
 the function is very small." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 90 "Now we carry out the Fourier transform whereby \+
the result is stored in the original array:" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 33 "FourierTransform(Z,inplace=true);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Z[1],Z[2],Z[N];" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 149 "Note the negligible imaginary parts. The discret
e FT recognized the symmetry of the function. We graph both real and i
maginary parts to confirm this:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 50 "plot([seq([pv[i],Re(Z[i])] ,i=1..N)],style=point);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot([seq([pv[i],Im(Z[i])] ,i=1..N)
],style=point);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 50 "How well does \+
the DFT agree with the known answer?" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 61 "fp:=int(fx*exp(I*p*x),x=-infinity..infinity) assuming
 p,real;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "evalf(limit(fp,
p=0));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "plot(fp,p=-10..10
);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "Z[1];" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 140 "Before we can compare in detail we have \+
to resolve some normalization issues. Neither our integral, nor the Ma
ple DFT contain the customary " }{TEXT 19 12 "1/sqrt(2*Pi)" }{TEXT -1 
32 " factor. Maple's DFT contains a " }{TEXT 19 9 "1/sqrt(N)" }{TEXT 
-1 26 " factor, but is missing a " }{TEXT 19 2 "dx" }{TEXT -1 103 " fa
ctor to form the Riemann sum (which connects the result to the integra
l). Thus, we correct for both:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 37 "f_DFT:=evalf(sqrt(N)*dx): f_DFT*Z[1];" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 60 "P1:=plot([seq([pv[i],Re(f_DFT*Z[i])] ,i=1..N)],s
tyle=point):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "P2:=plot(fp
,p=-10..10,color=blue): plots[display](P1,P2);" }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 27 "The agreement is very good." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 257 13 "Exercise 1.1:" }}{PARA 
0 "" 0 "" {TEXT -1 121 "Explore the question of agreement/disagreement
 in detail: calculate the percentage error and display it for all mome
nta (" }{TEXT 19 1 "p" }{TEXT -1 9 "-values)." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 258 13 "Exercise 1.2:" }}{PARA 
0 "" 0 "" {TEXT -1 48 "Explore the behaviour using different values of
 " }{TEXT 19 6 "N, dx:" }{TEXT -1 46 " can one obtain acceptable resul
ts with small " }{TEXT 19 1 "N" }{TEXT -1 1 "?" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 259 13 "Exercise 1.3:" }}{PARA 
0 "" 0 "" {TEXT -1 25 "Modify the function pair " }{TEXT 19 8 "\{fx, f
p\}" }{TEXT -1 22 " (you can replace the " }{TEXT 19 1 "1" }{TEXT -1 
141 " in the denominator by a constant), or change to a different pair
 (different power of the entire denominator?) and explore Exercises 1 \+
and 2." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 260 
13 "Mini-project:" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 129 "Exp
lore further issues with the FFT: calculate the second derivate of the
 function by going to transform space, multiplying with " }{TEXT 19 3 
"p^2" }{TEXT -1 30 ", and transforming back using " }{TEXT 19 23 "Inve
rseFourierTransform" }{TEXT -1 41 ", and compare with the analytical r
esult." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 70 "Now we try a modification of the mesh that excl
udes the origin in the " }{TEXT 19 1 "x" }{TEXT -1 80 "-domain. This i
s useful for cases when the original function is undefined there." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 20 "restart; Digits:=15:" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 25 "with(DiscreteTransforms);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 6 "N:=64;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Z
:=Array(1..N,datatype=complex[8]);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 8 "dx:=0.3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 
"xM:=Vector(N):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "for i fr
om 1 to N/2 do: xM[i]:=(i-1/2)*dx; xM[N-i+1]:=-xM[i]; od:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "xM[1],xM[N/2],xM[N/2+1],xM[N];" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dpx:=evalf(2*Pi/(N*dx));" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 71 "It will turn out that the transfo
rmed variable domain does contain the " }{TEXT 19 3 "p=0" }{TEXT -1 7 
" point:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "for i from 1 to
 N/2+1 do: pxM[i]:=(i-1)*dpx:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "if
 i>1 then pxM[N-i+2]:=-pxM[i]: fi: od: " }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 41 "pxM[1],pxM[2],pxM[N/2],pxM[N/2+1],pxM[N];" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 5 "a:=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "fx
:=1/(a^2+x^2)^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "fp:=int
(fx*exp(I*p*x),x=-infinity..infinity) assuming p,real;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "for i from 1 to N do: x:=xM[i]:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Z[i]:=fx: od: x:='x':" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "FourierTransform(Z,inplace=true);" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "Z[1],Z[2];" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 27 "Why is this complex-valued?" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 73 "Arrange the density in an array with p ra
nging from negative to positive." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 50 "plot([seq([pxM[i],Re(Z[i])],i=1..N)],style=point);" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot([seq([pxM[i],Im(Z[i])
],i=1..N)],style=point);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 
"f_DFT:=evalf(sqrt(N)*dx);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
60 "P1:=plot([seq([pxM[i],Re(f_DFT*Z[i])],i=1..N)],style=point):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "P2:=plot(fp,p=-10..10,color=
blue): plots[display](P1,P2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 219 
"As could have been guessed, the real part alone does not agree with t
he answer. If we take the absolute value, matters are corrected. Thus,
 for this mesh arrangement we obtain a DFT which contains a wrong comp
lex phase." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "P1:=plot([seq
([pxM[i],abs(f_DFT*Z[i])],i=1..N)],style=point):" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 57 "P2:=plot(fp,p=-10..10,color=blue): plots[disp
lay](P1,P2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 184 "Therefore, it mi
ght be better to work with the mesh that contains the origin, and to u
se a regulator in those instances where the function to be transformed
 has problems at the origin." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT 261 13 "Exercise 1.4:" }}{PARA 0 "" 0 "" {TEXT -1 50 "
Explore the effect of the phase error by reducing " }{TEXT 19 2 "dx" }
{TEXT -1 19 ". You can increase " }{TEXT 19 1 "N" }{TEXT -1 78 " or ke
ep it constant. What happens to the real and imaginary parts of the DF
T?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 265 13 "Ex
ercise 1.5:" }}{PARA 0 "" 0 "" {TEXT -1 132 "An interesting example to
 Fourier analyze is the Gaussian wavepacket which one can control thro
ugh various parameters such as width " }{TEXT 19 2 "Dx" }{TEXT -1 99 "
 (coordinate space parameter corresponds to inverse momentum width par
ameter), and center momentum " }{TEXT 19 2 "p0" }{TEXT -1 163 ". It is
 relatively straightforward to find parameter choices which exceed the
 mesh capabilities, either by exceeding the momentum resolution (very \+
broad packet in " }{TEXT 19 1 "x" }{TEXT -1 29 "-space = narrow struct
ure in " }{TEXT 19 1 "p" }{TEXT -1 23 "-space), or by putting " }
{TEXT 19 2 "p0" }{TEXT -1 70 " at the edge of the covered momentum ran
ge. Explore this on your own. " }{TEXT 19 35 "GWP=Norm*exp(-(x/Dx)^2)*
exp(I*p0*x)" }{TEXT -1 66 ". Find the analytic Fourier image and compa
re with the FFT result." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 262 22 "2. Two-dimensional FFT" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 192 "Now we p
roceed with the two-dimensional case. Comparison with analytical resul
ts will be more complicated. In contrast to the old FFT routines multi
-dimensional FFTs are calculated within the " }{TEXT 19 18 "DiscreteTr
ansforms" }{TEXT -1 9 " package." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 46 "restart; Digits:=15: with(DiscreteTransforms);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "N:=128;" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 40 "Z:=Array(1..N,1..N,datatype=complex[8]);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "dx:=0.3; dy:=dx:" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 62 "The simplest function to choose would be \+
one which factorizes." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "a:=
1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "fx:=1/(a^2+x^2)^3;" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "fpx:=int(fx*exp(I*px*x),x=
-infinity..infinity) assuming px,real;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 5 "b:=2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "fy
:=1/(b^2+y^2)^3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "fpy:=in
t(fy*exp(I*py*y),y=-infinity..infinity) assuming py,real;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 29 "xM:=Vector(N): yM:=Vector(N):" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 104 "for i from 1 to N/2 do: xM[i]:=(i-1/2)*dx;  x
M[N-i+1]:=-xM[i]; yM[i]:=(i-1/2)*dy; yM[N-i+1]:=-yM[i]; od:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "xM[1],xM[N/2],xM[N/2+1],xM[N
],yM[1],yM[N];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 125 "We chose a cub
e, but nothing prevented us from using a parallelepiped as a computati
onal domain. Simply use different dy, dz." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 63 "for i from 1 to N do: x:=xM[i]: for j from 1 to N d
o: y:=yM[j]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "r:=sqrt(x^2+y^2): Z
[i,j]:=fx*fy: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "DA
:=Array(1..N,1..N,1..3,datatype=float[8]):" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 119 "#for i from 1 to N do: x:=xM[i]: for j from 1 to N
 do: y:=yM[j]: DA[i,j,1]:=x: DA[i,j,2]:=y: DA[i,j,3]:=Z[i,j]; od: od:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 128 "for i from 1 to N do: \+
x:=xM[i]: for j from 1 to N do: y:=yM[j]: DA[i,j,1]:=x: DA[i,j,2]:=y: \+
DA[i,j,3]:=log(abs(Z[i,j])); od: od:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 4 "The " }{TEXT 19 8 "surfdata" }{TEXT -1 85 " procedure is terribl
e at plotting the original surface. That's why we graph the log." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plots[surfdata](DA,axes=boxe
d,style=patchcontour,shading=zhue);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 33 "FourierTransform(Z,inplace=true);" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 7 "Z[1,2];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
123 "Why is this complex-valued? The one-dimensional results imply tha
t the mesh distribution without an origin in the original " }{TEXT 19 
1 "r" }{TEXT -1 22 "-domain are to  blame." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 24 "dpx:=evalf(2*Pi/(N*dx));" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 24 "dpy:=evalf(2*Pi/(N*dy)):" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 63 "for i from 1 to N/2+1 do: pxM[i]:=(i-1)*dpx: pyM[i]
:=(i-1)*dpy:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "if i>1 then pxM[N-i
+2]:=-pxM[i]: pyM[N-i+2]:=-pyM[i]: fi: od: " }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 41 "pxM[1],pxM[2],pxM[N/2],pxM[N/2+1],pxM[N];" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 190 "At one point we were so dismayed \+
by the bad surface plot that we exported the data, and used Surfer 7 t
o graph them. The grid generation and contour plot in that program is \+
much superior to " }{TEXT 19 8 "surfdata" }{TEXT -1 76 ". We left the \+
commented commands to show how to export data via ASCII files." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "# fd := fopen(\"C:/marko/map
le9/testsrf.dat\",WRITE);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
44 "DA:=Array(1..N,1..N,1..3,datatype=float[8]):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 129 "for i from 1 to N do: px:=pxM[i]: for j from 1 to N \+
do: py:=pyM[j]: DA[i,j,1]:=px: DA[i,j,2]:=py: DA[i,j,3]:=abs(Z[i,j]); \+
od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "px:='px': py:='p
y':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "# fprintf(fd,\"%g %g
 %g\\n\",px,py,DA[i,j,3]); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
123 "for i from 1 to N do: for j from 1 to N do: if DA[i,j,3]>0 then D
A[i,j,3]:=log10(DA[i,j,3]): else DA[i,j,3]:=0: fi: od: od:" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 181 "We cheat here: the zero values that occu
r at the largest momentum value are transferred to the origin, where w
e enter the known value. Didn't know how else to get rid off the value
." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "for i from 1 to N do:
 for j from 1 to N do: if DA[i,j,3]=0 then DA[i,j,1]:=0: DA[i,j,2]:=0;
 DA[i,j,3]:=DA[1,1,3]; fi: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 13 "# fclose(fd);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 63 "plots[surfdata](DA,axes=boxed,style=patchcontour,shading=zhue)
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "Why the surface looks differ
ent viewed from the top versus the bottom is for the " }{TEXT 19 8 "su
rfdata" }{TEXT -1 25 " programmers to explain. " }}{PARA 0 "" 0 "" 
{TEXT -1 102 "There are cases when the graph looks bad from the top, b
ut has OK contours when looked at from below. " }{TEXT 19 11 "FT(exp(-
r))" }{TEXT -1 15 " is an example." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 139 "How do we know the calculated 2dimension
al FT is correct? We would have to use our known FT pair and plug in t
he correct conversion factor." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 37 "f_DFT2:=evalf(sqrt(N)*dx*sqrt(N)*dy);" }}}{EXCHG {PARA 0 "> " 0 
"" {MPLTEXT 1 0 14 "Z[1,1]*f_DFT2;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 31 "evalf(eval(fpx*fpy,\{x=0,y=0\}));" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 100 "plot3d(log10(abs(fpx*fpy/f_DFT2)),px=-10..
10,py=-10..10,axes=boxed,style=patchcontour,shading=zhue);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 38 "This looks OK what about the details?." }
}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 17 "n_x:=15; n_y:=10;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "
DA[n_x,n_y];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "evalf(eval(
log10(abs(fpx*fpy/f_DFT2)),\{px=pxM[n_x],py=pyM[n_y]\}));" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 59 "The agreement at intermediate momentum va
lues is excellent!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT 264 13 "Exercise 2.1:" }}{PARA 0 "" 0 "" {TEXT -1 95 "Find other
 FT function pairs (ideally non-factorized ones) and observe the accur
acy of the DFT." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT 266 24 "3. Three-dimensional FFT" }{TEXT 
-1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "restart; Digits:=
15: with(DiscreteTransforms);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 6 "N:=64;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "Z:=Array(1..
N,1..N,1..N,datatype=complex[8]);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 24 "dx:=0.3; dy:=dx: dz:=dy:" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 44 "xM:=Vector(N): yM:=Vector(N): zM:=Vector(N):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 143 "for i from 1 to N/2 do: xM[
i]:=(i-1/2)*dx;  xM[N-i+1]:=-xM[i]; yM[i]:=(i-1/2)*dy; yM[N-i+1]:=-yM[
i]; zM[i]:=(i-1/2)*dz;  zM[N-i+1]:=-zM[i]; od:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 36 "xM[1],xM[N],yM[1],yM[N],zM[1],zM[N];" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 118 "We chose a cube, but nothing prev
ented us from using a parallelepiped as a computational domain. Simply
 use different " }{TEXT 19 6 "dy, dz" }{TEXT -1 1 "." }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 5 "a:=1;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 95 "for i from 1 to N do: x:=xM[i]: for j from 1 to N do:
 y:=yM[j]: for k from 1 to N do: z:=zM[k]:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 54 "r:=sqrt(x^2+y^2+z^2): Z[i,j,k]:=exp(-a*r): od: od: od
:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "FourierTransform(Z,inp
lace=true);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "Z[1,1,2];" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Why is this complex-valued?" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "dpx:=evalf(2*Pi/(N*dx));" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "dpy:=evalf(2*Pi/(N*dy)): dp
z:=evalf(2*Pi/(N*dz)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "f
or i from 1 to N/2+1 do: pxM[i]:=(i-1)*dpx: pyM[i]:=(i-1)*dpy: pzM[i]:
=(i-1)*dpz:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "if i>1 then pxM[N-i+
2]:=-pxM[i]: pyM[N-i+2]:=-pyM[i]: pzM[N-i+2]:=-pzM[i]: fi: od: " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "pxM[1],pxM[2],pxM[N/2],pxM[N
/2+1],pxM[N];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "Now look at a cu
t. Arrange the density in an array with p ranging from negative to pos
itive." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "DA:=Array(1..N,1.
.N,1..3,datatype=float[8]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 198 "for
 i from 1 to N do: px:=pxM[i]: for j from 1 to N do: py:=pyM[j]: DA[i,
j,1]:=px: DA[i,j,2]:=py: DA[i,j,3]:=0:  for k from 1 to N do: pz:=pzM[
k]: DA[i,j,3]:=DA[i,j,3]+abs(Z[i,j,k])*dpz; od: od: od:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 123 "for i from 1 to N do: for j from 1
 to N do: if DA[i,j,3]>0 then DA[i,j,3]:=log10(DA[i,j,3]): else DA[i,j
,3]:=0: fi: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 141 "for
 i from 1 to N do: for j from 1 to N do: if DA[i,j,3]=0 then DA[i,j,3]
:=DA[1,1,3]: DA[i,j,1]:=DA[1,1,1]: DA[i,j,2]:=DA[1,1,2]: fi: od: od:" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plots[surfdata](DA,axes=b
oxed,style=patchcontour,shading=zhue);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 65 "This is a case that looks better from the bottom than fro
m above." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "f_DFT3:=evalf(s
qrt(N)*dx*sqrt(N)*dy*sqrt(N)*dz);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
170 "How can we verify the validity? We need to calculate the analytic
 FT for our function. This is done in spherical polar coordinates. The
 azimuthal angular integral yields " }{TEXT 19 4 "2*Pi" }{TEXT -1 87 "
. The polar angle integration is straightforward. The radial integral \+
is done in Maple:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "r:='r'
: fp:=simplify(4*Pi/p*int(exp(-a*r)*sin(p*r)*r,r=0..infinity)) assumin
g p>0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 126 "We have a spherically \+
symmetric answer. We should be able to verify the FT result by choosin
g a line along a beam starting at " }{TEXT 19 3 "p=0" }{TEXT -1 1 "." 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "ResN:=Vector(N): ResA:=Ve
ctor(N):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 151 "for i from 1 t
o N/2 do: px:=pxM[i]: py:=pyM[i]: pz:=pzM[i]: p:=sqrt(px^2+py^2+pz^2):
 ResA[i]:=[p,p*evalf(fp)]: ResN[i]:=[p,p*abs(f_DFT3*Z[i,i,i])]: od:" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "ResA[2],ResN[2];" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "P1:=plot([seq(ResN[i],i=1..N
/2)],style=point):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "P2:=p
lot([seq(ResA[i],i=1..N/2)],color=blue): plots[display](P1,P2);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "The FFT works remarkably well give
n that we have not resolved the maximum in the " }{TEXT 19 6 "p*f(p)" 
}{TEXT -1 64 " distribution, i.e., we are working on a relatively coar
se mesh." }}{PARA 0 "" 0 "" {TEXT -1 127 "Keep in mind that we used th
e absolute value of the transform. A mesh including the origin should \+
not suffer from this problem." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT 267 13 "Exercise 2.2:" }}{PARA 0 "" 0 "" {TEXT 
-1 21 "Change the parameter " }{TEXT 19 1 "a" }{TEXT -1 213 " in the w
avefunction and observe the change in the momentum distribution. Expla
in what it means physically. Find other function pairs (separable in C
artesians, and non-separable) and verify how well the FFT works." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 269 13 "Exercise
 2.3:" }}{PARA 0 "" 0 "" {TEXT -1 149 "For QM-literate students: consi
der the calculation of the kinetic energy expectation value in momentu
m space. You can calculate it by multiplying in " }{TEXT 19 1 "p" }
{TEXT -1 40 "-space with the multiplicative operator " }{TEXT 19 5 "p^
2/2" }{TEXT -1 15 " (in our units " }{TEXT 19 8 "m=hbar=1" }{TEXT -1 
81 ") and forming the average (and normalizing properly). This will em
phasize larger " }{TEXT 19 1 "p" }{TEXT -1 87 "-values. To check how a
ccurate this calculation will be you can probe the agreement of " }
{TEXT 19 8 "p^2*f(p)" }{TEXT -1 77 " in the above plot. What do you fi
nd? How can you improve upon the situation?" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT 270 13 "Exercise 2.4:" }}{PARA 0 "" 0 
"" {TEXT -1 88 "Change the mesh in the above code to include the origi
n also in the configuration space " }{TEXT 19 7 "(x,y,z)" }{TEXT -1 
95 " domain. Demonstrate that the problems with the incorrect complex \+
phase disappear in this case." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT 268 13 "Mini-project:" }}{PARA 0 "" 0 "" {TEXT 
-1 118 "Demonstrate the usefulness of the 3d FFT by making other cuts \+
for comparison between numerical and analytical results." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "146" 0 }{VIEWOPTS 1 1 
0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }