PROGRAM centralprojection;

{Set date in constant docdate}

{
The routines in this file are centered around the very simple
procedure  

     cproj(x);  (the argument x is a 3-dimensional vector x)

This procedure performs a central projection of the vector x from the
vector CENTER to the plane through the origin perpendicular to the vector
CENTER. The coordinates in the image plane are chosen so that the z-axis is
mapped to the second coordinate axis.  The procedure writes to the file

   cproj.tex

the two coordinates of the image point and, as a TeX comment, the three
coordinates of the original point.  The vector CENTER, and a few
corresponding constants entering in the central projection map, are set
by 
		   setcproj(c:vector);

   In addition, this file contains several routines for manipulating
with vectors.  The main part of this program is an example writing a
complete TeX file cproj.tex.  To tex the file, you need the spline
package: spline.sty.

}


CONST  pi=4*arctan(1);{3.1415926}
       epsilon=0.01;
       docdate='3 January 2005';

TYPE  vindex=1..3;
      vector=array[vindex] of REAL;

VAR   count, j, k:integer;
      x, x1, x2, NORMAL, TRANSVECTOR, CENTER : vector;
      r, l, r1, l1, 
      HEIGHT, ANGLE, cosPHI, sinPHI, cosTHETA, sinTHETA: real; 
      outputf: text;



FUNCTION vec(x1,x2,x3:real): vector; {return vector with 3 given coords}
begin
   vec[1]:=x1; vec[2]:=x2; vec[3]:=x3;
end{vec};


FUNCTION innerprod(x,y:vector):real; {return inner product of x and y}
var i:vindex;
begin
   innerprod:=0;
   for i:=1 to 3 do innerprod:=innerprod+x[i]*y[i];
end{innerprod};


FUNCTION veclength(x:vector): real; {return length of vector x}
begin
   veclength:=sqrt(innerprod(x,x));
end{veclength};


FUNCTION crossprod(x,y:vector):vector; {cross product x \times y}
begin
   crossprod[1]:=x[2]*y[3]-x[3]*y[2];
   crossprod[2]:=x[3]*y[1]-x[1]*y[3];
   crossprod[3]:=x[1]*y[2]-x[2]*y[1];
end{crossprod};


PROCEDURE  setnormal(c: vector); {sets the unitvector NORMAL}
var i:vindex; tmp: real;
begin
   tmp:=veclength(c);
   if tmp<epsilon then begin writeln('Error 4'); HALT; end; 
   for i:=1 to 3 do NORMAL[i]:=c[i]/tmp;
end{setnormal};


PROCEDURE rotate(var x:vector); {Uses unit vector NORMAL and angle ANGLE}
var i:vindex; ip: real; xhat: vector; 
begin
   ip:=innerprod(x,NORMAL);
   xhat:=crossprod(NORMAL,x);
   for i:=1 to 3 do 
       x[i]:=ip*NORMAL[i]
	   +cos(ANGLE)*(x[i]-ip*NORMAL[i])
	   +sin(ANGLE)*xhat[i];
end{rotate};


PROCEDURE vrotate(n:integer; var x:vector); {rotation 2pi/n around z-axis}
var x1:real;
begin
   x1:=x[1]; 
   x[1]:=cos(2*pi/n)*x1-sin(2*pi/n)*x[2];
   x[2]:=sin(2*pi/n)*x1+cos(2*pi/n)*x[2];
end{vrotate};


PROCEDURE rotateonefifth(var x:vector); {rotation around z-axis}
begin
   vrotate(5,x); 
end{rotateonefifth};



PROCEDURE rotateonefourth(var x:vector); {rotation around z-axis}
begin
   vrotate(4,x); 
end{rotateonefourth};


PROCEDURE settrans(t:vector); {set translation vector}
begin
  TRANSVECTOR:=t;
end{settrans};


FUNCTION vecplus(x:vector): vector; {return x + x/|x| }
var i:vindex; tmp: real;
begin
   tmp:=veclength(x);
   if tmp<epsilon 
	then begin writeln('Error 3'); HALT; end;
   tmp:=1+1/tmp;
   for i:=1 to 3 do vecplus[i]:=tmp*x[i];
end;


PROCEDURE writecoords(x:vector); {write 3 coordinates of vector x}
begin
   write(outputf,'(',x[1]:8:4,' , ',x[2]:8:4,' , ',x[3]:8:4,')');
end{writecoords};


