featpost (10oct14)

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

1 files changed, 249 insertions, 78 deletions

diff --git a/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
index 8a005bfd6d4..52a5f7789ed 100644
--- a/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
+++ b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
@@ -7,7 +7,7 @@
 % P. Jørgensen
 % S. Pakin
 % 
-% Copyright (C) 2013
+% Copyright (C) 2014
 
 % This set of macros adds a lot of features to
 % the MetaPost language and eases the production of 
@@ -23,7 +23,7 @@
 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 % GNU General Public License for more details.
 
-    message "Preloading FeatPost macros, version 0.8.7";
+    message "FeatPost-0.8.8";
 
     warningcheck := 0;
     background := 0.987white;
@@ -204,7 +204,6 @@
     enddef;
 
 % The normalized crossproduct of two vectors. 
-% Also check getangle below.
 
     def ncrossprod(expr A, B) = 
         N( ccrossprod( A, B ) )
@@ -290,6 +289,18 @@
       endgroup
     enddef;
 
+% Maybe you would like to calculate the angular arguments of kindofcube...
+    
+    def getanglepair( expr InVec ) =
+      begingroup
+	save alphaone, alphatwo;
+	numeric alphaone, alphatwo;
+	alphaone = angle( ( X(InVec), Y(InVec) ) );
+	alphatwo = angle( ( X(InVec) ++ Y(InVec), Z(InVec) ) );
+	( (alphaone,alphatwo) ) 
+      endgroup
+    enddef;
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%% Auxiliary:
 	    
@@ -366,6 +377,8 @@
         endgroup
     enddef;
 
+% The following two procedures are useful for getready.
+    
     vardef cstr( expr Cl ) =
       "(" &
 	decimal(X(Cl)) &
@@ -436,7 +449,8 @@
     enddef;
 
 % Much improved rigorous pseudo-projection algorithm that follows 
-% an idea from Cristian Barbarosie. This makes shadows.
+% an idea from Cristian Barbarosie.
+% This makes shadows caused by a light source point.
 
     def cb(expr R) =
       begingroup
@@ -477,6 +491,14 @@
       endgroup
     enddef;
 	
+% Vanishing point.     
+    
+    def vp( expr D ) =
+      begingroup
+	( rp( lineintersectplan( f, D, viewcentr, f-viewcentr) ) )
+      endgroup
+    enddef;
+    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%% Basic Functions:
 	    
@@ -1426,6 +1448,10 @@
 % And again wrong. The last is tdcircarrow.
 % Well, what can I say, really the last is ellipsoid.
 % Wait! It is torushadow!
+% Bullshit. Only death can stop me. Now it's revolparab.
+% Continuing with torus accessories...
+% And I forgot to mention intersectprolatespheroid!!!
+% The first of Red October, 2014, added vp.
     
     def spheroidshadow( expr CentrPoi, NorthPoleVec, Ray ) =
       begingroup
@@ -1528,7 +1554,8 @@
       endgroup
     enddef;
 
-% Another brute-force algorythm. It's advisable to use three orthogonal axes.
+% Another brute-force algorythm.
+% It's advisable to use three orthogonal axes.
 
     def ellipsoid( expr Centr, AxOne, AxTwo, AxThr ) =
       begingroup
@@ -1600,6 +1627,104 @@
       endgroup
     enddef;
   
