% --- start of displayed preamble in the book ---
% circulargradient as defined above
% --- end of displayed preamble in the book ---

% graphic converted to gray in book using 'color2gray'

def PARALLEL = 0 enddef;
def CIRCULAR   = 1 enddef;

def max(expr a,b)= (if a>b: a else: b fi) enddef;

def f(expr x)=((1+sind(2x*360))/2) enddef;
def g(expr x)=(x*x) enddef;

% f:   gradient function
% a,b: interval for f
% n:   steps
% ca: color 1
% cb: color 2
% p:   path
% c:   center of gradient (CIRCULAR) or direction (PARALLEL)
vardef gradient(text f)(expr a,b,n,ca,cb,type)(expr p,c)=
 save A,B,C,D,pic,rp,r;
 pair A,B,C,D;
 picture pic;
 path rp;
 if type=PARALLEL:rp=p rotated (-angle(c));else:rp=p;fi;
 A=llcorner(rp);
 C=urcorner(rp);
 B=(xpart C,ypart A);
 D=(xpart A,ypart C);
 if type=PARALLEL:
   A:=A rotated angle(c);
   B:=B rotated angle(c);
   C:=C rotated angle(c);
   D:=D rotated angle(c);
 fi;
 if type=CIRCULAR:
   r=max(arclength(c--A),arclength(c--B));
   r:=max(r,arclength(c--C));
   r:=max(r,arclength(c--D));
 fi;
 pic=nullpicture;
 if type=PARALLEL:
   for i:=0 upto n-1:
     addto pic
       contour ((i/n)[A,B]--((i+1)/n)[A,B]
                  --((i+1)/n)[D,C]
                  --(i/n)[D,C])--cycle
       withcolor ((f(i/n))[ca,cb]);
     endfor;
 elseif type=CIRCULAR:
   for i:=n step -1 until 1:
     addto pic
       contour fullcircle
               scaled (2*(i/n)*r) shifted c
       withcolor ((f(i/n))[ca,cb]);
     endfor;
 fi;
 clip pic to p;
 draw pic;
enddef;

def parallelgradient(text f)(expr a,b,n,ca,cb)(expr p,d)=
 gradient(f)(a,b,n,ca,cb,PARALLEL)(p,d);
enddef;

def circulargradient(text f)(expr a,b,n,ca,cb)(expr p,c)=
 gradient(f)(a,b,n,ca,cb,CIRCULAR)(p,c);
enddef;

defaultfont:="ptmr8r";
warningcheck:=0;
beginfig(1)
numeric u;
u=1cm;
path p;
p=fullcircle xscaled 4u
            yscaled 2u
            shifted (0,-4u);
circulargradient(g)(0,1,100,green,red)
               (p,center p+(2u,0));
endfig;
end;