% --- start of displayed preamble in the book ---

% --- end of displayed preamble in the book ---
defaultfont:="ptmr8r";
warningcheck:=0;
beginfig(1)
% Fraser's Spiral
% D. Roegel 29 April 2007

numeric u;
u=2.7cm;

vardef twisted_cord(expr i,twist,a,col)=
  save qa,qb,qc,P,ratio;
  path qa,qb,qc;
  pair P[];
  ratio=0.8;
  P1=A[i];
  P2=radius(i)*dir((1+i/2)*da);
  P3=radius(i)*dir((1+i/2-shift)*da);
  P5=B[i-1] rotated da;
  P6=radius(i)*dir((2+i/2)*da);
  P7=radius(i)*dir((2+i/2-shift)*da);
  P0=.5[P2,P7];
  P21=arclength((.5[ratio,1])[P1,P2]--origin)*unitvector(P0);
  P22=arclength((.5[ratio,1])[P5,P7]--origin)*unitvector(P0);

  P12=P2 rotatedaround(P0,twist);
  P13=P3 rotatedaround(P0,twist);
  P16=P6 rotatedaround(P0,twist);
  P17=P7 rotatedaround(P0,twist);
  P23=P21 rotatedaround(P0,twist);
  P24=P22 rotatedaround(P0,twist);


  P14=ratio[P1,P12];
  P18=ratio[P5,P17];

  qa=P1--P12--P13--cycle;
  qb=P5--P16--P17--cycle;
  qc=P14--P12...P24...P18--P17..P23..cycle;
  fill qa rotated a withcolor col;
  fill qb rotated a withcolor col;
  fill qc rotated a withcolor col;
enddef;

def circle(expr r)=
  fullcircle scaled 2r
enddef;

def radius(expr i)=
  (((i+10*i*i)/1000)*u)
enddef;

color cola,colb;
cola=.9(1,1,0);
colb=.8white;

cola:=.9(1,1,0);
colb:=.7green;

nc=16;
n=16;
da=360/n;
shift=0.4;
path p[];
pair A[],B[];

%beginfig(1);

  fill unitsquare shifted (-.5,-.5) scaled 3u withcolor cola;
  z0=origin ;
  if false:
  for i=1 upto nc:
    draw circle(radius(i));
  endfor;

  pickup pencircle scaled 5pt;
  for i=1 upto nc:
    for j=1 upto n:
      draw radius(i)*dir((j+i/2)*da) withcolor red;
      draw radius(i)*dir((j+i/2-shift)*da) withcolor green;
    endfor;
  endfor;

  pickup pencircle scaled .5;

  fi;


  p1=(radius(1)*dir((1+1/2-shift)*da))
  for i=2 upto nc:
    ..(radius(i)*dir((1+i/2-shift)*da))
  endfor;
  p2=p1 rotated (shift*da);
  p3=(radius(1)*dir((1+1/2-shift)*da))
  for i=2 upto nc:
    ..(radius(i)*dir((1-i/2-shift)*da))
  endfor;
  p4=p3 rotated (shift*da);

  if false:
  for j=1 upto n:
    draw p1 rotated (j*da);
    draw p2 rotated (j*da);
    draw p3 rotated (j*da);
    draw p4 rotated (j*da);
  endfor;
  fi;

  p5=p1--(reverse p2)--cycle;
  %fill p5 withcolor .8white;
  p6=p2--(reverse (p1 rotated da))--cycle;
  %fill p6 withcolor .8green;

  p7=p3--(reverse p4)--cycle;
  %fill p7 withcolor .8blue;
  p8=p4--(reverse (p3 rotated da))--cycle;
  %fill p8 withcolor .8red;

  for j=1 upto n:
    fill p7 rotated (j*da) withcolor colb;
  endfor;
  for j=1 upto n:
    fill p5 rotated (j*da) withcolor colb;
  endfor;


  pickup pencircle scaled 5pt;
  for i=1 upto nc-1:
    A[i]=p1 intersectionpoint (p4 rotated (i*da));
    B[i]=p2 intersectionpoint (p3 rotated (i*da+da));
    %draw A[i] withcolor (1,1,0);
    %draw B[i] withcolor (1,0,1);
  endfor;

  pickup pencircle scaled .5;
  for i=2 upto nc-1:
    for j=1 upto n:
      twisted_cord(i,3,j*da,(i/nc)[colb,black]);
    endfor;
  endfor;

  clip currentpicture to unitsquare shifted (-.5,-.5) scaled 3u;

  picture lab;
  % intentionally slightly shifted to the right
  lab=thelabel(btex Fraser's spiral etex,z0+.1u*right);
  unfill bbox lab;
  draw lab withcolor red;

%endfig;

%end

endfig;
end;