% featpost3Dplus2D.mp
% L. Nobre G., lnobreg@gmail.com, http://matagalatlante.org
% C. Barbarosie
% J. Schwaiger
% B. Jackowski
% P. J. Sebastião
% P. Jørgensen
% S. Pakin
% 
% Copyright (C) 2014

% This set of macros adds a lot of features to
% the MetaPost language and eases the production of 
% physics diagrams and animations.

% This is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 3
% of the License, or (at your option) any later version.

% This set of macros is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
% GNU General Public License for more details.

    message "FeatPost-0.8.8";

    warningcheck := 0;
    background := 0.987white;
    defaultscale := 0.75;
    defaultfont := "cmss17";   % This is used by cartaxes
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global Variables %%%%%%%%%%%

    boolean	 ParallelProj, SphericalDistortion, FCD[], ShadowOn;
    boolean      OverRidePolyhedricColor, MalcomX;
    numeric      Nobjects, RefDist[], HoriZon, RopeColorSeq[], PhotoMarks;
    numeric	 Spread, PrintStep, PageHeight, PageWidth, ActuC, Shifts;
    numeric      NL, npl[], NF, npf[], FC[], MaxFearLimit, TableColors;
    numeric      TDAtiplen, TDAhalftipbase, TDAhalfthick, RopeColors, NCL;
    pair	 OriginProjPagePos, ShiftV, PhotoPair[];
    path	 VGAborder, CLPath[];
    color        f, viewcentr, V[], L[]p[], F[]p[], TableC[]; 
    color        HigColor, SubColor, LightSource, PhotoPoint[];
    string       ostr[];
    pen          BackPen, ForePen;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Default Values %%%%%%%%%%%%%%%

            f := (3,5,4);   % This f is the point of view in 3D
	    
	    viewcentr := black; % This is the aim of the view
	    
            Spread := 140;   % Magnification

            Shifts := 105.00mm;
    
            ShiftV := Shifts*(1,1); % Central coordinates on paper

	    OriginProjPagePos := (Shifts,148.45mm); % This should be the
	                                            % page center. 

	    ParallelProj := false; % Kind of perspective %
				                         % Can't have both true
	    SphericalDistortion := false; % Kind of lens %

	    VGAborder := (182.05,210.00)--       % This definition assumes
			 (412.05,210.00)--       % ShiftV = 105.00mm(1,1)
			 (412.05,382.05)--       % Use: gs -r200 and you 
			 (182.05,382.05)--cycle; % get few extra pixels

	    PrintStep := 5;   % Coarseness, in resolvec

	    PageHeight := 9in;
	    PageWidth := 6in;   % And this is used by produce_auto_scale

            MaxFearLimit := 17;       % Valid Maximum Distance from Origin
	    
            HigColor := 0.85white;          % These two colors are used in
	    SubColor := 0.35white;          %        fillfacewithlight
	    LightSource := 10*(4,-3,4);     % This also
	    OverRidePolyhedricColor:=false; % And also this
	    ShadowOn := false; % Some objects may block the light and
	    HoriZon := 0;      % cast a shadow on a horizontal plane at this Z
	    	    	    
	    TableC0 := 0.85white;           % grey        %% G N U P L O T
	    TableC1 := red;                 % red         %% 
	    TableC2 := ( 0.2, 0.2, 1.0 );   % blue        %%    colors 
	    TableC3 := ( 1.0, 0.7, 0.0 );   % orange      %%
	    TableC4 := 0.85green;           % pale green  %%
	    TableC5 := 0.90*(red+blue);     % magenta     %%
	    TableC6 := 0.85*(green+blue);   % cyan        %%
	    TableC7 := 0.85*(red+green);    % yellow      %%
	    
	    TableColors := 7;
	    ActuC := 5;

	    RopeColorSeq0 := 3;             % 
	    RopeColorSeq1 := 3;             % 
	    RopeColorSeq2 := 1;             % 
	    RopeColorSeq3 := 3;             % ropepattern 
	    RopeColorSeq4 := 7;             % 
	    RopeColorSeq5 := 5;             % 
                                            %
	    RopeColors := 5;                % 

	    Nobjects := 0;                  % getready and doitnow

	    TDAtiplen := 0.05;              % tdarrow and tdcircarrow
	    TDAhalftipbase := 0.02;         % Three-Dimensional (Circular)
	    TDAhalfthick := 0.01;           % Arrow

	    NCL := 0;                       % closedline
	    ForePen := pencircle scaled 15pt;
	    BackPen := pencircle scaled 9pt;

	    MalcomX := false;
	    
            %%% The variables PhotoMarks, PhotoPair[], PhotoPoint[]
	    %%% and CLPath[] have NO default values.
	    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Part I:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Very basic:
	    
% Colors have three or four coordinates. Get one. 

    def X(expr A) = 
      if color A: redpart A elseif MalcomX: blackpart A else: cyanpart A fi
    enddef;
  
    def Y(expr A) =
      if color A: greenpart A else: magentapart A fi
    enddef;
  
    def Z(expr A) =  
      if color A: bluepart A else: yellowpart A fi
    enddef;
  
    def W(expr A) =   
      blackpart A
    enddef;
    
% The length of a vector.

    def conorm(expr A) = 
      ( X(A) ++ Y(A) ++ Z(A) )  
    enddef;

    def cmyknorm(expr A) = %% This is not good when MalcomX is true 
      ( X(A) ++ Y(A) ++ Z(A) ++ W(A) )  
    enddef;
  
    def makecmyk( expr A, B ) =
      ( ( X(A), Y(A), Z(A), B ) )
    enddef;
  
    def maketrio( expr A  ) =
      ( ( X(A), Y(A), Z(A) ) )
    enddef;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Vector Calculus:
	    
% Calculate the unit vector of a vector (or a point)

    def N(expr A) =
        begingroup
            save M, exitcolor;
            numeric M;
            color exitcolor;
            M = conorm( A );
            if M > 0:
                exitcolor = ( X(A)/M, Y(A)/M, Z(A)/M );
            else:
                exitcolor := black;
            fi;
            ( exitcolor )
        endgroup
    enddef;

    def cdotprod(expr A, B) = 
        ( X(A)*X(B) + Y(A)*Y(B) + Z(A)*Z(B) )
    enddef;

    def ccrossprod(expr A, B) = 
        ( Y(A)*Z(B) - Z(A)*Y(B), 
          Z(A)*X(B) - X(A)*Z(B), 
          X(A)*Y(B) - Y(A)*X(B) )
    enddef;

% The dotproduct of two normalized vectors is the cosine of the angle 
% they form.

    def ndotprod(expr A, B) = 
        begingroup
            save a, b;
            color a, b;
            a = N(A);
            b = N(B);
            ( ( X(a)*X(b) + Y(a)*Y(b) + Z(a)*Z(b) ) )
        endgroup
    enddef;

% The normalized crossproduct of two vectors. 

    def ncrossprod(expr A, B) = 
        N( ccrossprod( A, B ) )
    enddef;

% Haahaa! Trigonometry. 

    def getangle(expr A, B) =
        begingroup 
            save coss, sine;
            numeric coss, sine;
            coss := cdotprod( A, B );
            sine := conorm( ccrossprod( A, B ) );
            ( angle( coss, sine ) )
        endgroup
    enddef;

% Something I need for spatialhalfsfear.

    def getcossine( expr Center, Radius ) =
      begingroup
	save a, b;
	numeric a, b;
	a = conorm( f - Center );
	b = Radius/a;
	if abs(b) >= 1:
	  show "The point of view f is too close (getcossine).";
	  b := 2;                                                % DANGER!
	fi;
	( b )
      endgroup
    enddef;
    
% The following routine is used by circularsheet and may be used to
% rotate vectors elliptically. Also used by tdcircarrow.

    vardef planarrotation( expr VecX, VecY, TheAngle ) =
      ( VecX*cosd( TheAngle ) + VecY*sind( TheAngle ) )
    enddef;

% The following routine could be used by kindofcube and may be used to
% rotate polyhedra (must cycle through all Vs before calling makeface).

    def eulerrotation( expr AngA, AngB, AngC, Vec ) =
      begingroup
	save auxx, auxy, veca, vecb, vecc;
	color auxx, auxy, veca, vecb, vecc;
	veca = ( cosd(AngA)*cosd(AngB),
	  sind(AngA)*cosd(AngB),
	  sind(AngB) );
	auxx = ( cosd(AngA+90), sind(AngA+90), 0 );
	auxy = ccrossprod( veca, auxx );
	vecb = cosd(AngC)*auxx + sind(AngC)*auxy;
	vecc = cosd(AngC+90)*auxx + sind(AngC+90)*auxy;
	( X(Vec)*veca + Y(Vec)*vecb + Z(Vec)*vecc )
      endgroup
    enddef;

% Rotate a vector around another. Supposes all vectors share the same origin. 

    def rotvecaroundanother( expr Angle, RotVec, FixVec ) =
      begingroup
	save uf, cf, xr, yr;
	color uf, cf, xr, yr;
	uf = N( FixVec );
	yr = ccrossprod( uf, RotVec );
	cf = uf*cdotprod( uf, RotVec );
	xr = RotVec - cf;
	( cf + planarrotation( xr, yr, Angle ) )
      endgroup
    enddef;
    
% inplanarvolume is used by kindofcube.    
    
    def inplanarvolume( expr PointPerpA, PointPerpB, Point ) =
      begingroup
	save va, vb, vc;
	color va, vb, vc;
	va = Point - PointPerpA;
	vb = Point - PointPerpB;
	vc = PointPerpB - PointPerpA;
	( cdotprod(va,vc)*cdotprod(vb,vc) <= 0 )
      endgroup
    enddef;

% Maybe you would like to calculate the angular arguments of kindofcube...
    
    def getanglepair( expr InVec ) =
      begingroup
	save alphaone, alphatwo;
	numeric alphaone, alphatwo;
	alphaone = angle( ( X(InVec), Y(InVec) ) );
	alphatwo = angle( ( X(InVec) ++ Y(InVec), Z(InVec) ) );
	( (alphaone,alphatwo) ) 
      endgroup
    enddef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Auxiliary:
	    
% Projection Size. Meant for objects with size one.
% Used by signalvertex.

    def ps(expr A, Thicken_Factor) =
        Thicken_Factor/conorm(A-f)/3
    enddef;

% Rigorous Projection of a Point. Draws a circle with
% a diameter inversely proportional to the distance of
% that Point from the point of view.

    def signalvertex(expr A, TF, Col) =
      draw rp(A) withcolor Col withpen pencircle scaled (Spread*ps(A,TF))
    enddef;

    def signalshadowvertex(expr A, TF, Col) =
      begingroup
	save auxc, auxn;
	color auxc;
	numeric auxn;
	auxc := cb(A);
	auxn := TF*conorm(f-auxc)/conorm(LightSource-A);
	signalvertex( auxc, auxn, Col )
      endgroup
    enddef;

% Get the vector that projects onto the resolution

    def resolvec(expr A, B) =
        begingroup
            save sizel, returnvec;
            numeric sizel;
            color returnvec;
            sizel = abs( rp(A) - rp(B) );
            if sizel > 0:
                returnvec = PrintStep*(B-A)/sizel;
            else:
                returnvec = 0.3*(B-A);
            fi;
            ( returnvec )
        endgroup
    enddef;

% Movies need a constant frame

    def produce_vga_border =
      begingroup
	draw VGAborder withcolor background withpen pencircle scaled 0;
	clip currentpicture to VGAborder
      endgroup
    enddef;

% The following macro fits a figure in the page.
% Probably it is obsolete since MetaPost 1.000
% Should be the last command in a figure.

    def produce_auto_scale =
        begingroup
            picture storeall, scaleall;
	    numeric pwidth, pheight;
	    storeall = currentpicture shifted -(center currentpicture);
	    currentpicture := nullpicture;
	    pwidth = xpart ((lrcorner storeall)-(llcorner storeall));
	    pheight = ypart ((urcorner storeall)-(lrcorner storeall));
	    if PageHeight/PageWidth < pheight/pwidth:
	        scaleall = storeall scaled (PageHeight/pheight);
            else:
	        scaleall = storeall scaled (PageWidth/pwidth);
            fi;
	    draw scaleall shifted OriginProjPagePos
        endgroup
    enddef;

% The following two procedures are useful for getready.
    
    vardef cstr( expr Cl ) =
      "(" &
	decimal(X(Cl)) &
	"," &
	decimal(Y(Cl)) &
	"," &
	decimal(Z(Cl)) &
	")"
    enddef;

    vardef bstr( expr bv ) =
      save bstring; string bstring;
      if bv: bstring = "true"; else: bstring = "false"; fi;
      bstring
    enddef;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Fundamental:
	    
% Rigorous Projection. This the kernel of all these lines of code.
% It won't work if R belongs the plane that contains f and that is 
% ortogonal to vector f-viewcentr, unless SphericalDistortion is true.
% f-viewcentr must not be (anti-)parallel to zz.
 
    def rp(expr R) =
        begingroup

	    save projpoi;
            save v, u;
            save verti, horiz, eta, squarf, radio, ang, lenpl;
	    pair projpoi;
            color v, u;
            numeric verti, horiz, eta, squarf, radio, ang, lenpl;

            v = N( (-Y(f-viewcentr), X(f-viewcentr), 0) );
            u = ncrossprod( f-viewcentr, v );

	    horiz = cdotprod( R-viewcentr, v );
	    verti = cdotprod( R-viewcentr, u );

	    if SphericalDistortion:
		if ( horiz <> 0 ) or ( verti <> 0 ):
	    	    lenpl = ( horiz ++ verti )*20; %%%%%%%%%%%%%%% DANGER
	    	    ang = getangle( f-R, f-viewcentr );
		    horiz := ang*horiz/lenpl;
		    verti := ang*verti/lenpl;
		    projpoi = (horiz,verti);
		else:
		    projpoi = origin;
		fi;			
	    else:
	        if ParallelProj:
		    eta = 1;
	        else:
	            squarf = cdotprod( f-viewcentr, f-viewcentr );
		    radio = cdotprod( R-viewcentr, f-viewcentr );
		    eta = 1 - radio/squarf;
		    if abs((horiz,verti)) > MaxFearLimit*eta:  
		        eta := abs((horiz,verti))/MaxFearLimit;
		    fi;
	        fi;
	        projpoi = (horiz,verti)/eta;
	    fi;

	    ( projpoi*Spread + ShiftV )

        endgroup
    enddef;

% Much improved rigorous pseudo-projection algorithm that follows 
% an idea from Cristian Barbarosie.
% This makes shadows caused by a light source point.

    def cb(expr R) =
      begingroup
	save    ve, ho;
	numeric ve, ho;
	LightSource-ho*red-ve*green-HoriZon*blue=whatever*(LightSource-R);
	( ho*red + ve*green + HoriZon*blue ) 
      endgroup
    enddef;

% And this just projects points rigorously on some generic plane using
% LightSource as the point of convergence (focus).
    
    def projectpoint(expr ViewCentr, R) =
        begingroup
            save    verti, horiz;
            save    v, u, lray;
            numeric verti, horiz;
            color   v, u, lray;
	    lray = LightSource-ViewCentr;
            v = N( (-Y(lray), X(lray), 0) );
            u = ncrossprod( lray, v );
            lray - horiz*v - verti*u = whatever*( LightSource - R );
            ( horiz*v + verti*u + ViewCentr ) 
        endgroup
    enddef;

% And this is the way to calculate the intersection of some line with some
% plan.

    def lineintersectplan( expr LinePoi, LineDir, PlanPoi, PlanDir ) =
      begingroup
	save incognitus;
	color incognitus;
	cdotprod( incognitus-PlanPoi, PlanDir ) = 0;
	whatever*LineDir + LinePoi = incognitus;
	( incognitus )
      endgroup
    enddef;
	
% Vanishing point.     
    
    def vp( expr D ) =
      begingroup
	( rp( lineintersectplan( f, D, viewcentr, f-viewcentr) ) )
      endgroup
    enddef;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Basic Functions:
	    
% Get the 2D path of a straight line in beetween 3D points A and B.
% This would add rigor to rigorousdisc, if one would introduce the
% the concept of three-dimensional path. That is not possible now.
% Also this is only interesting when using SphericalDistortion:=true
    
    def pathofstraightline( expr A, B ) =
        begingroup
            save  k, i, mark, stepVec, returnp, pos; 
            numeric k, i;
            color mark, stepVec;
	    path returnp;
	    pair pos[];
            stepVec = resolvec(A,B);
	    pos0 = rp( A );
            k = 1;
            forever:
	        mark := A+(k*stepVec);
                exitif cdotprod(B-mark,stepVec) <= 0;
		pos[k] = rp( mark );
		k := incr(k);
	    endfor;
	    pos[k] = rp(B);
	    returnp = pos0 for i=1 upto k: ..pos[i] endfor;
	    ( returnp )
	endgroup
    enddef;

    def drawsegment( expr A, B ) =
        begingroup
	    if SphericalDistortion:
	      draw pathofstraightline( A, B );
	    else:
	      draw rp(A)--rp(B);
	    fi
	endgroup
    enddef;

% Cartesian axes with prescribed lengths.

    def cartaxes(expr axex, axey, axez) =
        begingroup
            save orig, axxc, ayyc, azzc;
            color orig, axxc, ayyc, azzc;
            orig = (0,0,0);
            axxc = (axex,0,0);
            ayyc = (0,axey,0);
            azzc = (0,0,axez);
            drawarrow rp(orig)..rp(axxc);
            drawarrow rp(orig)..rp(ayyc);
            drawarrow rp(orig)..rp(azzc);
            label.bot( "x" ,rp(axxc));       %%%%%%%%%%%%%%%%%%%%%%%%% 
            label.bot( "y" ,rp(ayyc));       %%   Some  Labels...   %%
            label.lft( "z" ,rp(azzc));       %%%%%%%%%%%%%%%%%%%%%%%%% 
        endgroup
    enddef;

% Orthogonal axes with prescribed lengths and labels.

    def orthaxes(expr axex, strx, axey, stry, axez, strz ) =
        begingroup
            save axxc, ayyc, azzc;
            color axxc, ayyc, azzc;
            axxc = (axex,0,0);
            ayyc = (0,axey,0);
            azzc = (0,0,axez);
            drawarrow rp(black)..rp(axxc);
            drawarrow rp(black)..rp(ayyc);
            drawarrow rp(black)..rp(azzc);
            label.bot( strx ,rp(axxc));        
            label.bot( stry ,rp(ayyc));       
            label.lft( strz ,rp(azzc));       
        endgroup
    enddef;

% This is it. Draw an arch beetween two straight lines with a
% common point (Or) in three-dimensional-euclidian-space and 
% place a label near the middle of the arch. Points A and B
% define the lines. The arch is at a distance W from Or. The
% label is S and the position is RelPos (rt,urt,top,ulft,lft,
% llft,bot,lrt). But arches must be smaller than 180 degrees.

    def angline(expr A, B, Or, W, S)(suffix RelPos) =
        begingroup
            save G, Dna, Dnb, al;
            numeric G;
            color Dna, Dnb;
            path al;
            G = conorm( W*( N(A-Or) - N(B-Or) ) )/2.5; %%%%%%% BIG DANGER!
            Dna = ncrossprod(ncrossprod(A-Or,B-Or),A-Or);
            Dnb = ncrossprod(ncrossprod(B-Or,A-Or),B-Or);
            al = rp(W*N(A-Or)+Or)..
                                  controls rp(W*N(A-Or)+Or+G*Dna) 
                                       and rp(W*N(B-Or)+Or+G*Dnb)..
                  rp(W*N(B-Or)+Or);
            draw al;
            label.RelPos( S, point 0.5*length al of al )
        endgroup
    enddef;

% As i don't know how to declare variables of type suffix,
% i provide a way to avoid the problem. This time RelPos may
% be 0,1,2,3,4,6,7 or anything else.

    def anglinen(expr A, B, Or, W, S, RelPos) =
        begingroup
            save G, Dna, Dnb, al, middlarc;
            numeric G;
            color Dna, Dnb;
            path al;
            pair middlarc;
            G = conorm( W*( N(A-Or) - N(B-Or) ) )/3;
            Dna = ncrossprod(ncrossprod(A-Or,B-Or),A-Or);
            Dnb = ncrossprod(ncrossprod(B-Or,A-Or),B-Or);
            al = rp(W*N(A-Or)+Or)..
                                  controls rp(W*N(A-Or)+Or+G*Dna) 
                                       and rp(W*N(B-Or)+Or+G*Dnb)..
                  rp(W*N(B-Or)+Or);
            draw al;
            middlarc = point 0.5*length al of al;
            if RelPos = 0:
                label.rt( S, middlarc );
            elseif RelPos =1:
                label.urt( S, middlarc );
            elseif RelPos =2:
                label.top( S, middlarc );
            elseif RelPos =3:
                label.ulft( S, middlarc );
            elseif RelPos =4:
                label.lft( S, middlarc );
            elseif RelPos =5:
                label.llft( S, middlarc );
            elseif RelPos =6:
                label.bot( S, middlarc );
            elseif RelPos =7:
                label.lrt( S, middlarc );
            else:
                label( S, middlarc );
            fi
        endgroup
    enddef;

