rayanim.c

#include <stdio.h>
#include <math.h>

#ifndef M_PI
#define M_PI 3.141592653589793238462643383279502884L
#endif

long double
sqrl (long double x)
{
  return x*x;
}

long double
fracl (long double x)
{
  return x - floor(x);
}

static long double sph;

void
ray (long double mh, long double ox, long double oy, long double oz,
     long double ux, long double uy, long double uz,
     long double *r, long double *g, long double *b)
{
  long double h,x,y,z,t,nx,ny,nz,vx,vy,vz;
  static short depth;
  if (depth>=300)
    {
      *r = 0;
      *g = 0;
      *b = 0;
      return;
    }
  depth++;
  h = mh;
  t = 0;
  x = ox;
  y = oy;
  z = oz;
  do
    {
      t = t+h;
      x = x+h*ux;
      y = y+h*uy;
      z = z+h*uz;
      if ( (x*x+y*y+sqrl(z-1-sph)-1)<=0 )
        {
          do
            {
              h = -(x*x+y*y+sqrl(z-1-sph)-1)/(2*x*ux+2*y*uy+2*(z-1-sph)*uz);
              t = t+h;
              x = x+ux*h;
              y = y+uy*h;
              z = z+uz*h;
            }
          while ( ! ( fabsl(x*x+y*y+sqrl(z-1-sph)-1)<=1E-10L ) );
          nx = 2*x;
          ny = 2*y;
          nz = 2*(z-1-sph);
          long double tmp = -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
          vx = 2*tmp*nx+ux;
          vy = 2*tmp*ny+uy;
          vz = 2*tmp*nz+uz;
          ray(mh,x,y,z,vx,vy,vz,r,g,b);
          *r = (*r)*0.9L;
          *g = (*g)*0.9L;
          *b = (*b)*0.9L;
          depth--;
          return;
        }
      if ( (sqrl(x+2)+y*y+sqrl(z-0.5L)-0.25L)<=0 )
        {
          do
            {
              h = -(sqrl(x+2)+y*y+sqrl(z-0.5L)-0.25L)/(2*(x+2)*ux+2*y*uy+2*(z-0.5L)*uz);
              t = t+h;
              x = x+ux*h;
              y = y+uy*h;
              z = z+uz*h;
            }
          while ( ! ( fabsl(sqrl(x+2)+y*y+sqrl(z-0.5L)-0.25L)<=1E-10L ) );
          nx = 2*(x+2);
          ny = 2*y;
          nz = 2*(z-0.5L);
          long double tmp = -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
          vx = 2*tmp*nx+ux;
          vy = 2*tmp*ny+uy;
          vz = 2*tmp*nz+uz;
          ray(mh,x,y,z,vx,vy,vz,r,g,b);
          *r = (*r)*0.9L;
          *g = (*g)*0.5L;
          *b = (*b)*0.5L;
          depth--;
          return;
        }
    if ( (sqrl(x-1)+sqrl(y+sqrtl(3))+sqrl(z-0.5L)-0.25L)<=0 )
      {
        do
          {
            h = -(sqrl(x-1)+sqrl(y+sqrtl(3))+sqrl(z-0.5L)-0.25L)/(2*(x-1)*ux+2*(y+sqrtl(3))*uy+2*(z-0.5L)*uz);
            t = t+h;
            x = x+ux*h;
            y = y+uy*h;
            z = z+uz*h;
          }
        while ( ! ( fabsl(sqrl(x-1)+sqrl(y+sqrtl(3))+sqrl(z-0.5L)-0.25L)<=1E-10L ) );
        nx = 2*(x-1);
        ny = 2*(y+sqrtl(3));
        nz = 2*(z-0.5L);
        long double tmp = -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
        vx = 2*tmp*nx+ux;
        vy = 2*tmp*ny+uy;
        vz = 2*tmp*nz+uz;
        ray(mh,x,y,z,vx,vy,vz,r,g,b);
        *r = (*r)*0.5L;
        *g = (*g)*0.9L;
        *b = (*b)*0.5L;
        depth--;
        return;
      }
    if ( (sqrl(x-1)+sqrl(y-sqrtl(3))+sqrl(z-0.5L)-0.25L)<=0 )
      {
        do
          {
            h = -(sqrl(x-1)+sqrl(y-sqrtl(3))+sqrl(z-0.5L)-0.25L)/(2*(x-1)*ux+2*(y-sqrtl(3))*uy+2*(z-0.5L)*uz);
            t = t+h;
            x = x+ux*h;
            y = y+uy*h;
            z = z+uz*h;
          }
        while ( ! ( fabsl(sqrl(x-1)+sqrl(y-sqrtl(3))+sqrl(z-0.5)-0.25)<=1E-10L ) );
        nx = 2*(x-1);
        ny = 2*(y-sqrtl(3));
        nz = 2*(z-0.5L);
        long double tmp = -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
        vx = 2*tmp*nx+ux;
        vy = 2*tmp*ny+uy;
        vz = 2*tmp*nz+uz;
        ray(mh,x,y,z,vx,vy,vz,r,g,b);
        *r = (*r)*0.5L;
        *g = (*g)*0.5L;
        *b = (*b)*0.9L;
        depth--;
        return;
      }
    if (x*x+y*y+z*z>100)
      {
        h = 0.001L*(x*x+y*y+z*z);
        if (z>0 && uz>=0)
          {
            *r = 0.5L;
            *g = 0.8L;
            *b = 1;
            depth--;
            return;
          }
      }
    if (z<=0)
      {
        do
          {
            h = -z/uz;
            t = t+h;
            x = x+ux*h;
            y = y+uy*h;
            z = z+uz*h;
          }
        while ( ! ( fabsl(z)<=1E-10 ) );
        if ( (fracl(x+100)<0.5L) != (fracl(y+100)<0.5L) )
          {
            *r = 100/(x*x+y*y+100);
            *g = *r;
            *b = *r;
          }
        else
          {
            *r = 0;
            *g = 0;
            *b = 0;
          }
        depth--;
        return;
      }
    }
  while ( 1 );
}

