%begin; !Tetrahedron demonstration

%include "inc:frame.imp"
%include "inc:maths.imp"

{%realfn minus(%real r)
{  %result=-r
{%end

%recordformat r3(%real x,y,z)
%recordformat r2(%real x,y)
%recordformat r3r3(%real xx,yx,zx,xy,yy,zy,xz,yz,zz)
%recordformat r3r2(%real xx,yx,xy,yy,xz,yz)
%recordformat point(%integer x,y)

%routine matmult(%record(r3r3)%name f,g,h); ! H(v) = G(F(v))
  h_xx=g_xx*f_xx+g_xy*f_yx+g_xz*f_zx
  h_xy=g_xx*f_xy+g_xy*f_yy+g_xz*f_zy
  h_xz=g_xx*f_xz+g_xy*f_yz+g_xz*f_zz
  h_yx=g_yx*f_xx+g_yy*f_yx+g_yz*f_zx
  h_yy=g_yx*f_xy+g_yy*f_yy+g_yz*f_zy
  h_yz=g_yx*f_xz+g_yy*f_yz+g_yz*f_zz
  h_zx=g_zx*f_xx+g_zy*f_yx+g_zz*f_zx
  h_zy=g_zx*f_xy+g_zy*f_yy+g_zz*f_zy
  h_zz=g_zx*f_xz+g_zy*f_yz+g_zz*f_zz
%end

%routine transform(%record(r3)%name in,out,%record(r3r3)%name f)
! 3d -> 3d transformation       out = f(in)
%real x,y,z
  x = f_xx*in_x + f_xy*in_y + f_xz*in_z
  y = f_yx*in_x + f_yy*in_y + f_yz*in_z
  z = f_zx*in_x + f_zy*in_y + f_zz*in_z
  out_x = x; out_y = y; out_z = z
%end

%routine project(%record(r3)%name in,%record(r2)%name out,
  %record(r3r2)%name f)
! 3d -> 2d transformation   out = f(in)
  out_x = f_xx*in_x + f_xy*in_y + f_xz*in_z
  out_y = f_yx*in_x + f_yy*in_y + f_yz*in_z
%end

! Perspective projection stuff

%record(r3)iris; !NB projection plane is parallel to xy, and 1.0 above iris

%routine perspex(%record(r3)%name p,%record(r2)%name q)
%real dx,dy,s
  dx = p_x-iris_x; dy = p_y-iris_y
  s = iris_z / (iris_z-p_z)
  q_x = s*dy; q_y = s*dx
%end

%routine p3(%record(r3)%name p,%real x,y,z)
  p_x=x; p_y=y; p_z=z
%end

%routine p2(%record(r2)%name p,%real x,y)
  p_x=x; p_y=y
%end

%routine add(%record(r3)%name this,to)
  to_x = to_x + this_x
  to_y = to_y + this_y
  to_z = to_z + this_z
%end

%routine form r3r3 matrix(%record(r3r3)%name m,%record(r3)%name i,j,k)
  m_xx = i_x; m_yx = i_y; m_zx = i_z
  m_xy = j_x; m_yy = j_y; m_zy = j_z
  m_xz = k_x; m_yz = k_y; m_zz = k_z
%end

%routine form r3r2 matrix(%record(r3r2)%name m,%record(r2)%name i,j,k)
  m_xx = i_x; m_yx = i_y
  m_xy = j_x; m_yy = j_y
  m_xz = k_x; m_yz = k_y
%end

%record(r3)c1,c2,c3,c4;    ! Four corners of tetrahedron
%record(r2)i1,i2,i3,i4;    ! Images thereof after projection
%record(point)j1,j2,j3,j4;    ! and after scaling
%record(point)k1=0,k2=0,k3=0,k4=0; !for erasing

%record(r3r3)rot;          ! Rotation matrix
%record(r3r2)proj;         ! Planar projection matrix
%integer page0=0,page1=512

%integerfn area(%record(point)%name a,b,c); !Vector product ABxAC
  %result = (a_y-b_y)*(c_x-a_x)+(b_x-a_x)*(c_y-a_y)
%end

%routine draw it
%owninteger xxl=0,xxr=0,yyb=0,yyt=0
%owninteger xl=0,xr=0,yb=0,yt=0,temp
%owninteger synctime = 0
  %routine note(%record(point)k)
    xl = k_x %if k_x<xl; xr = k_x %if k_x>xr
    yb = k_y %if k_y<yb; yt = k_y %if k_y>yt
  %end
  %routine myfill(%integer l,b,r,t)
    %if l<0 %start
      myfill(l&1023,b,1023,t); l = 0; printsymbol('L'); newline
    %finish
    %if r>1023 %start
      myfill(0,b,r&1023,t); r = 1023; printsymbol('R'); newline
    %finish
    %if b<0 %start
      myfill(l,b&1023,r,1023); b = 0; printsymbol('B'); newline
    %finish
    %if t>1023 %start
      myfill(l,0,r,t&1023); t = 1023; printsymbol('T'); newline
    %finish
    fill(l,b,r,t)
  %end
  %while cputime<synctime %cycle;%repeat
  setcolour(black)
  myfill(xxl,yyb+page1,xxr,yyt+page1)
  xxl = xl; xxr = xr; yyb = yb; yyt = yt
  temp = page1; page1 = page0; page0 = temp
  j1_y = j1_y+page0
  j2_y = j2_y+page0
  j3_y = j3_y+page0
  j4_y = j4_y+page0
  setcolour(red);    triangle(j1,j4,j3) %if area(j1,j4,j3)<0
  setcolour(green);  triangle(j1,j2,j4) %if area(j1,j2,j4)<0
  setcolour(yellow); triangle(j2,j3,j4) %if area(j2,j3,j4)<0
  setcolour(white);  triangle(j1,j3,j2) %if area(j1,j3,j2)<0
  offset(0,-page0)
  synctime = cputime+40
  j1_y = j1_y-page0
  j2_y = j2_y-page0
  j3_y = j3_y-page0
  j4_y = j4_y-page0
  xl = j1_x; xr = xl; yb = j1_y; yt = yb
  note(j2); note(j3); note(j4)
%end

%routine scale(%record(r2)%name i,%record(point)%name j)
%constinteger xf=150,yf=150,xm=350,ym=240
  j_x = int(i_x*xf)+xm
  j_y = int(i_y*yf)+ym
%end

%routine project it
  perspex(c1,i1); scale(i1,j1)
  perspex(c2,i2); scale(i2,j2)
  perspex(c3,i3); scale(i3,j3)
  perspex(c4,i4); scale(i4,j4)
%end

%routine rotate it
  transform(c1,c1,rot)
  transform(c2,c2,rot)
  transform(c3,c3,rot)
  transform(c4,c4,rot)
%end

%routine readvector(%record(r3)%name v)
  skipsymbol %andreturnif nextsymbol=nl
  read(v_x); read(v_y); read(v_z)
  skipsymbol
%end

%record(r3)centre;         ! of tetrahedron

%record(r3)i,j,k;          ! r3 spanning vectors
%record(r3)ii,jj,kk;       ! ditto after rotation
%record(r2)iii,jjj,kkk;    ! and after projection
%record(r3r3)peri,perj,perk,temp

%real r; r = 1.0;          ! length of one edge
%real c; c = 0.85;         ! cos(pi/6.0)
%real cc; cc = c*c

%real s; s = r*c
%real t; t = r*cc
%real hr,hs,ht
%real q,p,a
%integer t1,t2,t3,t4,t5,previous,now

  hr = r/2.0
  hs = s/2.0
  ht = t/2.0

p3(i,1.0,0,0)
p3(j,0,1.0,0)
p3(k,0,0,1.0)

skipsymbol
prompt("Centre of tetrahedron: ")
p3(centre,0,0,0); read vector(centre)

prompt("Iris:")
p3(iris,1.0,1.0,4.0); read vector(iris)

p3(c1,hs,-hr,-ht)
p3(c2,hs,hr,-ht)
p3(c3,-hs,0,-ht)
p3(c4,0,0,ht)

add(centre,c1)
add(centre,c2)
add(centre,c3)
add(centre,c4)

p2(iii,-0.4,-0.4)
p2(jjj,1.0,0)
p2(kkk,0,1.0)
form r3r2 matrix(proj,iii,jjj,kkk)

!! read rotation parameters

prompt("i:"); a=0.03; read(a)
p3(ii,1.0,0,0); p3(jj,0,cos(a),sin(a)); p3(kk,0,-sin(a),cos(a))
form r3r3 matrix(peri,ii,jj,kk)
prompt("j:"); a=-0.07; read(a)
p3(ii,cos(a),0,-sin(a)); p3(jj,0,1.0,0); p3(kk,sin(a),0,cos(a))
form r3r3 matrix(perj,ii,jj,kk)
prompt("k:"); a=0.05; read(a); skipsymbol
p3(ii,cos(a),sin(a),0); p3(jj,-sin(a),cos(a),0); p3(kk,0,0,1.0)
form r3r3 matrix(perk,ii,jj,kk)

matmult(peri,perj,temp)
matmult(temp,perk,rot)

! Draw coordinate axes in both pages

p3(ii,1.0,0,0); p3(jj,0,1.0,0); p3(kk,0,0,1.0)
perspex(ii,iii); perspex(jj,jjj); perspex(kk,kkk)
p3(ii,0,0,0)
perspex(ii,iii)

setup; clear

prompt("")
%cycle
  project it
  draw it
  rotate it
%repeat

%endofprogram