featpost

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

1 files changed, 234 insertions, 19 deletions

diff --git a/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp
index 091d867fc66..8a005bfd6d4 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) 2012
+% Copyright (C) 2013
 
 % 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.6";
+    message "Preloading FeatPost macros, version 0.8.7";
 
     warningcheck := 0;
     background := 0.987white;
@@ -64,9 +64,6 @@
 				                         % 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 
@@ -83,7 +80,9 @@
 	    SubColor := 0.35white;          %        fillfacewithlight
 	    LightSource := 10*(4,-3,4);     % This also
 	    OverRidePolyhedricColor:=false; % And also this
-	    	    
+	    ShadowOn := false; % Some objects may block the light and
+	    HoriZon := 0;      % cast a shadow on a horizontal plane at this Z
+	    	    	    
 	    TableC0 := 0.85white;           % grey        %% G N U P L O T
 	    TableC1 := red;                 % red         %% 
 	    TableC2 := ( 0.2, 0.2, 1.0 );   % blue        %%    colors 
@@ -107,8 +106,8 @@
 
 	    Nobjects := 0;                  % getready and doitnow
 
-	    TDAtiplen := 0.05;              % tdarrow
-	    TDAhalftipbase := 0.02;         % Three-Dimensional
+	    TDAtiplen := 0.05;              % tdarrow and tdcircarrow
+	    TDAhalftipbase := 0.02;         % Three-Dimensional (Circular)
 	    TDAhalfthick := 0.01;           % Arrow
 
 	    NCL := 0;                       % closedline
@@ -240,7 +239,7 @@
     enddef;
     
 % The following routine is used by circularsheet and may be used to
-% rotate vectors elliptically.
+% rotate vectors elliptically. Also used by tdcircarrow.
 
     vardef planarrotation( expr VecX, VecY, TheAngle ) =
       ( VecX*cosd( TheAngle ) + VecY*sind( TheAngle ) )
@@ -664,7 +663,7 @@
       endgroup
     enddef;
 
-% 3D arrow.
+% 3D straight arrow.
 
     def tdarrow(expr FromPos, ToTip ) =
       begingroup
@@ -695,6 +694,41 @@
       endgroup
     enddef;
 
+% 3D circular arrow.
+
+    def tdcircarrow(expr CenterPos, AngulMom, Ray, StartAngle, Amplituda ) =
+      begingroup
+	save veca, vecb, vecc, vecd, a, b, c, d, p, stepa, numa, anga, angb;
+	save signus, ca, da, i;
+	color veca, vecb, vecc, vecd;
+	pair a, b, c, d, ca, da, aa;
+	numeric stepa, numa, anga, angb, signus, i;
+	path p;
+	signus = Amplituda/abs(Amplituda);
+	stepa = 6signus; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER.
+	vecd = ncrossprod( CenterPos-f, AngulMom );
+	vecc = ncrossprod( vecd, AngulMom );
+	anga = signus*180*(TDAhalfthick/TDAhalftipbase)*TDAtiplen/(3.14159*Ray);
+	angb = signus*180*TDAtiplen/(3.14159*Ray);
+	numa = StartAngle+Amplituda;
+	a = rp(Ray*planarrotation(vecc,vecd,StartAngle+anga));
+	b = rp(Ray*planarrotation(vecc,vecd,numa+angb));
+	c = rp((Ray+TDAhalftipbase)*planarrotation(vecc,vecd,numa));
+	d = rp((Ray-TDAhalftipbase)*planarrotation(vecc,vecd,numa));
+	ca = rp((Ray+TDAhalfthick)*planarrotation(vecc,vecd,numa));
+	da = rp((Ray-TDAhalfthick)*planarrotation(vecc,vecd,numa));
+	aa = rp((Ray-TDAhalfthick)*planarrotation(vecc,vecd,StartAngle));
+	p = for i=StartAngle step stepa until numa:
+	      rp((Ray+TDAhalfthick)*planarrotation(vecc,vecd,i))..
+	    endfor ca--c--b--d--da..
+	    for i=numa-stepa step -stepa until StartAngle:
+	      rp((Ray-TDAhalfthick)*planarrotation(vecc,vecd,i))..
+	    endfor aa--a--cycle;
+	unfill p;
+	draw p
+      endgroup
+    enddef;
+
 % Draw lines with a better expression of three-dimensionality.
 
     def emptyline(expr JoinP,ThickenFactor,OutCol,InCol,theN,EmptyFrac,sN)(text LinFunc) = 
@@ -834,7 +868,7 @@
       begingroup
 	save ind;
 	numeric ind;
