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

KISS

-Safety, speed, usability, etc. all generally speaking originate from the same source, the system doing as little as possible in the simplest and understandable way possible while achieving the goal of doing a lot of the same thing it's good at.

-Transparency: the system must be transparent and open. It must not do anything behind the users back. The user must be able to view all that is going on in their system and understand it.

-Non-redundancy: in UI there ought to be only one way of doing things so as not to confuse the user or waste their time by making them learn many ways of doing the same thing.

-Limited abstraction: working close to the hardware while maintaining a level of readability and functionality. Unnecessary abstraction leads to waste of time, less speed, worse security and confusion.

-Stay out of the way of the user, never popup anything unless requested, never run anything in the background unless specifically requested.

Visit Russia, find subs, observe Militsiya, get a cold (after two years of health).

WM design ideas:

Applications are allowed to bring themselves on top only once (when they are started within certain time). After that the user must ALT+TAB to bring them on top.

OS/Compiler stuff:

1) Compile a static busybox and copy the files to a HDD image (after making the partition):
mount -o offset=32256 sda.bin tmp
2) mknod a few devices
3) Compile linux kernel

Tools these days...

4) qemu-system-x86_64 -m 512 -kernel bzImage -append "root=/dev/sda1" -enable-kvm -hda sda.bin

5) Modify tcc so it can be compiled into static binary (./configure --extra-cflags=-static), improve the assembler by adding syscall, use statically linked uClibc for startup code (plan: write the startup code entirely from scratch).

x86_64-asm.h: DEF_ASM_OP0(syscall, 0x0f05)
tccpp.c: get rid of ldexp and make your own d = d*((1<<((int)(exp_val-frac_bits))));
tccelf.c: get rid of all things related to libtcc1.a and crt?.o
tcc.c: add missing (dummy) function void __assert_fail (const char *assertion, const char *file, unsigned int line, const char *function) { abort(); }

tcc -DCONFIG_TCC_STATIC -o tcc.o -c tcc.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o libtcc.o -c libtcc.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccpp.o -c tccpp.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccgen.o -c tccgen.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccelf.o -c tccelf.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccasm.o -c tccasm.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o tccrun.o -c tccrun.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o x86_64-gen.o -c x86_64-gen.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o i386-asm.o -c i386-asm.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o libtcc1.o -c libtcc1.c -DTCC_TARGET_X86_64 -I./include -static
tcc -DCONFIG_TCC_STATIC -o alloca86_64.o -c alloca86_64.S -DTCC_TARGET_X86_64 -I./include -static

gcc -nostdlib -static -o tcc libtcc1.o alloca86_64.o tcc.o libtcc.o tccpp.o tccgen.o tccelf.o tccasm.o tccrun.o x86_64-gen.o i386-asm.o ../uClibc-0.9.33/lib/crt1.o ../uClibc-0.9.33/lib/crti.o ../uClibc-0.9.33/lib/crtn.o ../uClibc-0.9.33/lib/libc.a
strip tcc

6) End up with standalone 290 kB static C-compiler/assembler.

Boots in a second.

void* syscall(void* syscall_number, void* param1, void* param2, void* param3) {

  __asm__ __volatile__(
    "mov %0, %%rax\n"
    "mov %1, %%rdi\n"
    "mov %2, %%rsi\n"
    "mov %3, %%rdx\n"
    "syscall\n"
    :
    : "g" (syscall_number), "g" (param1), "g" (param2), "g" (param3)
    : "rax", "rdi", "rsi", "rdx");
}

typedef unsigned long int uintptr;
typedef long int intptr;

static intptr write(int fd, void const* data, uintptr nbytes) {
  syscall((void*)1, (void*)(intptr)fd, (void*)data, (void*)nbytes);
}

int main(int argc, char* argv[]) {
    write(1, "Hello, World!\n", 15);

    return 0;
}

_start() {
  __asm__ __volatile__(
    "call main\n"
    "mov $60, %rax\n"
    "mov 0, %rdi\n"
    "syscall\n");
}

7) Compile hello world on a few MB linux installation on a 64 bit CPU.

Belonging, purpose, transcendence and storytelling...

Navier-Stokes, good luck making the code shorter...

N  = 192; dt = 0.25; p = zeros(N, N);
u = zeros(N, N); v = zeros(N, N); ustar = zeros(N, N); vstar = zeros(N, N);
k = pi*[0:(N/2-1) (-N/2):(-1)]; [x y]  = meshgrid(k, k); ds = -(x.^2 + y.^2);
ds(1,1) = 1; px = 0.0; py = 0.0; gx = 0.0; gy = 0.0;
for t = 0:1000
    % constant velocity zone
    u(N/2-10:N/2+10, N/2-10:N/2+10) = 0.2;

    % solve intermediate velocity U*
    [ux uy] = gradient(u);
    [vx vy] = gradient(v);
    ustar = u + (u.*ux + v.*uy + del2(u)/20 - px/20 + gy)*dt;
    vstar = v + (u.*vx + v.*vy + del2(v)/20 - py/20 + gx)*dt;

    % solve Poisson equation: del2(p) = f
    fhat = fft2(divergence(ustar, vstar));
    p = N^2*real(ifft2(fhat./ds));

    % enforce del(U) = 0
    [px py] = gradient(p);
    u = ustar - px*dt;
    v = vstar - py*dt;

    cla; surface(sqrt(u.^2 + v.^2), 'edgecolor', 'none'); drawnow;
end

Espoo is full of furry animals.

Just a basic wave equation...

u  = zeros(192, 192);
du = zeros(192, 192);
[x y] = meshgrid(linspace(-1, 1, 192), linspace(-1, 1, 192));
for t = 0:0.002:0.4
    u = u + exp(-1000*((x-0.7*sin(7*t)).^2+(y-0.7*sin(11*t)-0.2).^2));
    u = u + 0.5*exp(-1000*((x+0.7*sin(5*t)).^2+(y+0.7).^2))*sin(3*31*t);

    du = du + del2(u);   
    u = u + du;
       
    cla; imagesc(u); drawnow;
end

...