void
anim (unsigned char img, FILE *f)
{
  int scx, scy;
  long double dst,theta,phi,p1,p2,ax,ay,az,bx,by,bz,cx,cy,cz,
    ex,ey,ez,ux,uy,uz,tmp;
  long double r,g,b;
  fprintf (f, "P3\n320 200\n255\n");
  dst = 4;
  theta = M_PI/4.L;
  phi = -5*M_PI/12.L;
  sph = img;
  sph = 1-sqrl((sph-16)/16.L);
  ax = sinl(theta)*cosl(phi);
  ay = sinl(theta)*sinl(phi);
  az = cosl(theta);
  bx = -sinl(phi);
  by = cosl(phi);
  bz = 0;
  cx = -cosl(theta)*cosl(phi);
  cy = -cosl(theta)*sinl(phi);
  cz = sinl(theta);
  ex = ax*dst;
  ey = ay*dst;
  ez = az*dst;
  for ( scy=0 ; scy<=199 ; scy++ )
    {
      for ( scx=0 ; scx<=319 ; scx++ )
        {
          p1 = (scx-160)/120.L;
          p2 = (100-scy)/100.L;
          ux = -ax+p1*bx+p2*cx;
          uy = -ay+p1*by+p2*cy;
          uz = -az+p1*bz+p2*cz;
          tmp = sqrtl(ux*ux+uy*uy+uz*uz);
          ux = ux/tmp;
          uy = uy/tmp;
          uz = uz/tmp;
          ray(0.1L,ex,ey,ez,ux,uy,uz,&r,&g,&b);
          unsigned char rr = floorl(r*6);  if (rr==6) rr=5;
          unsigned char gg = floorl(g*6);  if (gg==6) gg=5;
          unsigned char bb = floorl(b*6);  if (bb==6) bb=5;
          fprintf (f, "  %d %d %d", rr*48, gg*48, bb*48);
        }
      fprintf (f, "\n");
    }
}

int
main (void)
{
  unsigned char img;
  for ( img=0 ; img<=31 ; img++ )
    {
      char s[256];
      sprintf (s, "ima%02d.ppm", img);
      FILE *f = fopen(s, "w");
      anim (img, f);
      fclose (f);
    }
  return 0;
}

rayanim.pas

program rayanimation;

{$N+,E-}
{$M 65520,0,0}