-	( for ind=1 upto 36: 
+	( for ind=0 upto 35: 
 	    rp( CenterPos+planarrotation(OneAxe,OtherAxe,ind*10) )...
 	  endfor cycle ) 
       endgroup
@@ -1387,7 +1421,11 @@
     enddef;
 
 % Probably the last algorithm I'm going to write for featpost...
-% Wrong. The last is ultraimprovertex.    
+% Wrong. The last is ultraimprovertex. 
+% Wrong again. The last is necplusimprovertex.
+% And again wrong. The last is tdcircarrow.
+% Well, what can I say, really the last is ellipsoid.
+% Wait! It is torushadow!
     
     def spheroidshadow( expr CentrPoi, NorthPoleVec, Ray ) =
       begingroup
@@ -1490,6 +1528,78 @@
       endgroup
     enddef;
 
+% Another brute-force algorythm. It's advisable to use three orthogonal axes.
+
+    def ellipsoid( expr Centr, AxOne, AxTwo, AxThr ) =
+      begingroup
+	save count, i, j, axx, axy, cyc, cy, di, leng;
+	numeric count, i, j, leng;
+	color axx, axy, di[];
+	path cy, cyc[];
+	di1 = AxOne;
+	di2 = AxTwo;
+	di3 = AxThr;
+	count = 0;
+	for i=1 upto 3:
+	  if i=1:
+	    axx := di2;
+	    axy := di3;
+	  elseif i=2:
+	    axx := di1;
+	    axy := di3;
+	  else:
+	    axx := di1;
+	    axy := di2;
+	  fi;
+	  for j=5 step 10 until 175:
+	    cyc[incr(count)] = ellipticpath( Centr, di[i],
+		Centr + planarrotation( axx, axy, j ) );
+	  endfor;
+	endfor;
+	cy = cyc1;
+	for i=2 upto count-1:
+	  cy := twocyclestogether( cy, cyc[i] );
+	endfor;
+	leng = (length cy) - 1;
+	( point 0 of cy for i = 1 upto leng: ..point i of cy endfor ..cycle )
+      endgroup
+    enddef;
+  
+    def ellipsoidshadow( expr Centr, AxOne, AxTwo, AxThr ) =
+      begingroup
+	save count, i, j, axx, axy, cyc, cy, di, leng;
+	numeric count, i, j, leng;
+	color axx, axy, di[];
+	path cy, cyc[];
+	di1 = AxOne;
+	di2 = AxTwo;
+	di3 = AxThr;
+	count = 0;
+	for i=1 upto 3:
+	  if i=1:
+	    axx := di2;
+	    axy := di3;
+	  elseif i=2:
+	    axx := di1;
+	    axy := di3;
+	  else:
+	    axx := di1;
+	    axy := di2;
+	  fi;
+	  for j=5 step 10 until 175:
+	    cyc[incr(count)] = ellipticshadowpath( Centr, di[i],
+		Centr + planarrotation( axx, axy, j ) );
+	  endfor;
+	endfor;
+	cy = cyc1;
+	for i=2 upto count-1:
+	  cy := twocyclestogether( cy, cyc[i] );
+	endfor;
+	leng = (length cy) - 1;
+	( point 0 of cy for i = 1 upto leng: ..point i of cy endfor ..cycle )
+      endgroup
+    enddef;
+  
 % You can't see through this hole. f must not be on the hole axis.
 % Not yet documented because "buildcycle" doesn't work properly.
     
@@ -1761,6 +1871,9 @@
 	save tmoment;
 	color tmoment;
 	angstep= 4; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER!
+	if ShadowOn:
+	  torushadow( Tcenter, Tmoment, Bray, Sray );
+	fi;
 	viewline = f-Tcenter;
 	if cdotprod( viewline, Tmoment ) < 0:
 	  tmoment = -Tmoment;
@@ -1872,6 +1985,28 @@
       endgroup
     enddef;
 
+% The shadow of a torus
+    
+    def torushadow( expr Tcenter, Tmoment, Bray, Sray ) =
+      begingroup
+	save theplace, viewline, tmoment, refpair, sideaxe, nearaxe, i;
+	color theplace, viewline, tmoment, refpair, sideaxe, nearaxe;
+	numeric i;
+	viewline = f-Tcenter;
+	if cdotprod( viewline, Tmoment ) < 0:
+	  tmoment = -Tmoment;
+	else:
+	  tmoment = Tmoment;
+	fi;
+	sideaxe = Bray*ncrossprod( tmoment, viewline );
+	nearaxe = Bray*ncrossprod( sideaxe, tmoment );
+	for i=1 upto 60:  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER!
+	  theplace := Tcenter + planarrotation( sideaxe, nearaxe, 6i ); %%
+	  fill rigorousfearshadowpath( theplace, Sray );
+	endfor;
+      endgroup
+    enddef;
+    
 % Take a "quarter" of a donut
 % (no longer under construction but contains a bug). 
 
@@ -3615,37 +3750,117 @@
     enddef;
 
 % The approximation of the intersection of a plane, an infinite cylinder (tube)
-% and a spheroid   
+% and a spheroid (an oblate spheroid, a prolate spheroid or a perfect sphere)  
     
     def ultraimprovertex( expr PlanPoi, PlanDir, BaseCenter, Radius, LenVec,
 	CentrPoi, NorthPoleVec, Ray, IniV ) =
 	begingroup
 	save trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, factry, focn, focs;
-	save norpoldi, major;
+	save norpoldi, major, majortmp, nordi, tmpdi;
 	color trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, focn, focs;
-	color norpoldi;
-	numeric factry, major;
+	color norpoldi, nordi, tmpdi;
+	numeric factry, major, majortmp;
 	trypoi = IniV;
 	factry = 0.25;
 	plandi = N(PlanDir);
 	cyldi = N(LenVec);
 	norpoldi = N(NorthPoleVec);
 	major = conorm(NorthPoleVec);
-	focn = CentrPoi+(major+-+Ray)*norpoldi;
-	focs = CentrPoi-(major+-+Ray)*norpoldi;
+	if major>Ray:
+	  focn := CentrPoi+(major+-+Ray)*norpoldi;
+	  focs := CentrPoi-(major+-+Ray)*norpoldi;
+	fi;
 	for j=1 upto 50:
+	  if major<Ray:
+	    tmpdi := trypoi-CentrPoi;
+	    nordi := N(tmpdi-cdotprod(tmpdi,norpoldi)*norpoldi);
+	    focn := CentrPoi+(Ray+-+major)*nordi;
+	    focs := CentrPoi-(Ray+-+major)*nordi;
+	  elseif major=Ray:
+	    focn := CentrPoi;
+	    focs := CentrPoi;
+	  fi;  
 	  auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
 	  auxb := cyldi*cdotprod(trypoi-BaseCenter,cyldi)-trypoi+BaseCenter;
 	  auxb := N(auxb)*(conorm(auxb)-Radius);
 	  auxc := focn-trypoi;
 	  auxd := focs-trypoi;
-	  auxc := N(auxc+auxd)*(conorm(auxc)+conorm(auxd)-2*major);
+	  if major>Ray:
+	    auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*major);
+	  else:
+	    auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*Ray);
+	  fi;  
 	  trypoi := trypoi+factry*(auxa+auxb+auxc);
 	endfor;
 	( trypoi )
       endgroup
     enddef;
 
