Skip to content

Instantly share code, notes, and snippets.

@thedeemon thedeemon/path_int.d
Last active Sep 26, 2019

Embed
What would you like to do?
Toy quantum mechanics with path integrals
import std.complex, std.stdio, std.random, std.parallelism, std.math : hypot;
import std.getopt, std.datetime.stopwatch, std.array, std.algorithm : max, min, reduce;
import std.format, std.range : iota;
struct Wall { double y, x1, x2; } // x1 < x2
struct Point { double x, y; }
bool intersects(Point p1, Point p2, ref Wall wall) {
if (min(p1.y, p2.y) > wall.y) return false;
if (max(p1.y, p2.y) < wall.y) return false;
if (p1.y == p2.y) return p1.y == wall.y;
double x = p1.x + (p2.x - p1.x) * (wall.y - p1.y) / (p2.y - p1.y);
return x >= wall.x1 && x <= wall.x2;
}
bool canPass(Point p1, Point p2) {
foreach(ref w; walls)
if (intersects(p1, p2, w))
return false;
return true;
}
alias C = Complex!double;
enum MaxStops = 4;
__gshared {
real hc = 20; // 1/h
double Q = 0.01; // charge
uint NPaths = 50000;
}
__gshared walls = [Wall(0.5, 0.0, 0.4), Wall(0.5, 0.41, 0.59), Wall(0.5, 0.6, 1.0)];
C propagator(Point p1, Point p2) {
real len = hypot(p1.x - p2.x, p1.y - p2.y);
return C(expi(hc * len));
}
C amplitude(double endX, double endY) {
Point[MaxStops+1] points = void;
points[0] = Point(0.5, 0.0);
auto endPoint = Point(endX, endY);
C totalAmpl = canPass(points[0], endPoint) ? propagator(points[0], endPoint) : complex(0.0);
paths:
foreach(n; 0..NPaths) {
int npoints = 1 + uniform(0, MaxStops);
foreach(i; 1 .. 1 + npoints)
points[i] = Point(uniform01(), uniform01());
C ampl = complex(1.0, 0.0);
foreach(i; 0..npoints) {
if (!canPass(points[i], points[i+1])) continue paths;
ampl *= Q * propagator(points[i], points[i+1]);
}
if (!canPass(points[npoints], endPoint)) continue paths;
ampl *= propagator(points[npoints], endPoint);
totalAmpl += ampl;
}
return totalAmpl;
}
void main(string[] args) {
auto pic = new int[512*512];
auto prob = new double[512];
bool single = false;
auto opts = getopt(args, "hc", "1/h", &hc, "charge|q", &Q, "single|s", "only a single slit", &single,
"paths|p", "number of paths attempting to reach each end point", &NPaths);
if (opts.helpWanted)
return defaultGetoptPrinter("Toy QM.", opts.options);
if (single)
walls = [Wall(0.5, 0.0, 0.4), Wall(0.5, 0.41, 1.0)];
foreach(y; 0..512) {
writeln(y);
foreach(x; iota(512).parallel) {
C ampl = amplitude(x / 512.0, (512.0-y)/512.0);
prob[x] = sqAbs(ampl);
}
prob[] /= prob.reduce!max;
foreach(x; 0..512) {
int pr = cast(int)(prob[x] * 255);
pic[y*512+x] = pr*257;
}
}
foreach(w; walls) {
int y = 512 - cast(int)(w.y*512);
int x1 = cast(int)(w.x1*512);
int x2 = cast(int)(w.x2*512);
foreach(x; x1..x2)
pic[y*512+x] = 0x0000ff;
}
writeBmp(512, 512, cast(ubyte[])pic, format("paths_hc%03s_p%d_q%s%s.bmp", hc, NPaths, Q, single ? "s" : ""));
}
void writeBmp(int width, int height, ubyte[] data, string filename) {
import std.file : write;
auto app = appender!(byte[]);
void write4(uint what) { app.put( (cast(byte*) &what)[0..4]); }
void write2(ushort what){ app.put( (cast(byte*) &what)[0..2]); }
void write1(byte what) {app.put(what); }
uint fileSize = 54 + height * ((width * 3) + (!((width*3)%4) ? 0 : 4-((width*3)%4)));
write1('B');
write1('M');
write4(fileSize); // size of file in bytes
write4(0); // reserved
write4(54); // offset to the bitmap data
write4(40); // size of BITMAPINFOHEADER
write4(width); // width
write4(height); // height
write2(1); // planes
write2(24); // bpp
foreach(i; 0..6)
write4(0); // compression, size of uncompressed, x pels per meter, colors used, colors important
int offsetStart = cast(int) data.length;
for(int y = height; y > 0; y--) {
offsetStart -= width * 4;
int offset = offsetStart;
int b = 0;
foreach(x; 0 .. width) {
write1(data[offset + 2]); // blue
write1(data[offset + 1]); // green
write1(data[offset + 0]); // red
b += 3;
offset += 4;
}
int w = b%4;
if(w) for(int a = 0; a < 4-w; a++)
write1(0); // pad until divisible by four
}
write(filename, app.data);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.