+    def revolparab( expr BaseCenter, ParabTip, BaseRay ) =
+      begingroup
+	save bcpt, conetip, fakex, fakey, fakez, tipview, ellicenter, coneview;
+	save tanefe, majorvec, minorvec, cutvec, auxpoi, crux, auxx, baseview;
+	save xzero, yzero, xdelta, fakea, xpos, ypos, xneg, yneg, l, apertur;
+	save auxy, maxy, ellmaxang, ymin, xefe, auxray, auxcos, auxsin, ste, a;
+	save auxpath, tippath;
+	save conda, condb, condc;
+	color bcpt, conetip, fakex, fakey, fakez, tipview, ellicenter, coneview;
+	color tanefe, majorvec, minorvec, cutvec, auxpoi, crux, auxx, baseview;
+	numeric xzero, yzero, xdelta, fakea, xpos, ypos, xneg, yneg, l, apertur;
+	numeric auxy, maxy, ellmaxang, ymin, xefe, auxray, auxcos, auxsin, ste;
+	numeric a;
+	path auxpath, tippath;
+	boolean conda, condb, condc;
+	bcpt = BaseCenter-ParabTip;
+	conetip = BaseCenter-2*bcpt;
+	maxy = conorm(bcpt);
+	apertur = angle(2*maxy,BaseRay);
+	fakey = N(bcpt);
+	coneview = f-conetip;
+	tipview = f-ParabTip;
+	baseview = f-BaseCenter;
+	a = getangle(coneview,bcpt);
+	conda = a<=apertur;
+	condb = a>=180-apertur;
+	fakez = ncrossprod( tipview, bcpt );
+	fakex = ccrossprod( fakey, fakez );
+	xzero = cdotprod( fakex, tipview );
+	yzero = cdotprod( fakey, tipview );
+	fakea = maxy/(BaseRay**2);
+	condc = yzero>=fakea*(xzero**2);
+	if (conda and condc) or condb:
+	  auxpath = rigorouscircle( BaseCenter, bcpt, BaseRay );
+	  unfill auxpath;
+	  draw auxpath;
+	else:
+	  xdelta = sqrt(xzero**2 - yzero/fakea);
+	  xpos = xzero + xdelta;
+	  xneg = xzero - xdelta;
+	  ypos = fakea*(xpos**2);
+	  yneg = fakea*(xneg**2);
+	  auxy = 0.5[ypos,yneg];
+	  ellicenter = ParabTip+fakex*xzero+fakey*auxy;
+	  if yneg<ypos:
+	    ymin=yneg;
+	    xefe=xneg;
+	  else:
+	    ymin=ypos;
+	    xefe=xpos;
+	  fi;
+	  tanefe = ParabTip +fakex*xefe+fakey*ymin;
+	  minorvec = xdelta*fakez;
+	  majorvec = tanefe-ellicenter;
+	  if (yneg<maxy) and (ypos<maxy):
+	    draw ellipticpath( ellicenter, majorvec, minorvec );
+	  else:
+	    auxpoi = lineintersectplan(conetip,conetip-f,BaseCenter,bcpt);
+	    auxray = conorm( BaseCenter-auxpoi );
+	    auxcos = BaseRay/auxray;
+	    auxsin = 1 +-+ auxcos;
+	    if cdotprod(fakex,auxpoi-BaseCenter)>0:
+	      auxx = fakex;
+	    else:
+	      auxx = -fakex;
+	    fi;
+	    crux = BaseCenter+BaseRay*(auxx*auxcos +fakez*auxsin);
+	    cutvec = crux-ellicenter;
+	    auxcos := cdotprod(cutvec,N(majorvec))/conorm(majorvec);
+	    auxsin := cdotprod(cutvec,N(minorvec))/conorm(minorvec);      
+	    ellmaxang = angle(auxcos,auxsin);
+	    ste = ellmaxang/18; %%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER %%%%%%%%%%%%
+	    tippath =
+	      rp(ellicenter+planarrotation(majorvec,minorvec,-ellmaxang)) 
+	      for l=ste-ellmaxang step ste until ellmaxang:
+	      ..rp(ellicenter+planarrotation(majorvec,minorvec,l))
+	    endfor --rp(ellicenter+planarrotation(majorvec,minorvec,ellmaxang));
+	    if cdotprod( baseview, bcpt )<0:
+	      auxpath = rigorouscone(true,BaseCenter,bcpt,BaseRay,conetip);
+	      numeric len;
+	      len = length auxpath;
+	      auxpath := subpath (1,len-1) of auxpath;
+	      auxpath := tippath--(reverse auxpath)--cycle;
+	      unfill auxpath;
+	      draw auxpath;
+	    else:
+	      tippath := tippath--cycle;
+	      unfill tippath;
+	      draw tippath;
+	      auxpath = rigorouscircle( BaseCenter, bcpt, BaseRay );
+	      unfill auxpath;
+	      draw auxpath;
+	    fi;
+	  fi;
+	fi
+      endgroup
+    enddef;
+    
 % You can't see through this hole. f must not be on the hole axis.
 % Not yet documented because "buildcycle" doesn't work properly.
     
@@ -1748,18 +1873,6 @@
       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 ) =