+% The approximation of the intersection of a plane and two prolate spheroids
+
+    def necplusimprovertex( expr PlanPoi, PlanDir,
+	CentrPoiA, NorthPoleVecA, RayA,
+        CentrPoiB, NorthPoleVecB, RayB, IniV ) =
+        begingroup
+        save trypoi, auxa, auxb, auxc, auxd, plandi, factry, focni, focsi;
+        save norpoldi, norpoldj, maior, major, focnj, focsj, auxe;
+        color trypoi, auxa, auxb, auxc, auxd, plandi, focni, focsi, focnj, focsj;
+        color norpoldi, norpoldj, auxe;
+        numeric factry, maior, major;
+        trypoi = IniV;
+        factry = 0.25;
+        plandi = N(PlanDir);
+        norpoldi = N(NorthPoleVecA);
+        maior = conorm(NorthPoleVecA);        
+        norpoldj = N(NorthPoleVecB);
+        major = conorm(NorthPoleVecB);
+        focni = CentrPoiA+(maior+-+RayA)*norpoldi;
+        focsi = CentrPoiA-(maior+-+RayA)*norpoldi;
+        focnj = CentrPoiB+(major+-+RayB)*norpoldj;
+        focsj = CentrPoiB-(major+-+RayB)*norpoldj;
+        for j=1 upto 50:
+          auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
+          auxb := focni-trypoi;
+          auxe := focsi-trypoi;
+          auxc := focnj-trypoi;
+          auxd := focsj-trypoi;
+          auxb := (N(auxb)+N(auxe))*(conorm(auxb)+conorm(auxe)-2*maior);
+          auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*major);
+          trypoi := trypoi+factry*(auxa+auxb+auxc);
+        endfor;
+        ( trypoi )
+      endgroup
+    enddef;
+
+% The approximation of the intersection of a straight line and
+% one prolate spheroid
+
+    def intersectprolatespheroid( expr CentrPoi, NorthPoleVec, Ray,
+        LinePoi, LineDir ) =
+        begingroup
+	save trypoi, factry, linedi, norpol, focn, focs, maior, j;
+	save auxa, auxb, auxc, auxd;
+	color trypoi, linedi, norpol, focn, focs;
+	color auxa, auxb, auxc, auxd;
+	numeric factry, maior, j;
+        trypoi = LinePoi;
+        factry = 0.25;
+        linedi = N(LineDir);
+        norpol = N(NorthPoleVec);
+        maior = conorm(NorthPoleVec);
+        focn = CentrPoi+(maior+-+Ray)*norpol;
+        focs = CentrPoi-(maior+-+Ray)*norpol;
+        for j=1 upto 50:
+          auxb := focn-trypoi;
+          auxd := focs-trypoi;
+          auxa := (N(auxb)+N(auxd))*(conorm(auxb)+conorm(auxd)-2*maior);
+          auxc := linedi*cdotprod(auxa,linedi);
+          trypoi := trypoi+factry*auxc;
+        endfor;
+        ( trypoi )
+      endgroup
+    enddef;
+
 % 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)