uses
  dos,
  crt;

var
  screen : array[0..63999] of byte absolute $A000:$0000;
  regs : registers;
  pal : array[0..255] of record r,g,b : byte; end;
  f : file;
  s : string;

procedure putdot(x,y : word;c : byte);
begin
  regs.ah := $0C;
  regs.al := c;
  regs.bh := 0;
  regs.cx := x;
  regs.dx := y;
  intr($10,regs);
end;

var
  img : byte;
  i,j,k,l,m,n : longint;
  r,g,b : extended;
  scx,scy : longint;
  sph : extended;
  dst,theta,phi,p1,p2,ax,ay,az,bx,by,bz,cx,cy,cz,
  ex,ey,ez,ux,uy,uz,tmp,tmp2,costh : extended;
  depth : word;

procedure markdot;
var
  rr,gg,bb : byte;
begin
  rr := trunc(r*6); if rr=6 then rr := 5;
  gg := trunc(g*6); if gg=6 then gg := 5;
  bb := trunc(b*6); if bb=6 then bb := 5;
  putdot(scx,scy,rr*36 + gg*6 + bb);
end;

procedure ray(mh,ox,oy,oz,ux,uy,uz : extended;var r,g,b : extended);
var
  h,x,y,z,t,nx,ny,nz,vx,vy,vz : extended;
begin
  if depth>=300 then begin
    r := 0;
    g := 0;
    b := 0;
    exit;
  end;
  inc(depth);
  h := mh;
  t := 0;
  x := ox;
  y := oy;
  z := oz;
  repeat
    t := t+h;
    x := x+h*ux;
    y := y+h*uy;
    z := z+h*uz;
    if (x*x+y*y+sqr(z-1-sph)-1)<=0 then begin
      repeat
        h := -(x*x+y*y+sqr(z-1-sph)-1)/(2*x*ux+2*y*uy+2*(z-1-sph)*uz);
        t := t+h;
        x := x+ux*h;
        y := y+uy*h;
        z := z+uz*h;
      until abs(x*x+y*y+sqr(z-1-sph)-1)<=1E-10;
      nx := 2*x;
      ny := 2*y;
      nz := 2*(z-1-sph);
      tmp := -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
      vx := 2*tmp*nx+ux;
      vy := 2*tmp*ny+uy;
      vz := 2*tmp*nz+uz;
      ray(mh,x,y,z,vx,vy,vz,r,g,b);
      r := r*0.9;
      g := g*0.9;
      b := b*0.9;
      dec(depth);
      exit;
    end;
    if (sqr(x+2)+y*y+sqr(z-0.5)-0.25)<=0 then begin
      repeat
        h := -(sqr(x+2)+y*y+sqr(z-0.5)-0.25)/(2*(x+2)*ux+2*y*uy+2*(z-0.5)*uz);
        t := t+h;
        x := x+ux*h;
        y := y+uy*h;
        z := z+uz*h;
      until abs(sqr(x+2)+y*y+sqr(z-0.5)-0.25)<=1E-10;
      nx := 2*(x+2);
      ny := 2*y;
      nz := 2*(z-0.5);
      tmp := -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
      vx := 2*tmp*nx+ux;
      vy := 2*tmp*ny+uy;
      vz := 2*tmp*nz+uz;
      ray(mh,x,y,z,vx,vy,vz,r,g,b);
      r := r*0.9;
      g := g*0.5;
      b := b*0.5;
      dec(depth);
      exit;
    end;
    if (sqr(x-1)+sqr(y+sqrt(3))+sqr(z-0.5)-0.25)<=0 then begin
      repeat
        h := -(sqr(x-1)+sqr(y+sqrt(3))+sqr(z-0.5)-0.25)/(2*(x-1)*ux+2*(y+sqrt(3))*uy+2*(z-0.5)*uz);
        t := t+h;
        x := x+ux*h;
        y := y+uy*h;
        z := z+uz*h;
      until abs(sqr(x-1)+sqr(y+sqrt(3))+sqr(z-0.5)-0.25)<=1E-10;
      nx := 2*(x-1);
      ny := 2*(y+sqrt(3));
      nz := 2*(z-0.5);
      tmp := -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
      vx := 2*tmp*nx+ux;
      vy := 2*tmp*ny+uy;
      vz := 2*tmp*nz+uz;
      ray(mh,x,y,z,vx,vy,vz,r,g,b);
      r := r*0.5;
      g := g*0.9;
      b := b*0.5;
      dec(depth);
      exit;
    end;
    if (sqr(x-1)+sqr(y-sqrt(3))+sqr(z-0.5)-0.25)<=0 then begin
      repeat
        h := -(sqr(x-1)+sqr(y-sqrt(3))+sqr(z-0.5)-0.25)/(2*(x-1)*ux+2*(y-sqrt(3))*uy+2*(z-0.5)*uz);
        t := t+h;
        x := x+ux*h;
        y := y+uy*h;
        z := z+uz*h;
      until abs(sqr(x-1)+sqr(y-sqrt(3))+sqr(z-0.5)-0.25)<=1E-10;
      nx := 2*(x-1);
      ny := 2*(y-sqrt(3));
      nz := 2*(z-0.5);
      tmp := -(nx*ux+ny*uy+nz*uz)/(nx*nx+ny*ny+nz*nz);
      vx := 2*tmp*nx+ux;
      vy := 2*tmp*ny+uy;
      vz := 2*tmp*nz+uz;
      ray(mh,x,y,z,vx,vy,vz,r,g,b);
      r := r*0.5;
      g := g*0.5;
      b := b*0.9;
      dec(depth);
      exit;
    end;
    if (x*x+y*y+z*z>100) then begin
      h := 0.001*(x*x+y*y+z*z);
      if (z>0) and (uz>=0) then begin
        r := 0.5;
        g := 0.8;
        b := 1;
        dec(depth);
        exit;
      end;
    end;
    if z<=0 then begin
      repeat
        h := -z/uz;
        t := t+h;
        x := x+ux*h;
        y := y+uy*h;
        z := z+uz*h;
      until abs(z)<=1E-10;
      if (frac(x+100)<0.5) xor (frac(y+100)<0.5) then begin
        r := 100/(x*x+y*y+100);
        g := r;
        b := r;
      end else begin
        r := 0;
        g := 0;
        b := 0;
      end;
      dec(depth);
      exit;
    end;
  until false;
