trunk/Master/texmf-dist/metapost

git-svn-id: svn://tug.org/texlive/trunk@104 c570f23f-e606-0410-a88d-b1316a301751

2 files changed, 3607 insertions, 0 deletions

diff --git a/Master/texmf-dist/metapost/featpost/featpost.mp b/Master/texmf-dist/metapost/featpost/featpost.mp
new file mode 100644
index 00000000000..9c30b765fed
--- /dev/null
+++ b/Master/texmf-dist/metapost/featpost/featpost.mp
@@ -0,0 +1,36 @@
+% featpost.mp
+% L. Nobre G.
+% 2005
+%
+% Use "export TEX=latex" in your .bashrc.
+% Begin by pre-compiling this set of macros with "inimpost featpost.mp".
+% Produce your PostScript (PS) figures with "mpost -mem featpost file".
+% The output of this command (one or several files named "file.N")
+% must be in the current directory to use the bashscripts laproof,
+% lbproof and lcproof. Produce your EncapsulatedPostScript (EPS) figures 
+% with "bashscript/lXproof file N". laproof tranforms PS into EPS if
+% the figure is smaller than an A4 page. lbproof produces JPEG and EPS
+% that fits the width of an A4 portrait page and lcproof produces EPS 
+% that fits the height of an A4 portrait page independently of its 
+% original size.
+% N is the number of the figure (in file) that you want to encapsulate.
+% It is not necessary to use "export TEX=latex" nor "lXproof file N"
+% if the figures have no text. 
+% The resulting EPS figures are not insertable in LaTex documents.
+% The originals file.N are.
+
+input plain;
+input featpost3Dplus2D;
+
+dump;
+
+% It is possible to interactively experiment the effect
+% of figure parameters. This requires:
+% 1) MetaPost code written as "anglinerigorouscircle.mp".
+% 2) xcmd package (kindly provided by Pedro Sebastião);
+%	  it may be downloaded from http://lince.cii.fc.ul.pt/
+%	  in Debian, Red Hat or tar.gz package formats.
+% 3) Perl.
+% 4) gv.
+% 5) run command "xcmd/xmpost anglinerigorouscircle" or 
+%		 "xcmd/xmpost anymetapostwoextension".
diff --git a/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
new file mode 100644
index 00000000000..b6055afae58
--- /dev/null
+++ b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
@@ -0,0 +1,3571 @@
+% featpost3Dplus2D.mp
+% L. Nobre G., C. Barbarosie, J. Schwaiger and B. Jackowski
+% nobre@lince.cii.fc.ul.pt
+% http://matagalatlante.org
+% Copyright (C) 2005
+% see also featpost.mp
+
+% This set of macros extends the MetaPost language
+% to three dimensions and eases the production of 
+% physics diagrams.
+
+% 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 2
+% 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 "Preloading FeatPost macros, version 0.6.7";
+
+    warningcheck := 0;
+
+    background := 0.987white;
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global Variables %%%%%%%%%%%
+
+    boolean	 ParallelProj, SphericalDistortion, FCD[], ShadowOn;
+    boolean      OverRidePolyhedricColor;
+    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;
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Kept for backward compatibility
+
+            Shifts := 105.00mm;
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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
+
+            ShiftV := 105.00mm*(1,1); % Central coordinates on paper
+
+	    OriginProjPagePos := (105.00mm,148.45mm); % This should be the
+	                                              % page center. 
+
+	    ParallelProj := false; % Kind of perspective 
+				       % Can't have both true
+	    SphericalDistortion := false; % Kind of lens
+
+	    ShadowOn := false; % Some objects may block the light and
+
+	    HoriZon := 0; % cast a shadow on a horizontal plane at this Z
+	    
+	    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
+
+	    defaultscale := 0.75;
+	    defaultfont := "cmss17";   % This is used by cartaxes
+
+	    PageHeight := 9in;
+	    PageWidth := 6in;   % And this is used by produce_auto_scale
+
+            MaxFearLimit := 15;       % 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
+	    	    
+	    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
+	    TDAhalftipbase := 0.02;         % Three-Dimensional
+	    TDAhalfthick := 0.01;           % Arrow
+
+	    NCL := 0;                       % closedline
+	    ForePen := pencircle scaled 15pt;
+	    BackPen := pencircle scaled 9pt;
+	    
+            %%% The variables PhotoMarks, PhotoPair[], PhotoPoint[]
+	    %%% and CLPath[] have NO default values.
+	    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Part I:
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Very basic:
+	    
+% Colors have three coordinates. Get one. 
+
+    def X(expr A) =   
+        redpart A 
+    enddef;
+
+    def Y(expr A) = 
+        greenpart A 
+    enddef;
+
+    def Z(expr A) =  
+        bluepart A 
+    enddef;
+
+% The length of a vector.
+
+    def conorm(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. 
+% Also check getangle below.
+
+    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 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;
+
+% 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;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% 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;
+
+    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;
+
+    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, unless SphericalDistortion is true.
+% f must not be on a line parallel to zz and that contains the
+% viewcentr.
+ 
+    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 eta < 0.03:  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER
+		        eta := 0.03;
+		    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.
+
+    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;
+	
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% 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;
+
+% 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, Dna, Dnb, al;
+	numeric ind, G;
+	color vec[], 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(
+	    Dna:=ncrossprod(ncrossprod(vec[ind-1],vec[ind]),vec[ind-1]);
+	    Dnb:=ncrossprod(ncrossprod(vec[ind],vec[ind-1]),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 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;
+
+% 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
+	    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 cirath, ind, vec;
+	numeric ind;
+	color vec[];
+	path cirath;
+	for ind=1 upto 36:
+	  vec[ind] = CenterPos+OneAxe*cosd(ind*10)+OtherAxe*sind(ind*10);
+	endfor;
+	cirath = rp( vec1 ) 
+	  for ind=2 upto 36: 
+	    ...rp( vec[ind] )
+	  endfor
+	  ...cycle;
+	( cirath ) 
+      endgroup
+    enddef;
+
+% Shadow of an ellipse on the air.
+
+    def ellipticshadowpath(expr CenterPos, OneAxe, OtherAxe ) =
+      begingroup
+	save cirath, ind, vec;
+	numeric ind;
+	color vec[];
+	path cirath;
+	for ind=1 upto 36:
+	  vec[ind] = CenterPos+OneAxe*cosd(ind*10)+OtherAxe*sind(ind*10);
+	endfor;
+	cirath = rp( cb( vec1 ) ) 
+	  for ind=2 upto 36: 
+	    ...rp( cb( vec[ind] ) )
+	  endfor
+	  ...cycle;
+	( cirath ) 
+      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 )( text ConverterFunc )=
+      begingroup
+	save i, outpath;
+	numeric i;
+	path outpath;
+	outpath = for i=0.25 step 0.25 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;
+            numeric ind, decision;
+            color vecx, vecy, vec[], goodangulmom, view;
+            path cirath;
+	    view = f-CenterPos;
+	    decision = cdotprod( view, AngulMom );
+	    if decision < 0: 
+	        goodangulmom = -AngulMom;
+	    else:
+		goodangulmom = AngulMom;
+	    fi;
+            vecx = ncrossprod( view, goodangulmom );
+	    decision := getangle( view, goodangulmom );
+	    if decision > 0.5:                  %%%%%%%%%%%%%%% DANGER %%%
+                vecy = ncrossprod( goodangulmom, vecx );
+		for ind=1 upto 36:
+		    vec[ind] := vecx*cosd(ind*10) + vecy*sind(ind*10);
+		    vec[ind] := CenterPos + vec[ind]*Radius;
+                endfor;
+		cirath = rp( vec1 ) 
+		         for ind=2 upto 36: 
+			     ...rp( vec[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 cirath, vecx, vecy, ind, view;
+            numeric ind;
+            color vecx, vecy, vec[], view;
+            path cirath;
+	    view = f-Pos;
+            vecx = N( (-Y(view), X(view), 0) );
+            vecy = ncrossprod( view, vecx );
+            for ind=1 upto 36:
+                vec[ind] := vecx*cosd(ind*10) + vecy*sind(ind*10);
+		vec[ind] := Pos + vec[ind]*Radius; 
+            endfor;
+            cirath = rp( vec1 ) 
+	        for ind=2 upto 36: 
+		    ...rp( vec[ind] )
+		endfor
+		...cycle;
+            ( cirath ) 
+        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 auxil, auxih;
+	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 );
+	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, rect, 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, themargin, thelengthofc, thesubpath, fullpath;
+	save newlen, finalpath, i, auxpath, pa, pb, auxt, bigcirc;
+	save startt, endt;
+	path basepath, thesubpath, fullpath, finalpath, auxpath;
+	path bigcirc;
+	numeric themargin, newlen, i, auxt, startt, endt;
+	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;
+            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 auxil, auxih;
+            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;
+      
+% 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;
+	path patha, pathb, pathc;
+	numeric noption;
+	boolean hashole, poption;
+	color auxv, vv;
+	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:
+%	  draw patha withcolor green; show patha;
+%	  draw pathb withcolor green; show pathb;
+	  hashole := (-1,-1) <> ( patha intersectiontimes pathb );
+	  if hashole:
+	    pathc := buildcycle( patha, pathb ); % I don't get it!
+	    %fill pathc withcolor red;            % see fakehole.mp
+	  fi;
+	  noption = cdotprod( LenVec, vv );
+	  if  noption > (conorm(LenVec)**2):
+	    draw pathb;
+	    if hashole:
+	      draw pathc;
+	    fi;
+	  elseif noption < 0:
+	    draw patha;
+	    if hashole:
+	      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;
+    
+% 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;
+
+% 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;
+	      
+% 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!
+	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;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% 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.
+
+    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.
+
+    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;    
+
+% 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;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% 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, polyhedric vertex approximation):
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    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;
+
+    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;
+    
+% 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;
+	save jumplim, trypoi, halfpo, expfact, randcost, avercost;
+	numeric j, iternum, prevcost, factry, randcost, avercost, actfact;
+	numeric jumplim, expfact;
+	color actpoi, trypoi, halfpo;
+	iternum = 100;                                       % DANGER!
+	factry = 0.25;                                       % DANGER!
+	jumplim = 2.5;                                       % DANGER!
+	prevcost = calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,IniV);
+	trypoi = IniV;
+	show prevcost;
+	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;
+	show prevcost;
+	( trypoi )
+      endgroup;
+    enddef;
+	  
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Part V (strictly two-dimensional):
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 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 := sqrt( La**2 - (xa-x)**2 );
+                                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;
+
+% 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;
+      
+% EOF