% As a bigger avoidance, replace the arch by a paralellogram.

    def squareangline(expr A, B, Or, W) =
      begingroup
	save sal;
	path sal;
	sal = rp(Or)--rp(W*N(A-Or)+Or)-- 
	      rp(W*(N(B-Or)+N(A-Or))+Or)--rp(W*N(B-Or)+Or)--cycle;
	draw sal
      endgroup
    enddef;

% Just as we are here we can draw circles. (color,color,numeric)

    def rigorouscircle( expr CenterPos, AngulMom, Radius ) =
      begingroup
	save ind, G, Dn, Dna, Dnb, al, vec;
	numeric ind, G;
	color vec[], Dn, Dna, Dnb;
	path al;
	vec1 = ncrossprod( CenterPos-f, AngulMom);
	for ind=2 step 2 until 8:
	  vec[ind+1] = ncrossprod( vec[ind-1], AngulMom );
	  vec[ind] = N( vec[ind-1] + vec[ind+1] ); 
	endfor;
	G = conorm( Radius*( vec1 - vec2 ) )/3;
	al = rp(Radius*vec1+CenterPos)
	for ind=2 upto 8: 
	  hide(
	    Dn:=ncrossprod(vec[ind-1],vec[ind]);
	    Dna:=ncrossprod(Dn,vec[ind-1]);
	    Dnb:=ncrossprod(-Dn,vec[ind]) 
	    )
	  ..controls rp(Radius*vec[ind-1]+CenterPos+G*Dna) 
	  and rp(Radius*vec[ind]  +CenterPos+G*Dnb)
	  ..rp(Radius*vec[ind]  +CenterPos)
	endfor
	...cycle;
	( al ) 
      endgroup
    enddef;

% 3D straight arrow.

    def tdarrow(expr FromPos, ToTip ) =
      begingroup
	save basevec, longvec, a, b, c, d, e, g, h, len, p;
	color basevec, longvec, a, b, c, d, e, g, h;
	numeric len;
	path p;
	len = conorm( ToTip - FromPos );
	longvec := N( ToTip - FromPos );
	basevec := ncrossprod( FromPos-f, longvec );
	if len <= TDAtiplen:
	  b = basevec*TDAhalftipbase*len/TDAtiplen;
	  c = FromPos+b;
	  e = FromPos-b;
	  p = rp(ToTip)--rp(c)--rp(e)--cycle;
	else:
	  d = ToTip-longvec*TDAtiplen;
	  a = FromPos+basevec*TDAhalfthick;
	  h = FromPos-basevec*TDAhalfthick;
	  b = d+basevec*TDAhalfthick;
	  g = d-basevec*TDAhalfthick;
	  c = d+basevec*TDAhalftipbase;
	  e = d-basevec*TDAhalftipbase;
	  p = rp(a)--rp(b)--rp(c)--rp(ToTip)--rp(e)--rp(g)--rp(h)--cycle;
	fi;
	unfill p;
	draw p
      endgroup
    enddef;

% 3D circular arrow.

    def tdcircarrow(expr CenterPos, AngulMom, Ray, StartAngle, Amplituda ) =
      begingroup
	save veca, vecb, vecc, vecd, a, b, c, d, p, stepa, numa, anga, angb;
	save signus, ca, da, i;
	color veca, vecb, vecc, vecd;
	pair a, b, c, d, ca, da, aa;
	numeric stepa, numa, anga, angb, signus, i;
	path p;
	signus = Amplituda/abs(Amplituda);
	stepa = 6signus; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER.
	vecd = ncrossprod( CenterPos-f, AngulMom );
	vecc = ncrossprod( vecd, AngulMom );
	anga = signus*180*(TDAhalfthick/TDAhalftipbase)*TDAtiplen/(3.14159*Ray);
	angb = signus*180*TDAtiplen/(3.14159*Ray);
	numa = StartAngle+Amplituda;
	a = rp(Ray*planarrotation(vecc,vecd,StartAngle+anga));
	b = rp(Ray*planarrotation(vecc,vecd,numa+angb));
	c = rp((Ray+TDAhalftipbase)*planarrotation(vecc,vecd,numa));
	d = rp((Ray-TDAhalftipbase)*planarrotation(vecc,vecd,numa));
	ca = rp((Ray+TDAhalfthick)*planarrotation(vecc,vecd,numa));
	da = rp((Ray-TDAhalfthick)*planarrotation(vecc,vecd,numa));
	aa = rp((Ray-TDAhalfthick)*planarrotation(vecc,vecd,StartAngle));
	p = for i=StartAngle step stepa until numa:
	      rp((Ray+TDAhalfthick)*planarrotation(vecc,vecd,i))..
	    endfor ca--c--b--d--da..
	    for i=numa-stepa step -stepa until StartAngle:
	      rp((Ray-TDAhalfthick)*planarrotation(vecc,vecd,i))..
	    endfor aa--a--cycle;
	unfill p;
	draw p
      endgroup
    enddef;

% Draw lines with a better expression of three-dimensionality.

    def emptyline(expr JoinP,ThickenFactor,OutCol,InCol,theN,EmptyFrac,sN)(text LinFunc) = 
      begingroup
	save i, j, k; 
	numeric i, j, k;
	k = ThickenFactor*EmptyFrac;
	if ShadowOn:
	  for i = 0 upto theN:
	    signalshadowvertex( LinFunc(i/theN), ThickenFactor, black );
	  endfor;
	fi;
	for j = 0 upto sN-1:
	  signalvertex( LinFunc(j/theN), ThickenFactor, OutCol );
	endfor;
	if JoinP:
	  for j = -sN upto 0:
	    signalvertex( LinFunc(j/theN), k, InCol );
	  endfor;
	fi;
	for i = sN upto theN:
	  signalvertex( LinFunc( i/theN ), ThickenFactor, OutCol );
	  for j = sN downto 0: 
	    signalvertex( LinFunc( (i-j)/theN ), k, InCol );
	  endfor;
	endfor
      endgroup 
    enddef;

% Draw space-paths of possibly closed lines making use of "getready"

    def closedline( expr ThisIsClosed, theN, ForeFrac, BackFrac )( text LinFunc ) =
      begingroup
	save i, comm;
	numeric i;
	string comm;
	NCL := incr( NCL );
	if ThisIsClosed:
	  CLPath[NCL] :=
	    for i=1 upto theN:
	      rp(LinFunc(i/theN))..
	    endfor
	    cycle;
	  for i=1 upto theN:
	    comm:="draw subpath ("
	      & decimal(i-ForeFrac)
	      & ","
	      & decimal(i+ForeFrac)
	      & ") of CLPath"
	      & decimal(NCL)
	      & " withpen ForePen; undraw subpath ("
	      & decimal(i-BackFrac)
	      & ","
	      & decimal(i+BackFrac)
	      & ") of CLPath"
	      & decimal(NCL)
	      & " withpen BackPen;";
	    getready( comm, LinFunc(i/theN) );
	  endfor;
	else:
	  CLPath[NCL] := rp(LinFunc(0))
	    for i=1 upto theN: ..rp(LinFunc(i/theN)) endfor;
	  comm:="draw subpath (0,"
	    & decimal(ForeFrac)
	    & ") of CLPath"
	    & decimal(NCL)
	    & " withpen ForePen; undraw subpath (0,"
	    & decimal(BackFrac)
	    & ") of CLPath"
	    & decimal(NCL)
	    & " withpen BackPen;";
	  getready( comm, LinFunc(1/theN) );
	  for i=2 upto theN-1:
	    comm:="draw subpath ("
	      & decimal(i-ForeFrac)
	      & ","
	      & decimal(i+ForeFrac)
	      & ") of CLPath"
	      & decimal(NCL)
	      & " withpen ForePen; undraw subpath ("
	      & decimal(i-BackFrac)
	      & ","
	      & decimal(i+BackFrac)
	      & ") of CLPath"
	      & decimal(NCL)
	      & " withpen BackPen;";
	    getready( comm, LinFunc(i/theN) );
	  endfor;
	  comm:="draw subpath ("
	    & decimal(theN-ForeFrac)
	    & ","
	    & decimal(theN)
	    & ") of CLPath"
	    & decimal(NCL)
	    & " withpen ForePen; undraw subpath ("
	    & decimal(theN-BackFrac)
	    & ","
	    & decimal(theN)
	    & ") of CLPath"
	    & decimal(NCL)
	    & " withpen BackPen;";
	  getready( comm, LinFunc(1) );
	fi
      endgroup
    enddef;

% The next allows you to draw any solid that has no vertices and that has 
% two, exactly two, cyclic edges. In fact, it doesn't need to be a solid. 
% In order to complete the drawing of this solid you have to choose one of
% the edges to be drawn immediatly afterwards.    
    
    def twocyclestogether( expr CycleA, CycleB ) =
        begingroup
	    save TheLengthOfA, TheLengthOfB, TheMargin, Leng, i;
	    save SubPathA, SubPathB, PolygonPath, FinalPath;
	    numeric TheLengthOfA, TheLengthOfB, TheMargin, Leng, i;
	    path SubPathA, SubPathB, PolygonPath, FinalPath;
	    TheMargin = 0.02;
	    TheLengthOfA = ( length CycleA ) - TheMargin;
	    TheLengthOfB = ( length CycleB ) - TheMargin;
	    SubPathA = subpath ( 0, TheLengthOfA ) of CycleA;
	    SubPathB = subpath ( 0, TheLengthOfB ) of CycleB;
	    PolygonPath = makepath makepen ( SubPathA--SubPathB--cycle );
	    Leng = (length PolygonPath) - 1;
	    FinalPath = point 0 of PolygonPath
		        for i = 1 upto Leng:	  
			    --point i of PolygonPath
			endfor
			--cycle;
	    ( FinalPath )
        endgroup
    enddef;
    
% Ellipse on the air.

    def ellipticpath(expr CenterPos, OneAxe, OtherAxe ) =
      begingroup
	save ind;
	numeric ind;
	( for ind=0 upto 35: 
	    rp( CenterPos+planarrotation(OneAxe,OtherAxe,ind*10) )...
	  endfor cycle ) 
      endgroup
    enddef;

% Shadow of an ellipse on the air.

    def ellipticshadowpath(expr CenterPos, OneAxe, OtherAxe ) =
      begingroup
	save ind;
	numeric ind;
	( for ind=1 upto 36: 
	    rp( cb( CenterPos+planarrotation(OneAxe,OtherAxe,ind*10) ) )...
	  endfor cycle ) 
      endgroup
    enddef;

% It should be possible to attach some text to some plan.
% Unfortunately, this only works correctly when ParallelProj := true;
    
    def labelinspace(expr KeepRatio,RefPoi,BaseVec,UpVec)(text SomeString) =
      begingroup
	save labelpic, plak, lrc, ulc, llc, centerc, aratio, newbase;
	picture labelpic;
	pair lrc, ulc, llc;
	transform plak;
	color centerc, newbase;
	numeric aratio;
	labelpic = thelabel( SomeString, origin );
	lrc = lrcorner labelpic;
	ulc = ulcorner labelpic;
	llc = llcorner labelpic;
	aratio = (xpart  lrc - xpart llc)/(ypart ulc - ypart llc);
	if KeepRatio:
	  newbase = conorm(UpVec)*aratio*N(BaseVec);
	else:
	  newbase = BaseVec;
	fi;
	rp(RefPoi+newbase) = lrc transformed plak;
	rp(RefPoi+UpVec) = ulc transformed plak;
	centerc = RefPoi+0.5(newbase+UpVec);
	rp(RefPoi) = llc transformed plak;
	label( labelpic transformed plak, rp(centerc) )
      endgroup
    enddef;

% It should be possible to attach some path to some surface.
    
    def closedpathinspace( expr SomeTDPath, NDivide )( text ConverterFunc ) =
      begingroup
	save i, outpath, st;
	numeric i, st;
	path outpath;
	st = 1/NDivide;
	outpath = for i=st step st until (length SomeTDPath):
	  ConverterFunc( point i of SomeTDPath ) --
	endfor cycle;
	( outpath )
      endgroup
    enddef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Standard Objects:
	    
% And more precisely. The next routines spatialhalfcircle and 
% rigorousfear require circles drawn in a systematic and precise way.

    def goodcirclepath(expr CenterPos, AngulMom, Radius ) =
        begingroup
            save cirath, vecx, vecy, ind, goodangulmom, decision, view;
            numeric ind, decision;
            color vecx, vecy, goodangulmom, view;
            path cirath;
	    view = f-CenterPos;
	    decision = cdotprod( view, AngulMom );
	    if decision < 0: 
	        goodangulmom = -AngulMom;
	    else:
		goodangulmom = AngulMom;
	    fi;
            vecx = Radius*ncrossprod( view, goodangulmom );
	    decision := getangle( view, goodangulmom );
	    if decision > 0.5:                  %%%%%%%%%%%%%%% DANGER %%%
                vecy = Radius*ncrossprod( goodangulmom, vecx );
		cirath = for ind=10 step 10 until 360: 
			     rp( CenterPos + planarrotation(vecx,vecy,ind) )...
			 endfor cycle;
	    else:
	      cirath = head_on_circle( CenterPos, Radius );
	    fi;
            ( cirath ) 
        endgroup
    enddef;

% And its shadow. 

    def circleshadowpath(expr CenterPos, AngulMom, Radius ) =
      begingroup
	save cirath, vecx, vecy, view, decision;
	numeric decision;
	color vecx, vecy, view;
	path cirath;
	view = LightSource-CenterPos;
	vecx = ncrossprod( view, AngulMom );
	decision := getangle( view, AngulMom );
	if decision > 0.5:                  %%%%%%%%%%%%%%% DANGER %%%
	  vecy = ncrossprod( AngulMom, vecx );
	  cirath = ellipticshadowpath(CenterPos,vecx*Radius,vecy*Radius);
	else:
	  vecx := N( (-Y(view), X(view), 0) );
	  vecy = ncrossprod( view, vecx );
	  cirath = ellipticshadowpath(CenterPos,vecx*Radius,vecy*Radius);
	fi;
	( cirath ) 
      endgroup
    enddef;

% When there are numerical problems with the previous routine
% use the following alternative:

    def head_on_circle(expr Pos, Radius ) =
        begingroup
            save vecx, vecy, ind, view;
            numeric ind;
            color vecx, vecy, view;
	    view = f-Pos;
            vecx = Radius*N( (-Y(view), X(view), 0) );
            vecy = Radius*ncrossprod( view, vecx );
            ( for ind=10 step 10 until 360: 
		rp( Pos + planarrotation(vecx,vecy,ind) )...
	    endfor cycle )
        endgroup
    enddef;

% The nearest or the furthest part of a circle returned as a path.
% This function has been set to work for rigorousdisc (next).
% Very tough settings they were.

    def spatialhalfcircle(expr Center, AngulMom, Radius, ItsTheNearest ) =
      begingroup
	save va, vb, vc, cc, vd, ux, uy, pa, pb;
	save nr, cn, valx, valy, valr, choiceang;
	save auxil, auxih, fcirc, returnp;
	save choice;
	color va, vb, vc, cc, vd, ux, uy, pa, pb;
	numeric nr, cn, valx, valy, valr, choiceang;
	path auxil, auxih, fcirc, returnp;
	boolean choice;
	va := Center - f; 
	vb := N( AngulMom ); 
	vc := vb*( cdotprod( va, vb ) );
	cc := f + vc;
	vd := cc - Center;  % vd := va + vc;
	nr := conorm( vd );
	if Radius >= nr:
	  returnp := rp( cc ); % this single point will show up in spheroid
	else:
	  valr := Radius*Radius;
	  valx := valr/nr;
	  valy := sqrt( valr - valx*valx ); 
	  ux := N( vd );
	  choiceang := getangle( vc, va );                   %%%%%%%%%%%%%
	  choice := ( choiceang < 89 ) or ( choiceang > 91 );%%  DANGER  %
	  if choice:                                         %%%%%%%%%%%%%
	    uy := ncrossprod( vc, va );
	  else:
	    uy := ncrossprod( AngulMom, va );
	  fi;
	  pa := valx*ux + valy*uy + Center;
	  pb := pa - 2*valy*uy;
	  if choice:
	    auxil := rp(1.1[Center,pb])--rp(0.9[Center,pb]);
	    auxih := rp(1.1[Center,pa])--rp(0.9[Center,pa]);
	    fcirc := goodcirclepath( Center, AngulMom, Radius );
	    if ItsTheNearest:
	      returnp := (fcirc cutafter auxih) cutbefore auxil;
	    else:
	      returnp := (fcirc cutbefore auxih)..(fcirc cutafter auxil);
	    fi;
	  else:
	    if ItsTheNearest:
	      if cdotprod( va, AngulMom ) > 0:
		returnp := rp(pb)--rp(pa);
	      else:
		returnp := rp(pa)--rp(pb);
	      fi;
	    else:
	      if cdotprod( va, AngulMom ) < 0:
		returnp := rp(pb)--rp(pa);
	      else:
		returnp := rp(pa)--rp(pb);
	      fi;
	    fi;
	  fi;
	fi;
	( returnp )
      endgroup
    enddef;

% Cylinders or tubes ( numeric, boolean, color, numeric, color ).
% Great stuff. The "disc" in the name comes from the fact that
% when SphericalDistortion := true; the sides of cylinders are
% not drawn correctly (they are straight). And when it is a tube
% you should force the background to be white.
 
    def rigorousdisc(expr InRay, FullFill, BaseCenter, Radius, LenVec) =
      begingroup
        save va, vb, vc, cc, vd, base, holepic;
	save vA, cC, nr, vala, valb, hashole, istube;
	save auxil, auxih, halfl, halfh, thehole;
	save auxili, auxihi, rect, theshadow;

        color va, vb, vc, cc, vd, base;
	picture holepic;
	color vA, cC;
        numeric nr, vala, valb;
	boolean hashole, istube;
	path auxil, auxih, halfl, halfh, thehole;
	path auxili, auxihi, rect, theshadow;

	va := BaseCenter - f; 
        vb := N( LenVec ); 
        vc := vb*( cdotprod( va, vb ) );
        cc := f + vc;
	vd := cc - BaseCenter;
        nr := conorm( vd );
	base := BaseCenter + LenVec;
	vA := base - f;
	vala := conorm( va );
	valb := conorm( vA );
	if ShadowOn:
	  auxil := circleshadowpath( BaseCenter, LenVec, Radius );
	  auxih := circleshadowpath( base, LenVec, Radius );
	  fill twocyclestogether( auxil, auxih );
	fi;
	auxil := goodcirclepath( base, LenVec, Radius );
	auxih := goodcirclepath( BaseCenter, LenVec, Radius );
	istube := false;
	hashole := false;
	if InRay > 0:
	  istube := true;
	  auxili := goodcirclepath( base, LenVec, InRay );
	  auxihi := goodcirclepath( BaseCenter, LenVec, InRay );
	  hashole := (-1,-1) <> ( auxili intersectiontimes auxihi );
	  if hashole:
	    draw auxili;
	    draw auxihi;
	    holepic := currentpicture;
	    clip holepic to auxili;
	    clip holepic to auxihi;
	  fi;
	fi;
			 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if Radius >= nr: % THE CASE Radius > nr > InRay IS NOT SUPPORTED %
          if vala <= valb :
	    thehole := auxil;
	    auxil := auxih;
	    auxih := thehole;
	  fi;
	  if istube:
            if vala <= valb :
	      thehole := auxili;
	      auxili := auxihi;
	      auxihi := thehole;
	    fi;
	    holepic := currentpicture;
	    clip holepic to auxihi;
	  fi;
	  unfill auxil;
	  draw auxil;
	  if istube:
	    draw holepic;
	    draw auxihi;
	  fi;
	else:
	  cC := base + vd;
          if ( cdotprod( f - cc, f - cC ) <= 0 ) or ( not FullFill ):
	    halfl := spatialhalfcircle(BaseCenter,LenVec,Radius,true);
	    halfh := spatialhalfcircle(base,LenVec,Radius,true);
	    if FullFill:
	      rect := halfl--halfh--cycle;
	    else:
	      rect := halfl--(reverse halfh)--cycle;
	    fi;
	    unfill rect;
	    draw rect;
	  elseif vala > valb:
	    halfl := spatialhalfcircle(BaseCenter,LenVec,Radius,true);
	    halfh := spatialhalfcircle(base,LenVec,Radius,false);
	    rect := halfl--halfh--cycle;
	    unfill rect;
	    draw rect;
	    if istube:
	      if hashole:
	        draw holepic;
	      fi;
	      draw auxili;
	    fi;
	    draw auxil;
	  else:
	    halfl := spatialhalfcircle(BaseCenter,LenVec,Radius,false);
	    halfh := spatialhalfcircle(base,LenVec,Radius,true);
	    rect := halfl--halfh--cycle;
	    unfill rect;
	    draw rect;
	    if istube:
	      if hashole:
	        draw holepic;
	      fi;
	      draw auxihi;
	    fi;
	    draw auxih;
	  fi;
	fi
      endgroup
    enddef;

