Schrödinger sims, raytracing and Julia sets

Particle colliding with an obstacle.

N = 256;

[x y] = meshgrid(linspace(-1, 1, N), linspace(-1, 1, N));
psi = exp(-10*((x+0.5).^2 + (y+0.5).^2)).*exp(-j*(-50*x -50*y));
psi = psi/norm(psi);
V = zeros(N, N);
V(N/2-20:N/2-10, N/2:N/2+20) = 1;
V(N/2+10:N/2+20, N/2-10:N/2+10) = 1;
S = conv2(diff([0 0 1 0 0], 2), eye(N)); % del2 matrix form
for n = 1:N % H = V - del2
    Ex{n} = expm(-j*(diag(V(n, :)) - S(:, 2:end-1)));
    Ey{n} = expm(-j*(diag(V(:, n)) - S(:, 2:end-1)));
end
for t = 0:150 % time evolution psi(t+1) = psi(t)*exp(-j*H)
    for n = 1:N
        psi(:, n) = Ex{n}*psi(:, n);
    end
    for n = 1:N
        psi(n, :) = psi(n, :)*Ey{n};
    end
    A = abs(psi).^2/4.5e-4;
    cla; imagesc(A); caxis([0 1]); truesize; drawnow;
    imwrite(imresize(A, 0.5), ['png/' num2str(t+1, '%04.f') '.png']);
    % ffmpeg -r 30 -f image2 -i %04d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p test.mp4
    % convert -resize 50% -layers Optimize -delay 3 *.png test.gif
end

Three particles scattering, simplified approximate solution only.

N = 64;
[x y] = meshgrid(linspace(-1, 1, N), linspace(-1, 1, N));
A = 1.4; B = -0.2; C = 1.2;
psi{1} = exp(-50*((x-0.5*A-B).^2 + (y-0.4*A).^2)).*exp(-j*( 20*x*C + 10*y*C));
psi{2} = exp(-50*((x+0.2*A-B).^2 + (y+0.3*A).^2)).*exp(-j*(-20*x*C - 10*y*C));
psi{3} = exp(-50*((x-0.4*A-B).^2 + (y+0.3*A).^2)).*exp(-j*( 10*x*C - 20*y*C));
v = fft2(0.3./(x.^2+y.^2+0.01));
for x = 1:3
    psi{x} = psi{x}/norm(psi{x});
end
in = 1;
for t = 0:15000
    if mod(t, 150)==0
        B = 0;
        for x = 1:3
            A = abs(psi{x}).^2;
            B = B + A/max(max(A));
        end
        figure(1); cla; imagesc(imresize(B, 4)); truesize; drawnow;
        imwrite(B, ['png/' num2str(in, '%04.f') '.png']);
        in = in + 1
    end
    for x = 1:3 % fast 2D convolution with fft
        V{x} = fftshift(ifft2(v.*fft2(conj(psi{x}).*psi{x})));
    end
    for x = 1:3 % separate wave functions
        W = zeros(N, N);
        for y = 1:3
            if y~=x
                W = W + V{y};
            end
        end
        psi{x} = psi{x} + 1e-3*j*(4*del2(psi{x}) - W.*psi{x});
    end   
end

Traces ever pixel of the screen to z-direction until ray within the object and
then chooses color based on surface normal and distance.
The edges to infinity are funny due to "cutting corners" in the way the normal is computed.

/* gcc test.c -lGL -lSDL */
#include <GL/gl.h>
#include <SDL/SDL.h>

char *d =
  "float t = gl_Color.x*1000.;"
  "vec3 scene(vec2 p) {"
  "  float z = 5.;"
  "  mat3 X = mat3(1., 0., 0., 0., cos(2.*t), sin(2.*t), 0., -sin(2.*t), cos(2.*t));"
  "  mat3 Y = mat3(cos(3.*t), 0., -sin(3.*t), 0., 1., 0., sin(3.*t), 0., cos(3.*t));"
  "  vec3 c = X*Y*vec3(1.5*cos(.5*5.*t), 1.5*sin(.5*7.*t), 7.);"
  "  while(z<10.) {"
  "    vec3 r = X*Y*(vec3(sqrt(z)*p.x, sqrt(z)*p.y, z));"
  "    if(abs(r.x-c.x)<.8 && abs(r.y-c.y)<.8 && abs(r.z-c.z)<.8) return r;"
  "    if(distance(r, c)<1.) return r;"
  "    z += .05;"
  "  }"
  "  return vec3(10.);"
  "}"
  "void main() {"
  "  gl_FragColor = vec4(0.);"
  "  vec2 p = vec2(gl_FragCoord.x/256.-1., gl_FragCoord.y/256.-1.);"
  "  vec3 q = scene(p);"
  "  if(length(q)<10.) {"
  "    vec3 dx = scene(vec2(p.x + .05, p.y)) - scene(vec2(p.x - .05, p.y));"
  "    vec3 dy = scene(vec2(p.x, p.y + .05)) - scene(vec2(p.x, p.y - .05));"
  "    float a = 1.5*dot(normalize(cross(dx, dy)), normalize(q))/sqrt(length(q));"
  "    gl_FragColor = vec4(a);"
  "  }"
  "}";