{ Procedure CPROJ: 
  Project from c=center onto perpendicular plane with orthonormal
basis for the plane given by 
    f_2 \sim \bar e_3.

The projection is determined by  
     x \mapsto \bar x =  \lambda x + (1-\lambda)c,  
where  

   lambda=|c|^2 / (|c|^2-<x,c>)= h/(h-<x,n>).   h=|c|, n=c/h. 

Clearly,
    f_1  \sim   e_3 x c    = (-c_2,      c_1,           0)  
    f_2  \sim   c x  f_1   = (-c_1c3, -c_2c3, c_1^2+c_2^2)

Or, with  c_0=(c_1,c_2,0), h_0=|c_0|, n_0=c_0/h_0, 
    f_1  \sim  \hat c_0              
    f_2  \sim  -c_3c_0 + h_0^2 e_3   

    |f_1| = h_0,
    |f_2|^2 =c_3^2 |h_0|^2+h_0^4= h_0^2 h^2.

So, f_1  =  \hat n_0,
    f_2  =  -(c_3/h)n_0 + (h_0/h) e_3. 
and
    w_1=<x,f_1> = lambda <x,\hat n_0>
    w_2=<x,f_2> = -lambda (c_3/h)<x,n_0> + lambda (h_0/h)x_3


Better in polar coordinates: If
    c/h = cosTHETA e(PHI)+sinTHETA e_3,
then 
    lambda=h/ (h-x_1 cosTHETA cosPHI - x_2 cosTHETA sinPHI + x_3 sinTHETA)
    n_0 = e(PHI)=(cosPHI, sinPHI,0), and 
    h_0=h cosTHETA, c_3/h=sinTHETA
and
    w_1=lambda ( -x_1 sinPHI+x_2 cosPHI)
    w_2=lambda ( (-x_1 cosPHI -x_2 sinPHI)sinTHETA + x_3 cosTHETA )
}

PROCEDURE setcproj(c:vector); {set center and various parameters}
var height0: real;
begin 
   CENTER:=c;
   height0:=sqrt(c[1]*c[1]+c[2]*c[2]); 
   if height0<epsilon 
           then begin writeln('Error 5');HALT; end;
   cosPHI:=c[1]/height0; 
   sinPHI:=c[2]/height0;
   HEIGHT:=veclength(c);
   cosTHETA:=height0/height; 
   sinTHETA:=c[3]/height0;
end{setcproj};


PROCEDURE cproj(x:vector); {project x, and write coords to tex file}
                        {USES: HEIGHT, and COSines and SINes set by SETCPOJ,}
                        {and TRANSVECTOR set by translate}
var i:vindex; w1, w2, tmp, lambda: real;
begin
  for i:=1 to 3 do x[i]:=x[i]+TRANSVECTOR[i];
  tmp:= height - (x[1]*cosPHI+x[2]*sinPHI)*cosTHETA - x[3]*sinTHETA;
  if tmp < 1 
	then begin writeln('Error 6'); HALT; end;
  lambda:=height/tmp;
  w1:=lambda*(-x[1]*sinPHI+x[2]*cosPHI);
  w2:=lambda*(-(x[1]*cosPHI+x[2]*sinPHI)*sinTHETA + x[3]*cosTHETA );
  write(outputf, ' ', w1:8:4, ' ', w2:8:4, ' % ');
  for i:=1 to 3 do x[i]:=x[i]-TRANSVECTOR[i];
  writecoords(x);
  writeln(outputf);
end{cproj};


PROCEDURE w(t:string); {write string t to tex file}
begin
 writeln(outputf,t);
end;


PROCEDURE begintex; {Initialization, open tex file, start writing}
var i:vindex;
begin
  ANGLE:=0;
  for i:=1 to 3 do TRANSVECTOR[i]:=0;
  assign(outputf,'cproj.tex');
  rewrite(outputf);
  w('\input spline.sty');
  w('\PSunit=1cm');
end;


PROCEDURE endtex; {Close tex file}
begin
  w('\bye');
  close(outputf);
end;


PROCEDURE PSpicture; {begin a PSpicture}
begin
  w(' %');w(' %Picture generated by 3dspline.pas:');
  w('\PSpicture{15\PSunit}{15\PSunit}');
  w('\llcoords{-2.5 -13.5}');
end;


PROCEDURE endPSpicture; {end a PSpicutre}
begin
  w('\endPSpicture');w(' %');
end;