% And maybe a full cone border. The vertex may go anywhere.
% Choose the full cone border (UsualForm=true) or just the nearest
% part of the base edge (UsualForm=false).
% This is used by tropicalglobe as a generic spatialhalfcircle to
% draw only the in fact visible part of circular lines. Please, don't
% put the vertex too close to the base plan when UsualForm=false.
 
    def rigorouscone(expr UsualForm,CenterPos,AngulMom,Radius,VertexPos) =
      begingroup
	save basepath, thesubpath, fullpath, finalpath, auxpath, bigcirc;
	save themargin, newlen, i, auxt, startt, endt, thelengthofc;
	save pa, pb, pc, pd, pe;
	path basepath, thesubpath, fullpath, finalpath, auxpath;
	path bigcirc;
	numeric themargin, newlen, i, auxt, startt, endt, thelengthofc;
	pair pa, pb, pc, pd, pe;
	themargin = 0.02; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER
	basepath = goodcirclepath( CenterPos, AngulMom, Radius );
	thelengthofc = ( length basepath ) - themargin;
	thesubpath = subpath ( 0, thelengthofc ) of basepath;
	fullpath = makepath makepen ( rp(VertexPos)--thesubpath--cycle );
	pa = 0.995[rp(CenterPos),rp(VertexPos)];
	pb = 1.005[rp(CenterPos),rp(VertexPos)];
	auxpath = pa--pb;
	pc = auxpath intersectiontimes fullpath;
	if pc <> (-1,-1):
	  auxt = ypart pc;
	  newlen = length fullpath;
	  if UsualForm:
	    finalpath = point auxt of fullpath
	      --point auxt+1 of fullpath
	      for i = auxt+2 upto auxt+newlen-1:	  
	        ...point i of fullpath
	      endfor
	      --cycle;
	  else:
	    bigcirc = goodcirclepath( CenterPos, AngulMom, 1.005*Radius );
	    pd = bigcirc intersectiontimes fullpath;
	    pe = ( reverse bigcirc ) intersectiontimes fullpath;
	    startt = floor( xpart pd );
	    endt = ceiling( ( length bigcirc ) - ( xpart pe ) );
	    finalpath = subpath (startt,endt) of basepath;
	  fi;
	else:
	  finalpath = rp(VertexPos);
	fi;
	( finalpath )
      endgroup
    enddef;
    
    def verygoodcone(expr BackDash,CenterPos,AngulMom,Radius,VertexPos) =
      begingroup
	save thepath, lenpath, bonevec, sidevec, viewaxe, cipath;
	save thelengthofc, thesubpath, themargin, basepath;
	color bonevec, sidevec, viewaxe;
	path thepath, cipath, basepath, thesubpath;
	numeric lenpath, thelengthofc, themargin;
	themargin = 0.02;
	bonevec = VertexPos - CenterPos;
	if cdotprod( bonevec, AngulMom ) < 0:
	  sidevec = -N(AngulMom);
	else:
	  sidevec = N(AngulMom);
	fi;
	viewaxe = f-CenterPos;
	if ShadowOn:
	  basepath = circleshadowpath( CenterPos, AngulMom, Radius );
	  thelengthofc = ( length basepath ) - themargin;
	  thesubpath = subpath ( 0, thelengthofc ) of basepath;
	  fill makepath makepen ( rp(cb(VertexPos))--thesubpath--cycle );
	fi;
	thepath = rigorouscone(true,CenterPos,AngulMom,Radius,VertexPos);
	lenpath = length thepath;
	if lenpath<>0:
	  unfill thepath;
	  draw thepath;
	  if cdotprod( sidevec, viewaxe ) < 0:
	    draw goodcirclepath( CenterPos, AngulMom, Radius );
	  else:
	    if BackDash:
	      draw
	      goodcirclepath( CenterPos, AngulMom, Radius ) dashed evenly;
	    fi;
	  fi;
	else:
	  cipath = goodcirclepath( CenterPos, AngulMom, Radius );
	  unfill cipath;
	  draw cipath;
	  if cdotprod( sidevec, viewaxe ) > 0:
	    draw rp( VertexPos );
	  fi;	  
	fi
      endgroup
    enddef;

% Its a sphere, don't fear, but remember that the rigorous projection
% of a sphere is an ellipse. 

    def rigorousfearpath(expr Center, Radius ) =
        begingroup
	    save auxil, ux, uy, newcen, nr, valx, valy, valr;
            color ux, uy, newcen;
            numeric nr, valx, valy, valr;
	    path auxil;
            nr := conorm( Center - f );
	    valr := Radius**2;
	    valx := valr/nr;
	    valy := sqrt( valr - valx**2 );
	    newcen := valx*( f - Center )/nr;
	    auxil := head_on_circle( Center+newcen, valy );
	    ( auxil )
        endgroup
    enddef;

    def rigorousfearshadowpath(expr Center, Radius ) =
        begingroup
            save ux, uy, newcen;
            save nr, valx, valy, valr, lenr;
	    save auxil, auxih, fcirc, returnp;
	    save dcenter;
            color ux, uy, newcen;
            numeric nr, valx, valy, valr, lenr;
	    path auxil, auxih, fcirc, returnp;
	    pair dcenter;
            nr := conorm( Center - LightSource );
	    valr := Radius**2;
	    valx := valr/nr;
	    valy := sqrt( valr - valx**2 );
	    newcen := valx*( LightSource - Center )/nr;
	    auxil := circleshadowpath( Center+newcen, newcen, valy );
	    ( auxil )
        endgroup
    enddef;

% It's a globe (without land).

    def tropicalglobe( expr NumLats, TheCenter, Radius, AngulMom ) =
      begingroup
	save viewaxe, sinalfa, sinbeta, globaxe, aux, limicos, lc;
	save stepang, actang, newradius, foc, newcenter, cpath, i;
	save outerpath, conditiona, conditionb;
	color viewaxe, globaxe, foc, newcenter;
	numeric sinalfa, sinbeta, aux, limicos, stepang, actang;
	numeric newradius, lc, i;
	path cpath, outerpath;
	boolean conditiona, conditionb;
	if ShadowOn:
	  fill rigorousfearshadowpath( TheCenter, Radius );
	fi;
	viewaxe = f-TheCenter;
	sinalfa = Radius/conorm( viewaxe );
	aux = cdotprod( viewaxe, AngulMom );
	if aux < 0:
	  globaxe = -N(AngulMom);
	else:
	  globaxe = N(AngulMom);
	fi;
	sinbeta = cdotprod( globaxe, N(viewaxe) );
	aux := sqrt((1-sinalfa**2)*(1-sinbeta**2));
	limicos = aux - sinalfa*sinbeta;
	stepang = 180/NumLats;
	globaxe := globaxe*Radius;
	outerpath = rigorousfearpath(TheCenter,Radius);
	unfill outerpath;
	draw outerpath;
	for actang = 0.5*stepang step stepang until 179:
	  if cosd(actang) < limicos-0.005: %%%%%%%%%%%%%%%%%%%%%%%% DANGER
	    newradius := Radius*sind(actang);
	    newcenter := TheCenter - globaxe*cosd(actang);
	    conditiona:=(actang<94) and (actang>86); % DANGER % DANGER VV
	    conditionb:=abs(cdotprod(globaxe/Radius,N(f-newcenter)))<0.08;
	    if conditiona or conditionb:
	      draw spatialhalfcircle(newcenter,globaxe,newradius,true);
	    else:
	      foc := TheCenter - globaxe/cosd(actang);
	      lena := -Radius*cosd(actang);
	      lenb := cdotprod(viewaxe,globaxe/Radius);
	      if (actang <= 86) or ((lenb<lena) and (actang>=94)):
		cpath :=
		rigorouscone(false,newcenter,globaxe,newradius,foc);
		draw cpath;
	      else:
		cpath :=
		rigorouscone(true,newcenter,globaxe,newradius,foc);
		lc := length cpath;
		if lc <> 0:
		  draw subpath (1,lc-1) of cpath;
		else:
		  draw rigorouscircle( newcenter,globaxe,newradius );
		fi;
	      fi;
	    fi;
	  fi;
	endfor
      endgroup
    enddef;

% An elliptical frustum:
    
    def whatisthis(expr CenterPos,OneAxe,OtherAxe,CentersDist,TheFactor) =
      begingroup
	save patha, pathb, pathc, centersvec, noption;
	path patha, pathb, pathc;
	color centersvec;
	numeric noption;
	centersvec = CentersDist*ncrossprod( OneAxe, OtherAxe );
	if ShadowOn:
	  patha = ellipticshadowpath( CenterPos,
	    OneAxe,
	    OtherAxe );
	  pathb = ellipticshadowpath( CenterPos+centersvec,
	    TheFactor*OneAxe,
	    TheFactor*OtherAxe );
	  pathc = twocyclestogether( patha, pathb );
	  fill pathc;
	fi;
	patha := ellipticpath( CenterPos,
	  OneAxe,
	  OtherAxe );
	pathb := ellipticpath( CenterPos+centersvec,
	  TheFactor*OneAxe,
	  TheFactor*OtherAxe );
	pathc := twocyclestogether( patha, pathb );
	unfill pathc;
	draw pathc;
	noption = cdotprod( centersvec, f-CenterPos );
	if  noption > (CentersDist**2):
	  draw pathb;
	elseif noption < 0:
	  draw patha;
	fi
      endgroup
    enddef;

% Probably the last algorithm I'm going to write for featpost...
% Wrong. The last is ultraimprovertex. 
% Wrong again. The last is necplusimprovertex.
% And again wrong. The last is tdcircarrow.
% Well, what can I say, really the last is ellipsoid.
% Wait! It is torushadow!
% Bullshit. Only death can stop me. Now it's revolparab.
% Continuing with torus accessories...
% And I forgot to mention intersectprolatespheroid!!!
% The first of Red October, 2014, added vp.
    
    def spheroidshadow( expr CentrPoi, NorthPoleVec, Ray ) =
      begingroup
	save a, k, fx, fy, tmpa, tmpb, tmpc, ep, vax, wax, xax, yax, zax, cm, cp;
	save bdh, bdm, bdp, bdv, sm, sp, i, cac;
	numeric a, k, fx, fy, tmpa, tmpb, tmpc, cm, cp, sm, sp, i;
	path ep;
	color vax, wax, xax, yax, zax, cac;
	pair bdh, bdm, bdp, bdv;
	vax = LightSource-CentrPoi;
	if cdotprod(NorthPoleVec,vax)<0:
	  xax = -N(NorthPoleVec);
	else:
	  xax = N(NorthPoleVec);
	fi;
	a = conorm(NorthPoleVec);
	k = a/Ray;
	if getangle(xax,vax) > 0.5:               %%%%%%%%%%%%%%% DANGER %%%
	  zax = ncrossprod(xax,vax);
	else:
	  zax = N( ( 0, Z(vax), -Y(vax) ) );
	fi;
	yax = ncrossprod(zax,xax);
	fx = cdotprod(vax,xax);
	fy = cdotprod(vax,yax);
	tmpa = Ray*fx/k;
	tmpb = fy*(fy ++ (fx/k) +-+ Ray);
	tmpc = ((fx/k)**2)+(fy**2);
	cm = (tmpa-tmpb)/tmpc;
	cp = (tmpa+tmpb)/tmpc;
	sm = 1 +-+ cm;
	if fx<a:
	  sp = 1 +-+ cp;
	else:
	  sp = -( 1 +-+ cp );
	fi;
	bdm = (k*cm,sm)*Ray;
	bdp = (k*cp,sp)*Ray;
	bdh = 0.5[bdp,bdm];
	tmpc := Ray*( 1 +-+ ((xpart bdh)/a) );
	tmpb := tmpc +-+ (ypart bdh);
	bdv = bdm-bdp;
	wax = 0.5*( xax*( xpart (bdv) ) + yax*( ypart (bdv) ) );
	cac = CentrPoi+ xax*( xpart (bdh) ) + yax*( ypart (bdh) );
	fill ellipticshadowpath( cac, wax, zax*tmpb );
      endgroup
    enddef;

    def spheroid( expr CentrPoi, NorthPoleVec, Ray ) =
      begingroup
	save a, k, fx, fy, tmpa, tmpb, tmpc, ep;
	save vax, wax, xax, yax, zax, cm, cp;
	save bdh, bdm, bdp, bdv, sm, sp, i, cac;
	numeric a, k, fx, fy, tmpa, tmpb, tmpc, cm, cp, sm, sp, i;
	path ep;
	color vax, wax, xax, yax, zax, cac;
	pair bdh, bdm, bdp, bdv;
	if ShadowOn:
	  spheroidshadow( CentrPoi, NorthPoleVec, Ray );
	fi;
	vax = f-CentrPoi;
	if cdotprod(NorthPoleVec,vax)<0:
	  xax = -N(NorthPoleVec);
	else:
	  xax = N(NorthPoleVec);
	fi;
	a = conorm(NorthPoleVec);
	k = a/Ray;
	if getangle(xax,vax) > 0.5:               %%%%%%%%%%%%%%% DANGER %%%
	  zax = ncrossprod(xax,vax);
	else:
	  zax = N( (-Y(vax), X(vax), 0) );
	fi;
	yax = ncrossprod(zax,xax);
	fx = cdotprod(vax,xax);
	fy = cdotprod(vax,yax);
	tmpa = Ray*fx/k;
	tmpb = fy*(fy ++ (fx/k) +-+ Ray);
	tmpc = ((fx/k)**2)+(fy**2);
	cm = (tmpa-tmpb)/tmpc;
	cp = (tmpa+tmpb)/tmpc;
	sm = 1 +-+ cm;
	if fx<a:
	  sp = 1 +-+ cp;
	else:
	  sp = -( 1 +-+ cp );
	fi;
	bdm = (k*cm,sm)*Ray;
	bdp = (k*cp,sp)*Ray;
	bdh = 0.5[bdp,bdm];
	tmpc := Ray*( 1 +-+ ((xpart bdh)/a) );
	tmpb := tmpc +-+ (ypart bdh);
	bdv = bdm-bdp;
	wax = 0.5*( xax*( xpart (bdv) ) + yax*( ypart (bdv) ) );
	cac = CentrPoi+ xax*( xpart (bdh) ) + yax*( ypart (bdh) );
	ep = ellipticpath( cac, wax, zax*tmpb );
	unfill ep;
	draw ep;
	draw spatialhalfcircle( CentrPoi, NorthPoleVec, Ray, true );
      endgroup
    enddef;

% Another brute-force algorythm.
% It's advisable to use three orthogonal axes.

    def ellipsoid( expr Centr, AxOne, AxTwo, AxThr ) =
      begingroup
	save count, i, j, axx, axy, cyc, cy, di, leng;
	numeric count, i, j, leng;
	color axx, axy, di[];
	path cy, cyc[];
	di1 = AxOne;
	di2 = AxTwo;
	di3 = AxThr;
	count = 0;
	for i=1 upto 3:
	  if i=1:
	    axx := di2;
	    axy := di3;
	  elseif i=2:
	    axx := di1;
	    axy := di3;
	  else:
	    axx := di1;
	    axy := di2;
	  fi;
	  for j=5 step 10 until 175:
	    cyc[incr(count)] = ellipticpath( Centr, di[i],
		Centr + planarrotation( axx, axy, j ) );
	  endfor;
	endfor;
	cy = cyc1;
	for i=2 upto count-1:
	  cy := twocyclestogether( cy, cyc[i] );
	endfor;
	leng = (length cy) - 1;
	( point 0 of cy for i = 1 upto leng: ..point i of cy endfor ..cycle )
      endgroup
    enddef;
  
    def ellipsoidshadow( expr Centr, AxOne, AxTwo, AxThr ) =
      begingroup
	save count, i, j, axx, axy, cyc, cy, di, leng;
	numeric count, i, j, leng;
	color axx, axy, di[];
	path cy, cyc[];
	di1 = AxOne;
	di2 = AxTwo;
	di3 = AxThr;
	count = 0;
	for i=1 upto 3:
	  if i=1:
	    axx := di2;
	    axy := di3;
	  elseif i=2:
	    axx := di1;
	    axy := di3;
	  else:
	    axx := di1;
	    axy := di2;
	  fi;
	  for j=5 step 10 until 175:
	    cyc[incr(count)] = ellipticshadowpath( Centr, di[i],
		Centr + planarrotation( axx, axy, j ) );
	  endfor;
	endfor;
	cy = cyc1;
	for i=2 upto count-1:
	  cy := twocyclestogether( cy, cyc[i] );
	endfor;
	leng = (length cy) - 1;
	( point 0 of cy for i = 1 upto leng: ..point i of cy endfor ..cycle )
      endgroup
    enddef;
  
    def revolparab( expr BaseCenter, ParabTip, BaseRay ) =
      begingroup
	save bcpt, conetip, fakex, fakey, fakez, tipview, ellicenter, coneview;
	save tanefe, majorvec, minorvec, cutvec, auxpoi, crux, auxx, baseview;
	save xzero, yzero, xdelta, fakea, xpos, ypos, xneg, yneg, l, apertur;
	save auxy, maxy, ellmaxang, ymin, xefe, auxray, auxcos, auxsin, ste, a;
	save auxpath, tippath;
	save conda, condb, condc;
	color bcpt, conetip, fakex, fakey, fakez, tipview, ellicenter, coneview;
	color tanefe, majorvec, minorvec, cutvec, auxpoi, crux, auxx, baseview;
	numeric xzero, yzero, xdelta, fakea, xpos, ypos, xneg, yneg, l, apertur;
	numeric auxy, maxy, ellmaxang, ymin, xefe, auxray, auxcos, auxsin, ste;
	numeric a;
	path auxpath, tippath;
	boolean conda, condb, condc;
	bcpt = BaseCenter-ParabTip;
	conetip = BaseCenter-2*bcpt;
	maxy = conorm(bcpt);
	apertur = angle(2*maxy,BaseRay);
	fakey = N(bcpt);
	coneview = f-conetip;
	tipview = f-ParabTip;
	baseview = f-BaseCenter;
	a = getangle(coneview,bcpt);
	conda = a<=apertur;
	condb = a>=180-apertur;
	fakez = ncrossprod( tipview, bcpt );
	fakex = ccrossprod( fakey, fakez );
	xzero = cdotprod( fakex, tipview );
	yzero = cdotprod( fakey, tipview );
	fakea = maxy/(BaseRay**2);
	condc = yzero>=fakea*(xzero**2);
	if (conda and condc) or condb:
	  auxpath = rigorouscircle( BaseCenter, bcpt, BaseRay );
	  unfill auxpath;
	  draw auxpath;
	else:
	  xdelta = sqrt(xzero**2 - yzero/fakea);
	  xpos = xzero + xdelta;
	  xneg = xzero - xdelta;
	  ypos = fakea*(xpos**2);
	  yneg = fakea*(xneg**2);
	  auxy = 0.5[ypos,yneg];
	  ellicenter = ParabTip+fakex*xzero+fakey*auxy;
	  if yneg<ypos:
	    ymin=yneg;
	    xefe=xneg;
	  else:
	    ymin=ypos;
	    xefe=xpos;
	  fi;
	  tanefe = ParabTip +fakex*xefe+fakey*ymin;
	  minorvec = xdelta*fakez;
	  majorvec = tanefe-ellicenter;
	  if (yneg<maxy) and (ypos<maxy):
	    draw ellipticpath( ellicenter, majorvec, minorvec );
	  else:
	    auxpoi = lineintersectplan(conetip,conetip-f,BaseCenter,bcpt);
	    auxray = conorm( BaseCenter-auxpoi );
	    auxcos = BaseRay/auxray;
	    auxsin = 1 +-+ auxcos;
	    if cdotprod(fakex,auxpoi-BaseCenter)>0:
	      auxx = fakex;
	    else:
	      auxx = -fakex;
	    fi;
	    crux = BaseCenter+BaseRay*(auxx*auxcos +fakez*auxsin);
	    cutvec = crux-ellicenter;
	    auxcos := cdotprod(cutvec,N(majorvec))/conorm(majorvec);
	    auxsin := cdotprod(cutvec,N(minorvec))/conorm(minorvec);      
	    ellmaxang = angle(auxcos,auxsin);
	    ste = ellmaxang/18; %%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER %%%%%%%%%%%%
	    tippath =
	      rp(ellicenter+planarrotation(majorvec,minorvec,-ellmaxang)) 
	      for l=ste-ellmaxang step ste until ellmaxang:
	      ..rp(ellicenter+planarrotation(majorvec,minorvec,l))
	    endfor --rp(ellicenter+planarrotation(majorvec,minorvec,ellmaxang));
	    if cdotprod( baseview, bcpt )<0:
	      auxpath = rigorouscone(true,BaseCenter,bcpt,BaseRay,conetip);
	      numeric len;
	      len = length auxpath;
	      auxpath := subpath (1,len-1) of auxpath;
	      auxpath := tippath--(reverse auxpath)--cycle;
	      unfill auxpath;
	      draw auxpath;
	    else:
	      tippath := tippath--cycle;
	      unfill tippath;
	      draw tippath;
	      auxpath = rigorouscircle( BaseCenter, bcpt, BaseRay );
	      unfill auxpath;
	      draw auxpath;
	    fi;
	  fi;
	fi
      endgroup
    enddef;
    
