Toy quantum mechanics with path integrals

path_int.d

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);
}