From 45685e2c30714eae1fb1e722a26ca3d24b2f679e Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Sun, 4 Aug 2013 23:27:56 +0000 Subject: featpost git-svn-id: svn://tug.org/texlive/trunk@31346 c570f23f-e606-0410-a88d-b1316a301751 --- .../metapost/featpost/featpost3Dplus2D.mp | 253 +++++++++++++++++++-- 1 file changed, 234 insertions(+), 19 deletions(-) (limited to 'Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp') 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 majorRay: + 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