end;

procedure animg;
begin
  dst := 4;
  theta := pi/4;
  phi := -5*pi/12;
  sph := img;
  sph := 1-sqr((sph-16)/16);
  ax := sin(theta)*cos(phi);
  ay := sin(theta)*sin(phi);
  az := cos(theta);
  bx := -sin(phi);
  by := cos(phi);
  bz := 0;
  cx := -cos(theta)*cos(phi);
  cy := -cos(theta)*sin(phi);
  cz := sin(theta);
  ex := ax*dst;
  ey := ay*dst;
  ez := az*dst;
  for scy := 199 downto 0 do for scx := 0 to 319 do begin
    p1 := (scx-160)/120;
    p2 := (100-scy)/100;
    ux := -ax+p1*bx+p2*cx;
    uy := -ay+p1*by+p2*cy;
    uz := -az+p1*bz+p2*cz;
    tmp := sqrt(ux*ux+uy*uy+uz*uz);
    ux := ux/tmp;
    uy := uy/tmp;
    uz := uz/tmp;
    ray(0.1,ex,ey,ez,ux,uy,uz,r,g,b);
    markdot;
    if keypressed then halt;
  end;
end;

begin
  depth := 0;
  regs.ax := $0013;
  intr($10,regs);
  for i := 0 to 5 do for j := 0 to 5 do for k := 0 to 5 do begin
    pal[i*36+j*6+k].r := i*12;
    pal[i*36+j*6+k].g := j*12;
    pal[i*36+j*6+k].b := k*12;
  end;
  regs.ax := $1012;
  regs.bx := 0;
  regs.cx := 256;
  regs.es := seg(pal);
  regs.dx := ofs(pal);
  intr($10,regs);
  for img := 0 to 31 do begin
    s := 'IMA'+char(img div 10+48)+char(img mod 10+48)+'.SCR';
    assign(f,s);
    rewrite(f,1);
    animg;
    blockwrite(f,screen,64000);
    close(f);
  end;
end.

See https://twitter.com/gro_tsen/status/1564648994109890561 for context