% You can't see through this hole. f must not be on the hole axis.
% Not yet documented because "buildcycle" doesn't work properly.
    
    def fakehole( expr CenterPos, LenVec, Radius ) =
      begingroup
	save patha, pathb, pathc, noption, hashole, auxv, poption, vv;
	save ta, tb, taf, tbf, margint, stopair, pa, pb, testpath, isin;
	path patha, pathb, pathc, pa, pb, testpath;
	numeric noption, ta, tb, margint;
	boolean hashole, poption, isin;
	color auxv, vv;
	pair stopair;
	vv = f-CenterPos;
	patha := rigorouscircle( CenterPos, LenVec, Radius );
	pathb := rigorouscircle( CenterPos+LenVec, LenVec, Radius );
% 	patha := goodcirclepath( CenterPos, LenVec, Radius );
% 	pathb := goodcirclepath( CenterPos+LenVec, LenVec, Radius );
	auxv := ncrossprod( LenVec, ccrossprod( vv, LenVec ) );
	poption := abs( cdotprod( vv, auxv ) ) <= 1.05*Radius;% DANGER!
	if poption:
	  draw patha;
	  draw pathb;
	else:
	  noption = cdotprod( LenVec, vv );
	  if  noption > (conorm(LenVec)**2):
	    pa = patha;
	    pb = pathb;
	  elseif noption < 0:
	    pa = pathb;
	    pb = patha;
	  fi;
	  draw pb;
	  stopair = pa intersectiontimes pb;
	  hashole = (-1,-1) <> stopair;
	  if hashole:
 	    testpath = rp(CenterPos+0.5*LenVec)--(point 0 of pa);
 	    isin = (-1,-1) <> testpath intersectiontimes pb;
 	    if not isin:
 	      ta = xpart stopair;
 	      tb = ypart stopair;
 	      stopair := (reverse pa) intersectiontimes (reverse pb);
 	      taf = length pa - xpart stopair;
 	      tbf = length pb - ypart stopair;
 	      margint = 0.01;                                     % DANGER!
 	      draw (subpath (0,ta-margint) of pa)--
 	      (subpath (tb+margint,tbf-margint) of pb)--
 	      (subpath (taf+margint,length pa - margint) of pa)--
 	      cycle;
 	    else:
	      pathc := buildcycle( pa, pb ); % I don't get it!
	      % Why doesn't buildcycle work all the time??? See fakehole.mp
	      % When point 0 of pa is inside pb, builcycle doesn't work!!
	      draw pathc;
	    fi;
	  fi;
	fi
      endgroup
    enddef;
      
% It is time for a kind of cube. Don't use SphericalDistortion here.
    
    def kindofcube(expr WithDash, IsVertex, RefP, AngA, AngB, AngC, LenA, LenB, LenC ) =
      begingroup
	save star, pos, patw, patb, refv, near, centre, farv;
	save newa, newb, newc, veca, vecb, vecc, auxx, auxy, i;
	color star, pos[], refv, near, newa, newb, newc;
	color veca, vecb, vecc, auxx, auxy, centre, farv;
	path patw, patb;
	numeric i;
	veca = ( cosd(AngA)*cosd(AngB),
	  sind(AngA)*cosd(AngB),
	  sind(AngB) );
	auxx = ( cosd(AngA+90), sind(AngA+90), 0 );
	auxy = ccrossprod( veca, auxx );
	vecb = cosd(AngC)*auxx + sind(AngC)*auxy;
	vecc = cosd(AngC+90)*auxx + sind(AngC+90)*auxy;
	veca := LenA*veca;
	vecb := LenB*vecb;
	vecc := LenC*vecc;
	if IsVertex:
	  star = RefP;
	  centre = RefP + 0.5*( veca + vecb + vecc);
	else:
	  star = RefP - 0.5*( veca + vecb + vecc);
	  centre = RefP;
	fi;
	pos1 = star + veca;
	pos2 = pos1 + vecb;
	pos3 = pos2 + vecc;
	pos4 = pos3 - vecb;
	pos5 = pos4 - veca;
	pos6 = pos5 + vecb;
	pos7 = pos6 - vecc;
	if ShadowOn:
	  patw = rp(cb(star))--rp(cb(pos1))--rp(cb(pos2))
	    --rp(cb(pos3))--rp(cb(pos4))
	    --rp(cb(pos5))--rp(cb(pos6))--rp(cb(pos7))--cycle;
	  patb = makepath makepen patw;
	  fill patb;
	fi;  
	patw := rp(star)--rp(pos1)--rp(pos2)--rp(pos3)--rp(pos4)
	  --rp(pos5)--rp(pos6)--rp(pos7)--cycle;
	patb := makepath makepen patw;  
	unfill patb;
	draw patb;
	i = 0;
	if inplanarvolume( star, star+veca, f ): i := incr( i ); fi;
	if inplanarvolume( star, star+vecb, f ): i := incr( i ); fi;
	if inplanarvolume( star, star+vecc, f ): i := incr( i ); fi;
	if (i=2) and WithDash:
	  message "Unable to dash kindofcube " & cstr( RefP ) & ".";
	elseif i = 3:
	  message "f is inside kindofcube " & cstr( RefP ) & ".";
	else:  
	  refv = f - centre;
	  if cdotprod( refv, veca ) > 0:
	    newa = -veca;
	  else:
	    newa = veca;
	  fi;
	  if cdotprod( refv, vecb ) > 0:
	    newb = -vecb;
	  else:
	    newb = vecb;
	  fi;
	  if cdotprod( refv, vecc ) > 0:
	    newc = -vecc;
	  else:
	    newc = vecc;
	  fi;
	  near = centre - 0.5*( newa + newb + newc );
	  draw rp(near)--rp(near+newa);
	  draw rp(near)--rp(near+newb);
	  draw rp(near)--rp(near+newc);
	  if WithDash:
	    if i=1:
	      message "Unable to dash kindofcube " & cstr( RefP ) & ".";
	    else:
	      farv = centre + 0.5*( newa + newb + newc );
	      draw rp(farv)--rp(farv-newa) dashed evenly;
	      draw rp(farv)--rp(farv-newb) dashed evenly;
	      draw rp(farv)--rp(farv-newc) dashed evenly;
	    fi;
	  fi;
	fi
      endgroup
    enddef;
    
% It's a bit late now but the stage must be set.

    def setthestage( expr NumberOfSideSquares, SideSize ) =
      begingroup
	save i, j, squaresize, squarepath, ca, cb, cc, cd;
	numeric i, j, squaresize;
	path squarepath;
	color ca, cb, cc, cd;
	squaresize = SideSize/(2*NumberOfSideSquares-1);
	for i=-0.5*SideSize step 2*squaresize until 0.5*SideSize:
	  for j=-0.5*SideSize step 2*squaresize until 0.5*SideSize:
	    ca := (i,j,HoriZon);
	    cb := (i,j+squaresize,HoriZon);
	    cc := (i+squaresize,j+squaresize,HoriZon);
	    cd := (i+squaresize,j,HoriZon);
	    squarepath := rp(ca)--rp(cb)--rp(cc)--rp(cd)--cycle;
	    unfill squarepath;
	    draw squarepath;
	  endfor;
	endfor
      endgroup
    enddef;
    
    def setthearena( expr NumberOfDiameterCircles, ArenaDiameter ) =
      begingroup
	save i, j, circlesize, polar, currpos, phi, cpath;
	numeric i, j, circlesize, polar, phi;
	color currpos;
	path cpath;
 	circlesize = ArenaDiameter/NumberOfDiameterCircles;
	for i=0.5*ArenaDiameter step -circlesize until 0.4*circlesize:
	  polar := floor(6.28318*i/circlesize);
	  for j=1 upto polar:
	    phi := 360*j/polar;
	    currpos := i*(cosd(phi),sind(phi),0)+(0,0,HoriZon);
	    cpath := rigorouscircle( currpos, blue, 0.3*circlesize);
	    unfill cpath;
	    draw cpath;
	  endfor;
	endfor
      endgroup
    enddef;

% And a transparent dome. The angular momentum vector is supposed 
% to point from the concavity of the dome and into outer space.
% The pen can only be changed with a previous drawoptions().

    def spatialhalfsfear(expr Center, AngulMom, Radius ) =
      begingroup
	save spath, cpath, fpath, rpath, cutp;
	path spath, cpath, fpath, rpath, cutp;
	save ap, bp, cp, dp, cuti, cute, vp;    
	pair ap, bp, cp, dp, cuti, cute, vp;
	save auxcos, actcos, actsin, auxsin;
	numeric auxcos, actcos, actsin, auxsin;
	picture partoffear;
	save A, B;
	color A, B;
	spath = rigorousfearpath( Center, Radius );
	auxcos = getcossine( Center, Radius );
	actcos = cdotprod( N( f - Center ), N( AngulMom ) );
	actsin = sqrt(1-actcos**2);
	auxsin = sqrt(1-auxcos**2);
	if actsin <= auxcos:
	  if actcos >= 0:
	    cpath = goodcirclepath( Center, AngulMom, Radius );
	    draw cpath;
	  else:
	    draw spath;
	  fi;
	else:
	  fpath = spatialhalfcircle( Center, AngulMom, Radius, true );
	  rpath = spatialhalfcircle( Center, AngulMom, Radius, false );
	  cuti = point 0 of rpath;
	  cute = point ( length rpath ) of rpath;
	  ap = 1.1[cuti,cute];
	  bp = 1.1[cute,cuti];
	  partoffear = nullpicture;
	  addto partoffear doublepath spath;
	  A = ncrossprod( f-Center, ncrossprod( f-Center, AngulMom ) );
	  B = Center + 1.1*Radius*( auxcos*N( f-Center ) + auxsin*A );
	  vp = rp(B) - rp(Center);
	  cp = ap + vp;
	  dp = bp + vp;
	  cutp = ap--cp--dp--bp--cycle;
	  clip partoffear to cutp;
	  draw fpath;
	  draw partoffear;
	  if actcos >= 0:
	    draw rpath;
	  fi;
	fi
      endgroup
    enddef;
	      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Toroidal Stuff
    
% Take a donut. 

    def smoothtorus( expr Tcenter, Tmoment, Bray, Sray ) =
      begingroup
	save nearaxe, sideaxe, viewline, circlecenter, circlemoment;
	save ang, ind, i, anglim, angstep, cuspcond, holepic, coofrac;
	save cpath, apath, opath, ipath, wp, ep, refpair, distance, lr;
	color nearaxe, sideaxe, viewline, circlecenter, circlemoment;
	numeric ang, ind, i, anglim, angstep, distance, coofrac, lr;
	path cpath, apath, opath, ipath, wp, ep;
	pair outerp[], innerp[], refpair;
	boolean cuspcond;
	picture holepic;
	save tmoment;
	color tmoment;
	angstep= 4; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER!
	if ShadowOn:
	  torushadow( Tcenter, Tmoment, Bray, Sray );
	fi;
	viewline = f-Tcenter;
	if cdotprod( viewline, Tmoment ) < 0:
	  tmoment = -Tmoment;
	else:
	  tmoment = Tmoment;
	fi;
	refpair = unitvector( rp(Tcenter+tmoment)-rp(Tcenter) );
	sideaxe = Bray*ncrossprod( tmoment, viewline );
	nearaxe = Bray*ncrossprod( sideaxe, tmoment );
	coofrac = cdotprod( viewline, N( tmoment ) )/Sray;
	if (coofrac <= 1.04) and (coofrac >= 1.01): %%%%%%%%%%% DANGER!
	  ind = 360/angstep;
	  anglim = 0.5*angstep;
	  for i=1 upto ind:
	    ang := i*angstep-anglim-180.0;
	    circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter;
	    circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang);
	    cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true);
	    if i >= 0.5*ind+1:
	      outerp[i]=point 0 of cpath;
	    else:
	      outerp[i]=point (length cpath) of cpath;
	    fi;
	  endfor;
	  opath = for i=1 upto ind: outerp[i].. endfor cycle;
	  unfill opath;
	  draw opath;
	elseif coofrac < 1.01:  
	  distance = conorm( viewline );
	  lr = Bray + Sray*( 1 +-+ coofrac );
	  anglim = angle( ( lr, distance +-+ lr ) );
	  ind = 2*floor(anglim/angstep);
	  angstep := 2*anglim/(ind+1);
	  for i=0 upto 0.5*ind-1:
	    ang := i*angstep-anglim;
	    circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter;
	    circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang);
	    cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true);
	    innerp[i]=point 0 of cpath;
	    outerp[i]=point (length cpath) of cpath;
	  endfor;
	  for i=0.5*ind upto ind-2:
	    ang := (i+2)*angstep-anglim;
	    circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter;
	    circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang);
	    cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true);
	    outerp[i]=point 0 of cpath;
	    innerp[i]=point (length cpath) of cpath;
	  endfor;
	  if coofrac > 0.94:
	    apath = innerp0
	      for i=1 upto ind-2:
	        ..innerp[i]
	      endfor
	      --cycle;
	  else:
	    apath = innerp0 for i=2 upto ind-2: ..innerp[i] endfor
	    ..outerp[ind-2] for i=ind-3 downto 0: ..outerp[i] endfor
	    ..cycle;
	  fi;
	  unfill apath;
	  draw apath;
	else:
	  ind = 360/angstep;
	  anglim = 0.5*angstep;
	  for i=1 upto ind:
	    ang := i*angstep-anglim-180.0;
	    circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter;
	    circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang);
	    cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true);
	    if i >= 0.5*ind+1:
	      outerp[i]=point 0 of cpath;
	      innerp[i]=point (length cpath) of cpath;
	    else:
	      innerp[i]=point 0 of cpath;
	      outerp[i]=point (length cpath) of cpath;
	    fi;
	  endfor;
	  opath = for i=1 upto ind: outerp[i].. endfor cycle;
	  ipath = for i=1 upto ind: innerp[i].. endfor cycle;
	  holepic = currentpicture;
	  clip holepic to ipath;
	  unfill opath;
	  draw holepic;
	  draw opath;
	  draw ipath;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perhaps there is an analytic way of getting the angle of the cusp point?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	  i := ceiling(1+0.5*ind);
	  cuspcond = false;
	  forever:
	    i := incr( i );
	    exitif i > ind-1;
	    cuspcond :=
	      refpair dotprod innerp[i+1] < refpair dotprod innerp[i];
	    exitif cuspcond;
	  endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	  if cuspcond:
	    show "Torus shows cusp points.";
	    ep = outerp[ind-i+1]--innerp[ind-i+1];
	    wp = innerp[i]--outerp[i];
	    unfill buildcycle(reverse opath,ep,ipath,wp);
	    draw opath;
	    draw subpath (i-1,ind-i) of ipath;
	  fi;
	fi;
      endgroup
    enddef;

% The shadow of a torus
    
    def torushadow( expr Tcenter, Tmoment, Bray, Sray ) =
      begingroup
	save theplace, viewline, tmoment, refpair, sideaxe, nearaxe, i;
	color theplace, viewline, tmoment, refpair, sideaxe, nearaxe;
	numeric i;
	viewline = f-Tcenter;
	if cdotprod( viewline, Tmoment ) < 0:
	  tmoment = -Tmoment;
	else:
	  tmoment = Tmoment;
	fi;
	sideaxe = Bray*ncrossprod( tmoment, viewline );
	nearaxe = Bray*ncrossprod( sideaxe, tmoment );
	for i=1 upto 60:  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER!
	  theplace := Tcenter + planarrotation( sideaxe, nearaxe, 6i ); %%
	  fill rigorousfearshadowpath( theplace, Sray );
	endfor;
      endgroup
    enddef;
    
% Take a "quarter" of a donut
% (no longer under construction but contains a bug). 

    def quartertorus( expr Tcenter, Tstart, Tfinis, Sray ) =
      begingroup
	save sideaxe, viewline, circlecenter, circlemoment, amax;
	save i, angstep, tmoment, cpath, opath, ipath, fpath, finc;
	save vstart, vfinis, ostart, ofinis, cstart, cfinis, finis;
	color sideaxe, viewline, circlecenter, circlemoment;
	color tmoment, vstart, vfinis, finis, finc;
	numeric i, angstep, amax;
	path cpath, opath, ipath, cstart, cfinis, fpath;
	pair outerp[], innerp[];
	boolean ostart, ofinis;
	angstep = 4; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER!
	viewline = f-Tcenter;
	tmoment = ncrossprod( Tstart, Tfinis );
	vstart = ncrossprod( Tstart, tmoment );
	vfinis = ncrossprod( tmoment, Tfinis );
	finis = -vstart*conorm( Tstart );
	amax = getangle( Tstart, Tfinis );
	finc = N( Tfinis )*conorm( Tstart );
	if cdotprod( viewline, tmoment ) < 0:
	  tmoment := -tmoment;
	fi;
	sideaxe = ncrossprod( tmoment, viewline );
	for i=0 step angstep until amax:
	  circlecenter:= Tstart*cosd(i)+finis*sind(i);
	  circlemoment:= ccrossprod(circlecenter,tmoment);
	  cpath:=spatialhalfcircle(circlecenter+Tcenter,circlemoment,Sray,true);
	  if cdotprod( sideaxe, circlecenter ) < 0:
	    outerp[i/angstep]=point 0 of cpath;
	    innerp[i/angstep]=point (length cpath) of cpath;
	  else:
	    innerp[i/angstep]=point 0 of cpath;
	    outerp[i/angstep]=point (length cpath) of cpath;
	  fi;
	endfor;
	cstart = spatialhalfcircle(Tstart+Tcenter,vstart,Sray,true);
	if cdotprod( sideaxe, Tstart ) > 0:
	  cstart := reverse cstart;
	fi;
	cfinis = spatialhalfcircle(finc+Tcenter,vfinis,Sray,true);
	if cdotprod( sideaxe, finc ) < 0:
	  cfinis := reverse cfinis;
	fi;
	opath = outerp0 for i=angstep step angstep until amax:
	  ..outerp[i/angstep] endfor;
	ipath = innerp0 for i=angstep step angstep until amax:
	  ..innerp[i/angstep] endfor;
	fpath = ipath---cfinis---reverse opath---cstart---cycle;
	unfill fpath;
	draw fpath;
	ostart = cdotprod( viewline-Tstart, vstart ) > 0;
	if ostart:
	  cpath := rigorouscircle( Tcenter+Tstart, vstart, Sray );
	  unfill cpath;
	  draw cpath;
	fi;
	ofinis = cdotprod( viewline-finc, vfinis ) > 0;
	if ofinis:
	  cpath := rigorouscircle( Tcenter+finc, vfinis, Sray );
	  unfill cpath;
	  draw cpath;
	fi;
	%draw opath withcolor blue;
	%draw ipath withcolor red;
      endgroup
    enddef;

    def pointinsidetorus( expr Point, Tcenter, Tmoment, Bray, Sray ) =
      begingroup
	save outputboolean, height, dist, dray;
	boolean outputboolean;
	numeric height, dist, dray;
	height = cdotprod(N(Tmoment),Point-Tcenter);
	if abs(height)>=Sray:
	  outputboolean = false;
	else:
	  dist = conorm(Point-Tcenter-height*N(Tmoment));
	  dray = Sray +-+ abs(height);
	  if (dist<Bray+dray) and (dist>Bray-dray):
	    outputboolean = true;
	  else:
	    outputboolean = false;
	  fi;
	fi;
	( outputboolean )
      endgroup
    enddef;
	
    def pointrelativetotorus( expr Point, Tcenter, Tmoment, Bray, Sray ) =
      begingroup
	save height, dist;
	numeric height;
	color dist;
	height = cdotprod(N(Tmoment),Point-Tcenter);
	dist = N(Point-Tcenter-height*N(Tmoment));
	( conorm(Point-Bray*dist)-Sray )
      endgroup
    enddef;

    def intersectorus( expr Tcenter, Tmoment, Bray, Sray, LinePoi, LineDir ) =
      begingroup
	save trypoi, factry, linedi;
	color trypoi, linedi;
	numeric factry, j, auxd;
	trypoi = LinePoi;
        factry = 0.25;
        linedi = N(LineDir);
	for j= 1 upto 50:
	  auxd := pointrelativetotorus( trypoi, Tcenter, Tmoment, Bray, Sray );
	  trypoi := trypoi+factry*linedi*auxd;
	endfor;
	( trypoi )
      endgroup
    enddef;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Non-standard objects:
	    
    def positivecharge( expr InFactPositive, Center, BallRay ) =
      begingroup
	save auxc, axehorf, axeside, viewline, spath, pa, pb, pc, pd;
	color auxc, axehorf, axeside, viewline, pa, pb, pc, pd;
	path spath;
	viewline = f - Center;
	axehorf = N( ( X(viewline), Y(viewline), 0 ) );
	axeside = ccrossprod( axehorf, blue );
	if ShadowOn:
	  fill rigorousfearshadowpath( Center, BallRay );
	fi;
	spath = rigorousfearpath( Center, BallRay );
	unfill spath;
	draw spath;
	auxc = Center + sqrt(3)*axehorf;
	pa = auxc + axeside;
	pb = auxc - axeside;
	angline( pa, pb, Center, BallRay, "", top );
	if InFactPositive:
	  pc = auxc + blue;
	  pd = auxc - blue;
	  angline( pc, pd, Center, BallRay, "", top );
	fi
      endgroup
    enddef;

    def simplecar(expr RefP, AngCol, LenCol, FronWheelCol, RearWheelCol ) =
      begingroup
	save veca, vecb, vecc, anga, angb, angc, lena, lenb, lenc;
	save auxn, viewline, auxm, fl, fr, rl, rr, auxx, auxy;
	save fmar, fthi, fray, rmar, rthi, rray, inrefp;;
	color veca, auxx, auxy, vecb, vecc, viewline;
	color fl, fr, rl, rr, inrefp;
	numeric anga, angb, angc, lena, lenb, lenc, auxm, auxn;
	numeric fmar, fthi, fray, rmar, rthi, rray;
	anga = X( AngCol );
	angb = Y( AngCol );
	angc = Z( AngCol );
	lena = X( LenCol );
	lenb = Y( LenCol );
	lenc = Z( LenCol );
	fmar = X( FronWheelCol );
	fthi = Y( FronWheelCol );
	fray = Z( FronWheelCol );
	rmar = X( RearWheelCol );
	rthi = Y( RearWheelCol );
	rray = Z( RearWheelCol );
	veca = ( cosd(anga)*cosd(angb),
	  sind(anga)*cosd(angb),
	  sind(angb) );
	auxx = ( cosd(anga+90), sind(anga+90), 0 );
	auxy = ccrossprod( veca, auxx );
	vecb = cosd(angc)*auxx + sind(angc)*auxy;
	vecc = cosd(angc+90)*auxx + sind(angc+90)*auxy;
	viewline = f - RefP;
	auxm = cdotprod( viewline, veca );
	auxn = cdotprod( viewline, vecb );
	inrefp = RefP - 0.5*lenc*vecc;
	fl = inrefp + (0.5*lena-fmar-fray)*veca + 0.5*lenb*vecb;
	fr = inrefp + (0.5*lena-fmar-fray)*veca - 0.5*lenb*vecb;
	rl = inrefp - (0.5*lena-rmar-rray)*veca + 0.5*lenb*vecb;
	rr = inrefp - (0.5*lena-rmar-rray)*veca - 0.5*lenb*vecb;
	if auxn > 0.5*lenb:
	  if auxm > 0:
	    rigorousdisc( 0, true, rr, rray, -rthi*vecb );
	    rigorousdisc( 0, true, fr, fray, -fthi*vecb );
	    kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc);
	    rigorousdisc( 0, true, rl, rray, rthi*vecb );
	    rigorousdisc( 0, true, fl, fray, fthi*vecb );
	  else:
	    rigorousdisc( 0, true, fr, fray, -fthi*vecb );
	    rigorousdisc( 0, true, rr, rray, -rthi*vecb );
	    kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc);
	    rigorousdisc( 0, true, fl, fray, fthi*vecb );
	    rigorousdisc( 0, true, rl, rray, rthi*vecb );
	  fi;
	elseif auxn < -0.5*lenb:
	  if auxm > 0:
	    rigorousdisc( 0, true, rl, rray, rthi*vecb );
	    rigorousdisc( 0, true, fl, fray, fthi*vecb );
	    kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc);
	    rigorousdisc( 0, true, rr, rray, -rthi*vecb );
	    rigorousdisc( 0, true, fr, fray, -fthi*vecb );
	  else:
	    rigorousdisc( 0, true, fl, fray, fthi*vecb );
	    rigorousdisc( 0, true, rl, rray, rthi*vecb );
	    kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc);
	    rigorousdisc( 0, true, fr, fray, -fthi*vecb );
	    rigorousdisc( 0, true, rr, rray, -rthi*vecb );
	  fi;
	else:
	  if auxm > 0:
	    rigorousdisc( 0, true, rl, rray, rthi*vecb );
	    rigorousdisc( 0, true, fl, fray, fthi*vecb );
	    rigorousdisc( 0, true, rr, rray, -rthi*vecb );
	    rigorousdisc( 0, true, fr, fray, -fthi*vecb );
	    kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc);
	  else:
	    rigorousdisc( 0, true, fl, fray, fthi*vecb );
	    rigorousdisc( 0, true, rl, rray, rthi*vecb );
	    rigorousdisc( 0, true, fr, fray, -fthi*vecb );
	    rigorousdisc( 0, true, rr, rray, -rthi*vecb );
	    kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc);
	  fi;
	fi
      endgroup
    enddef;
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Differential Equations:
	    
