1 /* Quantum Minigolf, a computer game illustrating quantum mechanics
2 Copyright (C) 2007 Friedemann Reinhard <friedemann.reinhard@gmail.com>
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 */
18
19
20 #include "QuantumSimulator.h"
21
22 #include "quantumminigolf.h"
23
24 #define INTENS 120 // color intensity at maximal probability density
25
26 // constructor: setup the FFT engine and compute the Momentum Propagator
QuantumSimulator(int width,int height,float dt)27 QuantumSimulator::QuantumSimulator(int width, int height, float dt)
28 {
29 this->dt = dt;
30 this->width = width;
31 this->height = height;
32
33 psi = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*width*height);
34 prop = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*width*height);
35 xprop = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*width*height);
36
37 printf("Initializing FFT engine ... ");
38 fft = fftwf_plan_dft_2d(width, height,
39 psi, psi, FFTW_FORWARD, FFTW_MEASURE);
40 printf("done \n");
41
42 printf("Initializing inverse FFT engine ... ");
43 ifft = fftwf_plan_dft_2d(width, height,
44 psi, psi, FFTW_BACKWARD, FFTW_MEASURE);
45
46 BuildMomentumPropagator();
47
48 // construct a dummy position propagator. The right propagator
49 // is constructed from the track, once the user has made its choice
50 // where to play.
51 for(int x=0; x<width; x++){
52 for(int y=0; y<height; y++){
53 xprop[x*height+y][0] = 0;
54 xprop[x*height+y][1] = 0;
55 }
56 }
57
58 printf("done \n\n");
59
60 for(int i=0; i<width; i++)
61 for(int j=0; j<height; j++)
62 psi[i*height + j][0] = psi[i*height + j][1] = 0;
63 GaussNorm = 0;
64 }
65
~QuantumSimulator(void)66 QuantumSimulator::~QuantumSimulator(void)
67 {
68 fftwf_free(psi);
69 fftwf_free(prop);
70
71 fftwf_free(xprop);
72 }
73
BuildMomentumPropagator()74 void QuantumSimulator::BuildMomentumPropagator(){
75 int x, y;
76 float yscale = width/height*width/height; // scale factor to compensate for different
77 // k_0 in x and y direction due to different dimensions
78
79 for(x=0; x<width/2; x++){
80 for(y=0; y<height/2; y++){
81 prop[x*height+y][0] = cos(dt*(-x*x - yscale*y*y));
82 prop[x*height+y][1] = sin(dt*(-x*x - yscale*y*y));
83 }
84 for(y=height/2; y<height; y++){
85 prop[x*height+y][0] = cos(dt*(-x*x - yscale*(y-height)*(y-height)));
86 prop[x*height+y][1] = sin(dt*(-x*x - yscale*(y-height)*(y-height)));
87 }
88 }
89 for(x=width/2; x<width; x++){
90 for(y=0; y<height/2; y++){
91 prop[x*height+y][0] = cos(dt*(-(x-width)*(x-width) - yscale*y*y));
92 prop[x*height+y][1] = sin(dt*(-(x-width)*(x-width) - yscale*y*y));
93 }
94 for(y=height/2; y<height; y++){
95 prop[x*height+y][0] = cos(dt*(-(x-width)*(x-width) - yscale*(y-height)*(y-height)));
96 prop[x*height+y][1] = sin(dt*(-(x-width)*(x-width) - yscale*(y-height)*(y-height)));
97 }
98 }
99 }
100
101 // BuildPositionPropagator - build up the position from a bitmap
102 // with color-coded obstacle height
BuildPositionPropagator(SDL_Surface * V)103 void QuantumSimulator::BuildPositionPropagator(SDL_Surface *V){
104 int x, y;
105 Uint32 *V_dat = (Uint32 *)V->pixels;
106
107 // extract the potential
108 for(x=0; x<width; x++){
109 for(y=0; y<height; y++){
110 Uint8 dummy, red;
111
112 SDL_GetRGB(V_dat[x + y*width], V->format, &red, &dummy, &dummy);
113 xprop[x*height+y][0] = cos(-.5*(float)(red)*dt*30000/255);
114 xprop[x*height+y][1] = sin(-.5*(float)(red)*dt*30000/255);
115 if(red>250){
116 xprop[x*height+y][0] = 0;
117 xprop[x*height+y][1] = 0;
118 }
119 }
120 }
121 }
122
123 //PropagateMomentum -- FFT into k-space and apply the momentum propagator
124 // to the wave function
125 // effectively, this propagates the wavefunction by dt in a zero potential
PropagateMomentum(void)126 void QuantumSimulator::PropagateMomentum(void){
127 volatile float tre, tim, pre, pim; // swap register and propagator real and imaginary parts
128 volatile int x, y;
129
130 // propagate in momentum space
131 fftwf_execute(fft);
132
133 for(x=0; x<width; x++){
134 for(y=0; y<height; y++){
135 tre = psi[x*height+y][0];
136 tim = psi[x*height+y][1];
137 pre = prop[x*height+y][0];
138 pim = prop[x*height+y][1];
139
140 psi[x*height+y][0] = tre*pre - tim*pim;
141 psi[x*height+y][1] = tre*pim + tim*pre;
142 }
143 }
144
145 fftwf_execute(ifft);
146 }
147
148
149 //PropagatePosition -- propagate in position space
150 // and scale the wavefunction by a factor of quench
151 // note that this operation is not unitary due to the
152 // hard erase at infinite potentials
153 // return value: the new norm of the propagated wavefunction
PropagatePosition(float quench)154 float QuantumSimulator::PropagatePosition(float quench){
155 float norm = 0;
156 volatile float tre, tim, pre, pim, dnorm, mnorm=0;
157 volatile int x, y;
158
159 for(x=0; x<width; x++){
160 for(y=0; y<height; y++){
161 tre = psi[x*height+y][0];
162 tim = psi[x*height+y][1];
163
164 pre = xprop[x*height+y][0];
165 pim = xprop[x*height+y][1];
166
167 if(pre == 0 && pim == 0){
168 psi[x*height+y][0] = 0;
169 psi[x*height+y][1] = 0;
170 }
171 else{
172 //propagate the wavefunction.
173 // and correct for last time's shrink and the
174 // FFT's scaling
175 psi[x*height+y][0] = quench*(tre*pre - tim*pim);
176 psi[x*height+y][1] = quench*(tim*pre + tre*pim);
177 }
178
179 dnorm = (psi[x*height+y][0]*psi[x*height+y][0] +
180 psi[x*height+y][1]*psi[x*height+y][1]);
181
182 norm += dnorm;
183
184 }
185 }
186
187 norm /= GaussNorm*INTENS*INTENS;
188
189 return norm;
190 }
191
192
193 //PositionMeasurement
194 // performe a position measurement, i.e., randomly pick a point x, y
195 // according to the probability distribution defined by the wavefunction psi
PositionMeasurement(int * x,int * y)196 void QuantumSimulator::PositionMeasurement(int *x, int *y){
197 float norm = 0;
198 float psi2;
199 float cutoff = 1e-6;
200 float criterion;
201 float sucprob = 0;
202
203 int runs = 0;
204
205 int holex = 100, holey = 160;
206 int holer = 30;
207
208 for(int i=0; i<width; i++){
209 for(int j=0; j<height; j++){
210 float psi2 = psi[i*height+j][0]*psi[i*height+j][0] +
211 psi[i*height+j][1]*psi[i*height+j][1];
212 norm += psi2*psi2;
213 }
214 }
215
216 for(int i=holex-holer; i<holex+holer; i++){
217 for(int j=holey-holer; j<holey+holer; j++){
218 float psi2 = psi[i*height+j][0]*psi[i*height+j][0] +
219 psi[i*height+j][1]*psi[i*height+j][1];
220 sucprob += psi2*psi2;
221 }
222 }
223
224 sucprob /= norm;
225
226 do{
227 *x = rand()%width;
228 *y = rand()%height;
229 psi2 = psi[*x*height+*y][0]*psi[*x*height+*y][0] +
230 psi[*x*height+*y][1]*psi[*x*height+*y][1];
231 criterion = ((float)(rand())/RAND_MAX)+cutoff;
232 runs++;
233 }while(psi2*psi2/norm/2 < criterion);
234
235 }
236
237 // GenGauss
238 // generate a coherent state (i.e. a Gaussian wavepacket centered around
239 // cx, cy in position and around kx, ky in momentum space)
GenGauss(int cx,int cy,float kx,float ky,float w)240 void QuantumSimulator::GenGauss(int cx, int cy, float kx, float ky, float w ){
241 int x, y, xeff, yeff;
242 float r;
243
244 // commented out for uncertainty movie 070519
245 GaussNorm=0;
246
247 int xlower = (int)(cx-2.5*w); if(xlower < 0) xlower = 0;
248 int ylower = (int)(cy-2.5*w); if(ylower < 0) ylower = 0;
249 int xupper = (int)(cx+2.5*w); if(xupper > width)xupper =width;
250 int yupper = (int)(cy+2.5*w); if(yupper > height)yupper =height;
251
252 for(x=xlower; x<xupper; x++){
253 xeff = x-cx;
254 for(y=ylower; y<yupper; y++){
255 yeff = y-cy;
256 r=exp(-.25*(xeff*xeff + yeff*yeff)/w/w);
257 psi[height*x+y][0] = r*cos(kx*xeff + ky*yeff);
258 psi[height*x+y][1] = r*sin(kx*xeff + ky*yeff);
259 GaussNorm += psi[height*x+y][0]*psi[height*x+y][0] +
260 psi[height*x+y][1]*psi[height*x+y][1];
261 }
262 }
263
264 for(x=xlower; x<xupper; x++){
265 for(y=ylower; y<yupper; y++){
266 psi[height*x+y][0] *= INTENS;
267 psi[height*x+y][1] *= INTENS;
268 }
269 }
270 }
271
272 //ClearWave - initialize psi with zeros
ClearWave(void)273 void QuantumSimulator::ClearWave(void){
274 for(int x=0; x<width; x++){
275 for(int y=0; y<height; y++){
276 psi[height*x+y][0] = psi[height*x+y][1] = 0;
277 }
278 }
279 }
280