main(c) {
  SDL_Event e;
  SDL_SetVideoMode(512, 512, 0, SDL_OPENGL);
  glUseProgram(glCreateShaderProgramv(GL_FRAGMENT_SHADER, 1, &d));

/*
  c = 0;
  GLubyte **pixels = malloc(3*512*512);
  int cur;
  char filename[80];

*/

  while(e.type!=SDL_QUIT) {
    glColor3us(c++, 0, 0);
    glRecti(-1, -1, 1, 1);
    SDL_GL_SwapBuffers();
    SDL_PollEvent(&e);
/*
    sprintf(filename, "%04d.ppm", c);
    FILE *f = fopen(filename, "w");
    fprintf(f, "P3\n%d %d\n%d\n", 512, 512, 255);
    *pixels = realloc(*pixels, 3*512*512);
    glReadPixels(0, 0, 512, 512, GL_RGB, GL_UNSIGNED_BYTE, *pixels);
    for(int i=0; i<512; i++) {
        for(int j=0; j<512; j++) {
            cur = 3*((512-i-1)*512+j);
            fprintf(f, "%3d %3d %3d ", (*pixels)[cur], (*pixels)[cur+1], (*pixels)[cur+2]);
        }
        fprintf(f, "\n");

    }

*/

  }
}

Zoom into the Julia set (one of them).

N = 256; B = 0.75;

X = linspace(-1, 1, N);
Y = linspace(-1, 1, N);
while true
    B = B*sqrt(sqrt(0.75))
    y0 = 0.636;
    x0 = 0.107;
    X = linspace(x0-B, x0+B, N);
    Y = linspace(y0-B, y0+B, N);

    for y = 1:N

        for x = 1:N

            % z = 0.0;                  % Mandelbrot

            % z0 = X(x) + Y(y)*j;       % Mandelbrot

            z = X(x) + Y(y)*j;          % Julia

            z0 = 0.5 + 0.5*j;           % Julia

            iteration = 0;

            while abs(z) < 2 && iteration < 128

                z = z^2 + z0;

                iteration = iteration + 1;

            end

            s(x, y) = log(iteration);

        end

    end

    imagesc(s); drawnow;
end

#include <stdio.h>
#include <math.h>
main() {
  char filename[80];
  int N = 128, x, y, c=1;
  float f[128*128];  float p[128*128];
  float u[128*128];  float v[128*128];
  float ux[128*128]; float vx[128*128];
  float uy[128*128]; float vy[128*128];
  float u2[128*128]; float v2[128*128];
  float us[128*128]; float vs[128*128];
  float px[128*128]; float py[128*128];
  for(x=0; x<128*128; x++) {
    f[x] = 0.0;    p[x] = 0.0;
    u[x] = 0.0;    v[x] = 0.0;
    ux[x] = 0.0;   vx[x] = 0.0;
    uy[x] = 0.0;   vy[x] = 0.0;
    u2[x] = 0.0;   v2[x] = 0.0;
    us[x] = 0.0;   vs[x] = 0.0;
    px[x] = 0.0;   py[x] = 0.0;  }
  for(int t=0; t<4500; t++) {
    for(y=-5; y<5; y++)
      for(x=-5; x<5; x++) {
        u[(int)(0.75*N+x)+(N/2+y)*N] = 0.15;
        v[(int)(0.75*N+x)+(N/2+y)*N] = 0.0;
      }
    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        ux[x+y*N] = 0.50*(u[(x+1)+y*N] -                u[(x-1)+y*N]);
        vx[x+y*N] = 0.50*(v[(x+1)+y*N] -                v[(x-1)+y*N]);
        uy[x+y*N] = 0.50*(u[x+(y+1)*N] -                u[x+(y-1)*N]);
        vy[x+y*N] = 0.50*(v[x+(y+1)*N] -                v[x+(y-1)*N]);
        u2[x+y*N] = 0.25*(u[(x+1)+y*N] - 2.0*u[x+y*N] + u[(x-1)+y*N]);
        v2[x+y*N] = 0.25*(v[(x+1)+y*N] - 2.0*v[x+y*N] + v[(x-1)+y*N]);
        u2[x+y*N]+= 0.25*(u[x+(y+1)*N] - 2.0*u[x+y*N] + u[x+(y-1)*N]);
        v2[x+y*N]+= 0.25*(v[x+(y+1)*N] - 2.0*v[x+y*N] + v[x+(y-1)*N]);
      }
    for(x=0; x<N*N; x++) {
      us[x] = u[x] + u[x]*ux[x] + v[x]*uy[x] + u2[x]/20.0 - px[x]/20.0 + 0.00000;
      vs[x] = v[x] + u[x]*vx[x] + v[x]*vy[x] + v2[x]/20.0 - py[x]/20.0 - 0.00005;
    }
    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        f[x+y*N] = 0.5*(us[(x+1)+y*N] - us[(x-1)+y*N] + vs[x+(y+1)*N] - vs[x+(y-1)*N]);
        p[x+y*N] = 0.25*(p[(x+1)+y*N] + p[x+(y+1)*N] + p[(x-1)+y*N] + p[x+(y-1)*N] - f[x+y*N]);
    }
    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        px[x+y*N] = 0.5*(p[(x+1)+y*N] - p[(x-1)+y*N]);
        py[x+y*N] = 0.5*(p[x+(y+1)*N] - p[x+(y-1)*N]);
         u[x+y*N] = us[x+y*N] - px[x+y*N];
         v[x+y*N] = vs[x+y*N] - py[x+y*N];
      }
    if(t%40==0) {
      printf("%d\n", t);
      sprintf(filename, "ppm/%04d.pgm", c);
      FILE *f = fopen(filename, "w");
      fprintf(f, "P2\n%d %d\n%d\n", 128, 128, 255);
      c++;
      for(int j=0; j<128; j++) {
          for(int i=0; i<128; i++) {
             float val = 1500.0*sqrt( u[i+j*N]*u[i+j*N] + v[i+j*N]*v[i+j*N]);
             if(val>255.0) val=255.0;
             fprintf(f, "%3d ", (unsigned char)val);
          }
          fprintf(f, "\n");
      }
      close(f);
    }
  }
}