% Oh! Well... I couldn't do without differential equations.
% The point is that I want to draw vectorial field lines in space.
% Keep it simple: second-order Runge-Kutta method.
% This is for solving first order differential equations

    def fieldlinestep( expr Spos, Step )( text VecFunc ) =
      begingroup
	save kone, ktwo;
	color kone, ktwo;
	kone = Step*VecFunc( Spos );
	ktwo = Step*VecFunc( Spos+0.5*kone );
	( Spos+ktwo )
      endgroup
    enddef;

    def fieldlinepath( expr Numb, Spos, Step )( text VecFunc ) =
      begingroup
	save ind, flpath, prevpos, thispos;
	numeric ind;
	color prevpos, thispos;
	path flpath;
	prevpos = Spos;
	flpath = rp( Spos )
	for ind=1 upto Numb:
	  hide( thispos := fieldlinestep( prevpos, Step, VecFunc ) )
	  ..rp( thispos )
	  hide( prevpos := thispos )
	endfor;
	( flpath )
      endgroup
    enddef;
    
% Another point is that I want to draw trajectories in space.
% This is for solving second order differential equations

    def trajectorypath( expr Numb, Spos, Svel, Step )( text VecFunc ) =
      begingroup
	save ind, flpath, prevpos, thispos, prevvel, thisvel;
	save rone, rtwo, vone, vtwo;
	numeric ind;
	color prevpos, thispos, prevvel, thisvel;
	color rone, rtwo, vone, vtwo;
	path flpath;
	prevpos = Spos;
	prevvel = Svel;
	flpath = rp( Spos )
	for ind=1 upto Numb:
	  hide(
	    vone := Step*VecFunc( prevpos );
	    rone := Step*prevvel;
	    vtwo := Step*VecFunc( prevpos+0.5*rone );
	    rtwo := Step*( prevvel+0.5*vone );
	    thisvel := prevvel+vtwo;
	    thispos := prevpos+rtwo
	    )
	  ..rp( thispos )
	  hide(
	    prevpos := thispos;
	    prevvel := thisvel
	    )
	endfor;
	( flpath )
      endgroup
    enddef;    

% Another point is that I want to draw trajectories in space and
% dependent on velocity: VecFunc( position, velocity ).
% This time is fourth-order Runge-Kutta.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHANGES PrintStep!!!!
    
    def dragtrajectorypath( expr Spos, Svel, Step )( text VecFunc ) =
      begingroup
	save ind, flpath, prevpos, thispos, prevvel, thisvel;
	save rone, rtwo, rthr, rfou, vone, vtwo, vthr, vfou;
	numeric ind;
	color prevpos, thispos, prevvel, thisvel;
	color rone, rtwo, rthr, rfou, vone, vtwo, vthr, vfou;
	path flpath;
	prevpos = Spos;
	prevvel = Svel;
	flpath = rp( Spos );
	ind = 1;
	forever:
	  vone := Step*VecFunc( prevpos , prevvel );
	  rone := Step*prevvel;
	  vtwo := Step*VecFunc( prevpos+0.5*rone, prevvel+0.5*vone );
	  rtwo := Step*( prevvel+0.5*vone );
	  vthr := Step*VecFunc( prevpos+0.5*rtwo, prevvel+0.5*vtwo );
	  rthr := Step*( prevvel+0.5*vtwo );
	  vfou := Step*VecFunc( prevpos+rthr, prevvel+vthr );
	  rfou := Step*( prevvel+vthr );	    
	  thisvel := prevvel+(vtwo+vthr)/3+(vone+vfou)/6;
	  thispos := prevpos+(rtwo+rthr)/3+(rone+rfou)/6;
	  exitif Z( thispos ) < -0.0001;                     %%%%%%%%%% EDIT!
	  prevpos := thispos;
	  prevvel := thisvel;
	  flpath := flpath--rp( thispos );
	endfor;
	PrintStep := Y(thispos);
	( flpath )
      endgroup
    enddef;    

% And now i stop.

    def magnetictrajectorypath( expr Numb, Spos, Svel, Step )( text VecFunc ) =
      begingroup
	save ind, flpath, prevpos, thispos, prevvel, thisvel;
	save rone, rtwo, rthr, rfou, vone, vtwo, vthr, vfou;
	numeric ind;
	color prevpos, thispos, prevvel, thisvel;
	color rone, rtwo, rthr, rfou, vone, vtwo, vthr, vfou;
	path flpath;
	prevpos = Spos;
	prevvel = Svel;
	flpath = rp( Spos )
	for ind=1 upto Numb:
	  hide(
	    vone := Step*ccrossprod( VecFunc( prevpos ), prevvel );
	    rone := Step*prevvel;
	    vtwo :=
	      Step*ccrossprod(VecFunc(prevpos+0.5*rone),prevvel+0.5*vone);
	    rtwo := Step*( prevvel+0.5*vone );
	    vthr :=
	      Step*ccrossprod(VecFunc(prevpos+0.5*rtwo),prevvel+0.5*vtwo);
	    rthr := Step*( prevvel+0.5*vtwo );
	    vfou :=
	      Step*ccrossprod( VecFunc( prevpos+rthr ), prevvel+vthr );
	    rfou := Step*( prevvel+vthr );	    
	    thisvel := prevvel+(vtwo+vthr)/3+(vone+vfou)/6;
	    thispos := prevpos+(rtwo+rthr)/3+(rone+rfou)/6
	    )
	  ..rp( thispos )
	  hide(
	    prevpos := thispos;
	    prevvel := thisvel
	    )
	endfor;
	( flpath )
      endgroup
    enddef;    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Part II:
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Advanced 3D-Object Definition Functions %%%%%
% Please check the examples in planpht.mp or the default object below %%%%

    vardef makeline@#( text vertices ) = 
      save counter;
      numeric counter;
      counter = 0;
      for ind=vertices:
	counter := incr( counter );
	L@#p[counter] := V[ind];
      endfor;
      npl@# := counter;
      NL := @#
    enddef;

    vardef makeface@#( text vertices ) = 
      save counter;
      numeric counter;
      counter = 0;
      for ind=vertices:
	counter := incr( counter );
	F@#p[counter] := V[ind];
      endfor;
      npf@# := counter;
      NF := @#;
      FCD[NF] := false
    enddef;
    
    vardef getready( expr commstr, refpoi ) =
      Nobjects := incr( Nobjects );
      ostr[Nobjects] := commstr;
      RefDist[Nobjects] := conorm( f - refpoi )
    enddef;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Definition of a 3D-Object 
% define vertices
    V1 := (1,0,0);
    V2 := (0,0,0);
    V3 := (0,1,0);
    V4 := (-0.3,0.2,1);
    V5 := (1,0,1);
    V6 := (0,1,1);
    V7 := (0,0,2);
    V8 := (-0.5,0.6,1.2);
    V9 := (0.6,-0.5,1.2);
    makeline1(8,9);
    makeface1(1,2,7);
    makeface2(2,3,7);
    makeface3(5,4,6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% or the old way below %%%%%%%%%
% define lines
%    NL := 1;     % number of lines
%    npl1 := 2;   % number of vertices of the first line
%    L1p1 := V8;
%    L1p2 := V9;
% define faces
%    NF := 3;     % number of faces
%    npf1 := 3;   % number of vertices of the first face
%    F1p1 := V1;
%    F1p2 := V2;
%    F1p3 := V7;
%    npf2 := 3;   % number of vertices of the second face
%    F2p1 := V2;
%    F2p2 := V3;
%    F2p3 := V7;
%    npf3 := 3;   % number of vertices of the third face
%    F3p1 := V5;
%    F3p2 := V4;
%    F3p3 := V6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
% Flip first argument accordingly to the second

    def flipvector(expr A, B) =
        begingroup
            save  nv;
            color nv;
            if cdotprod( A, B) < 0 :
                nv = -A;
            else:
                nv = A;
            fi;
            ( nv )
        endgroup
    enddef;

% Frontside of a face given by three of its vertices

    def facevector(expr A, B, C) =
        begingroup
            save  nv;
            color nv;
            nv = ncrossprod( A-B, B-C );
            ( flipvector( nv, f-B ) )
        endgroup
    enddef;

% Center or inside of a face

    def masscenter(expr Nsides)(suffix Coords) =
        begingroup
            save mc, counter;
            numeric  counter;
            color mc;
            mc = (0,0,0);
            for counter=1 upto Nsides:
                mc := mc + Coords[counter];
            endfor;
            ( mc / Nsides )
        endgroup
    enddef;

% Direction of coverability. The trick is here.
% The condition for visibility is 
% that the angle beetween the vector that goes from the side of a face to
% the mark position and the covervector must be greater than 90 degrees.

    def covervector(expr A, B, MassCenter) =
        begingroup
            save  nv;
            color nv;
            nv = ncrossprod( A-f, B-f );
            ( flipvector( nv, MassCenter-B) )
        endgroup
    enddef;

% O.K., the following macro tests the visibility of a point

    def themarkisinview(expr Mark, OwnFace) =
        begingroup
            save  c, faceVec, centerPoint, coverVec, inview, l, m; 
            color c, faceVec, centerPoint, coverVec;
            boolean inview;
            numeric l, m;
            l = 0;
            forever:          % cycle over faces until the mark is covered
                l := incr(l);
                if l = OwnFace:
                    l := incr(l);
                fi;
                exitif l > NF;
                faceVec := facevector(F[l]p1,F[l]p2,F[l]p3);
                inview := true;
                if cdotprod(Mark-F[l]p1, faceVec) < 0:     
                    centerPoint := masscenter(npf[l], F[l]p);
                    m := 0;
                    forever:              % cycle over segments of a face 
                        m := incr(m);
                        exitif m > npf[l];
                        if m < npf[l]:
                            c := F[l]p[m+1];
                        else:
                            c := F[l]p1;
                        fi;
                        coverVec := covervector(F[l]p[m], c, centerPoint);
                        inview := cdotprod(Mark-c,coverVec) <= 0;
                        exitif inview;
                    endfor;
                fi;
                exitif not inview;
            endfor;
            ( inview )
        endgroup
    enddef;

% Check for possible intersection or crossing.

    def maycrossviewplan(expr Ea, Eb, La, Lb) =
        begingroup
            ( abs( cdotprod( ccrossprod(Ea-f,Eb-f), La-Lb ) ) > 0.001 )
        endgroup
    enddef;

% Calculate the intersection of two sides. This is very nice.

    def crossingpoint(expr Ea, Eb, La, Lb) =
        begingroup
            save  thecrossing, perpend, exten, aux;
            color thecrossing, perpend;
            numeric exten, aux;
            if ( Ea = Lb ) or ( Ea = La ):
                thecrossing = Ea;
            elseif ( Eb = Lb ) or ( Eb = La ):
                thecrossing = Eb;
            else:
                perpend = ccrossprod( Ea-f, Eb-f );
                if conorm( perpend ) = 0:
                    thecrossing = Eb;
                else: 
                    aux = cdotprod( perpend, f );
                    cdotprod( perpend, thecrossing ) = aux;
                    ( La-Lb )*exten = La-thecrossing;
                fi;
            fi; 
            ( thecrossing )
        endgroup
    enddef;

% Calculate the intersection of an edge and a face.

    def crossingpointf(expr Ea, Eb, Fen) =
        begingroup
            save  thecrossing, perpend, exten;
            color thecrossing, perpend;
            numeric exten;
            perpend = ccrossprod( F[Fen]p1-F[Fen]p2, F[Fen]p3-F[Fen]p2 );
            cdotprod(perpend,thecrossing) = cdotprod( perpend, F[Fen]p2 );
            ( Ea-Eb )*exten = Ea-thecrossing; 
            ( thecrossing )
        endgroup
    enddef;

% Check for possible intersection of an edge and a face.

    def maycrossviewplanf(expr Ea, Eb, Fen) =
        begingroup
            save perpend;
            color perpend;
            perpend = ccrossprod( F[Fen]p1-F[Fen]p2, F[Fen]p3-F[Fen]p2 );
            ( abs( cdotprod( perpend, Ea-Eb ) ) > 0.001 )
        endgroup
    enddef;

% The intersection point must be within the extremes of the segment.

    def insidedge(expr Point, Ea, Eb) =
        begingroup
            save fract;
            numeric fract;
            fract := cdotprod( Point-Ea, Point-Eb );
            ( fract < 0 )
        endgroup
    enddef;

% Skip edges that are too far away

    def insideviewsphere(expr Ea, Eb, La, Lb) =
        begingroup
            save furthestofedge, nearestofline, flag, exten;
            color nearestofline, furthestofedge;
            boolean flag;
            numeric exten;
            nearestofline = La+exten*(Lb-La);
            cdotprod( nearestofline-f, Lb-La ) = 0;
            if conorm(Ea-f) < conorm(Eb-f):
                furthestofedge := Eb;
            else:
                furthestofedge := Ea;
            fi;
            if conorm(nearestofline-f) < conorm(furthestofedge-f):
                flag := true;
            else:
                flag := false;
            fi;
            ( flag )
        endgroup
    enddef;

% The intersection point must be within the triangle defined by 
% three points. Really smart.

    def insidethistriangle(expr Point, A, B, C ) =
        begingroup
            save arep, area, areb, aret, flag;
            color arep, area, areb, aret;
            boolean flag;
            aret = ccrossprod( A-C, B-C );
            arep = flipvector( ccrossprod( C-Point, A-Point ), aret );
            area = flipvector( ccrossprod( A-Point, B-Point ), aret );
            areb = flipvector( ccrossprod( B-Point, C-Point ), aret );
            flag = ( conorm( arep + area + areb ) <= 2*conorm( aret ) ); 
            ( flag )
        endgroup
    enddef;

% The intersection point must be within the triangle defined by the
% point of view f and the extremes of some edge. 

    def insideviewtriangle(expr Point, Ea, Eb) =
        insidethistriangle( Point, Ea, Eb, f )
    enddef;

% The intersection point must be within the face

    def insidethisface(expr Point, FaN) =
        begingroup
            save flag, m, central;
            boolean flag;
            numeric m;
            color central;
            m = npf[FaN];
            central = masscenter( m, F[FaN]p );
            flag = insidethistriangle( Point, 
                              central, F[FaN]p[m], F[FaN]p[1] );
            for m=2 upto npf[FaN]:
                exitif flag;
                flag := insidethistriangle( Point, 
                              central, F[FaN]p[m-1], F[FaN]p[m] );
            endfor;
            ( flag )
        endgroup
    enddef;

% Draw the visible parts of a straight line in beetween points A and B
% changing the thickness of the line accordingly to the distance from f

    def coarse_line(expr A, B, Facen, Press, Col) =
        begingroup
            save  k, mark, stepVec; 
            numeric k;
            color mark, stepVec;
            stepVec := resolvec(A,B); 
            k := 0;
            forever:                    % cycle along a whole segment 
                mark := A+(k*stepVec);
                exitif cdotprod(B-mark,stepVec) < 0;
                if themarkisinview(mark,Facen):
                    signalvertex(mark, Press, Col);
                fi;
                k := incr(k);
            endfor
        endgroup
    enddef;

% Get the 2D rigorous projection path of a face.
% Don't use SphericalDistortion here.

    def facepath(expr Facen) =
        begingroup
            save thispath, counter;
            path thispath;
            numeric counter;
            thispath = rp(F[Facen]p[1])--
            for counter=2 upto (npf[Facen]):
                rp(F[Facen]p[counter])--
            endfor
            cycle;
            ( thispath )
        endgroup
    enddef;

    def faceshadowpath(expr Facen) =
        begingroup
            save thispath, counter;
            path thispath;
            numeric counter;
            thispath = rp(cb(F[Facen]p[1]))--
            for counter=2 upto (npf[Facen]):
                rp(cb(F[Facen]p[counter]))--
            endfor
            cycle;
            ( thispath )
        endgroup
    enddef;

% FillDraw a face

    def face_invisible( expr Facen )( text LineAtribs ) =
        begingroup
            save ghost;
            path ghost;
            ghost = facepath( Facen );
            unfill ghost;
            draw ghost LineAtribs
        endgroup
    enddef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Different kinds of renderers:

% Draw only the faces, calculating edge-intersections. 
% Mind-blogging kind of thing.
% Only two constraints: i) faces must be planar and
%                      ii) faces must be convex.
% Pay attention: depends on PrintStep (resolvec).
    
    def sharpraytrace