@@ -1855,6 +1968,9 @@
       endgroup
     enddef;
 	      
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Toroidal Stuff
+    
 % Take a donut. 
 
     def smoothtorus( expr Tcenter, Tmoment, Bray, Sray ) =
@@ -2077,6 +2193,54 @@
       endgroup
     enddef;
 
+    def pointinsidetorus( expr Point, Tcenter, Tmoment, Bray, Sray ) =
+      begingroup
+	save outputboolean, height, dist, dray;
+	boolean outputboolean;
+	numeric height, dist, dray;
+	height = cdotprod(N(Tmoment),Point-Tcenter);
+	if abs(height)>=Sray:
+	  outputboolean = false;
+	else:
+	  dist = conorm(Point-Tcenter-height*N(Tmoment));
+	  dray = Sray +-+ abs(height);
+	  if (dist<Bray+dray) and (dist>Bray-dray):
+	    outputboolean = true;
+	  else:
+	    outputboolean = false;
+	  fi;
+	fi;
+	( outputboolean )
+      endgroup
+    enddef;
+	
+    def pointrelativetotorus( expr Point, Tcenter, Tmoment, Bray, Sray ) =
+      begingroup
+	save height, dist;
+	numeric height;
+	color dist;
+	height = cdotprod(N(Tmoment),Point-Tcenter);
+	dist = N(Point-Tcenter-height*N(Tmoment));
+	( conorm(Point-Bray*dist)-Sray )
+      endgroup
+    enddef;
+
+    def intersectorus( expr Tcenter, Tmoment, Bray, Sray, LinePoi, LineDir ) =
+      begingroup
+	save trypoi, factry, linedi;
+	color trypoi, linedi;
+	numeric factry, j, auxd;
+	trypoi = LinePoi;
+        factry = 0.25;
+        linedi = N(LineDir);
+	for j= 1 upto 50:
+	  auxd := pointrelativetotorus( trypoi, Tcenter, Tmoment, Bray, Sray );
+	  trypoi := trypoi+factry*linedi*auxd;
+	endfor;
+	( trypoi )
+      endgroup
+    enddef;
+    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%% Non-standard objects:
 	    
@@ -2257,7 +2421,7 @@
     enddef;    
 
 % Another point is that I want to draw trajectories in space and
-% dependant on velocity: VecFunc( position, velocity ).
+% dependent on velocity: VecFunc( position, velocity ).
 % This time is fourth-order Runge-Kutta.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHANGES PrintStep!!!!
     
@@ -3537,7 +3701,7 @@
     enddef;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Part IV (automatic perspective tuning, polyhedric vertex approximation):
+% Part IV (automatic perspective tuning and minimization):
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     def randomfear =
@@ -3690,6 +3854,67 @@
       endgroup
     enddef;
 
+% Minimization routine for scalar functions like y=f(x) where an initial
+% triplet (x1,x2,x3) with x1<x2<x3 is given as a parabolic squeleton that
+% provides a way to search for the smallest value of y (if iterated)
+
+    def minimizestep( expr Abcisscolor )( text PlainFunc ) =
+      begingroup
+	save xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
+	save colout;
+	numeric xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
+	color colout;
+	xa = X( Abcisscolor );
+	xb = Y( Abcisscolor );
+	xc = Z( Abcisscolor );
+	ya = PlainFunc(xa); 
+	yb = PlainFunc(xb);
+	yc = PlainFunc(xc);
+	if ya = yb:
+	  colout = (-0.125[xa,xb],xb,xc);
+	elseif yb = yc:
+	  colout = (xa,xb,1.125[xb,xc]);
+	else:
+	  if (yb>ya) or (yb>yc):
+	    show  Abcisscolor;
+	    message "                             Unable to minimizestep!";
+	  fi;
+	  den = (xb-xc)*((xa**2)-(xb**2))-(xa-xb)*((xb**2)-(xc**2));
+	  if abs(den) < 0.0005:
+	    show den;
+	    message "                             Unable to minimizestep!";
+	  fi;
+	  coeb = ((yb-yc)*((xa**2)-(xb**2))-(ya-yb)*((xb**2)-(xc**2)))/den;
+	  coec = ((xb-xc)*(ya-yb)-(xa-xb)*(yb-yc))/den;
+	  xd = -0.5*coeb/coec;
+	  yd = PlainFunc( xd );
+	  if ((xa<xd) and (xd<xb)):
+	    if (yd<yb):
+	      colout = (xa,xd,xb);
+	    else:
+	      colout = (xd,xb,xc);
+	    fi;
+	  elseif ((xb<xd) and (xd<xc)):
+	    if (yd<yb):
+	      colout = (xb,xd,xc);
+	    else:
+	      colout = (xa,xb,xd);
+	    fi;
+	  else:
+	    aux := 0.125[xb,xc]-0.125[xb,xa];
+	    colout = (xa,0.125[xb,xa]+uniformdeviate(aux),xc);
+	  fi;
+	fi;
+	( colout )
+      endgroup
+    enddef;
+	  
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Part V:
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Intersections:
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% Check lineintersectplan above %%%%%%%%%%%%%%%
+    
     def calculatecostver(expr VerA,DisA,VerB,DisB,VerC,DisC,TryV) =
       begingroup
 	save a, b, c;
