%begin; !Tetrahedron demonstration

%endoflist
%include "level1:graphinc.imp"
%include "inc:maths.imp"
%external %routine %spec setframe %alias "FRED_GRAPHICS_FRAME" (%integer addr)
%external %routine %spec setup %alias "FRED_GRAPHICS_SETUP"

%integerfn int(%real r)
  %result = intpt(0.5+r) %unless r<0
  %result = -intpt(0.5-r)
%end

%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,half=0

%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 synctime = 0
%owninteger temp
  synctime=cputime %if synctime=0
  %while cputime<synctime %cycle; %repeat
  colour(black)
  half = 1-half
  half clear(half)
  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
! colour(white)
! line(j1_x, j1_y, j2_x, j2_y)
! line(j2_x, j2_y, j3_x, j3_y)
! line(j3_x, j3_y, j1_x, j1_y)
! line(j1_x, j1_y, j4_x, j4_y)
! line(j2_x, j2_y, j4_x, j4_y)
! line(j3_x, j3_y, j4_x, j4_y)
  colour(red);    triangle(j1_x,j1_y,j4_x,j4_y,j3_x,j3_y) %if area(j1,j4,j3)<0
  colour(green);  triangle(j1_x,j1_y,j2_x,j2_y,j4_x,j4_y) %if area(j1,j2,j4)<0
  colour(yellow); triangle(j2_x,j2_y,j3_x,j3_y,j4_x,j4_y) %if area(j2,j3,j4)<0
  colour(white);  triangle(j1_x,j1_y,j3_x,j3_y,j2_x,j2_y) %if area(j1,j3,j2)<0
  offset(0,-page0)
  synctime = synctime+60
  synctime = cputime+10 %and printsymbol(7) %if synctime<cputime
  j1_y = j1_y-page0
  j2_y = j2_y-page0
  j3_y = j3_y-page0
  j4_y = j4_y-page0
%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

%on %event 0 %start
   print symbol(27)
   print symbol('R')
   %stop
%finish
%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,minus(0.4),minus(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,minus(sin(a)),cos(a))
form r3r3 matrix(peri,ii,jj,kk)
prompt("j:"); a=minus(0.07); read(a)
p3(ii,cos(a),0,minus(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,minus(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)

setframe(16_E00000)
setup
clear

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

%endofprogram