PROCEDURE beginps(t:string); {start a PS-command with a left brace}
begin
  writeln(outputf,'\',t,'{');
end;


PROCEDURE endps; {closing brace for arguments of PS-command}
begin
  w('}');
end;


PROCEDURE spline;
begin
  beginps('spline');
end;


PROCEDURE cspline;
begin
  beginps('cspline');
end;


PROCEDURE PSmark;
begin
  beginps('PSmark');
end;


PROCEDURE PSline;
begin
  beginps('PSline');
end;


PROCEDURE PScline;
begin
  beginps('PScline');
end;


PROCEDURE thicklines;
begin
  w('\setPSlinewidth{0.8pt}');
end;


PROCEDURE thinlines;
begin
  w('\setPSlinewidth{0.2pt}');
end;


PROCEDURE dashedlines;
begin
  w('\setdash{5\PSlinewidth}');
end;


PROCEDURE normallines;
begin
    w('\setPSlinewidth{0.5pt}\nodashes');
end;


PROCEDURE PScircle(p:vector; r:real); {Draw (image of) circle with }
                  {center p, radius r, in a plane perpendicular to p;}
                  {Uses CENTER set by setcproj!}
var i:vindex; j:integer; ip, len:real; punit,qhat,q:vector;
begin
   len:=veclength(p);
   if len > epsilon then begin
	for i:=1 to 3 do punit[i]:=p[i]/len;
	ip:=innerprod(CENTER,punit);
	len:=sqrt(HEIGHT*HEIGHT-ip*ip);
	if len>epsilon then 
		for i:=1 to 3 do q[i]:=r*(ip*punit[i]-CENTER[i])/len;
   end;
   if len <= epsilon then begin 
		for i:=1 to 3 do punit[i]:=CENTER[i]/HEIGHT;
		q:=crossprod(punit,vec(0,0,1));
		len:=veclength(q);
		for i:=1 to 3 do q[i]:=r*q[i]/len;
   end;
{rotate: 60}
w(' %:');
w(' %Circle:');
   beginps('cspline');
   for j:=1 to 6 do begin
	cproj(q);
	ip:=innerprod(q,punit);
	qhat:=crossprod(punit,q);
	for i:=1 to 3 do q[i]:=ip*punit[i]
			+cos(2*pi/6)*(q[i]-ip*punit[i])
			+sin(2*pi/6)*qhat[i];
   end;
   endps;
end{PScircle};

PROCEDURE circlethrough(x:vector); {Draw circle through x,}
                                   {perpendicular to NORMAL}
var j:integer;
begin
  ANGLE:=2*pi/6;
  w(' %');w(' %Circlethrough:'); 
  beginps('cspline');
	for j:=1 to 6 do begin cproj(x);rotate(x); end;
  endps;
end{circlethrough};


PROCEDURE arcthrough(x:vector); {Draw arc, of size ANGLE, through x,}
                                   {perpendicular to NORMAL}
const N=4;
var i:vindex; j:integer; ip,len,saveangle:real; tangent, y: vector;
begin
  w(' %');w(' %Circlethrough:'); 
  saveangle:=ANGLE;
  ANGLE:=ANGLE/N;
  ip:=innerprod(x,NORMAL);
  for i:=1 to 3 do y[i]:=x[i]-ip*NORMAL[i];
  r:=veclength(y);
  beginps('dspline');
       tangent:=crossprod(NORMAL,x);
       len:=veclength(tangent);
       if len<epsilon then begin writeln('Error 9');HALT;end;
       for i:=1 to 3 do y[i]:=x[i]-1/3*angle/N*tangent[i]/len; 
    cproj(y);
    cproj(x);
    for j:=1 to N do begin rotate(x); 
    cproj(x); end;
       tangent:=crossprod(NORMAL,x);
       {veclength(tangent); should be unchanged}
       if abs(veclength(tangent)-len)>epsilon then
                           begin writeln('Error 10');HALT;end;
       for i:=1 to 3 do y[i]:=x[i]+1/3*angle/N*tangent[i]/len; 
    cproj(y);
  endps;
  ANGLE:=saveangle;
end{arcthrough};


PROCEDURE axis(v1,v2,v3:real); {Draw axis}
var v: vector;
begin
   w(' %');w(' %Axis:');
   v:=vec(v1,v2,v3);
   beginps('PSline'); 
      cproj(v); 
      cproj(vecplus(v));
   endps;
   v:=vec(-v1,-v2,-v3);
   beginps('PSline'); 
      cproj(v); 
      cproj(vecplus(v));
   endps;
   beginps('PSmark');
      cproj(vec(v1,v2,v3));cproj(vec(-v1,-v2,-v3));
   endps;
end{axis};


PROCEDURE hexahedron; {Draw hexahedron}
var i: integer; x1,x2: vector;
begin
   w(' %');w(' %Hexahedron:');
   x1:=vec(1,1,1); x2:=vec(1,-1,-1);
   for i:=1 to 4 do begin
	beginps('PSline');
	cproj(x2);rotateonefourth(x2); cproj(x2); cproj(x1);
	rotateonefourth(x1);cproj(x1);
	endps;
   end;
end{hexahedron};


PROCEDURE tetra1; {Draw tetrahedron (inscribed in hexahedron)}
begin
   w(' %');w(' %Tetra1:');
   beginps('PSline');
     cproj(vec(1,1,1));   
     cproj(vec(-1, 1,-1)); 
     cproj(vec(-1,-1, 1)); 
     cproj(vec(1,1,1)); 
     cproj(vec( 1,-1,-1)); 
     cproj(vec(-1, 1,-1));
   endps;
   beginps('PSline'); cproj(vec(-1,-1, 1)); cproj(vec( 1,-1,-1));
   endps;
end{tetra1};


{procedure tetra2; {Draw tetrahedron (inscribed in hexahedron)}
begin
   w(' %');w(' %Tetra2:');
   beginps('PSline');
     cproj(vec( 1,-1, 1)); 
     cproj(vec(-1, 1, 1)); 
     cproj(vec(-1,-1,-1)); 
     cproj(vec( 1,-1, 1)); 
     cproj(vec( 1, 1,-1)); 
     cproj(vec(-1,-1,-1));
   endps;
   beginps('PSline'); cproj(vec(-1, 1, 1)); cproj(vec( 1, 1,-1));
   endps;
end{tetra2};
}


PROCEDURE tetra2; {Second possibility: Draw tetrahedron (inscribed in hexahedron)}
var j:integer; p,x:vector;
begin
   w(' %');w(' %Tetra2:');
   p:=vec(1,-1,1);
   setnormal(p); ANGLE:=2*pi/3;
   x:=vec(-1, 1, 1); 
   for j:=1 to 3 do begin
	beginps('PSline');cproj(p);cproj(x);rotate(x);cproj(x);endps;
   end;
end{tetra2};


PROCEDURE octahedron; {Draw octahedron}
var count:integer; x: vector;
begin
   w(' %');w(' %Octahedron:');
   x:=vec(1,0,0);
   for count:=1 to 4 do begin
	beginps('PSline');
          cproj(vec(0,0,-1));cproj(x);rotateonefourth(x);
          cproj(x);cproj(vec(0,0,1));
	endps;
   end;
end{octahedron};


PROCEDURE icosahedron; {Draw icosahedron}
var count:integer; x,x1,x2: vector;
begin
   w(' %');w(' %Icosahedron:');
   x:=vec(0,0,1);
   x1:=vec(2/sqrt(5), 0, 1/sqrt(5));
   x2:=vec(2/sqrt(5)*cos(2*pi/10), 2/sqrt(5)*sin(2*pi/10), -1/sqrt(5));
   for count:=1 to 5 do begin
      beginps('PSline'); 
         cproj(x2);cproj(x1);
         rotateonefifth(x1);cproj(x1);
         cproj(x); rotateonefifth(x2);
      endps;
   end;
   x:=vec(0,0,-1);
   x1:=vec(-2/sqrt(5), 0, -1/sqrt(5));
   x2:=vec(-2/sqrt(5)*cos(2*pi/10), -2/sqrt(5)*sin(2*pi/10), 1/sqrt(5));
   for count:=1 to 5 do begin
     beginps('PSline'); 
       cproj(x2);cproj(x1); rotateonefifth(x1);
       cproj(x1); cproj(x); rotateonefifth(x2);
     endps;
   end;
end{icosahedron};



PROCEDURE dodecahedron; {Draw dodecahedron}
const  r=1/sqrt(3)/sin(2*pi/10);
       l=sqrt(1-r*r);
       r1=1/sqrt(3)*(sqrt(5)-1)/(2*sin(2*pi/10)); 
       l1=sqrt(1-r1*r1); 
var count:integer; x,x1,x2: vector; 
begin
   w(' %');w(' %Dedecahedron:');
   x1:=vec(r1,0,l1); x:=vec(r,0,l); 
   x2:=vec(r*cos(2*pi/10), r*sin(2*pi/10), -l);
   for count:=1 to 5 do begin
	beginps('PSline'); 
            cproj(x2);cproj(x);cproj(x1);
	    rotateonefifth(x1);cproj(x1);
	    rotateonefifth(x);rotateonefifth(x2);
	endps;
   end;
   x:=vec(-r,0,-l); x1:=vec(-r1,0,-l1); 
   x2:=vec(-r*cos(2*pi/10), -r*sin(2*pi/10), l);
   for count:=1 to 5 do begin
	beginps('PSline'); 
	    cproj(x2);cproj(x);cproj(x1);
	    rotateonefifth(x1);cproj(x1);
	    rotateonefifth(x);
	    rotateonefifth(x2);
	endps;
   end;
end{dodecahedron};


PROCEDURE sphere(r:real); {Draw sphere, radius r, with sevaral small circles}
begin
   w(' %');w(' %Sphere:');
   if r>HEIGHT-epsilon then begin writeln('Error 7');HALT end;
   beginps('PScircle');
       writeln(outputf, '0 0 ', r/sqrt(1-(r*r)/(HEIGHT*HEIGHT)):8:4); 
   endps;
   setnormal(vec(0,0,1)); 
      circlethrough(vec(3/sqrt(10),0, -1/sqrt(10)));
      circlethrough(vec(3/sqrt(10),0,  1/sqrt(10)));
   setnormal(vec(0,1,0));
      circlethrough(vec(1,0,0));
   setnormal(vec(-sin(2*pi/10), cos(2*pi/10), 0)); 
      circlethrough(vec(cos(2*pi/10), sin(2*pi/10),0));
end{sphere};


PROCEDURE spiral; {Draw spiral}
var k:integer; x2: real; x:vector; 
begin
   w(' %');w(' %Spiral:');
   beginps('dspline');
   cproj(vec(-2-1/(3*2*pi),-1/3,1));
   x:=vec(-2,0,1);
   for k:=1 to 24 do begin
     cproj(x); x[1]:=x[1]+0.25; x2:=x[2]; x[2]:=x[3]; x[3]:=-x2;
   end;
   cproj(x);
   cproj(vec(4+1/(3*2*pi),1/3,1));
   endps;
end{spiral};


FUNCTION improj(x:vector):vector; {Center for improj is (0,10,0)}
var trans, tmp: real; 
begin
   trans:=5;
   if x[2]> 10-epsilon then begin writeln('Error 8');HALT end;
   tmp:=(trans+10)/(10-x[2]);
   improj[1]:=tmp*x[1];improj[2]:=-trans; improj[3]:=tmp*x[3];
end{improj};


PROCEDURE image; {Draw a picture with a central projection}
var count: integer; 
begin
   setnormal(vec(0,1,5));
   x1:=vec(1, NORMAL[2]/sqrt(2), NORMAL[3]/sqrt(2));
   x2:=vec(1, -NORMAL[2]/sqrt(2), -NORMAL[3]/sqrt(2));
   angle:=-2*pi/10; rotate(x1); rotate(x2);
   angle:=2*pi/4;
   for count:=1 to 4 do begin
      w('\setdash{5\PSlinewidth}');
      beginps('PSline');cproj(x1);cproj(improj(x1));endps;
      beginps('PSline');cproj(x2);cproj(improj(x2));endps;
      w('\nodashes');
      beginps('PSline');cproj(x2);cproj(x1);endps;
      beginps('PSline');cproj(improj(x2));cproj(improj(x1));endps;
      rotate(x1);rotate(x2);
   end;
   beginps('PScline');
      for count:=1 to 4 do begin cproj(x1);rotate(x1);end;
   endps;
   beginps('PScline');
      for count:=1 to 4 do begin cproj(x2);rotate(x2);end;
   endps;
   beginps('PScline');
      for count:=1 to 4 do begin cproj(improj(x1));rotate(x1);end;
   endps;
   beginps('PScline');
      for count:=1 to 4 do begin cproj(improj(x2));rotate(x2);end;
   endps;
end{image};



begin{main body}


begintex;

w('\magnification=\magstep1 \hsize=15truecm');
w('\def\dfrac#1#2{{\displaystyle{#1\over#2}}}');
w('\def\display(#1)#2{\smallskip\line{(#1)\hfil#2\hfil}}');

w('');
w('\centerline{\bf Anders Thorup: Input for 3D splines}');
w('\medskip');
w('');

w('\noindent');
w(' The following pages contain examples of the use of a very simple');
w('Pascal program, {\tt 3dspline.pas}.');

w('The program is designed for generation of input data for curve');
w('drawings in \TeX, and it requires (knowledge of) the \TeX-package');
w('{\tt spline.sty}.  The routines in the Pascal program are centered');
w('around one procedure:');


w('\smallskip\centerline{\tt cproj.} \smallskip\noindent ');

w('To describe it, assume that the coordinates of a center in');
w('$3$-space is given.  Associate with it the central projection from');
w('the center onto the plane through the origin perpendicular to the');
w('center.');

w('Then the call {\tt cproj({\bf x})}, for a point {\bf x} in');
w('3-space, computes the image of {\bf x} under the projection, and');
w('writes the two coordinates of the image point to an output file.');
w('The purpose is to use the procedure to generate coordinates that');
w('may serve as input for various spline commands and other curve');
w('drawing commands from the spline package. The output coordinates');
w('depend on the choice of center, and the following pages show how');
w('the output varies under changes of the center. The source code of');
w('the Pascal program contains several routines for manipulating with');
w('vectors in $3$-space.  The body part of the program produces');
w('several examples.');

w('\bigskip\noindent');
w(' As mentioned above, the macros of the Pascal program are designed');
w('for the generation of input data for curve drawings in \TeX;');
w('however, the body part of the program as it is generates a');
w('complete \TeX\ file.  In fact, all the pages you see are in the');
w('source code of the program, and this whole document may be');
w('reproduced as follows: ');
w('');
w('(1) Compile the Pascal program {\tt 3dspline.pas} ');
w('');
w('(2) Run the resulting executable; it will produce a tex file, {\tt');
w('cproj.tex}, containing a lot of spline drawing commands.');
w('');
w('(3) Tex the file {\tt cproj.tex} (with Plain TeX or AmSTeX); the');
w('file {\tt spline.sty} is required. ');
w('');
w('(4) View the result on your favourite device; as always when the');
w('spline package is used it is required that the interpretation');
w('accepts imbedded PostScript code as described in DVIPS.');

w('\bigskip\bigskip \hfill Anders Thorup\par\hfill');
w(docdate);w('');



setcproj(vec(20,1,5));
w('\PSunit0.9truecm\PSpicture{12\PSunit}{3\PSunit}');
w('\llcoords{-6.5 -1}');
image; {sets x1 and x2}
w('\setPSlinewidth{0.3\PSlinewidth}');
for count:=1 to 4 do begin
   beginps('PSline');cproj(x1);cproj(vec(0,10,0));endps;
   beginps('PSline');cproj(x2);cproj(vec(0,10,0));endps;
   rotate(x1);rotate(x2);
end;
endpspicture;


w('\PSunit=1.15truecm');

for j:=1 to 21 do begin
  w('\vfill\eject\null');
  setcproj(vec(10*cos(j*pi/80), 10*sin(j*pi/80), 2.67));
  writeln(outputf, ' %page: ',j+1);
  PSpicture;
    hexahedron; 
    {vertices:}axis(1,1,1); axis(-1,1,1); axis(-1,-1,1); axis(1,-1,1);

    w('\llcoords{-6.5 -13.5}');
    hexahedron; 
    {edges:}axis(1,0,1); axis(0,1,1); axis(-1,0,1); axis(0,-1,1);
    axis(1,1,0); axis(-1,1,0);

    w('\llcoords{-10.5 -13.5}');
    hexahedron;
    {faces:} axis(1,0,0);axis(0,1,0); axis(0,0,1);

    w('\llcoords{-2 -9.5}');
    dodecahedron;

    writeln(outputf, '\llcoords{-5 -9.5}');
    icosahedron;

    writeln(outputf, '\llcoords{-8 -9.5}');
    sphere(1); pscircle(vec(1,-2,-3),1);

    writeln(outputf, '\llcoords{-11 -9.5}');
    angle:=pi; setnormal(vec(0,1,1));
    arcthrough(vec(1,0,0));setnormal(vec(0,5,1));
    arcthrough(vec(1,0,0));

    writeln(outputf, '\llcoords{-2.5 -5.5}');
    hexahedron; tetra1;

    writeln(outputf, '\llcoords{-5.7 -5.5}');
    hexahedron; tetra2;

    writeln(outputf, '\llcoords{-8.9 -5.5}');
    octahedron;hexahedron;

    writeln(outputf, '\llcoords{-11.6 -5.5}');
    w('\putPSmarks');
    PScircle(vec(0,1,0.5),1);w('\nomarks');

    writeln(outputf, '\llcoords{-5.3 -1.5}');
    spiral;

    writeln(outputf, '\llcoords{-12 -1.5}');
    image;
 endPSpicture;
end;

w('\vfill\eject');

w('\noindent{\bf Setup.} \ To obtain coordinates of the vertices of the');
w('regular polyhedra, consider a trapezoid inscribed in a unit circle, say');
w('with vertices of coordinates $(\pm r,l)$ and $(\pm r'',l'')$. If $c$');
w('is the common length of the two non-horizontal sides then           ');
w(' $$r^2+l^2=1,\quad r^{\prime2}+l^{\prime2}=1,\quad(r''-r)^2+(l''-l)^2=c^2.');
w(' $$');
w(' Assume that ${l''}>l$ and ${l''}>0$. Then, in terms of the two ratios');
w('$\mu =r/c$ and $\mu''=r''/c$, the coordinates are given by the formulas');
w('(where $L:=1-(\mu-\mu'')^2$ and $N:=1+4\mu\mu''$):');

w('\display(1''){$ c^2=\dfrac{4L}N,\quad r^2=\dfrac{4\mu^2L}N,\quad');
w('r^{\prime\,2}=\dfrac{4\mu^{\prime\,2}L}N$,}');

w('\smallskip');
w('\display(1''''){${l''}-l=\dfrac{2L}{\sqrt N},\quad ');
w('l=\dfrac{-1+2\mu (\mu -\mu'') }{\sqrt N},\quad');
w('{l''}=\dfrac{1+2{\mu''} (\mu -{\mu''}) }{\sqrt N}\,$.}');

w('\bigskip');

w(' Note the degenerate case of a triangle: ${r''}:=0$, ${l''}=1$. Then,');
w('$\mu''=0$, and ');

w('\smallskip');
  
w('\display(2){$ c^2=4(1-\mu^2),\quad r^2=4\mu^2(1-\mu^2),');
w('\quad l=-1+2\mu^2.$}');

w('\bigskip');

w('If $r$ and $r''$ are the radii of the circles circumscribing the base and');
w('the top of a regular truncated pyramid, based on a regular $n$-gon, and');
w('inscribed in a unit sphere, then the side lengths are $b=2r\sin \pi/n$');
w('and $b''=2r''\sin \pi/n$. Therefore,');

w('\display(3){$\mu =\dfrac{b/c}{2\sin \pi/n}, \quad \mu''');
w('=\dfrac{b''/c}{2\sin \pi/n}.$}');

w('\bigskip');



{Set center}
setcproj(vec(30, 30, 7));


w('\PSmarkradius=3\PSlinewidth\PSunit=2truecm');
w('\PSpicture{8\PSunit}{2\PSunit}');
  w('\llcoords{-2 -1}');
  {Outline of sphere:}
  {r:=r/sqrt(1-(r*r)/(HEIGHT*HEIGHT));}
  beginps('PScircle');
  writeln(outputf, '0 0 ', 1/sqrt(1-1/(HEIGHT*HEIGHT)):8:4);
  endps;
  {Draw truncated pyramid}
  l:=-1/sqrt(10);r:=3/sqrt(10);
  r1:=1/sqrt(5); l1:=2/sqrt(5);
  x:=vec(r,0,l); x1:=vec(r1,0,l1);
  beginps('PSmark'); cproj(x);cproj(x1);
  endps;
  beginps('TX{\hss$(r,0,l)$\hskip5truemm}');cproj(x);
  endps;
  beginps('TX{\hss$(r'',0,l'')$\hskip7truemm}');cproj(x1);
  endps;
  beginps('TX{$\,\,c$\hss}');cproj(vec((r+r1)/2,0,(l+l1)/2));endps;
  beginps('TX{\hskip9truemm$b$\hss}');cproj(vec(r,0,l-0.15));endps;
  beginps('TX{\hskip3truemm$b''$\hss}');cproj(vec(r1,0,l1-0.15));endps;
  beginps('cspline');
    for count:=1 to 5 do begin cproj(x); rotateonefifth(x);end;
  endps;
  beginps('cspline');
    for count:=1 to 5 do begin cproj(x1); rotateonefifth(x1);end;
  endps;
  rotateonefifth(x1);
  for count:=1 to 5 do begin
  	     beginps('PSline'); 
                  cproj(x);
  		rotateonefifth(x);cproj(x); 
  		cproj(x1); rotateonefifth(x1);cproj(x1);
               endps;
  end;
    
  w('\llcoords{-6 -1}');
  {Outline of sphere:}
  {r:=r/sqrt(1-(r*r)/(HEIGHT*HEIGHT));}
  beginps('PScircle');
  writeln(outputf, '0 0 ', 1/sqrt(1-1/(HEIGHT*HEIGHT)):8:4);
  endps;
  {Draw pyramid}
  l:=1/sqrt(17);r:=4/sqrt(17);
  x:=vec(r,0,l); 
  beginps('PSmark');cproj(x);endps;
  beginps('TX{\hss$(r,0,l)$\hskip7.5truemm}');cproj(x);endps;
  beginps('cspline');
    for count:=1 to 5 do begin cproj(x); rotateonefifth(x);end;
  endps;
  for count:=1 to 5 do begin
  	     beginps('PSline'); 
                  cproj(vec(0,0,1));
  		cproj(x); rotateonefifth(x);cproj(x); 
               endps;
  end;
endPSpicture;


w('\bigskip\bigskip\noindent');
w(' {\bf The side lengths of the regular polyhedra.} \ For example, by');
w('using (2) it is easy to determine the side length $c$ of the regular');
w('polyhedra inscribed in a unit sphere, with a vertex at the North Pole:');
w('');
w('For the tetrahedron, $n=3$ and $b=c$. Hence $\mu^2=1/(2\sin\pi/3)=1/3$');
w('and $c^2=8/3$.');
w('');
w('For the hexahedron, $n=3$ and $b=\sqrt 2c$. Hence $\mu ^2=2/(2\sin');
w('\pi/3)^2=2/3$, and $c^2=4/3$.');
w('');
w('For the octahedron, $n=4$ and $b=c$. Hence $\mu^2=1/(2\sin');
w('\pi/4)^2=1/2$, ');
w('and $c^2=2$.');
w('');
w('For the icosahedron, $n=5$ and $b=c$. Hence $\mu^2=1/(2\sin \pi/5)^2$.');
w('As $2\sin^2\pi/5=1-\cos 2\pi/5=(5-\sqrt5)/4$, it follows that');
w('$\mu^2=2/(5-\sqrt5)=(5+\sqrt 5)/10$, and $c^2=2(5-\sqrt 5)/5$. ');
w('');
w('For the dodecahedron, $n=3$, and $b/c$ is the ratio between the diagonal');
w('and the side of regular pentagon, and hence equal to $(2\sin');
w('2\pi/5)/(2\sin \pi /5)=2\cos \pi/5$. So $\mu =\cos(\pi/5)/\sin(\pi/3)$');
w('and $c^2=4(3-4\cos^2(\pi/5))/3=4(3-2(\cos(2\pi/5))/3 =(4\cos');
w('2\pi/5)^2/3$.');
w('');
w('\bigskip\noindent');
w(' {\bf Coordinates.} \ For the icosahedron (one vertex on the');
w('North Pole), use (2) to determine the coordinates of the vertices: One');
w('vertex is $(r,0,l)$, where $r,l$ are given by (2): $r=2/\sqrt 5$,');
w('$l=1/\sqrt 5$. ');
w('');

w('For the dodecahedron (top face horizontal), use (1):');
w(' $b/c=\sin(2\pi /5)/\sin(\pi /5)$ $=2\cos \pi /5$ and $b''/c=1$. Hence');
w('$\mu =\cot \pi/5$ and $\mu''=1/(2\sin \pi/5)$.  It is fastest to use');
w('$c$ determined above: $r=\mu c=(4/\sqrt 3)(\cos 2\pi/5)(\cot \pi/5)$,');
w('and $r''=\mu''c=(2/\sqrt 3)(\cos 2\pi/5)/(\sin \pi/5)$.  ');
w('');
w('\bigskip');


w('\PSmarkradius=5\PSlinewidth\PSunit=2truecm');
w('\PSpicture{5\PSunit}{3\PSunit}');
 w('\llcoords{-3 -2}');
 dodecahedron;
 r1:=sqrt(3)/3*(sqrt(5)-1)/(2*sin(2*pi/10)); 
 l1:=sqrt(1-r1*r1); 
 r:=2/sqrt(3)/(2*sin(2*pi/10));
 l:=sqrt(1-r*r);
 x:=vec(r,0,l); 
 x1:=vec(r1,0,l1);
 x2:=vec(r*cos(2*pi/10), r*sin(2*pi/10), -l);
 beginps('PSmark');cproj(x1);cproj(x);cproj(x2);
 endps;
 beginps('TX{\hss{\tt x}$=(r,0,l)$\hskip5truemm}');cproj(x);
 endps;
 beginps('TX{\hss{\tt x1}$=(r'',0,l'')$\hskip5truemm}');cproj(x1);
 endps;
 beginps('TX{\hskip22truemm$(r\cos2\pi/10,r\sin2\pi/10,-l)=${\tt x2}\hss}');
    cproj(x2);
 endps;
 w(' \setPSlinewidth{3\PSlinewidth}');
 beginps('PSline');
   cproj(x2);cproj(x);cproj(x1);rotateonefifth(x1);cproj(x1);
 endps;
 w(' \setPSlinewidth{0.5pt}');
 w('\RTX{\tt r1:=1/sqrt(3)*(sqrt(5)-1)/(2*sin(2*pi/10))}{-1 -1.05}');
 w('\RTX{\tt l1:=sqrt(1-r1*r1)}{-1 -1.25}'); 
 w('\RTX{\tt r:=1/sqrt(3)/sin(2*pi/10)}{-1 -1.45}');
 w('\RTX{\tt l:=sqrt(1-r*r)}{-1 -1.65}');
endpspicture;


w('\noindent');
w(' The path marked above on the dodecahedron begins at the point x2');
w('and ends at the point obtained by rotating x1 one fifth around');
w('the z-axis. So the path may be drawn by a PSline, produced by');
w('the following block of code:');
w('\smallskip');
w('{\tt beginps(''PSline'');\par');
w('\ \ cproj(x2);\par');
w('\ \ cproj(x);\par');
w('\ \ cproj(x1);\par');
w('\ \ rotateonefifth(x1);cproj(x1);\par');
w('endps;\par');
w('rotateonefifth(x2);rotateonfifth(x);\par}');
w('\smallskip');
w('\noindent After the first `rotateonefifth'', the point x1 becomes');
w('the end point of the path. ');
w('The next two `rotateonefifth''s do not contribute to the PSline,');
w('but when the code has been completed, all three points x2, x, x1 have');
w('been rotated. So, the top part of the dodecahedron may be drawn by');
w('repeating the code 5 times. The bottom part may be drawn by reflecting');
w('the top part in the origin. ');
w('');



endtex;

writeln('*********************Done!*********************');;

end.