@@ -3700,7 +3925,8 @@
 	( ( a ++ b ++ c )**2 )
       endgroup
     enddef;
-    
+
+% Approximation of the intersection of three spheres.    
 % Be aware of the next three danger parameters.
     
     def improvertex( expr VerA, DisA, VerB, DisB, VerC, DisC, IniV ) =
@@ -3861,64 +4087,10 @@
       endgroup
     enddef;
 
-% Minimization routine for scalar functions like y=f(x) where an initial
-% triplet (x1,x2,x3) with x1<x2<x3 is given as a parabolic squeleton that
-% provides a way to search for the smallest value of y (if iterated)
-
-    def minimizestep( expr Abcisscolor )( text PlainFunc ) =
-      begingroup
-	save xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
-	save colout;
-	numeric xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
-	color colout;
-	xa = X( Abcisscolor );
-	xb = Y( Abcisscolor );
-	xc = Z( Abcisscolor );
-	ya = PlainFunc(xa); 
-	yb = PlainFunc(xb);
-	yc = PlainFunc(xc);
-	if ya = yb:
-	  colout = (-0.125[xa,xb],xb,xc);
-	elseif yb = yc:
-	  colout = (xa,xb,1.125[xb,xc]);
-	else:
-	  if (yb>ya) or (yb>yc):
-	    show  Abcisscolor;
-	    message "                             Unable to minimizestep!";
-	  fi;
-	  den = (xb-xc)*((xa**2)-(xb**2))-(xa-xb)*((xb**2)-(xc**2));
-	  if abs(den) < 0.0005:
-	    show den;
-	    message "                             Unable to minimizestep!";
-	  fi;
-	  coeb = ((yb-yc)*((xa**2)-(xb**2))-(ya-yb)*((xb**2)-(xc**2)))/den;
-	  coec = ((xb-xc)*(ya-yb)-(xa-xb)*(yb-yc))/den;
-	  xd = -0.5*coeb/coec;
-	  yd = PlainFunc( xd );
-	  if ((xa<xd) and (xd<xb)):
-	    if (yd<yb):
-	      colout = (xa,xd,xb);
-	    else:
-	      colout = (xd,xb,xc);
-	    fi;
-	  elseif ((xb<xd) and (xd<xc)):
-	    if (yd<yb):
-	      colout = (xb,xd,xc);
-	    else:
-	      colout = (xa,xb,xd);
-	    fi;
-	  else:
-	    aux := 0.125[xb,xc]-0.125[xb,xa];
-	    colout = (xa,0.125[xb,xa]+uniformdeviate(aux),xc);
-	  fi;
-	fi;
-	( colout )
-      endgroup
-    enddef;
-	  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%% Part V (strictly two-dimensional):
+%%%% Part VI (strictly two-dimensional):
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 % Verify if a path is cyclic (written by Scott Pakin)
     
     def is_cyclic expr cpath =
@@ -4246,8 +4418,7 @@
     enddef;
 
 % Shrink or swell a cyclic path without cusp points and without
-% coinciding pre and post control points. This algorithm should
-% be improved to add circular arcs on the outside of convex corners.
+% coinciding pre and post control points. 
 
     def lasermachine( expr DefinedPath, Beam, CosLimit ) =
       begingroup