%%                     ( expr LabelCrossPoints ) 
						  =
       begingroup
          save i, j, k, l, counter, a, b, c, d, currcross; 
          save flag, swapc, swapn, somepoint, frac, exten;
          save trythis, refpath, otherpath, intertimes;
          save counter, infolabel;
          numeric i, j, k, l, counter, swapn;
          color a, b, c, d, currcross, swapc;
          boolean flag, trythis;
          path refpath, otherpath;
          pair intertimes;
          string infolabel;
          color crosspoin[];
          numeric sortangle[];
          for i=1 upto NF:                   % scan all faces
             for j=1 upto npf[i]:            % scan all edges
                a := F[i]p[j];
                if j <> npf[i]:
                   b := F[i]p[j+1];
                else:
                   b := F[i]p1;
                fi;
                crosspoin[1] := a;
                counter := 2;
                refpath := rp(a)--rp(b); % The limits a and b of one
		                         % side of one face 
                for k=1 upto NF:
                   otherpath := facepath(k);
                   intertimes := refpath intersectiontimes otherpath;
                   trythis := xpart intertimes <> 0;
                   if trythis and (xpart intertimes <> 1) and (k <> i):
                      for l=1 upto npf[k]:    
                         c := F[k]p[l];
                         if l < npf[k]:
                            d := F[k]p[l+1];
                         else:
                            d := F[k]p1;
                         fi; 
                         if insideviewsphere( a, b, c, d ):
                           if maycrossviewplan( a, b, c, d ):
                             currcross := crossingpoint( a, b, c, d );
                             if insideviewtriangle( currcross, a, b ):
                               if insidedge( currcross, c, d ):
                                 swapc := ccrossprod( a-b, f-currcross);
                                 swapc := ccrossprod(swapc,f-currcross);
                                 color somepo;
                                 numeric fract;
                                 (b-a)*fract = somepo-a;
                                 cdotprod(swapc,somepo) =cdotprod(swapc,f);
                                 if (fract>0) and (fract<1):
                                   crosspoin[counter] := somepo;
                                   counter := incr(counter);
                                 fi;
                               fi;
                             fi;
                           fi;
                         fi;
                      endfor;
                      if maycrossviewplanf( a, b, k ):
                         currcross := crossingpointf( a, b, k );
                         if insidethisface( currcross, k ):
                            if insidedge( currcross, a, b ):
                               crosspoin[counter] := currcross;
                               counter := incr(counter);
                            fi;
                         fi;
                      fi;
                   fi;
                endfor;
                crosspoin[counter] := b;
                sortangle[1] := 0;
                for k=2 upto counter:
                   sortangle[k] := conorm(crosspoin[k]-a);
                endfor;
                forever:
                   flag := true;
                   for k=2 upto counter:
                      if sortangle[k] < sortangle[k-1]:
                         swapn := sortangle[k-1];
                         sortangle[k-1] := sortangle[k];
                         sortangle[k] := swapn;
                         swapc := crosspoin[k-1];
                         crosspoin[k-1] := crosspoin[k];
                         crosspoin[k] := swapc;
                         flag := false;
                      fi;
                   endfor;
                   exitif flag;
                endfor;
                for k=2 upto counter:
                   swapc := resolvec(crosspoin[k-1],crosspoin[k]);
                   flag := themarkisinview( crosspoin[k-1]+swapc, i );
                   if flag and themarkisinview( crosspoin[k]-swapc, i ):
                      draw rp(crosspoin[k-1])--rp(crosspoin[k]);
                   fi;
                endfor;
%                if LabelCrossPoints:
%                  for k=1 upto counter:
%                    infolabel:=decimal(i)&","&decimal(j)&","&decimal(k); 
%                    infolabel := "0";
%                    dotlabelrand(infolabel,rp(crosspoin[k]));
%                  endfor;
%                fi;
             endfor;
          endfor
       endgroup
    enddef;

% Draw three-dimensional lines checking visibility. 

    def lineraytrace(expr Press, Col) =
        begingroup
            save  i, j, a, b; 
            numeric i, j;
            color a, b;
            for i=1 upto NL:                    % scan all lines
                for j=1 upto npl[i]-1:
                    a := L[i]p[j];
                    b := L[i]p[j+1]; 
                    coarse_line( a, b, 0, Press, Col);
                endfor;    
            endfor
        endgroup
    enddef;

% Draw only the faces, rigorously projecting the edges. 

    def faceraytrace(expr Press, Col) =
        begingroup
            save  i, j, a, b; 
            numeric i, j;
            color a, b;
            for i=1 upto NF:                    % scan all faces
                for j=1 upto npf[i]:
                    a := F[i]p[j];
                    if j <> npf[i]:
                        b := F[i]p[j+1];
                    else:
                        b := F[i]p1;
                    fi; 
                    coarse_line( a, b, i, Press, Col);
                endfor;    
            endfor
        endgroup
    enddef;

% Fast test for your three-dimensional object

    def draw_all_test( expr AlsoDrawLines ) =
      begingroup
	save  i, j, a, b; 
	numeric i, j;
	color a, b;
	if ShadowOn:
	  for i=1 upto NF:                    
	    fill faceshadowpath( i );
	  endfor;
	  if AlsoDrawLines:
	    for i=1 upto NL:                    % scan all lines
	      for j=1 upto npl[i]-1:
		a := L[i]p[j];
		b := L[i]p[j+1]; 
		drawsegment( cb(a), cb(b) );
	      endfor;    
	    endfor;
	  fi;
	fi;
	for i=1 upto NF:                    % scan all faces
	  for j=1 upto npf[i]:
	    a := F[i]p[j];
	    if j <> npf[i]:
	      b := F[i]p[j+1];
	    else:
	      b := F[i]p1;
	    fi; 
	    drawsegment( a, b );
	  endfor;    
	endfor;
	if AlsoDrawLines:
	  for i=1 upto NL:                    % scan all lines
	    for j=1 upto npl[i]-1:
	      a := L[i]p[j];
	      b := L[i]p[j+1]; 
	      drawsegment( a, b );
	    endfor;    
	  endfor;
	fi
      endgroup
    enddef;

% Don't use SphericalDistortion here.

    def fill_faces( text LineAtribs ) =
        begingroup
            save  i; 
            numeric i;
	    if ShadowOn:
	      for i=1 upto NF:                    
                fill faceshadowpath( i );
	      endfor;
	    fi;
            for i=1 upto NF:                    
                face_invisible( i, LineAtribs );
            endfor
        endgroup
    enddef;

    def doitnow =
      begingroup
	save i, j, farone;
	numeric i, j, farone[];
	for i=1 upto Nobjects:
	  farone[i] := i;
	endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shell's Method of sorting %%%%%%%%
	save inc, v, vv;
	numeric inc, v, vv;
	inc = 1;
	forever:
	  inc := 3*inc+1;
	  exitif inc > Nobjects;
	endfor;
	forever:
	  inc := round(inc/3);
	  for i=inc+1 upto Nobjects:
	    v := RefDist[i];
	    vv:= farone[i];
	    j := i;
	    forever:
	      exitunless RefDist[j-inc] > v;
	      RefDist[j] := RefDist[j-inc];
	      farone[j] := farone[j-inc];
	      j := j-inc;
	      exitif j <= inc;
	    endfor;
	    RefDist[j] := v;
	    farone[j] := vv;
	  endfor;
	  exitunless inc > 1;
	endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	for i=Nobjects downto 1:
	  j := farone[i];
	  scantokens ostr[j];
	endfor
      endgroup
    enddef;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Nematic Liquid Crystal wise:

    def generatedirline(expr Lin, Phi, Theta, Long, Currpos ) =
        begingroup
            save longvec;
            color longvec;
            npl[Lin] := 2;
            longvec := Long*( cosd(Phi)*cosd(Theta),
                              sind(Phi)*cosd(Theta),
                              sind(Theta) );
            L[Lin]p1 := Currpos-0.5*longvec;
            L[Lin]p2 := Currpos+0.5*longvec
        endgroup
    enddef;

    def generatedirface(expr Fen, Phi, Theta, Long, Base, Currpos ) =
        begingroup
            save basevec, longvec;
            color basevec, longvec;
            npf[Fen] := 3;
            longvec := Long*( cosd(Phi)*cosd(Theta),
                              sind(Phi)*cosd(Theta),
                              sind(Theta) );
            basevec := Base*ncrossprod( Currpos-f, longvec );
            F[Fen]p1 := Currpos-0.5*(longvec+basevec);
            F[Fen]p2 := Currpos+0.5*longvec;
            F[Fen]p3 := Currpos-0.5*(longvec-basevec)
        endgroup
    enddef;

    def generateonebiax(expr Lin, Phi, Theta, Long, SndDirAngl, Base, Currpos ) =
        begingroup
            save basevec, longvec, u, v;
            color basevec, longvec, u, v;
            npl[Lin] := 4;
            longvec := Long*( cosd(Phi)*cosd(Theta),
                              sind(Phi)*cosd(Theta),
                              sind(Theta) );
            v = (-sind(Phi), cosd(Phi), 0);
            u = ( cosd(Phi)*cosd(Theta+90),
                  sind(Phi)*cosd(Theta+90),
                  sind(Theta+90) );
            basevec := Base*( v*cosd(SndDirAngl)+u*sind(SndDirAngl) );
            L[Lin]p1 := Currpos-0.5*longvec;
            L[Lin]p2 := Currpos+0.5*basevec;
            L[Lin]p3 := Currpos+0.5*longvec;
            L[Lin]p4 := Currpos-0.5*basevec
        endgroup
    enddef;

    def director_invisible( expr SortEmAll, ThickenFactor, CyclicLines ) =
      begingroup
	save i, j, k, farone, thisfar;
	save outerr, innerr, direc, ounum; 
	numeric i, j, k, farone[], dist[], thisfar, ounum;
	pen actualpen, outerr, innerr;
	path direc;
	actualpen = currentpen;
	if SortEmAll:
	  for i=1 upto NL:                    % scan all lines
	    dist[i] := conorm( masscenter( npl[i], L[i]p ) - f );
	    farone[i] := i;
	  endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shell's Method of sorting %%%%%%%%
	  save inc, v, vv;
	  numeric inc, v, vv;
	  inc = 1;
	  forever:
	    inc := 3*inc+1;
	    exitif inc > NL;
	  endfor;
	  forever:
	    inc := round(inc/3);
	    for i=inc+1 upto NL:
	      v := dist[i];
	      vv:= farone[i];
	      j := i;
	      forever:
		exitunless dist[j-inc] > v;
		dist[j] := dist[j-inc];
		farone[j] := farone[j-inc];
		j := j-inc;
		exitif j <= inc;
	      endfor;
	      dist[j] := v;
	      farone[j] := vv;
	    endfor;
	    exitunless inc > 1;
	  endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	else:
	  for i=1 upto NL:
	    farone[i] := i;
	  endfor;
	fi;
	for i=NL downto 1:                    % draw all pathes
	  j := farone[i];
	  direc := rp( L[j]p1 )
	  for k=2 upto npl[j]:
	    --rp( L[j]p[k] )
	  endfor;
	  if CyclicLines:
	    direc := direc--cycle;
	  fi;
	  ounum := Spread*ps( masscenter(npl[j],L[j]p), ThickenFactor );
	  outerr := pencircle scaled ounum;
	  innerr := pencircle scaled (0.8*ounum); %% DANGER %%
	  pickup outerr;
	  draw direc withcolor black;
	  pickup innerr;
	  draw direc withcolor background;
	endfor;
	pickup actualpen
      endgroup
    enddef;

% Now two routines to draw "bananas", well, sort of...
% Initially coded by Pedro J. Sebastião

    def circularsheet( expr CenterP, Rad, VecX, VecY, StartA, FinisA, Width ) =
      begingroup
	save vecz, ind;
	color vecz;
	numeric ind;
	vecz = ncrossprod( VecX, VecY );
	( rp( CenterP + Width*vecz + Rad*planarrotation( VecX, VecY, StartA ) )
	  for ind=StartA+1 upto FinisA:
	    --rp( CenterP + Width*vecz + Rad*planarrotation( VecX, VecY, ind ))
	  endfor
	  for ind=FinisA downto StartA:
	    --rp( CenterP + Rad*planarrotation( VecX, VecY, ind ) )
	  endfor
	  --cycle )
      endgroup
    enddef;
  
    def banana( expr CenterPos,  Radius, AngleColor, Wid, Amp  ) =
      begingroup
	save ind, sinbeta, cosbeta, aux, delta, angfvbx, angpos, angneg, au;
	save bx, by, bz, fv, beta, gamma, alfa;
	save outpath, outpathb, outpathc, visneg, vispos;
	numeric ind, sinbeta, cosbeta, aux, delta, angfvbx, angpos, angneg, au;
	numeric beta, gamma, alfa;
	color bx, by, bz, fv;
	path outpath, outpathb, outpathc;
	boolean visneg, vispos;
	
	alfa = X(AngleColor);
	beta = Y(AngleColor);
	gamma= Z(AngleColor);
	bx = eulerrotation( alfa, beta, gamma, red );
	by = eulerrotation( alfa, beta, gamma, green );
	bz = eulerrotation( alfa, beta, gamma, blue );
	
	au = cdotprod( f-CenterPos, by );
	if 0 > au:
	  by := -by;
	fi;
	fv = cdotprod( f, bx )*bx + cdotprod( f, by )*by;
	aux = conorm( fv-CenterPos );
	if aux > Radius:
	  cosbeta = Radius/aux;
	  sinbeta = ( aux +-+ Radius )/aux;
	  delta = angle( cosbeta, sinbeta );
	  angfvbx = getangle( fv - CenterPos, bx );
	  if 0 > cdotprod( fv - CenterPos, by ):
	    angfvbx := -angfvbx;
	  fi;
	  angpos = delta + angfvbx;
	  angneg = angfvbx - delta;
	  if ( angneg > -Amp ) and ( Amp > angneg ):
	    visneg = true;
	  else:
	    visneg = false;
	  fi;
	  if ( angpos > -Amp ) and ( Amp > angpos ):
	    vispos = true;
	  else:
	    vispos = false;
	  fi;
	  if visneg and not vispos:
	    outpath = circularsheet(CenterPos,Radius,bx,by,-Amp,angneg,Wid);
	    unfill outpath;
	    draw outpath;
	    outpathb = circularsheet(CenterPos,Radius,bx,by,angneg,Amp,Wid);
	    unfill outpathb;
	    draw outpathb
	  fi;
	  if vispos and not visneg:
	    outpath = circularsheet(CenterPos,Radius,bx,by,-Amp,angpos,Wid);
	    unfill outpath;
	    draw outpath;
	    outpathb = circularsheet(CenterPos,Radius,bx,by,angpos,Amp,Wid);
	    unfill outpathb;
	    draw outpathb
	  fi;
	  if (not vispos) and (not visneg):
	    outpath = circularsheet(CenterPos,Radius,bx,by,-Amp,Amp,Wid);
	    unfill outpath;
	    draw outpath;
	  fi;
	  if vispos and visneg:
	    if 0 > cdotprod( f-CenterPos, bx ):
	      outpath=circularsheet(CenterPos,Radius,bx,by,angneg,angpos,Wid);
	      unfill outpath;
	      draw outpath;
	      outpathb = circularsheet(CenterPos,Radius,bx,by,-Amp,angneg,Wid);
	      unfill outpathb;
	      draw outpathb;
	      outpathc = circularsheet(CenterPos,Radius,bx,by,angpos,Amp,Wid);
	      unfill outpathc;
	      draw outpathc;
	    else:
	      outpathb = circularsheet(CenterPos,Radius,bx,by,-Amp,angneg,Wid);
	      unfill outpathb;
	      draw outpathb;
	      outpathc = circularsheet(CenterPos,Radius,bx,by,angpos,Amp,Wid);
	      unfill outpathc;
	      draw outpathc;
	      outpath=circularsheet(CenterPos,Radius,bx,by,angneg,angpos,Wid);
	      unfill outpath;
	      draw outpath;	  
	    fi;
	  fi;
	else:
	  outpath = circularsheet( CenterPos, Radius, bx, by, -Amp, Amp, Wid );
	  unfill outpath;
	  draw outpath;
	fi;
	draw rp(CenterPos+Radius*bx)--rp(CenterPos+Radius*bx+Wid*bz);
      endgroup
    enddef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Plotting:

% Define and draw surfaces with a triangular mesh.
% On a hexagonal or triangular area. Without sorting (no need).
    
    def hexagonaltrimesh( expr BeHexa,theN,SideSize )( text SurFunc ) = 
        begingroup
	    save i, j, posx, posy, posz, higx, higy,
	    counter, stepx, stepy, poi, lowx, lowy,
	    newn, bola, bolb, bolc;
	    numeric i, j, posx, posy, posz, higx, higy,
                    counter, stepx, stepy, lowx, lowy, newn;
	    color poi[][];
	    boolean bola, bolb, bolc;
	    npf0 := 3;
	    FCD0 := true; % this is used in the calls to fillfacewithlight
	    ActuC := incr( ActuC );
	    if ActuC > TableColors:
	      ActuC := 1;
	    fi;
	    FC0  := ActuC;    %%
	    counter = 0;
            stepy = SideSize/theN;
            stepx = 0.5*stepy*sqrt(3);
	    lowy = -0.5*SideSize;
	    lowx = -sqrt(3)*SideSize/6;
	    higy = -lowy;
	    higx = sqrt(3)*SideSize/3;
	    for i=0 upto theN:
	        for j=0 upto theN-i:
                    posx := lowx + i*stepx;
                    posy := lowy + i*stepx/sqrt(3) + j*stepy;
	    	    posz := SurFunc( posx, posy );
	    	    poi[i][j] := ( posx, posy, posz );
	        endfor;
	    endfor;
	    if BeHexa:
	        newn = round((theN+1)/3)+1;
	    else:
	        newn = 1;
	    fi;
	    for j=newn upto theN-newn+1:
	        F0p1 := poi[0][j-1];
	        F0p2 := poi[0][j];
	        F0p3 := poi[1][j-1];
	        fillfacewithlight( 0 );   % see below
	    endfor;
	    for i=1 upto theN-1:
	        for j=1 upto theN-i:
		    bola := ( i < newn ) and ( j < newn-i );
		    bolb := ( i < newn ) and ( j > theN-newn+1 );
		    bolc := ( i > theN-newn );
		    if not ( bola or bolb or bolc ):
	                F0p1 := poi[i-1][j];
			F0p2 := poi[i][j-1];
			F0p3 := poi[i][j];
			fillfacewithlight( 0 );
			F0p1 := poi[i+1][j-1];
			fillfacewithlight( 0 );
		    fi;
	        endfor;
	    endfor;
	    i := theN-newn+1;
	    for j=1 upto newn-1:
	        F0p1 := poi[i-1][j];
		F0p2 := poi[i][j-1];
		F0p3 := poi[i][j];
		fillfacewithlight( 0 );
	    endfor;
        endgroup 
    enddef;
    
    def fillfacewithlight( expr FaceN ) =
      begingroup
	save perpvec, reflectio, viewvec, inciden, refpos, projincid;
	save shiftv, fcol, lcol, theangle, ghost, pa, pb, pc, j, lowcolor;
	color perpvec, reflectio, viewvec, inciden, refpos, projincid;
	color shiftv, fcol, lcol, pa, pb, pc, lowcolor;
	numeric theangle, j;
	path ghost;
	ghost := rp( F[FaceN]p1 )
	  for j=2 upto npf[FaceN]:
	    --rp( F[FaceN]p[j] )
	  endfor
	  --cycle;
	if OverRidePolyhedricColor:
	  unfill ghost;
	else:
	  refpos = masscenter( npf[FaceN], F[FaceN]p );
	  pa = F[FaceN]p1;
	  pb = F[FaceN]p2;
	  pc = F[FaceN]p3;
	  inciden = LightSource - refpos;
	  viewvec = f - refpos;
	  perpvec = ncrossprod( pa-pb, pb-pc );
	  if cdotprod( perpvec, blue ) < 0:
	    perpvec := -perpvec;
	  fi;
	  projincid = perpvec*cdotprod( perpvec, inciden );
	  shiftv = inciden - projincid;
	  reflectio = projincid - shiftv;
	  theangle = getangle( reflectio, viewvec );
	  if FCD[FaceN]:
	    lowcolor = TableC[FC[FaceN]];
	  else:
	    lowcolor = TableC0;
	  fi;
	  lcol = (cosd( theangle ))[lowcolor,HigColor];
	  if cdotprod( viewvec, perpvec ) < 0:
	    fcol = lcol - SubColor;
	  else:
	    fcol = lcol;
	  fi;
	  fill ghost withcolor fcol;
	fi;
	draw ghost;
      endgroup
    enddef;
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Part III (parametric plots and another renderer):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%                                         %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%   Kindly contributed by Jens Schwaiger  %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%                                         %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dmin_ is the minimal distance of all faces to f;
% dmax_ is the maximal distance of all faces to f;
% (both values are determined in "draw_invisible")
% Facen is the number of the face to be filled;
% ColAtrib is the color used for filling;
% ColAtribone the color used for drawing;
% Colordensity depends on distance of the face from f

    def face_drawfill( expr Facen, dmin_, dmax_ ,ColAtrib, ColAtribone ) =
      begingroup
	save j, ptmp, colfac_, coltmp_;
	path ghost;
	numeric j, colfac_;
	color ptmp, coltmp_;
	ghost := rp( F[Facen]p1 );
	for j=2 upto npf[Facen]:
	  ghost := ghost--rp( F[Facen]p[j] );
	endfor;
	ghost := ghost--cycle;
	ptmp:= masscenter( npf[Facen], F[Facen]p ) - f;
        % 0<=colfac_<=1
	colfac_ := ( conorm(ptmp)-dmin_ )/( dmax_ - dmin_ );
        % color should be brighter, if distance to f is smaller
	colfac_ := 1 - colfac_;
        % color should be not to dark, i.e., >=0.1 (and <=1)
	colfac_ := 0.9colfac_ + 0.1; 
        % now filling and drawing with appropriate color;
	fill ghost withcolor colfac_*ColAtrib;
	draw ghost withcolor colfac_*ColAtribone;
      endgroup
    enddef;

