featpost (16apr12)

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

1 files changed, 106 insertions, 70 deletions

diff --git a/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
index 85417a17271..091d867fc66 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) 2011
+% Copyright (C) 2012
 
 % 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.5";
+    message "Preloading FeatPost macros, version 0.8.6";
 
     warningcheck := 0;
     background := 0.987white;
@@ -60,13 +60,12 @@
 	    OriginProjPagePos := (Shifts,148.45mm); % This should be the
 	                                            % page center. 
 
-	    ParallelProj := false; % Kind of perspective 
-				       % Can't have both true
-	    SphericalDistortion := false; % Kind of lens
+	    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
+	    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)
@@ -243,9 +242,9 @@
 % The following routine is used by circularsheet and may be used to
 % rotate vectors elliptically.
 
- vardef planarrotation( expr VecX, VecY, TheAngle ) =
-   ( VecX*cosd( TheAngle ) + VecY*sind( TheAngle ) )
- enddef;
+    vardef planarrotation( expr VecX, VecY, TheAngle ) =
+      ( VecX*cosd( TheAngle ) + VecY*sind( TheAngle ) )
+    enddef;
 
 % The following routine could be used by kindofcube and may be used to
 % rotate polyhedra (must cycle through all Vs before calling makeface).
@@ -368,21 +367,21 @@
         endgroup
     enddef;
 
- vardef cstr( expr Cl ) =
-   "(" &
-     decimal(X(Cl)) &
-     "," &
-     decimal(Y(Cl)) &
-     "," &
-     decimal(Z(Cl)) &
-     ")"
- 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;
+    vardef bstr( expr bv ) =
+      save bstring; string bstring;
+      if bv: bstring = "true"; else: bstring = "false"; fi;
+      bstring
+    enddef;
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%% Fundamental:
@@ -1388,6 +1387,7 @@
     enddef;
 
 % Probably the last algorithm I'm going to write for featpost...
+% Wrong. The last is ultraimprovertex.    
     
     def spheroidshadow( expr CentrPoi, NorthPoleVec, Ray ) =
       begingroup
@@ -2203,36 +2203,36 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 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;
+    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
@@ -3576,7 +3576,7 @@
 	numeric jumplim, expfact;
 	color actpoi, trypoi, halfpo, auxa, auxb, auxc;
 	trypoi = IniV;
-	iternum = 28;                                        % DANGER!
+	iternum = 24;                                        % DANGER!
 	factry = 0.25;                                       % DANGER!
 	jumplim = 2.5;                                       % DANGER!
 	prevcost = calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
@@ -3594,26 +3594,62 @@
 	prevcost := calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
 	message "phase one cost: " & decimal(prevcost)
 	                           & " at point " & cstr(trypoi);
-	for j=0 upto iternum:
-	  expfact := mexp(-j*200/iternum);                    % DANGER!
-	  actpoi := trypoi+factry*expfact*randomfear;
-	  randcost :=
-	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,actpoi);
-	  halfpo := 0.5[trypoi,actpoi];
-	  avercost :=
-	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,halfpo);
-	  actfact := calculatejump(avercost,prevcost,randcost,jumplim);
-	  trypoi := actfact[trypoi,actpoi];
-	  prevcost :=
-	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
-%	  show prevcost;
-	endfor;
-	message "final cost: " & decimal(prevcost)
-	                       & " at point " & cstr(trypoi);
+ 	for j=0 upto iternum:
+ 	  expfact := mexp(-j*200/iternum);                    % DANGER!
+ 	  actpoi := trypoi+factry*expfact*randomfear;
+ 	  randcost :=
+ 	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,actpoi);
+ 	  halfpo := 0.5[trypoi,actpoi];
+ 	  avercost :=
+ 	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,halfpo);
+ 	  actfact := calculatejump(avercost,prevcost,randcost,jumplim);
+ 	  trypoi := actfact[trypoi,actpoi];
+ 	  prevcost :=
+ 	    calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
+ %	  show prevcost;
+ 	endfor;
+ 	message "final cost: " & decimal(prevcost)
+ 	                       & " at point " & cstr(trypoi);
 	( trypoi )
       endgroup;
     enddef;
 
+% The approximation of the intersection of a plane, an infinite cylinder (tube)
+% and a spheroid   
+    
+    def ultraimprovertex( expr PlanPoi, PlanDir, BaseCenter, Radius, LenVec,
+	CentrPoi, NorthPoleVec, Ray, IniV ) =
+	begingroup
+	save trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, factry, focn, focs;
+	save norpoldi, major;
+	color trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, focn, focs;
+	color norpoldi;
+	numeric factry, major;
+	trypoi = IniV;
+	factry = 0.25;
+	plandi = N(PlanDir);
+	cyldi = N(LenVec);
+	norpoldi = N(NorthPoleVec);
+	major = conorm(NorthPoleVec);
+	focn = CentrPoi+(major+-+Ray)*norpoldi;
+	focs = CentrPoi-(major+-+Ray)*norpoldi;
+	for j=1 upto 50:
+	  auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
+	  auxb := cyldi*cdotprod(trypoi-BaseCenter,cyldi)-trypoi+BaseCenter;
+	  auxb := N(auxb)*(conorm(auxb)-Radius);
+	  auxc := focn-trypoi;
+	  auxd := focs-trypoi;
+	  auxc := N(auxc+auxd)*(conorm(auxc)+conorm(auxd)-2*major);
+	  trypoi := trypoi+factry*(auxa+auxb+auxc);
+	endfor;
+	( trypoi )
+      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;