From 97e4684a37e4899e8cdd42ab670b6e25dd319d1a Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Sat, 3 Dec 2011 16:47:39 +0000 Subject: featpost (17nov11) git-svn-id: svn://tug.org/texlive/trunk@24738 c570f23f-e606-0410-a88d-b1316a301751 --- Master/texmf-dist/metapost/featpost/featpost.mp | 36 - .../metapost/featpost/featpost3Dplus2D.mp | 1525 +++++++++++++++----- 2 files changed, 1132 insertions(+), 429 deletions(-) delete mode 100644 Master/texmf-dist/metapost/featpost/featpost.mp (limited to 'Master/texmf-dist/metapost') diff --git a/Master/texmf-dist/metapost/featpost/featpost.mp b/Master/texmf-dist/metapost/featpost/featpost.mp deleted file mode 100644 index 9c30b765fed..00000000000 --- a/Master/texmf-dist/metapost/featpost/featpost.mp +++ /dev/null @@ -1,36 +0,0 @@ -% 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 index b6055afae58..3c1c809eabd 100644 --- a/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp +++ b/Master/texmf-dist/metapost/featpost/featpost3Dplus2D.mp @@ -1,17 +1,21 @@ % 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. +% L. Nobre G., lnobreg@gmail.com, http://matagalatlante.org +% C. Barbarosie +% J. Schwaiger +% B. Jackowski +% P. J. Sebastião +% P. Jørgensen +% S. Pakin +% +% Copyright (C) 2011 + +% This set of macros adds a lot of features to +% the MetaPost language and eases the production of +% physics diagrams and animations. % 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 +% as published by the Free Software Foundation; either version 3 % of the License, or (at your option) any later version. % This set of macros is distributed in the hope that it will be useful, @@ -19,16 +23,17 @@ % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. - message "Preloading FeatPost macros, version 0.6.7"; + message "Preloading FeatPost macros, version 0.8.2"; warningcheck := 0; - background := 0.987white; + defaultscale := 0.75; + defaultfont := "cmss17"; % This is used by cartaxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global Variables %%%%%%%%%%% boolean ParallelProj, SphericalDistortion, FCD[], ShadowOn; - boolean OverRidePolyhedricColor; + boolean OverRidePolyhedricColor, MalcomX; numeric Nobjects, RefDist[], HoriZon, RopeColorSeq[], PhotoMarks; numeric Spread, PrintStep, PageHeight, PageWidth, ActuC, Shifts; numeric NL, npl[], NF, npf[], FC[], MaxFearLimit, TableColors; @@ -40,10 +45,6 @@ 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 @@ -52,10 +53,12 @@ Spread := 140; % Magnification - ShiftV := 105.00mm*(1,1); % Central coordinates on paper + Shifts := 105.00mm; + + ShiftV := Shifts*(1,1); % Central coordinates on paper - OriginProjPagePos := (105.00mm,148.45mm); % This should be the - % page center. + OriginProjPagePos := (Shifts,148.45mm); % This should be the + % page center. ParallelProj := false; % Kind of perspective % Can't have both true @@ -72,13 +75,10 @@ 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 + MaxFearLimit := 17; % Valid Maximum Distance from Origin HigColor := 0.85white; % These two colors are used in SubColor := 0.35white; % fillfacewithlight @@ -115,6 +115,8 @@ NCL := 0; % closedline ForePen := pencircle scaled 15pt; BackPen := pencircle scaled 9pt; + + MalcomX := false; %%% The variables PhotoMarks, PhotoPair[], PhotoPoint[] %%% and CLPath[] have NO default values. @@ -124,26 +126,42 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Very basic: -% Colors have three coordinates. Get one. +% Colors have three or four coordinates. Get one. - def X(expr A) = - redpart A + def X(expr A) = + if color A: redpart A elseif MalcomX: blackpart A else: cyanpart A fi enddef; - - def Y(expr A) = - greenpart A + + def Y(expr A) = + if color A: greenpart A else: magentapart A fi enddef; - + def Z(expr A) = - bluepart A + if color A: bluepart A else: yellowpart A fi enddef; - + + def W(expr A) = + blackpart A + enddef; + % The length of a vector. def conorm(expr A) = - ( X(A) ++ Y(A) ++ Z(A) ) + ( X(A) ++ Y(A) ++ Z(A) ) enddef; + def cmyknorm(expr A) = %% This is not good when MalcomX is true + ( X(A) ++ Y(A) ++ Z(A) ++ W(A) ) + enddef; + + def makecmyk( expr A, B ) = + ( ( X(A), Y(A), Z(A), B ) ) + enddef; + + def maketrio( expr A ) = + ( ( X(A), Y(A), Z(A) ) ) + enddef; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Vector Calculus: @@ -222,6 +240,13 @@ endgroup enddef; +% 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; + % The following routine could be used by kindofcube and may be used to % rotate polyhedra (must cycle through all Vs before calling makeface). @@ -240,6 +265,20 @@ endgroup enddef; +% Rotate a vector around another. Supposes all vectors share the same origin. + + def rotvecaroundanother( expr Angle, RotVec, FixVec ) = + begingroup + save uf, cf, xr, yr; + color uf, cf, xr, yr; + uf = N( FixVec ); + yr = ccrossprod( uf, RotVec ); + cf = uf*cdotprod( uf, RotVec ); + xr = RotVec - cf; + ( cf + planarrotation( xr, yr, Angle ) ) + endgroup + enddef; + % inplanarvolume is used by kindofcube. def inplanarvolume( expr PointPerpA, PointPerpB, Point ) = @@ -308,6 +347,10 @@ endgroup enddef; +% The following macro fits a figure in the page. +% Probably it is obsolete since MetaPost 1.000 +% Should be the last command in a figure. + def produce_auto_scale = begingroup picture storeall, scaleall; @@ -325,30 +368,29 @@ 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; + 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. +% ortogonal to vector f-viewcentr, unless SphericalDistortion is true. +% f-viewcentr must not be (anti-)parallel to zz. def rp(expr R) = begingroup @@ -383,8 +425,8 @@ squarf = cdotprod( f-viewcentr, f-viewcentr ); radio = cdotprod( R-viewcentr, f-viewcentr ); eta = 1 - radio/squarf; - if eta < 0.03: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER - eta := 0.03; + if abs((horiz,verti)) > MaxFearLimit*eta: + eta := abs((horiz,verti))/MaxFearLimit; fi; fi; projpoi = (horiz,verti)/eta; @@ -467,7 +509,7 @@ endgroup enddef; - def drawsegment( expr A, B )= + def drawsegment( expr A, B ) = begingroup if SphericalDistortion: draw pathofstraightline( A, B ); @@ -496,6 +538,24 @@ endgroup enddef; +% Orthogonal axes with prescribed lengths and labels. + + def orthaxes(expr axex, strx, axey, stry, axez, strz ) = + begingroup + save axxc, ayyc, azzc; + color axxc, ayyc, azzc; + axxc = (axex,0,0); + ayyc = (0,axey,0); + azzc = (0,0,axez); + drawarrow rp(black)..rp(axxc); + drawarrow rp(black)..rp(ayyc); + drawarrow rp(black)..rp(azzc); + label.bot( strx ,rp(axxc)); + label.bot( stry ,rp(ayyc)); + label.lft( strz ,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 @@ -579,9 +639,9 @@ def rigorouscircle( expr CenterPos, AngulMom, Radius ) = begingroup - save ind, G, Dna, Dnb, al; + save ind, G, Dn, Dna, Dnb, al, vec; numeric ind, G; - color vec[], Dna, Dnb; + color vec[], Dn, Dna, Dnb; path al; vec1 = ncrossprod( CenterPos-f, AngulMom); for ind=2 step 2 until 8: @@ -592,8 +652,9 @@ 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]) + Dn:=ncrossprod(vec[ind-1],vec[ind]); + Dna:=ncrossprod(Dn,vec[ind-1]); + Dnb:=ncrossprod(-Dn,vec[ind]) ) ..controls rp(Radius*vec[ind-1]+CenterPos+G*Dna) and rp(Radius*vec[ind] +CenterPos+G*Dnb) @@ -637,8 +698,7 @@ % Draw lines with a better expression of three-dimensionality. - def emptyline(expr JoinP,ThickenFactor,OutCol,InCol,theN,EmptyFrac,sN) - (text LinFunc) = + def emptyline(expr JoinP,ThickenFactor,OutCol,InCol,theN,EmptyFrac,sN)(text LinFunc) = begingroup save i, j, k; numeric i, j, k; @@ -667,8 +727,7 @@ % Draw space-paths of possibly closed lines making use of "getready" - def closedline( expr ThisIsClosed, theN, ForeFrac, BackFrac ) - ( text LinFunc ) = + def closedline( expr ThisIsClosed, theN, ForeFrac, BackFrac )( text LinFunc ) = begingroup save i, comm; numeric i; @@ -748,8 +807,10 @@ % 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 )= + def twocyclestogether( expr CycleA, CycleB ) = begingroup + save TheLengthOfA, TheLengthOfB, TheMargin, Leng, i; + save SubPathA, SubPathB, PolygonPath, FinalPath; numeric TheLengthOfA, TheLengthOfB, TheMargin, Leng, i; path SubPathA, SubPathB, PolygonPath, FinalPath; TheMargin = 0.02; @@ -772,19 +833,11 @@ def ellipticpath(expr CenterPos, OneAxe, OtherAxe ) = begingroup - save cirath, ind, vec; + save ind; 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 ) + ( for ind=1 upto 36: + rp( CenterPos+planarrotation(OneAxe,OtherAxe,ind*10) )... + endfor cycle ) endgroup enddef; @@ -792,27 +845,18 @@ def ellipticshadowpath(expr CenterPos, OneAxe, OtherAxe ) = begingroup - save cirath, ind, vec; + save ind; 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 ) + ( for ind=1 upto 36: + rp( cb( CenterPos+planarrotation(OneAxe,OtherAxe,ind*10) ) )... + endfor cycle ) 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)= + def labelinspace(expr KeepRatio,RefPoi,BaseVec,UpVec)(text SomeString) = begingroup save labelpic, plak, lrc, ulc, llc, centerc, aratio, newbase; picture labelpic; @@ -840,12 +884,13 @@ % It should be possible to attach some path to some surface. - def closedpathinspace( expr SomeTDPath )( text ConverterFunc )= + def closedpathinspace( expr SomeTDPath, NDivide )( text ConverterFunc ) = begingroup - save i, outpath; - numeric i; + save i, outpath, st; + numeric i, st; path outpath; - outpath = for i=0.25 step 0.25 until (length SomeTDPath): + st = 1/NDivide; + outpath = for i=st step st until (length SomeTDPath): ConverterFunc( point i of SomeTDPath ) -- endfor cycle; ( outpath ) @@ -860,9 +905,9 @@ def goodcirclepath(expr CenterPos, AngulMom, Radius ) = begingroup - save cirath, vecx, vecy, ind, goodangulmom, decision; + save cirath, vecx, vecy, ind, goodangulmom, decision, view; numeric ind, decision; - color vecx, vecy, vec[], goodangulmom, view; + color vecx, vecy, goodangulmom, view; path cirath; view = f-CenterPos; decision = cdotprod( view, AngulMom ); @@ -871,19 +916,13 @@ else: goodangulmom = AngulMom; fi; - vecx = ncrossprod( view, goodangulmom ); + vecx = Radius*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; + vecy = Radius*ncrossprod( goodangulmom, vecx ); + cirath = for ind=10 step 10 until 360: + rp( CenterPos + planarrotation(vecx,vecy,ind) )... + endfor cycle; else: cirath = head_on_circle( CenterPos, Radius ); fi; @@ -919,23 +958,15 @@ def head_on_circle(expr Pos, Radius ) = begingroup - save cirath, vecx, vecy, ind, view; + save vecx, vecy, ind, view; numeric ind; - color vecx, vecy, vec[], view; - path cirath; + color vecx, vecy, view; 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 ) + vecx = Radius*N( (-Y(view), X(view), 0) ); + vecy = Radius*ncrossprod( view, vecx ); + ( for ind=10 step 10 until 360: + rp( Pos + planarrotation(vecx,vecy,ind) )... + endfor cycle ) endgroup enddef; @@ -945,7 +976,10 @@ def spatialhalfcircle(expr Center, AngulMom, Radius, ItsTheNearest ) = begingroup - save auxil, auxih; + save va, vb, vc, cc, vd, ux, uy, pa, pb; + save nr, cn, valx, valy, valr, choiceang; + save auxil, auxih, fcirc, returnp; + save choice; color va, vb, vc, cc, vd, ux, uy, pa, pb; numeric nr, cn, valx, valy, valr, choiceang; path auxil, auxih, fcirc, returnp; @@ -957,7 +991,7 @@ vd := cc - Center; % vd := va + vc; nr := conorm( vd ); if Radius >= nr: - returnp := rp( cc ); + returnp := rp( cc ); % this single point will show up in spheroid else: valr := Radius*Radius; valx := valr/nr; @@ -1011,7 +1045,7 @@ 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 auxil, auxih, halfl, halfh, thehole; save auxili, auxihi, rect, theshadow; color va, vb, vc, cc, vd, base; @@ -1126,14 +1160,14 @@ % 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)= + 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; + save basepath, thesubpath, fullpath, finalpath, auxpath, bigcirc; + save themargin, newlen, i, auxt, startt, endt, thelengthofc; + save pa, pb, pc, pd, pe; path basepath, thesubpath, fullpath, finalpath, auxpath; path bigcirc; - numeric themargin, newlen, i, auxt, startt, endt; + numeric themargin, newlen, i, auxt, startt, endt, thelengthofc; pair pa, pb, pc, pd, pe; themargin = 0.02; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER basepath = goodcirclepath( CenterPos, AngulMom, Radius ); @@ -1169,7 +1203,7 @@ endgroup enddef; - def verygoodcone(expr BackDash,CenterPos,AngulMom,Radius,VertexPos)= + def verygoodcone(expr BackDash,CenterPos,AngulMom,Radius,VertexPos) = begingroup save thepath, lenpath, bonevec, sidevec, viewaxe, cipath; save thelengthofc, thesubpath, themargin, basepath; @@ -1219,7 +1253,7 @@ def rigorousfearpath(expr Center, Radius ) = begingroup - save auxil; + save auxil, ux, uy, newcen, nr, valx, valy, valr; color ux, uy, newcen; numeric nr, valx, valy, valr; path auxil; @@ -1235,7 +1269,10 @@ def rigorousfearshadowpath(expr Center, Radius ) = begingroup - save auxil, auxih; + save ux, uy, newcen; + save nr, valx, valy, valr, lenr; + save auxil, auxih, fcirc, returnp; + save dcenter; color ux, uy, newcen; numeric nr, valx, valy, valr, lenr; path auxil, auxih, fcirc, returnp; @@ -1252,7 +1289,7 @@ % It's a globe (without land). - def tropicalglobe( expr NumLats, TheCenter, Radius, AngulMom )= + def tropicalglobe( expr NumLats, TheCenter, Radius, AngulMom ) = begingroup save viewaxe, sinalfa, sinbeta, globaxe, aux, limicos, lc; save stepang, actang, newradius, foc, newcenter, cpath, i; @@ -1315,7 +1352,7 @@ % An elliptical frustum: - def whatisthis(expr CenterPos,OneAxe,OtherAxe,CentersDist,TheFactor)= + def whatisthis(expr CenterPos,OneAxe,OtherAxe,CentersDist,TheFactor) = begingroup save patha, pathb, pathc, centersvec, noption; path patha, pathb, pathc; @@ -1349,17 +1386,122 @@ fi endgroup enddef; - + +% Probably the last algorithm I'm going to write for featpost... + + def spheroidshadow( expr CentrPoi, NorthPoleVec, Ray ) = + begingroup + save a, k, fx, fy, tmpa, tmpb, tmpc, ep, vax, wax, xax, yax, zax, cm, cp; + save bdh, bdm, bdp, bdv, sm, sp, i, cac; + numeric a, k, fx, fy, tmpa, tmpb, tmpc, cm, cp, sm, sp, i; + path ep; + color vax, wax, xax, yax, zax, cac; + pair bdh, bdm, bdp, bdv; + vax = LightSource-CentrPoi; + if cdotprod(NorthPoleVec,vax)<0: + xax = -N(NorthPoleVec); + else: + xax = N(NorthPoleVec); + fi; + a = conorm(NorthPoleVec); + k = a/Ray; + if getangle(xax,vax) > 0.5: %%%%%%%%%%%%%%% DANGER %%% + zax = ncrossprod(xax,vax); + else: + zax = N( ( 0, Z(vax), -Y(vax) ) ); + fi; + yax = ncrossprod(zax,xax); + fx = cdotprod(vax,xax); + fy = cdotprod(vax,yax); + tmpa = Ray*fx/k; + tmpb = fy*(fy ++ (fx/k) +-+ Ray); + tmpc = ((fx/k)**2)+(fy**2); + cm = (tmpa-tmpb)/tmpc; + cp = (tmpa+tmpb)/tmpc; + sm = 1 +-+ cm; + if fx 0.5: %%%%%%%%%%%%%%% DANGER %%% + zax = ncrossprod(xax,vax); + else: + zax = N( (-Y(vax), X(vax), 0) ); + fi; + yax = ncrossprod(zax,xax); + fx = cdotprod(vax,xax); + fy = cdotprod(vax,yax); + tmpa = Ray*fx/k; + tmpb = fy*(fy ++ (fx/k) +-+ Ray); + tmpc = ((fx/k)**2)+(fy**2); + cm = (tmpa-tmpb)/tmpc; + cp = (tmpa+tmpb)/tmpc; + sm = 1 +-+ cm; + if fx ( 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; + pa = patha; + pb = pathb; elseif noption < 0: - draw patha; - if hashole: + pa = pathb; + pb = patha; + fi; + draw pb; + stopair = pa intersectiontimes pb; + hashole = (-1,-1) <> stopair; + if hashole: + testpath = rp(CenterPos+0.5*LenVec)--(point 0 of pa); + isin = (-1,-1) <> testpath intersectiontimes pb; + if not isin: + ta = xpart stopair; + tb = ypart stopair; + stopair := (reverse pa) intersectiontimes (reverse pb); + taf = length pa - xpart stopair; + tbf = length pb - ypart stopair; + margint = 0.01; % DANGER! + draw (subpath (0,ta-margint) of pa)-- + (subpath (tb+margint,tbf-margint) of pb)-- + (subpath (taf+margint,length pa - margint) of pa)-- + cycle; + else: + pathc := buildcycle( pa, pb ); % I don't get it! + % Why doesn't buildcycle work all the time??? See fakehole.mp + % When point 0 of pa is inside pb, builcycle doesn't work!! draw pathc; fi; fi; @@ -1396,8 +1551,7 @@ % 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 ) = + 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; @@ -1498,7 +1652,7 @@ % It's a bit late now but the stage must be set. - def setthestage( expr NumberOfSideSquares, SideSize )= + def setthestage( expr NumberOfSideSquares, SideSize ) = begingroup save i, j, squaresize, squarepath, ca, cb, cc, cd; numeric i, j, squaresize; @@ -1519,7 +1673,7 @@ endgroup enddef; - def setthearena( expr NumberOfDiameterCircles, ArenaDiameter )= + def setthearena( expr NumberOfDiameterCircles, ArenaDiameter ) = begingroup save i, j, circlesize, polar, currpos, phi, cpath; numeric i, j, circlesize, polar, phi; @@ -1718,6 +1872,88 @@ endgroup enddef; +% Take a "quarter" of a "stretched" donut (under construction) + + def quartertorus( expr Tcenter, Tstart, Tfinis, Sray ) = + begingroup + save sideaxe, viewline, circlecenter, circlemoment; + save i, angstep, cuspcond, coofrac, tmoment, tstart, tfinis; + save cpath, opath, ipath, wp, ep, refpair, fpath; + save vstart, vfinis, ostart, ofinis, cstart, cfinis; + color sideaxe, viewline, circlecenter, circlemoment; + color tmoment, vstart, vfinis, tstart, tfinis; + numeric i, angstep, coofrac; + path cpath, opath, ipath, wp, ep, cstart, cfinis, fpath; + pair outerp[], innerp[], refpair; + boolean cuspcond, ostart, ofinis; + angstep = 6; + viewline = f-Tcenter; +% if cdotprod( viewline, Tstart ) < cdotprod( viewline, Tfinis ): +% tstart = Tfinis; +% tfinis = Tstart; +% else: + tstart = Tstart; + tfinis = Tfinis; +% fi; + tmoment = ncrossprod( tstart, tfinis ); + vstart = ncrossprod( tstart, tmoment ); + vfinis = ncrossprod( tmoment, tfinis ); + ostart = cdotprod( viewline-tstart, vstart ) > 0; + ofinis = cdotprod( viewline-tfinis, vfinis ) > 0; + cstart = spatialhalfcircle(tstart+Tcenter,vstart,Sray,true); + cfinis = spatialhalfcircle(tfinis+Tcenter,vfinis,Sray,true); + if cdotprod( viewline, tmoment ) < 0: + tmoment := -tmoment; + fi; + refpair = unitvector( rp(Tcenter+tmoment)-rp(Tcenter) ); + sideaxe = ncrossprod( tmoment, viewline ); + coofrac = cdotprod( viewline, tmoment )/Sray; + + for i=0 step angstep until 90: + circlecenter:= tstart*cosd(i)+tfinis*sind(i); + circlemoment:= ccrossprod(circlecenter,tmoment); + cpath:=spatialhalfcircle(circlecenter+Tcenter,circlemoment,Sray,true); + if cdotprod( sideaxe, circlecenter ) < 0: + outerp[i/angstep]=point 0 of cpath; + innerp[i/angstep]=point (length cpath) of cpath; + else: + innerp[i/angstep]=point 0 of cpath; + outerp[i/angstep]=point (length cpath) of cpath; + fi; + endfor; + opath = outerp0 for i=angstep step angstep until 90: + ..outerp[i/angstep] endfor; + ipath = innerp0 for i=angstep step angstep until 90: + ..innerp[i/angstep] endfor; + fpath = ipath---cfinis---reverse opath---cstart---cycle; + unfill fpath; + draw fpath; +% draw cstart---ipath---cfinis; + if ostart: + cpath := rigorouscircle( Tcenter+tstart, vstart, Sray ); + unfill cpath; + draw cpath; + fi; + if ofinis: + cpath := rigorouscircle( Tcenter+tfinis, vfinis, Sray ); + unfill cpath; + draw cpath; + fi; + +% i := 0; +% cuspcond = false; +% forever: +% exitif i > (90/angstep)-2; +% i := incr( i ); +% cuspcond := +% refpair dotprod outerp[i+1] < +% refpair dotprod outerp[i]; +% exitif cuspcond; +% endfor; +% undraw subpath (i,90/angstep) of opath; + endgroup + enddef; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Non-standard objects: @@ -1747,7 +1983,7 @@ endgroup enddef; - def simplecar(expr RefP, AngCol, LenCol, FronWheelCol, RearWheelCol )= + 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; @@ -1835,8 +2071,9 @@ % 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. +% This is for solving first order differential equations - def fieldlinestep( expr Spos, Step )( text VecFunc )= + def fieldlinestep( expr Spos, Step )( text VecFunc ) = begingroup save kone, ktwo; color kone, ktwo; @@ -1846,7 +2083,7 @@ endgroup enddef; - def fieldlinepath( expr Numb, Spos, Step )( text VecFunc )= + def fieldlinepath( expr Numb, Spos, Step )( text VecFunc ) = begingroup save ind, flpath, prevpos, thispos; numeric ind; @@ -1864,8 +2101,9 @@ enddef; % Another point is that I want to draw trajectories in space. +% This is for solving second order differential equations - def trajectorypath( expr Numb, Spos, Svel, Step )( text VecFunc )= + def trajectorypath( expr Numb, Spos, Svel, Step )( text VecFunc ) = begingroup save ind, flpath, prevpos, thispos, prevvel, thisvel; save rone, rtwo, vone, vtwo; @@ -1895,10 +2133,47 @@ endgroup enddef; +% Another point is that I want to draw trajectories in space and +% dependant on velocity: VecFunc( position, velocity ). +% This time is fourth-order Runge-Kutta. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHANGES PrintStep!!!! + + def dragtrajectorypath( expr 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 ); + ind = 1; + forever: + vone := Step*VecFunc( prevpos , prevvel ); + rone := Step*prevvel; + vtwo := Step*VecFunc( prevpos+0.5*rone, prevvel+0.5*vone ); + rtwo := Step*( prevvel+0.5*vone ); + vthr := Step*VecFunc( prevpos+0.5*rtwo, prevvel+0.5*vtwo ); + rthr := Step*( prevvel+0.5*vtwo ); + vfou := Step*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; + exitif Z( thispos ) < -0.0001; %%%%%%%%%% EDIT! + prevpos := thispos; + prevvel := thisvel; + flpath := flpath--rp( thispos ); + endfor; + PrintStep := Y(thispos); + ( flpath ) + endgroup + enddef; + % And now i stop. - def magnetictrajectorypath( expr Numb, Spos, Svel, Step ) - ( text VecFunc )= + 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; @@ -1940,36 +2215,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 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 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 @@ -2294,7 +2569,7 @@ % FillDraw a face - def face_invisible( expr Facen )( text LineAtribs )= + def face_invisible( expr Facen )( text LineAtribs ) = begingroup save ghost; path ghost; @@ -2363,7 +2638,7 @@ color somepo; numeric fract; (b-a)*fract = somepo-a; - cdotprod(swapc,somepo)=cdotprod(swapc,f); + cdotprod(swapc,somepo) =cdotprod(swapc,f); if (fract>0) and (fract<1): crosspoin[counter] := somepo; counter := incr(counter); @@ -2442,7 +2717,7 @@ % Draw only the faces, rigorously projecting the edges. - def faceraytrace(expr Press, Col)= + def faceraytrace(expr Press, Col) = begingroup save i, j, a, b; numeric i, j; @@ -2463,7 +2738,7 @@ % Fast test for your three-dimensional object - def draw_all_test( expr AlsoDrawLines )= + def draw_all_test( expr AlsoDrawLines ) = begingroup save i, j, a, b; numeric i, j; @@ -2507,7 +2782,7 @@ % Don't use SphericalDistortion here. - def fill_faces( text LineAtribs )= + def fill_faces( text LineAtribs ) = begingroup save i; numeric i; @@ -2594,8 +2869,7 @@ endgroup enddef; - def generateonebiax(expr Lin, Phi, Theta, Long, - SndDirAngl, Base, Currpos ) = + def generateonebiax(expr Lin, Phi, Theta, Long, SndDirAngl, Base, Currpos ) = begingroup save basevec, longvec, u, v; color basevec, longvec, u, v; @@ -2615,7 +2889,7 @@ endgroup enddef; - def director_invisible( expr SortEmAll, ThickenFactor, CyclicLines )= + def director_invisible( expr SortEmAll, ThickenFactor, CyclicLines ) = begingroup save i, j, k, farone, thisfar; save outerr, innerr, direc, ounum; @@ -2681,6 +2955,123 @@ endgroup enddef; +% Now two routines to draw "bananas", well, sort of... +% Initially coded by Pedro J. Sebastião + + def circularsheet( expr CenterP, Rad, VecX, VecY, StartA, FinisA, Width ) = + begingroup + save vecz, ind; + color vecz; + numeric ind; + vecz = ncrossprod( VecX, VecY ); + ( rp( CenterP + Width*vecz + Rad*planarrotation( VecX, VecY, StartA ) ) + for ind=StartA+1 upto FinisA: + --rp( CenterP + Width*vecz + Rad*planarrotation( VecX, VecY, ind )) + endfor + for ind=FinisA downto StartA: + --rp( CenterP + Rad*planarrotation( VecX, VecY, ind ) ) + endfor + --cycle ) + endgroup + enddef; + + def banana( expr CenterPos, Radius, AngleColor, Wid, Amp ) = + begingroup + save ind, sinbeta, cosbeta, aux, delta, angfvbx, angpos, angneg, au; + save bx, by, bz, fv, beta, gamma, alfa; + save outpath, outpathb, outpathc, visneg, vispos; + numeric ind, sinbeta, cosbeta, aux, delta, angfvbx, angpos, angneg, au; + numeric beta, gamma, alfa; + color bx, by, bz, fv; + path outpath, outpathb, outpathc; + boolean visneg, vispos; + + alfa = X(AngleColor); + beta = Y(AngleColor); + gamma= Z(AngleColor); + bx = eulerrotation( alfa, beta, gamma, red ); + by = eulerrotation( alfa, beta, gamma, green ); + bz = eulerrotation( alfa, beta, gamma, blue ); + + au = cdotprod( f-CenterPos, by ); + if 0 > au: + by := -by; + fi; + fv = cdotprod( f, bx )*bx + cdotprod( f, by )*by; + aux = conorm( fv-CenterPos ); + if aux > Radius: + cosbeta = Radius/aux; + sinbeta = ( aux +-+ Radius )/aux; + delta = angle( cosbeta, sinbeta ); + angfvbx = getangle( fv - CenterPos, bx ); + if 0 > cdotprod( fv - CenterPos, by ): + angfvbx := -angfvbx; + fi; + angpos = delta + angfvbx; + angneg = angfvbx - delta; + if ( angneg > -Amp ) and ( Amp > angneg ): + visneg = true; + else: + visneg = false; + fi; + if ( angpos > -Amp ) and ( Amp > angpos ): + vispos = true; + else: + vispos = false; + fi; + if visneg and not vispos: + outpath = circularsheet(CenterPos,Radius,bx,by,-Amp,angneg,Wid); + unfill outpath; + draw outpath; + outpathb = circularsheet(CenterPos,Radius,bx,by,angneg,Amp,Wid); + unfill outpathb; + draw outpathb + fi; + if vispos and not visneg: + outpath = circularsheet(CenterPos,Radius,bx,by,-Amp,angpos,Wid); + unfill outpath; + draw outpath; + outpathb = circularsheet(CenterPos,Radius,bx,by,angpos,Amp,Wid); + unfill outpathb; + draw outpathb + fi; + if (not vispos) and (not visneg): + outpath = circularsheet(CenterPos,Radius,bx,by,-Amp,Amp,Wid); + unfill outpath; + draw outpath; + fi; + if vispos and visneg: + if 0 > cdotprod( f-CenterPos, bx ): + outpath=circularsheet(CenterPos,Radius,bx,by,angneg,angpos,Wid); + unfill outpath; + draw outpath; + outpathb = circularsheet(CenterPos,Radius,bx,by,-Amp,angneg,Wid); + unfill outpathb; + draw outpathb; + outpathc = circularsheet(CenterPos,Radius,bx,by,angpos,Amp,Wid); + unfill outpathc; + draw outpathc; + else: + outpathb = circularsheet(CenterPos,Radius,bx,by,-Amp,angneg,Wid); + unfill outpathb; + draw outpathb; + outpathc = circularsheet(CenterPos,Radius,bx,by,angpos,Amp,Wid); + unfill outpathc; + draw outpathc; + outpath=circularsheet(CenterPos,Radius,bx,by,angneg,angpos,Wid); + unfill outpath; + draw outpath; + fi; + fi; + else: + outpath = circularsheet( CenterPos, Radius, bx, by, -Amp, Amp, Wid ); + unfill outpath; + draw outpath; + fi; + draw rp(CenterPos+Radius*bx)--rp(CenterPos+Radius*bx+Wid*bz); + endgroup + enddef; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Plotting: @@ -2817,7 +3208,7 @@ % ColAtribone the color used for drawing; % Colordensity depends on distance of the face from f - def face_drawfill( expr Facen, dmin_, dmax_ ,ColAtrib, ColAtribone )= + def face_drawfill( expr Facen, dmin_, dmax_ ,ColAtrib, ColAtribone ) = begingroup save j, ptmp, colfac_, coltmp_; path ghost; @@ -2854,7 +3245,7 @@ % ColAtrib=color for filling faces % ColAtribone=color for drawing edges - def draw_invisible( expr Option, DoJS, ColAtrib, ColAtribone )= + 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_; @@ -2945,11 +3336,7 @@ % 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 ) = + 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; @@ -3224,13 +3611,70 @@ ( trypoi ) endgroup; enddef; + + 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 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; + 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 := La +-+ (xa-x); + 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: - aux = Index; + selectedpoint = Bs; 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; + 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: + 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 ) + 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; + +% 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. + + def lasermachine( expr DefinedPath, Beam, CosLimit ) = + begingroup + save patlen, j; + save apoi, bpoi, cpoi, dpoi, epoi; + save apat, bpat, cpat, dpat, a, b, c; + save anew, bnew, cnew; + save aang, bang, cang, dang; + save newp, pairvector, cou, val; + numeric patlen, j; + pair apoi, bpoi, cpoi, dpoi, epoi; + pair apat, bpat, cpat, dpat, a, b, c; + pair anew, bnew, cnew, pairvector[]; + numeric aang, bang, cang, dang, cou, val; + path newp; + patlen = length DefinedPath; + cou = 0; + for j=0 upto patlen-1: + apoi := precontrol j of DefinedPath; + bpoi := point j of DefinedPath; + cpoi := postcontrol j of DefinedPath; + dpoi := precontrol j+1 of DefinedPath; + epoi := point j+1 of DefinedPath; + apat := apoi-bpoi; + bpat := bpoi-cpoi; + cpat := cpoi-dpoi; + dpat := dpoi-epoi; + aang := angle( apat ); + bang := angle( bpat ); + cang := angle( cpat ); + dang := angle( dpat ); + val := cosd( 0.5*(aang-bang) ); + bnew := cpoi+Beam*dir( 90+0.5*(bang+cang) )/cosd( 0.5*(bang-cang) ); + cnew := dpoi+Beam*dir( 90+0.5*(cang+dang) )/cosd( 0.5*(cang-dang) ); + if ( val > 0 ) and + ( val < CosLimit ) and + ( Beam*sind( 0.5*(aang-bang) ) < 0 ): + + a := bpoi+Beam*dir(90+aang); + b := a-Beam*dir(aang); + anew := bpoi+Beam*dir(90+bang); + c := anew+Beam*dir(bang); + pairvector[3*cou] = a; + pairvector[3*cou+1] = b; + pairvector[3*cou+2] = c; + cou := incr(cou); + pairvector[3*cou] = anew; + pairvector[3*cou+1] = bnew; + pairvector[3*cou+2] = cnew; + cou := incr(cou); + else: + anew := bpoi+Beam*dir( 90+0.5*(aang+bang) )/val; + pairvector[3*cou] = anew; + pairvector[3*cou+1] = bnew; + pairvector[3*cou+2] = cnew; + cou := incr(cou); + fi; endfor; - cutp = auxp intersectionpoint ( origin--( 0, length Path ) ); - ( point ( ypart cutp ) of Path ) - endgroup - enddef; + newp = for j=0 upto cou-1: + pairvector[3*j]..controls pairvector[3*j+1] and pairvector[3*j+2].. + endfor cycle; + ( newp ) + endgroup + enddef; + +% Move the starting point of a cyclic path along that path + def startahead( expr DefinedPath, JumpTime ) = + begingroup + save patlen, j; + save apoi, bpoi, cpoi; + save newp; + numeric patlen, j; + pair apoi, bpoi, cpoi; + path newp; + patlen = length DefinedPath; + newp = for j=0 upto patlen-1: + hide( + apoi := point JumpTime+j of DefinedPath; + bpoi := postcontrol JumpTime+j of DefinedPath; + cpoi := precontrol JumpTime+j+1 of DefinedPath; + ) + apoi..controls bpoi and cpoi.. + endfor + cycle; + ( newp ) + endgroup + enddef; + +% In order to use a "lasermachine" one needs a single full outline. +% One may have two somewhat concentric cyclic paths intersecting in +% several points. +% The next routine may help but first the paths must be adapted with +% "startahead" and/or "reverse" so that they both rotate in the same +% direction and they start on consecutive "lobes" (hard to explain). +% Now pay attention: given the direction of rotation (clockwise or +% counter-clockwise) the SecondPath must start BEFORE the FirstPath. +% And another problem: there must be at least four intersection points. +% Very nasty routine. All because of finispoi... + + def crossingline( expr FirstPath, SecondPath, TimeTolerance ) = + begingroup + save m; + save its, finispoi; + save increm, fo, ma, tmpp, mastarter; + numeric m; + pair its, finispoi; + path increm, fo, ma, tmpp, mastarter; + m = TimeTolerance; + fo = FirstPath; + ma = SecondPath; + its := fo intersectiontimes ma; + increm := subpath (m, (xpart its) - m ) of fo; + mastarter = subpath (m, (ypart its) - m ) of ma; + finispoi = reverse ma intersectionpoint reverse fo; + + forever: + + fo := subpath ( (xpart its)+m, length fo ) of fo; + ma := subpath ( (ypart its)+m, length ma ) of ma; + its := ma intersectiontimes fo; + tmpp := subpath ( m, (xpart its)-m ) of ma; + increm := increm...tmpp; + + fo := subpath ( (ypart its)+m, length fo ) of fo; + ma := subpath ( (xpart its)+m, length ma ) of ma; + its := fo intersectiontimes ma; + tmpp := subpath ( m, (xpart its)-m ) of fo; + increm := increm...tmpp; + + exitif abs( point (xpart its) of fo - finispoi ) < m; + endfor; + + fo := subpath ( (xpart its)+m, length fo ) of fo; + ma := (subpath ( (ypart its)+m, (length ma)-m ) of ma)...mastarter; + its := ma intersectiontimes fo; + tmpp := subpath ( m, (xpart its)-m ) of ma; + increm := increm...tmpp; + + tmpp := subpath ( (ypart its)+m, (length fo)-m ) of fo; + ( increm...tmpp...cycle ) + 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 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; + +% Mark bidimensional angles (contributed by Palle Jørgensen +% to the metapost mailing list) + + vardef archangle@#( expr _p, _q, _s, archwidth ) text _t = + begingroup; + save _a, _b, _w, _arch, _halfangle, _label_origin; + ( _a, _b ) = _p intersectiontimes _q; + pair _w; + _w = whatever[ + point _a of _p + + archwidth * unitvector direction _a of _p, + point _a of _p + + archwidth * unitvector direction _a of _p + + (ypart.direction _a of _p, -xpart.direction _a of _p) + ]; + _w = whatever[point _b of _q, point _b of _q + direction _b of _q]; + path _arch; + _arch = point _a of _p + + archwidth * unitvector direction _a of _p{ + (if direction _a of _p dotprod unitvector direction _b of _q > 0: + 1 + else: + -1 + fi) * + ( _w - (point _a of _p + archwidth * unitvector direction _a of _p) ) + }..point _b of _q + + archwidth * unitvector direction _b of _q; + draw _arch _t; + path _halfangle; + _halfangle = point _a of _p - 2*archwidth* + unitvector( direction _a of _p + direction _b of _q )--point _a of _p + + 2*archwidth*unitvector( direction _a of _p + direction _b of _q ); + pair _label_origin; + _label_origin = _halfangle intersectionpoint _arch; + label@#( _s, _label_origin ) _t; + endgroup; +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 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% Part VI (strictly two-dimensional and related to planifications): +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Default Values %%%%%%%%%%%% + pen grossa, usual, fina; + grossa= pencircle scaled 0.9mm; + usual = pencircle scaled 0.6mm; + fina = pencircle scaled 0.3mm; + numeric u, arestaref, defaultaresta; + u = 0.5mm; + defaultaresta = 30mm; + arestaref = defaultaresta; + picture pequeno, grande, nada; + linecap := squared; %%%%%%%%%%%%%%%%%%%%% + pequeno = dashpattern( off 3u on 1u ); + grande = dashpattern( off 8u on 6u ); + nada = dashpattern( on 1pt ); + linecap := rounded; %%%%%%%%%%%%%%%%%%%%% + pair urCorner,ulCorner,llCorner,lrCorner,PageCenter; + urCorner := (8.13in,11.533in); + llCorner := (0.118in,35.5bp); + ulCorner := (xpart llCorner,ypart urCorner); + lrCorner := (xpart urCorner,ypart llCorner); + PageCenter := (urCorner+ulCorner+llCorner+lrCorner)/4; + path thebigframe; + thebigframe = llCorner--lrCorner--urCorner--ulCorner--cycle; + boolean EPSmode, Athreemode; + EPSmode = false; + Athreemode = false; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + def centerof( text t) = + begingroup + save momentum, counter; + numeric counter; + pair momentum; + counter = 0; + momentum = (0,0); + forsuffixes $=t: + counter := incr(counter); + momentum := momentum + z$; + endfor; + ( momentum/counter ) + endgroup enddef; - -% EOF + + def truefyepsmode = + begingroup + EPSmode := true; + arestaref := arestaref/2; + endgroup + enddef; + + def makearestarefathree = + begingroup + if not EPSmode: + Athreemode := true; + arestaref := arestaref*1.414; + fi; + endgroup + enddef; + + def terminar( expr ScaleFactor, RotAngle ) = + begingroup + picture inicial,final; + inicial = currentpicture; + currentpicture := nullpicture; + final = inicial scaled ScaleFactor; + if EPSmode: + final := final scaled 2 + fi; + final := final rotated RotAngle; + if Athreemode: + draw final shifted (PageCenter*1.414); + arestaref := defaultaresta; + Athreemode := false; + else: + draw final shifted PageCenter; + if EPSmode: + pickup grossa; + draw thebigframe; + inicial := currentpicture; + currentpicture := nullpicture; + final := inicial rotated -90; + draw final + fi; + fi; + endgroup + enddef; + + def desvia( text zes ) = + begingroup + picture inicial,final; + inicial = currentpicture; + currentpicture := nullpicture; + final = inicial shifted -centerof( zes ); + draw final + endgroup + enddef; + + def desvec( expr vect ) = + begingroup + picture inicial,final; + inicial = currentpicture; + currentpicture := nullpicture; + final = inicial shifted -vect; + draw final + endgroup + enddef; + +% EOF \ No newline at end of file -- cgit v1.2.3