% Now a much faster faces-only-ray-tracer based upon the unfill
% command and the constraint of non-intersecting faces of similar
% sizes. Faces are sorted by distance from nearest vertex or
% masscenter to the point of view. This routine may be used to 
% draw graphs of 3D surfaces (use masscenters in this case).
%    
% Option=true: test all vertices
% Option=false: test masscenters of faces
% DoJS=true: use face_drawfill by J. Schwaiger
% DoJS=false: use fillfacewithlight by L. Nobre G.
% ColAtrib=color for filling faces
% ColAtribone=color for drawing edges

    def draw_invisible( expr Option, DoJS, ColAtrib, ColAtribone ) =
      begingroup
	save  i, j, a, b, thisfar, ptmp, farone;
	numeric i, j, farone[], dist[], thisfar, distmin_, distmax_;
	color a, b, ptmp;
	for i=1 upto NF:                         % scan all faces
	  if Option:                             % for distances of
	    dist[i] = conorm( F[i]p1 - f );      % nearest vertices
	    if i=1: 
	      distmin_ := dist1;                 % initialisation of
	      distmax_ := dist1;                 % dmin_ and dmax_
	    fi;
	    distmin_ := min( distmin_, dist[i] );
	    distmax_ := max( distmax_, dist[i] );
	    for j=2 upto npf[i]:
	      thisfar := conorm( F[i]p[j] - f );
	      distmin_ := min( distmin_, thisfar );
	      distmax_ := max( distmax_, thisfar );
	      if thisfar < dist[i]:
		dist[i] := thisfar;
	      fi;
	    endfor;
	  else:                        % for distances of centers of mass
	    dist[i] := conorm( masscenter( npf[i], F[i]p ) - f );
	    if i=1:
	      distmin_ := dist1;       % initialisation of
	      distmax_ := dist1;       % dmin_ and dmax_ in this case 
	    fi;
	    distmin_ := min( distmin_, dist[i] );
	    distmax_ := max( distmax_, dist[i] );
	  fi;
	  farone[i] = i;
	endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shell's Method of sorting %%%%%%%%
	save inc, v, vv;
	numeric inc, v, vv;
	inc = 1;
	forever:
	  inc := 3*inc+1;
	  exitif inc > NF;
	endfor;
	forever:
	  inc := round(inc/3);
	  for i=inc+1 upto NF:
	    v := dist[i];
	    vv:= farone[i];
	    j := i;
	    forever:
	      exitunless dist[j-inc] > v;
	      dist[j] := dist[j-inc];
	      farone[j] := farone[j-inc];
	      j := j-inc;
	      exitif j <= inc;
	    endfor;
	    dist[j] := v;
	    farone[j] := vv;
	  endfor;
	  exitunless inc > 1;
	endfor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	for i=NF downto 1:                     % draw and fill all pathes
	  j := farone[i];
	  if DoJS:
	    face_drawfill( j, distmin_, distmax_, ColAtrib, ColAtribone );
	  else:
	    fillfacewithlight( j );
	  fi;
	endfor;
      endgroup
    enddef;

% Move to the good range (-1,1).

    def bracket( expr low, poi, hig ) = 
        begingroup
            save zout;
            numeric zout;
            zout = (2*poi-hig-low)/(hig-low);
            if zout > 1:
                zout := 1;
            fi;
            if zout < -1:
                zout := -1;
            fi;
            ( zout )
        endgroup 
    enddef;

% Define parametric surfaces with a triangular mesh... unless a
% quadrangular mesh can do a fine, rigorous job just as well.

    def partrimesh( expr nt,ns,lowt,higt,lows,higs,lowx,higx,lowy,higy,lowz,higz,facz)( text parSurFunc ) =
      begingroup
	save i, j, k, posx, posy, posz;
	save counter, stept, steps, poss, post, tmpaux;
	save veca, vecb, vecc, vecd;
	numeric i, j, k, posx, posy, posz, counter, stept, steps;
	color poi[][], tmpaux, veca, vecb, vecc, vecd;
	counter := NF;                     % <-- NF must be initialized!
	ActuC := incr( ActuC );
	if ActuC > TableColors:
	  ActuC := 1;
	fi;
	steps = ( higs - lows )/ns;
	stept = ( higt - lowt )/nt;
	for i=0 upto ns: 
	  for j=0 upto nt:
	    poss := lows + i*steps;
	    post := lowt + j*stept;
	    tmpaux := parSurFunc( poss, post );
	    posz := Z(tmpaux); 
	    posx := X(tmpaux);
	    posy := Y(tmpaux);
	    posx := bracket(lowx,posx,higx);      
	    posy := bracket(lowy,posy,higy);      
	    posz := bracket(lowz,posz,higz)/facz; 
	    poi[i][j] := ( posx, posy, posz );
	  endfor;
	endfor;
	for i=1 upto ns:
	  for j=1 step 1 until nt:
	    veca := poi[i][j]-poi[i-1][j];
	    vecb := poi[i][j]-poi[i-1][j-1];
	    vecc := poi[i][j]-poi[i][j-1];
	    if abs(cdotprod(ccrossprod(veca,vecb),vecc))<0.00005: %DANGER!
	      counter := incr(counter);
	      npf[counter] := 4;
	      F[counter]p1 := poi[i-1][j-1];
	      F[counter]p2 := poi[i-1][j];
	      F[counter]p3 := poi[i][j];
	      F[counter]p4 := poi[i][j-1];
	      FC[counter] := ActuC;
	      FCD[counter] := true;
	    else:
	      tmpaux:=
	        0.25*(poi[i-1][j-1]+poi[i-1][j]+poi[i][j]+poi[i][j-1]);
	      veca := poi[i-1][j-1]-tmpaux;
	      vecb := poi[i-1][j]-tmpaux;
	      vecc := poi[i][j]-tmpaux;
	      vecd := poi[i][j-1]-tmpaux;
	      if getangle(vecb,vecd)>getangle(veca,vecc):
		for k=-1 upto 0:
		  counter := incr(counter);
		  npf[counter] := 3;
		  F[counter]p1 := poi[i-1][j-1];
		  F[counter]p2 := poi[i+k][j-1-k];
		  F[counter]p3 := poi[i][j];
		  FC[counter] := ActuC;
		  FCD[counter] := true;
		endfor;
	      else:
		for k=-1 upto 0:
		  counter := incr(counter);
		  npf[counter] := 3;
		  F[counter]p1 := poi[i+k][j+k];
		  F[counter]p2 := poi[i][j-1];
		  F[counter]p3 := poi[i-1][j];
		  FC[counter] := ActuC;
		  FCD[counter] := true;
		endfor;
	      fi;
	    fi;
	  endfor;
	endfor;
	NF := counter;
      endgroup
    enddef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Part IV (automatic perspective tuning and minimization):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    def randomfear =
      begingroup
	save a, b, c, i, h;
	numeric a, b, c, i, h;
	i = uniformdeviate( 360 );
	h = uniformdeviate( 180 );
	a = cosd( h )*cosd( i );
	b = sind( h )*cosd( i );
	c = sind( i );
	( (a,b,c) )
      endgroup
    enddef;

    def renormalizevc( expr inF, inVC ) =
      begingroup
	save a;
	color a;
	cdotprod( inF, a ) = 0;
	inF - a = whatever*( inF - inVC );
	( a )
      endgroup
    enddef;

    def calculatecost( expr TryF, TryVc, TrySp, TrySh ) =
      begingroup
	save sumsquares, i, difpair, xx, yy;
	numeric sumsquares, i, xx, yy;
	pair difpair;
	sumsquares = 0;
	f := TryF;
	viewcentr := TryVc;
	Spread := TrySp;
	ShiftV := TrySh;
	for i=1 upto PhotoMarks:
	  difpair := 0.05*(rp(PhotoPoint[i]) - PhotoPair[i]); % DANGER
	  xx := ((xpart difpair)**2)/PhotoMarks;
	  yy := ((ypart difpair)**2)/PhotoMarks;
	  sumsquares := sumsquares + xx + yy;
	endfor;
	( sumsquares )
      endgroup
    enddef;
	
    def forcepointinsidefear( text A ) =
      begingroup
	if abs( X(A) ) > 0.5*MaxFearLimit :
	  A := 0.5*A*MaxFearLimit/abs( X(A) );
	fi;
	if abs( Y(A) ) > 0.5*MaxFearLimit :
	  A := 0.5*A*MaxFearLimit/abs( Y(A) );
	fi;
	if abs( Z(A) ) > 0.5*MaxFearLimit :
	  A := 0.5*A*MaxFearLimit/abs( Z(A) );
	fi;
	( A )
      endgroup
    enddef;

    def forcepairinsidepage( text A ) =
      begingroup
	if xpart A > 2*(xpart OriginProjPagePos):
	  A := ( 2*(xpart OriginProjPagePos), ypart A );
	fi;
	if ypart A > 2*(ypart OriginProjPagePos):
	  A := ( xpart A, 2*(ypart OriginProjPagePos) );
	fi;
	if xpart A < 0:
	  A := ( 0, ypart A );
	fi;
	if ypart A < 0:
	  A := ( xpart A, 0 );
	fi;
	( A )
      endgroup
    enddef;

    def calculatejump( expr AverCost, PrevCost, RandCost, JumpLimit ) =
      begingroup
	save funfact, numer, denom;
	numeric funfact, numer, denom;
	if RandCost+PrevCost > 2*AverCost:
	  numer := 3*PrevCost - 4*AverCost + RandCost;
	  denom := 4*( RandCost + PrevCost - 2*AverCost );
	  if abs( denom )*JumpLimit < abs( numer ):
	    funfact := 0;
	  else:
	    funfact := numer/denom;
	  fi;
	else:
	  funfact := 0;
	fi;
	( funfact )
      endgroup
    enddef;
    	
    def photoreverse( expr IterNum, ExpTao, JumpFact ) =
      begingroup
	save j, auxvc, actfact, auxfac;
	save prevf, randf, averf, prevvc, randvc, avervc;
	save prevsp, randsp, aversp, prevsh, randsh, aversh;
	save prevcost, randcost, avercost, spreadlimit, expfact;
	numeric j, prevsp, randsp, aversp, actfact;
	numeric prevcost, randcost, avercost;
	numeric spreadlimit, expfact;
	color prevf, randf, averf, prevvc, randvc, avervc, auxvc;
	pair prevsh, randsh, aversh;
	spreadlimit = 20;                                         % DANGER
	prevf = f;
	prevvc= viewcentr;
	prevsp= Spread;
	prevsh= ShiftV;
	prevcost = calculatecost( prevf, prevvc, prevsp, prevsh );
	show prevcost;
	auxfac = -1280.0/ExpTao/IterNum;
	for j=0 upto IterNum:
	  expfact := mexp(auxfac*j);
	  randf := prevf + expfact*randomfear;
	  randcost := calculatecost( randf, prevvc, prevsp, prevsh );
	  averf := 0.5[prevf,randf];
	  avercost := calculatecost( averf, prevvc, prevsp, prevsh );
	  actfact := calculatejump(avercost,prevcost,randcost,JumpFact);
	  prevf := actfact[prevf,randf];	  
	  auxvc := prevvc + expfact*randomfear;
	  randvc:= renormalizevc( randf, auxvc );
	  randcost := calculatecost( prevf, randvc, prevsp, prevsh );
	  auxvc := 0.5[prevvc,randvc];
	  avervc:= renormalizevc( averf, auxvc );
	  avercost := calculatecost( prevf, avervc, prevsp, prevsh );
	  actfact := calculatejump(avercost,prevcost,randcost,JumpFact);
	  auxvc := actfact[prevvc,randvc];
	  prevvc:= renormalizevc( prevf, auxvc );
	  randsp:= prevsp+expfact*spreadlimit*(uniformdeviate( 1 ) - 0.5);
	  randcost := calculatecost( prevf, prevvc, randsp, prevsh );
	  aversp:= 0.5[prevsp,randsp];
	  avercost := calculatecost( prevf, prevvc, aversp, prevsh );
	  actfact := calculatejump(avercost,prevcost,randcost,JumpFact);
	  prevsp:= actfact[prevsp,randsp];
	  randsh:= prevsh + expfact*dir( uniformdeviate( 360 ) ); % DANGER
	  randcost := calculatecost( prevf, prevvc, prevsp, randsh );
	  aversh:= 0.5[prevsh,randsh];
	  avercost := calculatecost( prevf, prevvc, prevsp, aversh );
	  actfact := calculatejump(avercost,prevcost,randcost,JumpFact);
	  prevsh:= actfact[prevsh,randsh];
	  prevcost := calculatecost( prevf, prevvc, prevsp, prevsh );
	  %show (prevcost,expfact);
	endfor;
	show prevcost;
      endgroup
    enddef;

% Minimization routine for scalar functions like y=f(x) where an initial
% triplet (x1,x2,x3) with x1<x2<x3 is given as a parabolic squeleton that
% provides a way to search for the smallest value of y (if iterated)

    def minimizestep( expr Abcisscolor )( text PlainFunc ) =
      begingroup
	save xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
	save colout;
	numeric xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
	color colout;
	xa = X( Abcisscolor );
	xb = Y( Abcisscolor );
	xc = Z( Abcisscolor );
	ya = PlainFunc(xa); 
	yb = PlainFunc(xb);
	yc = PlainFunc(xc);
	if ya = yb:
	  colout = (-0.125[xa,xb],xb,xc);
	elseif yb = yc:
	  colout = (xa,xb,1.125[xb,xc]);
	else:
	  if (yb>ya) or (yb>yc):
	    show  Abcisscolor;
	    message "                             Unable to minimizestep!";
	  fi;
	  den = (xb-xc)*((xa**2)-(xb**2))-(xa-xb)*((xb**2)-(xc**2));
	  if abs(den) < 0.0005:
	    show den;
	    message "                             Unable to minimizestep!";
	  fi;
	  coeb = ((yb-yc)*((xa**2)-(xb**2))-(ya-yb)*((xb**2)-(xc**2)))/den;
	  coec = ((xb-xc)*(ya-yb)-(xa-xb)*(yb-yc))/den;
	  xd = -0.5*coeb/coec;
	  yd = PlainFunc( xd );
	  if ((xa<xd) and (xd<xb)):
	    if (yd<yb):
	      colout = (xa,xd,xb);
	    else:
	      colout = (xd,xb,xc);
	    fi;
	  elseif ((xb<xd) and (xd<xc)):
	    if (yd<yb):
	      colout = (xb,xd,xc);
	    else:
	      colout = (xa,xb,xd);
	    fi;
	  else:
	    aux := 0.125[xb,xc]-0.125[xb,xa];
	    colout = (xa,0.125[xb,xa]+uniformdeviate(aux),xc);
	  fi;
	fi;
	( colout )
      endgroup
    enddef;
	  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Part V:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Intersections:
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Check lineintersectplan above %%%%%%%%%%%%%%%
    
    def calculatecostver(expr VerA,DisA,VerB,DisB,VerC,DisC,TryV) =
      begingroup
	save a, b, c;
	numeric a, b, c;
	a = conorm( TryV - VerA ) - DisA;
	b = conorm( TryV - VerB ) - DisB;
	c = conorm( TryV - VerC ) - DisC;
	( ( a ++ b ++ c )**2 )
      endgroup
    enddef;

% Approximation of the intersection of three spheres.    
% Be aware of the next three danger parameters.
    
    def improvertex( expr VerA, DisA, VerB, DisB, VerC, DisC, IniV ) =
      begingroup
	save j, iternum, prevcost, factry, actpoi, actfact, auxa, auxb;
	save jumplim, trypoi, halfpo, expfact, randcost, avercost, auxc;
	numeric j, iternum, prevcost, factry, randcost, avercost, actfact;
	numeric jumplim, expfact;
	color actpoi, trypoi, halfpo, auxa, auxb, auxc;
	trypoi = IniV;
	iternum = 24;                                        % DANGER!
	factry = 0.25;                                       % DANGER!
	jumplim = 2.5;                                       % DANGER!
	prevcost = calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
	message "improvertex initial cost: " & decimal(prevcost)
	                                     & " at point " & cstr(trypoi);
	for j=1 upto 5:
	  auxa := VerA-trypoi;
	  auxa := N(auxa)*(conorm(auxa)-DisA);
	  auxb := VerB-trypoi;
	  auxb := N(auxb)*(conorm(auxb)-DisB);
	  auxc := VerC-trypoi;
	  auxc := N(auxc)*(conorm(auxc)-DisC);
	  trypoi := trypoi+factry*(auxa+auxb+auxc);
	endfor;
	prevcost := calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
	message "phase one cost: " & decimal(prevcost)
	                           & " at point " & cstr(trypoi);
 	for j=0 upto iternum:
 	  expfact := mexp(-j*200/iternum);                    % DANGER!
 	  actpoi := trypoi+factry*expfact*randomfear;
 	  randcost :=
 	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,actpoi);
 	  halfpo := 0.5[trypoi,actpoi];
 	  avercost :=
 	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,halfpo);
 	  actfact := calculatejump(avercost,prevcost,randcost,jumplim);
 	  trypoi := actfact[trypoi,actpoi];
 	  prevcost :=
 	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
 %	  show prevcost;
 	endfor;
 	message "final cost: " & decimal(prevcost)
 	                       & " at point " & cstr(trypoi);
	( trypoi )
      endgroup;
    enddef;

% The approximation of the intersection of a plane, an infinite cylinder (tube)
% and a spheroid (an oblate spheroid, a prolate spheroid or a perfect sphere)  
    
    def ultraimprovertex( expr PlanPoi, PlanDir, BaseCenter, Radius, LenVec,
	CentrPoi, NorthPoleVec, Ray, IniV ) =
	begingroup
	save trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, factry, focn, focs;
	save norpoldi, major, majortmp, nordi, tmpdi;
	color trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, focn, focs;
	color norpoldi, nordi, tmpdi;
	numeric factry, major, majortmp;
	trypoi = IniV;
	factry = 0.25;
	plandi = N(PlanDir);
	cyldi = N(LenVec);
	norpoldi = N(NorthPoleVec);
	major = conorm(NorthPoleVec);
	if major>Ray:
	  focn := CentrPoi+(major+-+Ray)*norpoldi;
	  focs := CentrPoi-(major+-+Ray)*norpoldi;
	fi;
	for j=1 upto 50:
	  if major<Ray:
	    tmpdi := trypoi-CentrPoi;
	    nordi := N(tmpdi-cdotprod(tmpdi,norpoldi)*norpoldi);
	    focn := CentrPoi+(Ray+-+major)*nordi;
	    focs := CentrPoi-(Ray+-+major)*nordi;
	  elseif major=Ray:
	    focn := CentrPoi;
	    focs := CentrPoi;
	  fi;  
	  auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
	  auxb := cyldi*cdotprod(trypoi-BaseCenter,cyldi)-trypoi+BaseCenter;
	  auxb := N(auxb)*(conorm(auxb)-Radius);
	  auxc := focn-trypoi;
	  auxd := focs-trypoi;
	  if major>Ray:
	    auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*major);
	  else:
	    auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*Ray);
	  fi;  
	  trypoi := trypoi+factry*(auxa+auxb+auxc);
	endfor;
	( trypoi )
      endgroup
    enddef;

% The approximation of the intersection of a plane and two prolate spheroids

    def necplusimprovertex( expr PlanPoi, PlanDir,
	CentrPoiA, NorthPoleVecA, RayA,
        CentrPoiB, NorthPoleVecB, RayB, IniV ) =
        begingroup
        save trypoi, auxa, auxb, auxc, auxd, plandi, factry, focni, focsi;
        save norpoldi, norpoldj, maior, major, focnj, focsj, auxe;
        color trypoi, auxa, auxb, auxc, auxd, plandi, focni, focsi, focnj, focsj;
        color norpoldi, norpoldj, auxe;
        numeric factry, maior, major;
        trypoi = IniV;
        factry = 0.25;
        plandi = N(PlanDir);
        norpoldi = N(NorthPoleVecA);
        maior = conorm(NorthPoleVecA);        
        norpoldj = N(NorthPoleVecB);
        major = conorm(NorthPoleVecB);
        focni = CentrPoiA+(maior+-+RayA)*norpoldi;
        focsi = CentrPoiA-(maior+-+RayA)*norpoldi;
        focnj = CentrPoiB+(major+-+RayB)*norpoldj;
        focsj = CentrPoiB-(major+-+RayB)*norpoldj;
        for j=1 upto 50:
          auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
          auxb := focni-trypoi;
          auxe := focsi-trypoi;
          auxc := focnj-trypoi;
          auxd := focsj-trypoi;
          auxb := (N(auxb)+N(auxe))*(conorm(auxb)+conorm(auxe)-2*maior);
          auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*major);
          trypoi := trypoi+factry*(auxa+auxb+auxc);
        endfor;
        ( trypoi )
      endgroup
    enddef;

% The approximation of the intersection of a straight line and
% one prolate spheroid

    def intersectprolatespheroid( expr CentrPoi, NorthPoleVec, Ray,
        LinePoi, LineDir ) =
        begingroup
	save trypoi, factry, linedi, norpol, focn, focs, maior, j;
	save auxa, auxb, auxc, auxd;
	color trypoi, linedi, norpol, focn, focs;
	color auxa, auxb, auxc, auxd;
	numeric factry, maior, j;
        trypoi = LinePoi;
        factry = 0.25;
        linedi = N(LineDir);
        norpol = N(NorthPoleVec);
        maior = conorm(NorthPoleVec);
        focn = CentrPoi+(maior+-+Ray)*norpol;
        focs = CentrPoi-(maior+-+Ray)*norpol;
        for j=1 upto 50:
          auxb := focn-trypoi;
          auxd := focs-trypoi;
          auxa := (N(auxb)+N(auxd))*(conorm(auxb)+conorm(auxd)-2*maior);
          auxc := linedi*cdotprod(auxa,linedi);
          trypoi := trypoi+factry*auxc;
        endfor;
        ( trypoi )
      endgroup
    enddef;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Part VI (strictly two-dimensional):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Verify if a path is cyclic (written by Scott Pakin)
    
    def is_cyclic expr cpath =
      (point 0 of cpath = point (length cpath) of cpath)
    enddef;

% Produce the schematics of a spring.

    def springpath( expr begp, endp, piturnum, piturnproj, spgfrac ) =  
	begingroup
	    boolean leftside;
	    numeric counter, springwidth;
	    pair stepdir, leftdir, rightdir, auxil, lastp, endvec;
	    path thespring;
	    leftside = true;
	    stepdir = spgfrac*(endp-begp)/piturnum;
	    endvec = (1-spgfrac)*(endp-begp)/2;
	    springwidth = piturnproj +-+ abs( stepdir );
	    auxil = ( -ypart stepdir, xpart stepdir );
	    leftdir = springwidth*unitvector( auxil );
	    auxil := ( ypart stepdir, -xpart stepdir );
	    rightdir = springwidth*unitvector( auxil );
	    leftdir := (stepdir + leftdir)/2;
	    rightdir := (stepdir + rightdir)/2;
	    lastp = begp+endvec;
	    thespring := begp--lastp;
	    for counter=1 upto piturnum:
		if leftside:
		    auxil := lastp + leftdir;
		else:
		    auxil := lastp + rightdir;
		fi;
		lastp := lastp + stepdir;
		thespring := thespring--auxil--lastp;
		leftside := not leftside;
	    endfor;
	    thespring := thespring--endp;
	    ( thespring )
	endgroup
    enddef;

% Summarize a great length in a zig-zag frontier line

    def zigzagfrontier( expr begp, endp, nzigs, dev, zthick, tthick, fthick, excol, incol ) =  
	begingroup
	    interim linecap := squared;
	    interim linejoin := mitered;
	    boolean leftside;
	    numeric counter, tmpval;
	    pair stepdir, leftdir, rightdir, auxil, lastp;
	    path thefrontier;
	    leftside = true;
	    stepdir = (endp-begp)/nzigs;
	    leftdir = unitvector( ( -ypart stepdir, xpart stepdir ) );
	    rightdir = unitvector( ( ypart stepdir, -xpart stepdir ) );
	    thefrontier = begp;
	    lastp = begp;
	    for counter=1 upto nzigs-1:
		lastp := lastp + stepdir;
		tmpval := zthick+normaldeviate*dev;
		if leftside:
		    auxil := lastp + leftdir*tmpval;
		else:
		    auxil := lastp + rightdir*tmpval;
		fi;
		thefrontier := thefrontier--auxil;
		leftside := not leftside;
	    endfor;
	    thefrontier := thefrontier--endp;
	    draw
	      thefrontier withcolor excol withpen pencircle scaled tthick;
	    draw
	      thefrontier withcolor incol withpen pencircle scaled fthick;
	    ( thefrontier )
	endgroup
    enddef;

% The name says it all.

    def randomcirc( expr radi, stddev, numpois ) =  
	begingroup
	    numeric i, astep;
	    path ranc;
	    astep = 360/numpois;
	    ranc = (radi+stddev*normaldeviate)*dir(-180);
	    for i= -180+astep step astep until 180:
		ranc := ranc--((radi+stddev*normaldeviate)*dir(i));
	    endfor;
	    ranc := ranc--cycle;
	    ( ranc )
	endgroup
    enddef;

% Just label some 2D-place in a way similar to anglinen.

    def labeln(expr S, Pos, RelPos) =
        begingroup
            if RelPos = 0:
                label.rt( S, Pos );
            elseif RelPos =1:
                label.urt( S, Pos );
            elseif RelPos =2:
                label.top( S, Pos );
            elseif RelPos =3:
                label.ulft( S, Pos );
            elseif RelPos =4:
                label.lft( S, Pos );
            elseif RelPos =5:
                label.llft( S, Pos );
            elseif RelPos =6:
                label.bot( S, Pos );
            elseif RelPos =7:
                label.lrt( S, Pos );
            else:
                label( S, Pos );
            fi
        endgroup
    enddef;

% There must be a sort of planar cross-product. The z coordinate. 

    def paircrossprod(expr A, B) = 
        ( (xpart A)*(ypart B) - (xpart B)*(ypart A) )
    enddef;

% Just dot-label some 2D-place in a random way.

    def dotlabelrand(expr S, Pos ) =
        begingroup
            save RelPos;
            numeric RelPos;
            RelPos = floor( uniformdeviate(9) );
            if RelPos = 0:
                dotlabel.rt( S, Pos );
            elseif RelPos =1:
                dotlabel.urt( S, Pos );
            elseif RelPos =2:
                dotlabel.top( S, Pos );
            elseif RelPos =3:
                dotlabel.ulft( S, Pos );
            elseif RelPos =4:
                dotlabel.lft( S, Pos );
            elseif RelPos =5:
                dotlabel.llft( S, Pos );
            elseif RelPos =6:
                dotlabel.bot( S, Pos );
            elseif RelPos =7:
                dotlabel.lrt( S, Pos );
            else:
                dotlabel( S, Pos );
            fi
        endgroup
    enddef;

    def radialcross( expr A, la, B, lb, GoUp) =
      begingroup
	numeric x, y, xa, xb, ya, yb, YM, YA, La, Lb;
	numeric AA, BB, CC, auxil, na, nb, norm; 
	pair As, Bs, selectedpoint;
	na = abs(A);
	nb = abs(B);
	norm := 0;
	for t = na, nb, la, lb:
	  if norm < t:
	    norm := t;
	  fi;
	endfor;
	xa = xpart A/norm;
	xb = xpart B/norm;
	ya = ypart A/norm;
	yb = ypart B/norm;
	La = la/norm;
	Lb = lb/norm;
	if abs( ya - yb ) < 0.005 :
	  x := La**2 - Lb**2 + xb**2 - xa**2;
	  x := 0.5*x/( xb - xa );
	  auxil := La +-+ (xa-x);
	  As  = ( x, ya + auxil );
	  Bs = ( x, ya - auxil );
	else:
	  YM := (xb-xa)/(ya-yb);
	  YA := Lb**2 - La**2 + xa**2 - xb**2;
	  YA := 0.5*( YA - (ya-yb)**2 )/(ya-yb);
	  AA := 1 + YM**2;
	  BB := 2*( YM*YA - xa );
	  CC := xa**2 - La**2 + YA**2;
	  CC := sqrt( BB**2 - 4*AA*CC );
	  x  := -0.5*( BB + CC )/AA;
	  y  := YA + ya + YM*x;
	  Bs = ( x, y );
	  x  := -0.5*( BB - CC )/AA;
	  y  := YA + ya + YM*x;
	  As = ( x, y );
	fi;
	if ypart As > ypart Bs:
	  if GoUp:
	    selectedpoint = As;
	  else:
	    selectedpoint = Bs;
	  fi;
	elseif ypart As = ypart Bs:
	  if xpart As > xpart Bs:
	    if GoUp:
	      selectedpoint = As;
	    else:
	      selectedpoint = Bs;
	    fi;
	  else:
	    if GoUp:
	      selectedpoint = Bs;
	    else:
	      selectedpoint = As;
	    fi;
	  fi;
	else:
	  if GoUp:
	    selectedpoint = Bs;
	  else:
	    selectedpoint = As;
	  fi;
	fi;                        
	( norm*selectedpoint )
      endgroup
    enddef;                        

    def ropethread( expr Index ) =
      begingroup
	save aux;
	numeric aux;
	if Index > RopeColors:
	  aux = 0;
	else:
	  aux = Index;
	fi;
	( aux )
      endgroup
    enddef;
    
    def ropepattern( expr BasePath, RopeWidth, Nturns ) =
      begingroup
	save indturns, nmoves, indthread, movelen, turnlen, totlen;
	numeric indturns, nmoves, indthread, movelen, turnlen, totlen;
	save lenpos, timar, steplen, indstep, startdownc, startupcol;
	numeric lenpos, timar, steplen, startdownc, indstep;
	save actuc, actdc, stepwidth;
	numeric actuc, actdc, stepwidth, startupcol;
	save p;
	pair p[];
	save actcolor;
	color actcolor;
	nmoves = 2*(RopeColors+1);
	totlen = arclength BasePath;
	turnlen = totlen/Nturns;
	movelen = turnlen/nmoves;
	steplen = movelen/2;
	startdownc = 0;
	startupcol = RopeColors;
	stepwidth = RopeWidth/RopeColors;
	for indturns=0 upto Nturns-1:
	  for indmove=0 upto nmoves-1:
	    for indstep=0 upto 3:
	      lenpos :=
		indturns*turnlen+indmove*movelen+indstep*steplen;
	      timar := arctime lenpos of BasePath;
	      p[indstep] := direction timar of BasePath rotated 90;
	      p[indstep] := unitvector( p[indstep] );
	      p[indstep+4] := point timar of BasePath; 
	    endfor;
	    actdc := startdownc;
	    for indthread=0 upto RopeColors:
	      p8 := p5-p1*(0.5*RopeWidth-(indthread-0.5)*stepwidth);
	      p9 := p4-p0*(0.5*RopeWidth-indthread*stepwidth);
	      p10:= p5-p1*(0.5*RopeWidth-(indthread+0.5)*stepwidth);
	      p11:= p6-p2*(0.5*RopeWidth-indthread*stepwidth);
	      actcolor := TableC[RopeColorSeq[actdc]];
	      fill p8--p9--p10--p11--cycle withcolor actcolor;
	      actdc := ropethread( incr( actdc ) );
	    endfor;
	    startdownc := ropethread( incr( startdownc ) );
	    actuc := startupcol;
	    p9 := p5+p1*0.5*(RopeWidth+stepwidth);
	    p10:= p6+p2*0.5*RopeWidth;
	    p11:= p7+p3*0.5*(RopeWidth+stepwidth);
	    actcolor := TableC[RopeColorSeq[actuc]];
	    fill p9--p10--p11--cycle withcolor actcolor;
	    actuc := ropethread( incr( actuc ) );
	    for indthread=0 upto RopeColors-1:
	      p8 := p6+p2*(0.5*RopeWidth-indthread*stepwidth);
	      p9 := p5+p1*(0.5*RopeWidth-(indthread+0.5)*stepwidth);
	      p10:= p6+p2*(0.5*RopeWidth-(indthread+1)*stepwidth);
	      p11:= p7+p3*(0.5*RopeWidth-(indthread+0.5)*stepwidth);
	      actcolor := TableC[RopeColorSeq[actuc]];
	      fill p8--p9--p10--p11--cycle withcolor actcolor;
	      actuc := ropethread( incr( actuc ) );
	    endfor;
	    p8 := p6-p2*0.5*RopeWidth;
	    p9 := p5-0.5*p1*(RopeWidth+stepwidth);
	    p11:= p7-0.5*p3*(RopeWidth+stepwidth);
	    actcolor := TableC[RopeColorSeq[actuc]];
	    fill p8--p9--p11--cycle withcolor actcolor;
	    startupcol := ropethread( incr( startupcol ) );
	  endfor;
	endfor
      endgroup
    enddef;
	
    def firsttangencypoint( expr Path, Point, ResolvN ) =
      begingroup
	save auxp, i, cutp, va, vb;
	path auxp;
	numeric i;
	pair cutp, va, vb;
	auxp =
	  hide( va := unitvector( point 0 of Path - Point );
	  vb := unitvector( direction 0 of Path ); )
	  ( paircrossprod( va, vb ), 0 )
	  for i=1/ResolvN step 1/ResolvN until length Path:
	  hide( va := unitvector( point i of Path - Point );
	    vb := unitvector( direction i of Path ); )
	    ...( paircrossprod( va, vb ), i )
	endfor;
	cutp = auxp intersectionpoint ( origin--( 0, length Path ) );
	( point ( ypart cutp ) of Path )
      endgroup
    enddef;

% Shrink or swell a cyclic path without cusp points and without
% coinciding pre and post control points. 

    def lasermachine( expr DefinedPath, Beam, CosLimit ) =
      begingroup
	save patlen, j;
	save apoi, bpoi, cpoi, dpoi, epoi;
	save apat, bpat, cpat, dpat, a, b, c;
	save anew, bnew, cnew;
	save aang, bang, cang, dang;
	save newp, pairvector, cou, val;
	numeric patlen, j;
	pair apoi, bpoi, cpoi, dpoi, epoi;
	pair apat, bpat, cpat, dpat, a, b, c;
	pair anew, bnew, cnew, pairvector[];
	numeric aang, bang, cang, dang, cou, val;
	path newp;
	patlen = length DefinedPath;
	cou = 0;	
	for j=0 upto patlen-1:
	  apoi := precontrol j of DefinedPath;
	  bpoi := point j of DefinedPath;
	  cpoi := postcontrol j of DefinedPath;
	  dpoi := precontrol j+1 of DefinedPath;
	  epoi := point j+1 of DefinedPath;
	  apat := apoi-bpoi;
	  bpat := bpoi-cpoi;
	  cpat := cpoi-dpoi;
	  dpat := dpoi-epoi;
	  aang := angle( apat );
	  bang := angle( bpat );
	  cang := angle( cpat );
	  dang := angle( dpat );
	  val := cosd( 0.5*(aang-bang) );
	  bnew := cpoi+Beam*dir( 90+0.5*(bang+cang) )/cosd( 0.5*(bang-cang) );
	  cnew := dpoi+Beam*dir( 90+0.5*(cang+dang) )/cosd( 0.5*(cang-dang) );
	  if ( val > 0 ) and
	    ( val < CosLimit ) and
	      ( Beam*sind( 0.5*(aang-bang) ) < 0 ):
	      
	      a := bpoi+Beam*dir(90+aang);
	    b := a+0.5*Beam*dir(aang); %%% DANGER %%%%%0.5%%%%%%%
	    %draw b withpen pencircle scaled 3pt withcolor green;
	    anew := bpoi+Beam*dir(90+bang);
	    c := anew-0.5*Beam*dir(bang); %%% DANGER %%0.5%%%%%%%
	    %draw c withpen pencircle scaled 3pt withcolor blue;
	    pairvector[3*cou] = a;
	    pairvector[3*cou+1] = b;
	    pairvector[3*cou+2] = c;
	    cou := incr(cou);
	    pairvector[3*cou] = anew;
	    pairvector[3*cou+1] = bnew;
	    pairvector[3*cou+2] = cnew;
	    cou := incr(cou);
	  else:
	    anew := bpoi+Beam*dir( 90+0.5*(aang+bang) )/val;
	    pairvector[3*cou] = anew;
	    pairvector[3*cou+1] = bnew;
	    pairvector[3*cou+2] = cnew;
	    cou := incr(cou);
	  fi;	
	endfor;
	newp = for j=0 upto cou-1:
	    pairvector[3*j]..controls pairvector[3*j+1] and pairvector[3*j+2]..
	endfor cycle;
	( newp )
      endgroup
    enddef;  

% Move the starting point of a cyclic path along that path

    def startahead( expr DefinedPath, JumpTime ) =
      begingroup
	save patlen, j;
	save apoi, bpoi, cpoi;
	save newp;
	numeric patlen, j;
	pair apoi, bpoi, cpoi;
	path newp;
	patlen = length DefinedPath;
	newp = for j=0 upto patlen-1:
	    hide( 
	      apoi := point JumpTime+j of DefinedPath;
	    bpoi := postcontrol JumpTime+j of DefinedPath;
	    cpoi := precontrol JumpTime+j+1 of DefinedPath;
	    )
	    apoi..controls bpoi and cpoi..
        endfor
	cycle;
	( newp )
      endgroup
    enddef;

% In order to use a "lasermachine" one needs a single full outline.
% One may have two somewhat concentric cyclic paths intersecting in
% several points.
% The next routine may help but first the paths must be adapted with
% "startahead" and/or "reverse" so that they both rotate in the same
% direction and they start on consecutive "lobes" (hard to explain).
% Now pay attention: given the direction of rotation (clockwise or
% counter-clockwise) the SecondPath must start BEFORE the FirstPath.
% And another problem: there must be at least four intersection points.
% Very nasty routine. All because of finispoi...
  
    def crossingline( expr FirstPath, SecondPath, TimeTolerance ) =
      begingroup
	save m;
	save its, finispoi;
	save increm, fo, ma, tmpp, mastarter;
	numeric m;
	pair its, finispoi;
	path increm, fo, ma, tmpp, mastarter;
	m = TimeTolerance;
	fo = FirstPath;
	ma = SecondPath;
	its := fo intersectiontimes ma;
	increm := subpath (m, (xpart its) - m ) of fo;
	mastarter = subpath (m, (ypart its) - m ) of ma;
	finispoi = reverse ma intersectionpoint reverse fo;
	
	forever:
	  
	  fo := subpath ( (xpart its)+m, length fo ) of fo;
	  ma := subpath ( (ypart its)+m, length ma ) of ma;
	  its := ma intersectiontimes fo;
	  tmpp := subpath ( m, (xpart its)-m ) of ma;
	  increm := increm...tmpp;
	  
	  fo := subpath ( (ypart its)+m, length fo ) of fo;
	  ma := subpath ( (xpart its)+m, length ma ) of ma;
	  its := fo intersectiontimes ma;
	  tmpp := subpath ( m, (xpart its)-m ) of fo;
	  increm := increm...tmpp;
	  
	  exitif abs( point (xpart its) of fo - finispoi ) < m;
	endfor;
	
	fo := subpath ( (xpart its)+m, length fo ) of fo;
	ma := (subpath ( (ypart its)+m, (length ma)-m ) of ma)...mastarter;
	its := ma intersectiontimes fo;
	tmpp := subpath ( m, (xpart its)-m ) of ma;
	increm := increm...tmpp;
	
	tmpp := subpath ( (ypart its)+m, (length fo)-m ) of fo;
	( increm...tmpp...cycle )
      endgroup
    enddef;
    
% Calculate path areas (contributed by Boguslaw Jackowski
% to the metapost mailing list)

 vardef segmentarea( expr Ps ) = 
   save xa, xb, xc, xd, ya, yb, yc, yd;
   ( xa, 20ya ) = point 0 of Ps;
   ( xb, 20yb ) = postcontrol 0 of Ps;
   ( xc, 20yc ) = precontrol 1 of Ps;
   ( xd, 20yd ) = point 1 of Ps;
   ( xb - xa )*( 10ya + 6yb + 3yc +   yd )
     + ( xc - xb )*( 4ya + 6yb + 6yc +  4yd )
     + ( xd - xc )*( ya + 3yb + 6yc + 10yd )
 enddef;
 
 vardef cyclicpatharea( expr P ) = % result = area of the interior
   segmentarea(subpath (0,1) of P)
     for t=1 upto length(P)-1: + segmentarea(subpath (t,t+1) of P) endfor
 enddef;
 
% Mark bidimensional angles (contributed by Palle Jørgensen 
% to the metapost mailing list)

 vardef archangle@#( expr _p, _q, _s, archwidth ) text _t =
   begingroup;
     save _a, _b, _w, _arch, _halfangle, _label_origin;
     ( _a, _b ) = _p intersectiontimes _q;
     pair _w;
     _w = whatever[
       point _a of _p +
       archwidth * unitvector direction _a of _p,
       point _a of _p +
       archwidth * unitvector direction _a of _p +
       (ypart.direction _a of _p, -xpart.direction _a of _p)
       ];
     _w  = whatever[point _b of _q, point _b of _q + direction _b of _q];
     path _arch;
     _arch = point _a of _p +
       archwidth * unitvector direction _a of _p{
       (if direction _a of _p dotprod unitvector direction _b of _q > 0:
	 1
     else:
       -1
     fi) *
   ( _w - (point _a of _p + archwidth * unitvector direction _a of _p) )
     }..point _b of _q +
     archwidth * unitvector direction _b of _q;
   draw _arch _t;
   path _halfangle;
   _halfangle = point _a of _p - 2*archwidth*
     unitvector( direction _a of _p + direction _b of _q )--point _a of _p +
     2*archwidth*unitvector( direction _a of _p + direction _b of _q );
   pair _label_origin;
   _label_origin = _halfangle intersectionpoint _arch;
   label@#( _s, _label_origin ) _t;
 endgroup; 
enddef;

% EOF