% featpost3Dplus2D.mp % 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) 2014 % 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 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, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. message "FeatPost-0.8.8"; warningcheck := 0; background := 0.987white; defaultscale := 0.75; defaultfont := "cmss17"; % This is used by cartaxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global Variables %%%%%%%%%%% boolean ParallelProj, SphericalDistortion, FCD[], ShadowOn; boolean OverRidePolyhedricColor, MalcomX; numeric Nobjects, RefDist[], HoriZon, RopeColorSeq[], PhotoMarks; numeric Spread, PrintStep, PageHeight, PageWidth, ActuC, Shifts; numeric NL, npl[], NF, npf[], FC[], MaxFearLimit, TableColors; numeric TDAtiplen, TDAhalftipbase, TDAhalfthick, RopeColors, NCL; pair OriginProjPagePos, ShiftV, PhotoPair[]; path VGAborder, CLPath[]; color f, viewcentr, V[], L[]p[], F[]p[], TableC[]; color HigColor, SubColor, LightSource, PhotoPoint[]; string ostr[]; pen BackPen, ForePen; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Default Values %%%%%%%%%%%%%%% f := (3,5,4); % This f is the point of view in 3D viewcentr := black; % This is the aim of the view Spread := 140; % Magnification Shifts := 105.00mm; ShiftV := Shifts*(1,1); % Central coordinates on paper 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 % 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 (182.05,382.05)--cycle; % get few extra pixels PrintStep := 5; % Coarseness, in resolvec PageHeight := 9in; PageWidth := 6in; % And this is used by produce_auto_scale MaxFearLimit := 17; % Valid Maximum Distance from Origin HigColor := 0.85white; % These two colors are used in 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 TableC3 := ( 1.0, 0.7, 0.0 ); % orange %% TableC4 := 0.85green; % pale green %% TableC5 := 0.90*(red+blue); % magenta %% TableC6 := 0.85*(green+blue); % cyan %% TableC7 := 0.85*(red+green); % yellow %% TableColors := 7; ActuC := 5; RopeColorSeq0 := 3; % RopeColorSeq1 := 3; % RopeColorSeq2 := 1; % RopeColorSeq3 := 3; % ropepattern RopeColorSeq4 := 7; % RopeColorSeq5 := 5; % % RopeColors := 5; % Nobjects := 0; % getready and doitnow TDAtiplen := 0.05; % tdarrow and tdcircarrow TDAhalftipbase := 0.02; % Three-Dimensional (Circular) TDAhalfthick := 0.01; % Arrow NCL := 0; % closedline ForePen := pencircle scaled 15pt; BackPen := pencircle scaled 9pt; MalcomX := false; %%% The variables PhotoMarks, PhotoPair[], PhotoPoint[] %%% and CLPath[] have NO default values. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Part I: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Very basic: % Colors have three or four coordinates. Get one. def X(expr A) = if color A: redpart A elseif MalcomX: blackpart A else: cyanpart A fi enddef; def Y(expr A) = if color A: greenpart A else: magentapart A fi enddef; def Z(expr 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) ) 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: % Calculate the unit vector of a vector (or a point) def N(expr A) = begingroup save M, exitcolor; numeric M; color exitcolor; M = conorm( A ); if M > 0: exitcolor = ( X(A)/M, Y(A)/M, Z(A)/M ); else: exitcolor := black; fi; ( exitcolor ) endgroup enddef; def cdotprod(expr A, B) = ( X(A)*X(B) + Y(A)*Y(B) + Z(A)*Z(B) ) enddef; def ccrossprod(expr A, B) = ( Y(A)*Z(B) - Z(A)*Y(B), Z(A)*X(B) - X(A)*Z(B), X(A)*Y(B) - Y(A)*X(B) ) enddef; % The dotproduct of two normalized vectors is the cosine of the angle % they form. def ndotprod(expr A, B) = begingroup save a, b; color a, b; a = N(A); b = N(B); ( ( X(a)*X(b) + Y(a)*Y(b) + Z(a)*Z(b) ) ) endgroup enddef; % The normalized crossproduct of two vectors. def ncrossprod(expr A, B) = N( ccrossprod( A, B ) ) enddef; % Haahaa! Trigonometry. def getangle(expr A, B) = begingroup save coss, sine; numeric coss, sine; coss := cdotprod( A, B ); sine := conorm( ccrossprod( A, B ) ); ( angle( coss, sine ) ) endgroup enddef; % Something I need for spatialhalfsfear. def getcossine( expr Center, Radius ) = begingroup save a, b; numeric a, b; a = conorm( f - Center ); b = Radius/a; if abs(b) >= 1: show "The point of view f is too close (getcossine)."; b := 2; % DANGER! fi; ( b ) endgroup enddef; % The following routine is used by circularsheet and may be used to % rotate vectors elliptically. Also used by tdcircarrow. 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). def eulerrotation( expr AngA, AngB, AngC, Vec ) = begingroup save auxx, auxy, veca, vecb, vecc; color auxx, auxy, veca, vecb, vecc; veca = ( cosd(AngA)*cosd(AngB), sind(AngA)*cosd(AngB), sind(AngB) ); auxx = ( cosd(AngA+90), sind(AngA+90), 0 ); auxy = ccrossprod( veca, auxx ); vecb = cosd(AngC)*auxx + sind(AngC)*auxy; vecc = cosd(AngC+90)*auxx + sind(AngC+90)*auxy; ( X(Vec)*veca + Y(Vec)*vecb + Z(Vec)*vecc ) 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 ) = begingroup save va, vb, vc; color va, vb, vc; va = Point - PointPerpA; vb = Point - PointPerpB; vc = PointPerpB - PointPerpA; ( cdotprod(va,vc)*cdotprod(vb,vc) <= 0 ) 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: % Projection Size. Meant for objects with size one. % Used by signalvertex. def ps(expr A, Thicken_Factor) = Thicken_Factor/conorm(A-f)/3 enddef; % Rigorous Projection of a Point. Draws a circle with % a diameter inversely proportional to the distance of % that Point from the point of view. def signalvertex(expr A, TF, Col) = draw rp(A) withcolor Col withpen pencircle scaled (Spread*ps(A,TF)) enddef; def signalshadowvertex(expr A, TF, Col) = begingroup save auxc, auxn; color auxc; numeric auxn; auxc := cb(A); auxn := TF*conorm(f-auxc)/conorm(LightSource-A); signalvertex( auxc, auxn, Col ) endgroup enddef; % Get the vector that projects onto the resolution def resolvec(expr A, B) = begingroup save sizel, returnvec; numeric sizel; color returnvec; sizel = abs( rp(A) - rp(B) ); if sizel > 0: returnvec = PrintStep*(B-A)/sizel; else: returnvec = 0.3*(B-A); fi; ( returnvec ) endgroup enddef; % Movies need a constant frame def produce_vga_border = begingroup draw VGAborder withcolor background withpen pencircle scaled 0; clip currentpicture to VGAborder 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; numeric pwidth, pheight; storeall = currentpicture shifted -(center currentpicture); currentpicture := nullpicture; pwidth = xpart ((lrcorner storeall)-(llcorner storeall)); pheight = ypart ((urcorner storeall)-(lrcorner storeall)); if PageHeight/PageWidth < pheight/pwidth: scaleall = storeall scaled (PageHeight/pheight); else: scaleall = storeall scaled (PageWidth/pwidth); fi; draw scaleall shifted OriginProjPagePos endgroup enddef; % The following two procedures are useful for getready. 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-viewcentr, unless SphericalDistortion is true. % f-viewcentr must not be (anti-)parallel to zz. def rp(expr R) = begingroup save projpoi; save v, u; save verti, horiz, eta, squarf, radio, ang, lenpl; pair projpoi; color v, u; numeric verti, horiz, eta, squarf, radio, ang, lenpl; v = N( (-Y(f-viewcentr), X(f-viewcentr), 0) ); u = ncrossprod( f-viewcentr, v ); horiz = cdotprod( R-viewcentr, v ); verti = cdotprod( R-viewcentr, u ); if SphericalDistortion: if ( horiz <> 0 ) or ( verti <> 0 ): lenpl = ( horiz ++ verti )*20; %%%%%%%%%%%%%%% DANGER ang = getangle( f-R, f-viewcentr ); horiz := ang*horiz/lenpl; verti := ang*verti/lenpl; projpoi = (horiz,verti); else: projpoi = origin; fi; else: if ParallelProj: eta = 1; else: squarf = cdotprod( f-viewcentr, f-viewcentr ); radio = cdotprod( R-viewcentr, f-viewcentr ); eta = 1 - radio/squarf; if abs((horiz,verti)) > MaxFearLimit*eta: eta := abs((horiz,verti))/MaxFearLimit; fi; fi; projpoi = (horiz,verti)/eta; fi; ( projpoi*Spread + ShiftV ) endgroup enddef; % Much improved rigorous pseudo-projection algorithm that follows % an idea from Cristian Barbarosie. % This makes shadows caused by a light source point. def cb(expr R) = begingroup save ve, ho; numeric ve, ho; LightSource-ho*red-ve*green-HoriZon*blue=whatever*(LightSource-R); ( ho*red + ve*green + HoriZon*blue ) endgroup enddef; % And this just projects points rigorously on some generic plane using % LightSource as the point of convergence (focus). def projectpoint(expr ViewCentr, R) = begingroup save verti, horiz; save v, u, lray; numeric verti, horiz; color v, u, lray; lray = LightSource-ViewCentr; v = N( (-Y(lray), X(lray), 0) ); u = ncrossprod( lray, v ); lray - horiz*v - verti*u = whatever*( LightSource - R ); ( horiz*v + verti*u + ViewCentr ) endgroup enddef; % And this is the way to calculate the intersection of some line with some % plan. def lineintersectplan( expr LinePoi, LineDir, PlanPoi, PlanDir ) = begingroup save incognitus; color incognitus; cdotprod( incognitus-PlanPoi, PlanDir ) = 0; whatever*LineDir + LinePoi = incognitus; ( incognitus ) endgroup enddef; % Vanishing point. def vp( expr D ) = begingroup ( rp( lineintersectplan( f, D, viewcentr, f-viewcentr) ) ) endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Basic Functions: % Get the 2D path of a straight line in beetween 3D points A and B. % This would add rigor to rigorousdisc, if one would introduce the % the concept of three-dimensional path. That is not possible now. % Also this is only interesting when using SphericalDistortion:=true def pathofstraightline( expr A, B ) = begingroup save k, i, mark, stepVec, returnp, pos; numeric k, i; color mark, stepVec; path returnp; pair pos[]; stepVec = resolvec(A,B); pos0 = rp( A ); k = 1; forever: mark := A+(k*stepVec); exitif cdotprod(B-mark,stepVec) <= 0; pos[k] = rp( mark ); k := incr(k); endfor; pos[k] = rp(B); returnp = pos0 for i=1 upto k: ..pos[i] endfor; ( returnp ) endgroup enddef; def drawsegment( expr A, B ) = begingroup if SphericalDistortion: draw pathofstraightline( A, B ); else: draw rp(A)--rp(B); fi endgroup enddef; % Cartesian axes with prescribed lengths. def cartaxes(expr axex, axey, axez) = begingroup save orig, axxc, ayyc, azzc; color orig, axxc, ayyc, azzc; orig = (0,0,0); axxc = (axex,0,0); ayyc = (0,axey,0); azzc = (0,0,axez); drawarrow rp(orig)..rp(axxc); drawarrow rp(orig)..rp(ayyc); drawarrow rp(orig)..rp(azzc); label.bot( "x" ,rp(axxc)); %%%%%%%%%%%%%%%%%%%%%%%%% label.bot( "y" ,rp(ayyc)); %% Some Labels... %% label.lft( "z" ,rp(azzc)); %%%%%%%%%%%%%%%%%%%%%%%%% 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 % define the lines. The arch is at a distance W from Or. The % label is S and the position is RelPos (rt,urt,top,ulft,lft, % llft,bot,lrt). But arches must be smaller than 180 degrees. def angline(expr A, B, Or, W, S)(suffix RelPos) = begingroup save G, Dna, Dnb, al; numeric G; color Dna, Dnb; path al; G = conorm( W*( N(A-Or) - N(B-Or) ) )/2.5; %%%%%%% BIG DANGER! Dna = ncrossprod(ncrossprod(A-Or,B-Or),A-Or); Dnb = ncrossprod(ncrossprod(B-Or,A-Or),B-Or); al = rp(W*N(A-Or)+Or).. controls rp(W*N(A-Or)+Or+G*Dna) and rp(W*N(B-Or)+Or+G*Dnb).. rp(W*N(B-Or)+Or); draw al; label.RelPos( S, point 0.5*length al of al ) endgroup enddef; % As i don't know how to declare variables of type suffix, % i provide a way to avoid the problem. This time RelPos may % be 0,1,2,3,4,6,7 or anything else. def anglinen(expr A, B, Or, W, S, RelPos) = begingroup save G, Dna, Dnb, al, middlarc; numeric G; color Dna, Dnb; path al; pair middlarc; G = conorm( W*( N(A-Or) - N(B-Or) ) )/3; Dna = ncrossprod(ncrossprod(A-Or,B-Or),A-Or); Dnb = ncrossprod(ncrossprod(B-Or,A-Or),B-Or); al = rp(W*N(A-Or)+Or).. controls rp(W*N(A-Or)+Or+G*Dna) and rp(W*N(B-Or)+Or+G*Dnb).. rp(W*N(B-Or)+Or); draw al; middlarc = point 0.5*length al of al; if RelPos = 0: label.rt( S, middlarc ); elseif RelPos =1: label.urt( S, middlarc ); elseif RelPos =2: label.top( S, middlarc ); elseif RelPos =3: label.ulft( S, middlarc ); elseif RelPos =4: label.lft( S, middlarc ); elseif RelPos =5: label.llft( S, middlarc ); elseif RelPos =6: label.bot( S, middlarc ); elseif RelPos =7: label.lrt( S, middlarc ); else: label( S, middlarc ); fi endgroup enddef; % As a bigger avoidance, replace the arch by a paralellogram. def squareangline(expr A, B, Or, W) = begingroup save sal; path sal; sal = rp(Or)--rp(W*N(A-Or)+Or)-- rp(W*(N(B-Or)+N(A-Or))+Or)--rp(W*N(B-Or)+Or)--cycle; draw sal endgroup enddef; % Just as we are here we can draw circles. (color,color,numeric) def rigorouscircle( expr CenterPos, AngulMom, Radius ) = begingroup save ind, G, Dn, Dna, Dnb, al, vec; numeric ind, G; color vec[], Dn, Dna, Dnb; path al; vec1 = ncrossprod( CenterPos-f, AngulMom); for ind=2 step 2 until 8: vec[ind+1] = ncrossprod( vec[ind-1], AngulMom ); vec[ind] = N( vec[ind-1] + vec[ind+1] ); endfor; G = conorm( Radius*( vec1 - vec2 ) )/3; al = rp(Radius*vec1+CenterPos) for ind=2 upto 8: hide( 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) ..rp(Radius*vec[ind] +CenterPos) endfor ...cycle; ( al ) endgroup enddef; % 3D straight arrow. def tdarrow(expr FromPos, ToTip ) = begingroup save basevec, longvec, a, b, c, d, e, g, h, len, p; color basevec, longvec, a, b, c, d, e, g, h; numeric len; path p; len = conorm( ToTip - FromPos ); longvec := N( ToTip - FromPos ); basevec := ncrossprod( FromPos-f, longvec ); if len <= TDAtiplen: b = basevec*TDAhalftipbase*len/TDAtiplen; c = FromPos+b; e = FromPos-b; p = rp(ToTip)--rp(c)--rp(e)--cycle; else: d = ToTip-longvec*TDAtiplen; a = FromPos+basevec*TDAhalfthick; h = FromPos-basevec*TDAhalfthick; b = d+basevec*TDAhalfthick; g = d-basevec*TDAhalfthick; c = d+basevec*TDAhalftipbase; e = d-basevec*TDAhalftipbase; p = rp(a)--rp(b)--rp(c)--rp(ToTip)--rp(e)--rp(g)--rp(h)--cycle; fi; unfill p; draw p 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) = begingroup save i, j, k; numeric i, j, k; k = ThickenFactor*EmptyFrac; if ShadowOn: for i = 0 upto theN: signalshadowvertex( LinFunc(i/theN), ThickenFactor, black ); endfor; fi; for j = 0 upto sN-1: signalvertex( LinFunc(j/theN), ThickenFactor, OutCol ); endfor; if JoinP: for j = -sN upto 0: signalvertex( LinFunc(j/theN), k, InCol ); endfor; fi; for i = sN upto theN: signalvertex( LinFunc( i/theN ), ThickenFactor, OutCol ); for j = sN downto 0: signalvertex( LinFunc( (i-j)/theN ), k, InCol ); endfor; endfor endgroup enddef; % Draw space-paths of possibly closed lines making use of "getready" def closedline( expr ThisIsClosed, theN, ForeFrac, BackFrac )( text LinFunc ) = begingroup save i, comm; numeric i; string comm; NCL := incr( NCL ); if ThisIsClosed: CLPath[NCL] := for i=1 upto theN: rp(LinFunc(i/theN)).. endfor cycle; for i=1 upto theN: comm:="draw subpath (" & decimal(i-ForeFrac) & "," & decimal(i+ForeFrac) & ") of CLPath" & decimal(NCL) & " withpen ForePen; undraw subpath (" & decimal(i-BackFrac) & "," & decimal(i+BackFrac) & ") of CLPath" & decimal(NCL) & " withpen BackPen;"; getready( comm, LinFunc(i/theN) ); endfor; else: CLPath[NCL] := rp(LinFunc(0)) for i=1 upto theN: ..rp(LinFunc(i/theN)) endfor; comm:="draw subpath (0," & decimal(ForeFrac) & ") of CLPath" & decimal(NCL) & " withpen ForePen; undraw subpath (0," & decimal(BackFrac) & ") of CLPath" & decimal(NCL) & " withpen BackPen;"; getready( comm, LinFunc(1/theN) ); for i=2 upto theN-1: comm:="draw subpath (" & decimal(i-ForeFrac) & "," & decimal(i+ForeFrac) & ") of CLPath" & decimal(NCL) & " withpen ForePen; undraw subpath (" & decimal(i-BackFrac) & "," & decimal(i+BackFrac) & ") of CLPath" & decimal(NCL) & " withpen BackPen;"; getready( comm, LinFunc(i/theN) ); endfor; comm:="draw subpath (" & decimal(theN-ForeFrac) & "," & decimal(theN) & ") of CLPath" & decimal(NCL) & " withpen ForePen; undraw subpath (" & decimal(theN-BackFrac) & "," & decimal(theN) & ") of CLPath" & decimal(NCL) & " withpen BackPen;"; getready( comm, LinFunc(1) ); fi endgroup enddef; % The next allows you to draw any solid that has no vertices and that has % two, exactly two, cyclic edges. In fact, it doesn't need to be a solid. % 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 ) = 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; TheLengthOfA = ( length CycleA ) - TheMargin; TheLengthOfB = ( length CycleB ) - TheMargin; SubPathA = subpath ( 0, TheLengthOfA ) of CycleA; SubPathB = subpath ( 0, TheLengthOfB ) of CycleB; PolygonPath = makepath makepen ( SubPathA--SubPathB--cycle ); Leng = (length PolygonPath) - 1; FinalPath = point 0 of PolygonPath for i = 1 upto Leng: --point i of PolygonPath endfor --cycle; ( FinalPath ) endgroup enddef; % Ellipse on the air. def ellipticpath(expr CenterPos, OneAxe, OtherAxe ) = begingroup save ind; numeric ind; ( for ind=0 upto 35: rp( CenterPos+planarrotation(OneAxe,OtherAxe,ind*10) )... endfor cycle ) endgroup enddef; % Shadow of an ellipse on the air. def ellipticshadowpath(expr CenterPos, OneAxe, OtherAxe ) = begingroup save ind; numeric ind; ( 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) = begingroup save labelpic, plak, lrc, ulc, llc, centerc, aratio, newbase; picture labelpic; pair lrc, ulc, llc; transform plak; color centerc, newbase; numeric aratio; labelpic = thelabel( SomeString, origin ); lrc = lrcorner labelpic; ulc = ulcorner labelpic; llc = llcorner labelpic; aratio = (xpart lrc - xpart llc)/(ypart ulc - ypart llc); if KeepRatio: newbase = conorm(UpVec)*aratio*N(BaseVec); else: newbase = BaseVec; fi; rp(RefPoi+newbase) = lrc transformed plak; rp(RefPoi+UpVec) = ulc transformed plak; centerc = RefPoi+0.5(newbase+UpVec); rp(RefPoi) = llc transformed plak; label( labelpic transformed plak, rp(centerc) ) endgroup enddef; % It should be possible to attach some path to some surface. def closedpathinspace( expr SomeTDPath, NDivide )( text ConverterFunc ) = begingroup save i, outpath, st; numeric i, st; path outpath; st = 1/NDivide; outpath = for i=st step st until (length SomeTDPath): ConverterFunc( point i of SomeTDPath ) -- endfor cycle; ( outpath ) endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Standard Objects: % And more precisely. The next routines spatialhalfcircle and % rigorousfear require circles drawn in a systematic and precise way. def goodcirclepath(expr CenterPos, AngulMom, Radius ) = begingroup save cirath, vecx, vecy, ind, goodangulmom, decision, view; numeric ind, decision; color vecx, vecy, goodangulmom, view; path cirath; view = f-CenterPos; decision = cdotprod( view, AngulMom ); if decision < 0: goodangulmom = -AngulMom; else: goodangulmom = AngulMom; fi; vecx = Radius*ncrossprod( view, goodangulmom ); decision := getangle( view, goodangulmom ); if decision > 0.5: %%%%%%%%%%%%%%% DANGER %%% 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; ( cirath ) endgroup enddef; % And its shadow. def circleshadowpath(expr CenterPos, AngulMom, Radius ) = begingroup save cirath, vecx, vecy, view, decision; numeric decision; color vecx, vecy, view; path cirath; view = LightSource-CenterPos; vecx = ncrossprod( view, AngulMom ); decision := getangle( view, AngulMom ); if decision > 0.5: %%%%%%%%%%%%%%% DANGER %%% vecy = ncrossprod( AngulMom, vecx ); cirath = ellipticshadowpath(CenterPos,vecx*Radius,vecy*Radius); else: vecx := N( (-Y(view), X(view), 0) ); vecy = ncrossprod( view, vecx ); cirath = ellipticshadowpath(CenterPos,vecx*Radius,vecy*Radius); fi; ( cirath ) endgroup enddef; % When there are numerical problems with the previous routine % use the following alternative: def head_on_circle(expr Pos, Radius ) = begingroup save vecx, vecy, ind, view; numeric ind; color vecx, vecy, view; view = f-Pos; 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; % The nearest or the furthest part of a circle returned as a path. % This function has been set to work for rigorousdisc (next). % Very tough settings they were. def spatialhalfcircle(expr Center, AngulMom, Radius, ItsTheNearest ) = begingroup 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; boolean choice; va := Center - f; vb := N( AngulMom ); vc := vb*( cdotprod( va, vb ) ); cc := f + vc; vd := cc - Center; % vd := va + vc; nr := conorm( vd ); if Radius >= nr: returnp := rp( cc ); % this single point will show up in spheroid else: valr := Radius*Radius; valx := valr/nr; valy := sqrt( valr - valx*valx ); ux := N( vd ); choiceang := getangle( vc, va ); %%%%%%%%%%%%% choice := ( choiceang < 89 ) or ( choiceang > 91 );%% DANGER % if choice: %%%%%%%%%%%%% uy := ncrossprod( vc, va ); else: uy := ncrossprod( AngulMom, va ); fi; pa := valx*ux + valy*uy + Center; pb := pa - 2*valy*uy; if choice: auxil := rp(1.1[Center,pb])--rp(0.9[Center,pb]); auxih := rp(1.1[Center,pa])--rp(0.9[Center,pa]); fcirc := goodcirclepath( Center, AngulMom, Radius ); if ItsTheNearest: returnp := (fcirc cutafter auxih) cutbefore auxil; else: returnp := (fcirc cutbefore auxih)..(fcirc cutafter auxil); fi; else: if ItsTheNearest: if cdotprod( va, AngulMom ) > 0: returnp := rp(pb)--rp(pa); else: returnp := rp(pa)--rp(pb); fi; else: if cdotprod( va, AngulMom ) < 0: returnp := rp(pb)--rp(pa); else: returnp := rp(pa)--rp(pb); fi; fi; fi; fi; ( returnp ) endgroup enddef; % Cylinders or tubes ( numeric, boolean, color, numeric, color ). % Great stuff. The "disc" in the name comes from the fact that % when SphericalDistortion := true; the sides of cylinders are % not drawn correctly (they are straight). And when it is a tube % you should force the background to be white. def rigorousdisc(expr InRay, FullFill, BaseCenter, Radius, LenVec) = begingroup save va, vb, vc, cc, vd, base, holepic; save vA, cC, nr, vala, valb, hashole, istube; save auxil, auxih, halfl, halfh, thehole; save auxili, auxihi, rect, theshadow; color va, vb, vc, cc, vd, base; picture holepic; color vA, cC; numeric nr, vala, valb; boolean hashole, istube; path auxil, auxih, halfl, halfh, thehole; path auxili, auxihi, rect, theshadow; va := BaseCenter - f; vb := N( LenVec ); vc := vb*( cdotprod( va, vb ) ); cc := f + vc; vd := cc - BaseCenter; nr := conorm( vd ); base := BaseCenter + LenVec; vA := base - f; vala := conorm( va ); valb := conorm( vA ); if ShadowOn: auxil := circleshadowpath( BaseCenter, LenVec, Radius ); auxih := circleshadowpath( base, LenVec, Radius ); fill twocyclestogether( auxil, auxih ); fi; auxil := goodcirclepath( base, LenVec, Radius ); auxih := goodcirclepath( BaseCenter, LenVec, Radius ); istube := false; hashole := false; if InRay > 0: istube := true; auxili := goodcirclepath( base, LenVec, InRay ); auxihi := goodcirclepath( BaseCenter, LenVec, InRay ); hashole := (-1,-1) <> ( auxili intersectiontimes auxihi ); if hashole: draw auxili; draw auxihi; holepic := currentpicture; clip holepic to auxili; clip holepic to auxihi; fi; fi; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if Radius >= nr: % THE CASE Radius > nr > InRay IS NOT SUPPORTED % if vala <= valb : thehole := auxil; auxil := auxih; auxih := thehole; fi; if istube: if vala <= valb : thehole := auxili; auxili := auxihi; auxihi := thehole; fi; holepic := currentpicture; clip holepic to auxihi; fi; unfill auxil; draw auxil; if istube: draw holepic; draw auxihi; fi; else: cC := base + vd; if ( cdotprod( f - cc, f - cC ) <= 0 ) or ( not FullFill ): halfl := spatialhalfcircle(BaseCenter,LenVec,Radius,true); halfh := spatialhalfcircle(base,LenVec,Radius,true); if FullFill: rect := halfl--halfh--cycle; else: rect := halfl--(reverse halfh)--cycle; fi; unfill rect; draw rect; elseif vala > valb: halfl := spatialhalfcircle(BaseCenter,LenVec,Radius,true); halfh := spatialhalfcircle(base,LenVec,Radius,false); rect := halfl--halfh--cycle; unfill rect; draw rect; if istube: if hashole: draw holepic; fi; draw auxili; fi; draw auxil; else: halfl := spatialhalfcircle(BaseCenter,LenVec,Radius,false); halfh := spatialhalfcircle(base,LenVec,Radius,true); rect := halfl--halfh--cycle; unfill rect; draw rect; if istube: if hashole: draw holepic; fi; draw auxihi; fi; draw auxih; fi; fi endgroup enddef; % And maybe a full cone border. The vertex may go anywhere. % Choose the full cone border (UsualForm=true) or just the nearest % part of the base edge (UsualForm=false). % This is used by tropicalglobe as a generic spatialhalfcircle to % 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) = begingroup 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, thelengthofc; pair pa, pb, pc, pd, pe; themargin = 0.02; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER basepath = goodcirclepath( CenterPos, AngulMom, Radius ); thelengthofc = ( length basepath ) - themargin; thesubpath = subpath ( 0, thelengthofc ) of basepath; fullpath = makepath makepen ( rp(VertexPos)--thesubpath--cycle ); pa = 0.995[rp(CenterPos),rp(VertexPos)]; pb = 1.005[rp(CenterPos),rp(VertexPos)]; auxpath = pa--pb; pc = auxpath intersectiontimes fullpath; if pc <> (-1,-1): auxt = ypart pc; newlen = length fullpath; if UsualForm: finalpath = point auxt of fullpath --point auxt+1 of fullpath for i = auxt+2 upto auxt+newlen-1: ...point i of fullpath endfor --cycle; else: bigcirc = goodcirclepath( CenterPos, AngulMom, 1.005*Radius ); pd = bigcirc intersectiontimes fullpath; pe = ( reverse bigcirc ) intersectiontimes fullpath; startt = floor( xpart pd ); endt = ceiling( ( length bigcirc ) - ( xpart pe ) ); finalpath = subpath (startt,endt) of basepath; fi; else: finalpath = rp(VertexPos); fi; ( finalpath ) endgroup enddef; def verygoodcone(expr BackDash,CenterPos,AngulMom,Radius,VertexPos) = begingroup save thepath, lenpath, bonevec, sidevec, viewaxe, cipath; save thelengthofc, thesubpath, themargin, basepath; color bonevec, sidevec, viewaxe; path thepath, cipath, basepath, thesubpath; numeric lenpath, thelengthofc, themargin; themargin = 0.02; bonevec = VertexPos - CenterPos; if cdotprod( bonevec, AngulMom ) < 0: sidevec = -N(AngulMom); else: sidevec = N(AngulMom); fi; viewaxe = f-CenterPos; if ShadowOn: basepath = circleshadowpath( CenterPos, AngulMom, Radius ); thelengthofc = ( length basepath ) - themargin; thesubpath = subpath ( 0, thelengthofc ) of basepath; fill makepath makepen ( rp(cb(VertexPos))--thesubpath--cycle ); fi; thepath = rigorouscone(true,CenterPos,AngulMom,Radius,VertexPos); lenpath = length thepath; if lenpath<>0: unfill thepath; draw thepath; if cdotprod( sidevec, viewaxe ) < 0: draw goodcirclepath( CenterPos, AngulMom, Radius ); else: if BackDash: draw goodcirclepath( CenterPos, AngulMom, Radius ) dashed evenly; fi; fi; else: cipath = goodcirclepath( CenterPos, AngulMom, Radius ); unfill cipath; draw cipath; if cdotprod( sidevec, viewaxe ) > 0: draw rp( VertexPos ); fi; fi endgroup enddef; % Its a sphere, don't fear, but remember that the rigorous projection % of a sphere is an ellipse. def rigorousfearpath(expr Center, Radius ) = begingroup save auxil, ux, uy, newcen, nr, valx, valy, valr; color ux, uy, newcen; numeric nr, valx, valy, valr; path auxil; nr := conorm( Center - f ); valr := Radius**2; valx := valr/nr; valy := sqrt( valr - valx**2 ); newcen := valx*( f - Center )/nr; auxil := head_on_circle( Center+newcen, valy ); ( auxil ) endgroup enddef; def rigorousfearshadowpath(expr Center, Radius ) = begingroup 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; pair dcenter; nr := conorm( Center - LightSource ); valr := Radius**2; valx := valr/nr; valy := sqrt( valr - valx**2 ); newcen := valx*( LightSource - Center )/nr; auxil := circleshadowpath( Center+newcen, newcen, valy ); ( auxil ) endgroup enddef; % It's a globe (without land). def tropicalglobe( expr NumLats, TheCenter, Radius, AngulMom ) = begingroup save viewaxe, sinalfa, sinbeta, globaxe, aux, limicos, lc; save stepang, actang, newradius, foc, newcenter, cpath, i; save outerpath, conditiona, conditionb; color viewaxe, globaxe, foc, newcenter; numeric sinalfa, sinbeta, aux, limicos, stepang, actang; numeric newradius, lc, i; path cpath, outerpath; boolean conditiona, conditionb; if ShadowOn: fill rigorousfearshadowpath( TheCenter, Radius ); fi; viewaxe = f-TheCenter; sinalfa = Radius/conorm( viewaxe ); aux = cdotprod( viewaxe, AngulMom ); if aux < 0: globaxe = -N(AngulMom); else: globaxe = N(AngulMom); fi; sinbeta = cdotprod( globaxe, N(viewaxe) ); aux := sqrt((1-sinalfa**2)*(1-sinbeta**2)); limicos = aux - sinalfa*sinbeta; stepang = 180/NumLats; globaxe := globaxe*Radius; outerpath = rigorousfearpath(TheCenter,Radius); unfill outerpath; draw outerpath; for actang = 0.5*stepang step stepang until 179: if cosd(actang) < limicos-0.005: %%%%%%%%%%%%%%%%%%%%%%%% DANGER newradius := Radius*sind(actang); newcenter := TheCenter - globaxe*cosd(actang); conditiona:=(actang<94) and (actang>86); % DANGER % DANGER VV conditionb:=abs(cdotprod(globaxe/Radius,N(f-newcenter)))<0.08; if conditiona or conditionb: draw spatialhalfcircle(newcenter,globaxe,newradius,true); else: foc := TheCenter - globaxe/cosd(actang); lena := -Radius*cosd(actang); lenb := cdotprod(viewaxe,globaxe/Radius); if (actang <= 86) or ((lenb=94)): cpath := rigorouscone(false,newcenter,globaxe,newradius,foc); draw cpath; else: cpath := rigorouscone(true,newcenter,globaxe,newradius,foc); lc := length cpath; if lc <> 0: draw subpath (1,lc-1) of cpath; else: draw rigorouscircle( newcenter,globaxe,newradius ); fi; fi; fi; fi; endfor endgroup enddef; % An elliptical frustum: def whatisthis(expr CenterPos,OneAxe,OtherAxe,CentersDist,TheFactor) = begingroup save patha, pathb, pathc, centersvec, noption; path patha, pathb, pathc; color centersvec; numeric noption; centersvec = CentersDist*ncrossprod( OneAxe, OtherAxe ); if ShadowOn: patha = ellipticshadowpath( CenterPos, OneAxe, OtherAxe ); pathb = ellipticshadowpath( CenterPos+centersvec, TheFactor*OneAxe, TheFactor*OtherAxe ); pathc = twocyclestogether( patha, pathb ); fill pathc; fi; patha := ellipticpath( CenterPos, OneAxe, OtherAxe ); pathb := ellipticpath( CenterPos+centersvec, TheFactor*OneAxe, TheFactor*OtherAxe ); pathc := twocyclestogether( patha, pathb ); unfill pathc; draw pathc; noption = cdotprod( centersvec, f-CenterPos ); if noption > (CentersDist**2): draw pathb; elseif noption < 0: draw patha; fi endgroup enddef; % Probably the last algorithm I'm going to write for featpost... % 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! % 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 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=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 yneg0: 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. def fakehole( expr CenterPos, LenVec, Radius ) = begingroup save patha, pathb, pathc, noption, hashole, auxv, poption, vv; save ta, tb, taf, tbf, margint, stopair, pa, pb, testpath, isin; path patha, pathb, pathc, pa, pb, testpath; numeric noption, ta, tb, margint; boolean hashole, poption, isin; color auxv, vv; pair stopair; vv = f-CenterPos; patha := rigorouscircle( CenterPos, LenVec, Radius ); pathb := rigorouscircle( CenterPos+LenVec, LenVec, Radius ); % patha := goodcirclepath( CenterPos, LenVec, Radius ); % pathb := goodcirclepath( CenterPos+LenVec, LenVec, Radius ); auxv := ncrossprod( LenVec, ccrossprod( vv, LenVec ) ); poption := abs( cdotprod( vv, auxv ) ) <= 1.05*Radius;% DANGER! if poption: draw patha; draw pathb; else: noption = cdotprod( LenVec, vv ); if noption > (conorm(LenVec)**2): pa = patha; pb = pathb; elseif noption < 0: 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; fi endgroup enddef; % 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 ) = begingroup save star, pos, patw, patb, refv, near, centre, farv; save newa, newb, newc, veca, vecb, vecc, auxx, auxy, i; color star, pos[], refv, near, newa, newb, newc; color veca, vecb, vecc, auxx, auxy, centre, farv; path patw, patb; numeric i; veca = ( cosd(AngA)*cosd(AngB), sind(AngA)*cosd(AngB), sind(AngB) ); auxx = ( cosd(AngA+90), sind(AngA+90), 0 ); auxy = ccrossprod( veca, auxx ); vecb = cosd(AngC)*auxx + sind(AngC)*auxy; vecc = cosd(AngC+90)*auxx + sind(AngC+90)*auxy; veca := LenA*veca; vecb := LenB*vecb; vecc := LenC*vecc; if IsVertex: star = RefP; centre = RefP + 0.5*( veca + vecb + vecc); else: star = RefP - 0.5*( veca + vecb + vecc); centre = RefP; fi; pos1 = star + veca; pos2 = pos1 + vecb; pos3 = pos2 + vecc; pos4 = pos3 - vecb; pos5 = pos4 - veca; pos6 = pos5 + vecb; pos7 = pos6 - vecc; if ShadowOn: patw = rp(cb(star))--rp(cb(pos1))--rp(cb(pos2)) --rp(cb(pos3))--rp(cb(pos4)) --rp(cb(pos5))--rp(cb(pos6))--rp(cb(pos7))--cycle; patb = makepath makepen patw; fill patb; fi; patw := rp(star)--rp(pos1)--rp(pos2)--rp(pos3)--rp(pos4) --rp(pos5)--rp(pos6)--rp(pos7)--cycle; patb := makepath makepen patw; unfill patb; draw patb; i = 0; if inplanarvolume( star, star+veca, f ): i := incr( i ); fi; if inplanarvolume( star, star+vecb, f ): i := incr( i ); fi; if inplanarvolume( star, star+vecc, f ): i := incr( i ); fi; if (i=2) and WithDash: message "Unable to dash kindofcube " & cstr( RefP ) & "."; elseif i = 3: message "f is inside kindofcube " & cstr( RefP ) & "."; else: refv = f - centre; if cdotprod( refv, veca ) > 0: newa = -veca; else: newa = veca; fi; if cdotprod( refv, vecb ) > 0: newb = -vecb; else: newb = vecb; fi; if cdotprod( refv, vecc ) > 0: newc = -vecc; else: newc = vecc; fi; near = centre - 0.5*( newa + newb + newc ); draw rp(near)--rp(near+newa); draw rp(near)--rp(near+newb); draw rp(near)--rp(near+newc); if WithDash: if i=1: message "Unable to dash kindofcube " & cstr( RefP ) & "."; else: farv = centre + 0.5*( newa + newb + newc ); draw rp(farv)--rp(farv-newa) dashed evenly; draw rp(farv)--rp(farv-newb) dashed evenly; draw rp(farv)--rp(farv-newc) dashed evenly; fi; fi; fi endgroup enddef; % It's a bit late now but the stage must be set. def setthestage( expr NumberOfSideSquares, SideSize ) = begingroup save i, j, squaresize, squarepath, ca, cb, cc, cd; numeric i, j, squaresize; path squarepath; color ca, cb, cc, cd; squaresize = SideSize/(2*NumberOfSideSquares-1); for i=-0.5*SideSize step 2*squaresize until 0.5*SideSize: for j=-0.5*SideSize step 2*squaresize until 0.5*SideSize: ca := (i,j,HoriZon); cb := (i,j+squaresize,HoriZon); cc := (i+squaresize,j+squaresize,HoriZon); cd := (i+squaresize,j,HoriZon); squarepath := rp(ca)--rp(cb)--rp(cc)--rp(cd)--cycle; unfill squarepath; draw squarepath; endfor; endfor endgroup enddef; def setthearena( expr NumberOfDiameterCircles, ArenaDiameter ) = begingroup save i, j, circlesize, polar, currpos, phi, cpath; numeric i, j, circlesize, polar, phi; color currpos; path cpath; circlesize = ArenaDiameter/NumberOfDiameterCircles; for i=0.5*ArenaDiameter step -circlesize until 0.4*circlesize: polar := floor(6.28318*i/circlesize); for j=1 upto polar: phi := 360*j/polar; currpos := i*(cosd(phi),sind(phi),0)+(0,0,HoriZon); cpath := rigorouscircle( currpos, blue, 0.3*circlesize); unfill cpath; draw cpath; endfor; endfor endgroup enddef; % And a transparent dome. The angular momentum vector is supposed % to point from the concavity of the dome and into outer space. % The pen can only be changed with a previous drawoptions(). def spatialhalfsfear(expr Center, AngulMom, Radius ) = begingroup save spath, cpath, fpath, rpath, cutp; path spath, cpath, fpath, rpath, cutp; save ap, bp, cp, dp, cuti, cute, vp; pair ap, bp, cp, dp, cuti, cute, vp; save auxcos, actcos, actsin, auxsin; numeric auxcos, actcos, actsin, auxsin; picture partoffear; save A, B; color A, B; spath = rigorousfearpath( Center, Radius ); auxcos = getcossine( Center, Radius ); actcos = cdotprod( N( f - Center ), N( AngulMom ) ); actsin = sqrt(1-actcos**2); auxsin = sqrt(1-auxcos**2); if actsin <= auxcos: if actcos >= 0: cpath = goodcirclepath( Center, AngulMom, Radius ); draw cpath; else: draw spath; fi; else: fpath = spatialhalfcircle( Center, AngulMom, Radius, true ); rpath = spatialhalfcircle( Center, AngulMom, Radius, false ); cuti = point 0 of rpath; cute = point ( length rpath ) of rpath; ap = 1.1[cuti,cute]; bp = 1.1[cute,cuti]; partoffear = nullpicture; addto partoffear doublepath spath; A = ncrossprod( f-Center, ncrossprod( f-Center, AngulMom ) ); B = Center + 1.1*Radius*( auxcos*N( f-Center ) + auxsin*A ); vp = rp(B) - rp(Center); cp = ap + vp; dp = bp + vp; cutp = ap--cp--dp--bp--cycle; clip partoffear to cutp; draw fpath; draw partoffear; if actcos >= 0: draw rpath; fi; fi endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Toroidal Stuff % Take a donut. def smoothtorus( expr Tcenter, Tmoment, Bray, Sray ) = begingroup save nearaxe, sideaxe, viewline, circlecenter, circlemoment; save ang, ind, i, anglim, angstep, cuspcond, holepic, coofrac; save cpath, apath, opath, ipath, wp, ep, refpair, distance, lr; color nearaxe, sideaxe, viewline, circlecenter, circlemoment; numeric ang, ind, i, anglim, angstep, distance, coofrac, lr; path cpath, apath, opath, ipath, wp, ep; pair outerp[], innerp[], refpair; boolean cuspcond; picture holepic; 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; else: tmoment = Tmoment; fi; refpair = unitvector( rp(Tcenter+tmoment)-rp(Tcenter) ); sideaxe = Bray*ncrossprod( tmoment, viewline ); nearaxe = Bray*ncrossprod( sideaxe, tmoment ); coofrac = cdotprod( viewline, N( tmoment ) )/Sray; if (coofrac <= 1.04) and (coofrac >= 1.01): %%%%%%%%%%% DANGER! ind = 360/angstep; anglim = 0.5*angstep; for i=1 upto ind: ang := i*angstep-anglim-180.0; circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter; circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang); cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true); if i >= 0.5*ind+1: outerp[i]=point 0 of cpath; else: outerp[i]=point (length cpath) of cpath; fi; endfor; opath = for i=1 upto ind: outerp[i].. endfor cycle; unfill opath; draw opath; elseif coofrac < 1.01: distance = conorm( viewline ); lr = Bray + Sray*( 1 +-+ coofrac ); anglim = angle( ( lr, distance +-+ lr ) ); ind = 2*floor(anglim/angstep); angstep := 2*anglim/(ind+1); for i=0 upto 0.5*ind-1: ang := i*angstep-anglim; circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter; circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang); cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true); innerp[i]=point 0 of cpath; outerp[i]=point (length cpath) of cpath; endfor; for i=0.5*ind upto ind-2: ang := (i+2)*angstep-anglim; circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter; circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang); cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true); outerp[i]=point 0 of cpath; innerp[i]=point (length cpath) of cpath; endfor; if coofrac > 0.94: apath = innerp0 for i=1 upto ind-2: ..innerp[i] endfor --cycle; else: apath = innerp0 for i=2 upto ind-2: ..innerp[i] endfor ..outerp[ind-2] for i=ind-3 downto 0: ..outerp[i] endfor ..cycle; fi; unfill apath; draw apath; else: ind = 360/angstep; anglim = 0.5*angstep; for i=1 upto ind: ang := i*angstep-anglim-180.0; circlecenter:= nearaxe*cosd(ang)+sideaxe*sind(ang)+Tcenter; circlemoment:=-nearaxe*sind(ang)+sideaxe*cosd(ang); cpath:=spatialhalfcircle(circlecenter,circlemoment,Sray,true); if i >= 0.5*ind+1: outerp[i]=point 0 of cpath; innerp[i]=point (length cpath) of cpath; else: innerp[i]=point 0 of cpath; outerp[i]=point (length cpath) of cpath; fi; endfor; opath = for i=1 upto ind: outerp[i].. endfor cycle; ipath = for i=1 upto ind: innerp[i].. endfor cycle; holepic = currentpicture; clip holepic to ipath; unfill opath; draw holepic; draw opath; draw ipath; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Perhaps there is an analytic way of getting the angle of the cusp point? %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% i := ceiling(1+0.5*ind); cuspcond = false; forever: i := incr( i ); exitif i > ind-1; cuspcond := refpair dotprod innerp[i+1] < refpair dotprod innerp[i]; exitif cuspcond; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if cuspcond: show "Torus shows cusp points."; ep = outerp[ind-i+1]--innerp[ind-i+1]; wp = innerp[i]--outerp[i]; unfill buildcycle(reverse opath,ep,ipath,wp); draw opath; draw subpath (i-1,ind-i) of ipath; fi; fi; 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). def quartertorus( expr Tcenter, Tstart, Tfinis, Sray ) = begingroup save sideaxe, viewline, circlecenter, circlemoment, amax; save i, angstep, tmoment, cpath, opath, ipath, fpath, finc; save vstart, vfinis, ostart, ofinis, cstart, cfinis, finis; color sideaxe, viewline, circlecenter, circlemoment; color tmoment, vstart, vfinis, finis, finc; numeric i, angstep, amax; path cpath, opath, ipath, cstart, cfinis, fpath; pair outerp[], innerp[]; boolean ostart, ofinis; angstep = 4; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DANGER! viewline = f-Tcenter; tmoment = ncrossprod( Tstart, Tfinis ); vstart = ncrossprod( Tstart, tmoment ); vfinis = ncrossprod( tmoment, Tfinis ); finis = -vstart*conorm( Tstart ); amax = getangle( Tstart, Tfinis ); finc = N( Tfinis )*conorm( Tstart ); if cdotprod( viewline, tmoment ) < 0: tmoment := -tmoment; fi; sideaxe = ncrossprod( tmoment, viewline ); for i=0 step angstep until amax: circlecenter:= Tstart*cosd(i)+finis*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; cstart = spatialhalfcircle(Tstart+Tcenter,vstart,Sray,true); if cdotprod( sideaxe, Tstart ) > 0: cstart := reverse cstart; fi; cfinis = spatialhalfcircle(finc+Tcenter,vfinis,Sray,true); if cdotprod( sideaxe, finc ) < 0: cfinis := reverse cfinis; fi; opath = outerp0 for i=angstep step angstep until amax: ..outerp[i/angstep] endfor; ipath = innerp0 for i=angstep step angstep until amax: ..innerp[i/angstep] endfor; fpath = ipath---cfinis---reverse opath---cstart---cycle; unfill fpath; draw fpath; ostart = cdotprod( viewline-Tstart, vstart ) > 0; if ostart: cpath := rigorouscircle( Tcenter+Tstart, vstart, Sray ); unfill cpath; draw cpath; fi; ofinis = cdotprod( viewline-finc, vfinis ) > 0; if ofinis: cpath := rigorouscircle( Tcenter+finc, vfinis, Sray ); unfill cpath; draw cpath; fi; %draw opath withcolor blue; %draw ipath withcolor red; 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 (distBray-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: def positivecharge( expr InFactPositive, Center, BallRay ) = begingroup save auxc, axehorf, axeside, viewline, spath, pa, pb, pc, pd; color auxc, axehorf, axeside, viewline, pa, pb, pc, pd; path spath; viewline = f - Center; axehorf = N( ( X(viewline), Y(viewline), 0 ) ); axeside = ccrossprod( axehorf, blue ); if ShadowOn: fill rigorousfearshadowpath( Center, BallRay ); fi; spath = rigorousfearpath( Center, BallRay ); unfill spath; draw spath; auxc = Center + sqrt(3)*axehorf; pa = auxc + axeside; pb = auxc - axeside; angline( pa, pb, Center, BallRay, "", top ); if InFactPositive: pc = auxc + blue; pd = auxc - blue; angline( pc, pd, Center, BallRay, "", top ); fi endgroup enddef; 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; save fmar, fthi, fray, rmar, rthi, rray, inrefp;; color veca, auxx, auxy, vecb, vecc, viewline; color fl, fr, rl, rr, inrefp; numeric anga, angb, angc, lena, lenb, lenc, auxm, auxn; numeric fmar, fthi, fray, rmar, rthi, rray; anga = X( AngCol ); angb = Y( AngCol ); angc = Z( AngCol ); lena = X( LenCol ); lenb = Y( LenCol ); lenc = Z( LenCol ); fmar = X( FronWheelCol ); fthi = Y( FronWheelCol ); fray = Z( FronWheelCol ); rmar = X( RearWheelCol ); rthi = Y( RearWheelCol ); rray = Z( RearWheelCol ); veca = ( cosd(anga)*cosd(angb), sind(anga)*cosd(angb), sind(angb) ); auxx = ( cosd(anga+90), sind(anga+90), 0 ); auxy = ccrossprod( veca, auxx ); vecb = cosd(angc)*auxx + sind(angc)*auxy; vecc = cosd(angc+90)*auxx + sind(angc+90)*auxy; viewline = f - RefP; auxm = cdotprod( viewline, veca ); auxn = cdotprod( viewline, vecb ); inrefp = RefP - 0.5*lenc*vecc; fl = inrefp + (0.5*lena-fmar-fray)*veca + 0.5*lenb*vecb; fr = inrefp + (0.5*lena-fmar-fray)*veca - 0.5*lenb*vecb; rl = inrefp - (0.5*lena-rmar-rray)*veca + 0.5*lenb*vecb; rr = inrefp - (0.5*lena-rmar-rray)*veca - 0.5*lenb*vecb; if auxn > 0.5*lenb: if auxm > 0: rigorousdisc( 0, true, rr, rray, -rthi*vecb ); rigorousdisc( 0, true, fr, fray, -fthi*vecb ); kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc); rigorousdisc( 0, true, rl, rray, rthi*vecb ); rigorousdisc( 0, true, fl, fray, fthi*vecb ); else: rigorousdisc( 0, true, fr, fray, -fthi*vecb ); rigorousdisc( 0, true, rr, rray, -rthi*vecb ); kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc); rigorousdisc( 0, true, fl, fray, fthi*vecb ); rigorousdisc( 0, true, rl, rray, rthi*vecb ); fi; elseif auxn < -0.5*lenb: if auxm > 0: rigorousdisc( 0, true, rl, rray, rthi*vecb ); rigorousdisc( 0, true, fl, fray, fthi*vecb ); kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc); rigorousdisc( 0, true, rr, rray, -rthi*vecb ); rigorousdisc( 0, true, fr, fray, -fthi*vecb ); else: rigorousdisc( 0, true, fl, fray, fthi*vecb ); rigorousdisc( 0, true, rl, rray, rthi*vecb ); kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc); rigorousdisc( 0, true, fr, fray, -fthi*vecb ); rigorousdisc( 0, true, rr, rray, -rthi*vecb ); fi; else: if auxm > 0: rigorousdisc( 0, true, rl, rray, rthi*vecb ); rigorousdisc( 0, true, fl, fray, fthi*vecb ); rigorousdisc( 0, true, rr, rray, -rthi*vecb ); rigorousdisc( 0, true, fr, fray, -fthi*vecb ); kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc); else: rigorousdisc( 0, true, fl, fray, fthi*vecb ); rigorousdisc( 0, true, rl, rray, rthi*vecb ); rigorousdisc( 0, true, fr, fray, -fthi*vecb ); rigorousdisc( 0, true, rr, rray, -rthi*vecb ); kindofcube(false,false,RefP,anga,angb,angc,lena,lenb,lenc); fi; fi endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Differential Equations: % 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 ) = begingroup save kone, ktwo; color kone, ktwo; kone = Step*VecFunc( Spos ); ktwo = Step*VecFunc( Spos+0.5*kone ); ( Spos+ktwo ) endgroup enddef; def fieldlinepath( expr Numb, Spos, Step )( text VecFunc ) = begingroup save ind, flpath, prevpos, thispos; numeric ind; color prevpos, thispos; path flpath; prevpos = Spos; flpath = rp( Spos ) for ind=1 upto Numb: hide( thispos := fieldlinestep( prevpos, Step, VecFunc ) ) ..rp( thispos ) hide( prevpos := thispos ) endfor; ( flpath ) endgroup 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 ) = begingroup save ind, flpath, prevpos, thispos, prevvel, thisvel; save rone, rtwo, vone, vtwo; numeric ind; color prevpos, thispos, prevvel, thisvel; color rone, rtwo, vone, vtwo; path flpath; prevpos = Spos; prevvel = Svel; flpath = rp( Spos ) for ind=1 upto Numb: hide( vone := Step*VecFunc( prevpos ); rone := Step*prevvel; vtwo := Step*VecFunc( prevpos+0.5*rone ); rtwo := Step*( prevvel+0.5*vone ); thisvel := prevvel+vtwo; thispos := prevpos+rtwo ) ..rp( thispos ) hide( prevpos := thispos; prevvel := thisvel ) endfor; ( flpath ) endgroup enddef; % Another point is that I want to draw trajectories in space and % dependent 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 ) = 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 ) for ind=1 upto Numb: hide( vone := Step*ccrossprod( VecFunc( prevpos ), prevvel ); rone := Step*prevvel; vtwo := Step*ccrossprod(VecFunc(prevpos+0.5*rone),prevvel+0.5*vone); rtwo := Step*( prevvel+0.5*vone ); vthr := Step*ccrossprod(VecFunc(prevpos+0.5*rtwo),prevvel+0.5*vtwo); rthr := Step*( prevvel+0.5*vtwo ); vfou := Step*ccrossprod( 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 ) ..rp( thispos ) hide( prevpos := thispos; prevvel := thisvel ) endfor; ( flpath ) endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Part II: %%%%%%%%%%%%%%%%%%%%%%%%%%%% 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Definition of a 3D-Object % define vertices V1 := (1,0,0); V2 := (0,0,0); V3 := (0,1,0); V4 := (-0.3,0.2,1); V5 := (1,0,1); V6 := (0,1,1); V7 := (0,0,2); V8 := (-0.5,0.6,1.2); V9 := (0.6,-0.5,1.2); makeline1(8,9); makeface1(1,2,7); makeface2(2,3,7); makeface3(5,4,6); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% or the old way below %%%%%%%%% % define lines % NL := 1; % number of lines % npl1 := 2; % number of vertices of the first line % L1p1 := V8; % L1p2 := V9; % define faces % NF := 3; % number of faces % npf1 := 3; % number of vertices of the first face % F1p1 := V1; % F1p2 := V2; % F1p3 := V7; % npf2 := 3; % number of vertices of the second face % F2p1 := V2; % F2p2 := V3; % F2p3 := V7; % npf3 := 3; % number of vertices of the third face % F3p1 := V5; % F3p2 := V4; % F3p3 := V6; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Flip first argument accordingly to the second def flipvector(expr A, B) = begingroup save nv; color nv; if cdotprod( A, B) < 0 : nv = -A; else: nv = A; fi; ( nv ) endgroup enddef; % Frontside of a face given by three of its vertices def facevector(expr A, B, C) = begingroup save nv; color nv; nv = ncrossprod( A-B, B-C ); ( flipvector( nv, f-B ) ) endgroup enddef; % Center or inside of a face def masscenter(expr Nsides)(suffix Coords) = begingroup save mc, counter; numeric counter; color mc; mc = (0,0,0); for counter=1 upto Nsides: mc := mc + Coords[counter]; endfor; ( mc / Nsides ) endgroup enddef; % Direction of coverability. The trick is here. % The condition for visibility is % that the angle beetween the vector that goes from the side of a face to % the mark position and the covervector must be greater than 90 degrees. def covervector(expr A, B, MassCenter) = begingroup save nv; color nv; nv = ncrossprod( A-f, B-f ); ( flipvector( nv, MassCenter-B) ) endgroup enddef; % O.K., the following macro tests the visibility of a point def themarkisinview(expr Mark, OwnFace) = begingroup save c, faceVec, centerPoint, coverVec, inview, l, m; color c, faceVec, centerPoint, coverVec; boolean inview; numeric l, m; l = 0; forever: % cycle over faces until the mark is covered l := incr(l); if l = OwnFace: l := incr(l); fi; exitif l > NF; faceVec := facevector(F[l]p1,F[l]p2,F[l]p3); inview := true; if cdotprod(Mark-F[l]p1, faceVec) < 0: centerPoint := masscenter(npf[l], F[l]p); m := 0; forever: % cycle over segments of a face m := incr(m); exitif m > npf[l]; if m < npf[l]: c := F[l]p[m+1]; else: c := F[l]p1; fi; coverVec := covervector(F[l]p[m], c, centerPoint); inview := cdotprod(Mark-c,coverVec) <= 0; exitif inview; endfor; fi; exitif not inview; endfor; ( inview ) endgroup enddef; % Check for possible intersection or crossing. def maycrossviewplan(expr Ea, Eb, La, Lb) = begingroup ( abs( cdotprod( ccrossprod(Ea-f,Eb-f), La-Lb ) ) > 0.001 ) endgroup enddef; % Calculate the intersection of two sides. This is very nice. def crossingpoint(expr Ea, Eb, La, Lb) = begingroup save thecrossing, perpend, exten, aux; color thecrossing, perpend; numeric exten, aux; if ( Ea = Lb ) or ( Ea = La ): thecrossing = Ea; elseif ( Eb = Lb ) or ( Eb = La ): thecrossing = Eb; else: perpend = ccrossprod( Ea-f, Eb-f ); if conorm( perpend ) = 0: thecrossing = Eb; else: aux = cdotprod( perpend, f ); cdotprod( perpend, thecrossing ) = aux; ( La-Lb )*exten = La-thecrossing; fi; fi; ( thecrossing ) endgroup enddef; % Calculate the intersection of an edge and a face. def crossingpointf(expr Ea, Eb, Fen) = begingroup save thecrossing, perpend, exten; color thecrossing, perpend; numeric exten; perpend = ccrossprod( F[Fen]p1-F[Fen]p2, F[Fen]p3-F[Fen]p2 ); cdotprod(perpend,thecrossing) = cdotprod( perpend, F[Fen]p2 ); ( Ea-Eb )*exten = Ea-thecrossing; ( thecrossing ) endgroup enddef; % Check for possible intersection of an edge and a face. def maycrossviewplanf(expr Ea, Eb, Fen) = begingroup save perpend; color perpend; perpend = ccrossprod( F[Fen]p1-F[Fen]p2, F[Fen]p3-F[Fen]p2 ); ( abs( cdotprod( perpend, Ea-Eb ) ) > 0.001 ) endgroup enddef; % The intersection point must be within the extremes of the segment. def insidedge(expr Point, Ea, Eb) = begingroup save fract; numeric fract; fract := cdotprod( Point-Ea, Point-Eb ); ( fract < 0 ) endgroup enddef; % Skip edges that are too far away def insideviewsphere(expr Ea, Eb, La, Lb) = begingroup save furthestofedge, nearestofline, flag, exten; color nearestofline, furthestofedge; boolean flag; numeric exten; nearestofline = La+exten*(Lb-La); cdotprod( nearestofline-f, Lb-La ) = 0; if conorm(Ea-f) < conorm(Eb-f): furthestofedge := Eb; else: furthestofedge := Ea; fi; if conorm(nearestofline-f) < conorm(furthestofedge-f): flag := true; else: flag := false; fi; ( flag ) endgroup enddef; % The intersection point must be within the triangle defined by % three points. Really smart. def insidethistriangle(expr Point, A, B, C ) = begingroup save arep, area, areb, aret, flag; color arep, area, areb, aret; boolean flag; aret = ccrossprod( A-C, B-C ); arep = flipvector( ccrossprod( C-Point, A-Point ), aret ); area = flipvector( ccrossprod( A-Point, B-Point ), aret ); areb = flipvector( ccrossprod( B-Point, C-Point ), aret ); flag = ( conorm( arep + area + areb ) <= 2*conorm( aret ) ); ( flag ) endgroup enddef; % The intersection point must be within the triangle defined by the % point of view f and the extremes of some edge. def insideviewtriangle(expr Point, Ea, Eb) = insidethistriangle( Point, Ea, Eb, f ) enddef; % The intersection point must be within the face def insidethisface(expr Point, FaN) = begingroup save flag, m, central; boolean flag; numeric m; color central; m = npf[FaN]; central = masscenter( m, F[FaN]p ); flag = insidethistriangle( Point, central, F[FaN]p[m], F[FaN]p[1] ); for m=2 upto npf[FaN]: exitif flag; flag := insidethistriangle( Point, central, F[FaN]p[m-1], F[FaN]p[m] ); endfor; ( flag ) endgroup enddef; % Draw the visible parts of a straight line in beetween points A and B % changing the thickness of the line accordingly to the distance from f def coarse_line(expr A, B, Facen, Press, Col) = begingroup save k, mark, stepVec; numeric k; color mark, stepVec; stepVec := resolvec(A,B); k := 0; forever: % cycle along a whole segment mark := A+(k*stepVec); exitif cdotprod(B-mark,stepVec) < 0; if themarkisinview(mark,Facen): signalvertex(mark, Press, Col); fi; k := incr(k); endfor endgroup enddef; % Get the 2D rigorous projection path of a face. % Don't use SphericalDistortion here. def facepath(expr Facen) = begingroup save thispath, counter; path thispath; numeric counter; thispath = rp(F[Facen]p[1])-- for counter=2 upto (npf[Facen]): rp(F[Facen]p[counter])-- endfor cycle; ( thispath ) endgroup enddef; def faceshadowpath(expr Facen) = begingroup save thispath, counter; path thispath; numeric counter; thispath = rp(cb(F[Facen]p[1]))-- for counter=2 upto (npf[Facen]): rp(cb(F[Facen]p[counter]))-- endfor cycle; ( thispath ) endgroup enddef; % FillDraw a face def face_invisible( expr Facen )( text LineAtribs ) = begingroup save ghost; path ghost; ghost = facepath( Facen ); unfill ghost; draw ghost LineAtribs endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Different kinds of renderers: % Draw only the faces, calculating edge-intersections. % Mind-blogging kind of thing. % Only two constraints: i) faces must be planar and % ii) faces must be convex. % Pay attention: depends on PrintStep (resolvec). def sharpraytrace %% ( expr LabelCrossPoints ) = begingroup save i, j, k, l, counter, a, b, c, d, currcross; save flag, swapc, swapn, somepoint, frac, exten; save trythis, refpath, otherpath, intertimes; save counter, infolabel; numeric i, j, k, l, counter, swapn; color a, b, c, d, currcross, swapc; boolean flag, trythis; path refpath, otherpath; pair intertimes; string infolabel; color crosspoin[]; numeric sortangle[]; for i=1 upto NF: % scan all faces for j=1 upto npf[i]: % scan all edges a := F[i]p[j]; if j <> npf[i]: b := F[i]p[j+1]; else: b := F[i]p1; fi; crosspoin[1] := a; counter := 2; refpath := rp(a)--rp(b); % The limits a and b of one % side of one face for k=1 upto NF: otherpath := facepath(k); intertimes := refpath intersectiontimes otherpath; trythis := xpart intertimes <> 0; if trythis and (xpart intertimes <> 1) and (k <> i): for l=1 upto npf[k]: c := F[k]p[l]; if l < npf[k]: d := F[k]p[l+1]; else: d := F[k]p1; fi; if insideviewsphere( a, b, c, d ): if maycrossviewplan( a, b, c, d ): currcross := crossingpoint( a, b, c, d ); if insideviewtriangle( currcross, a, b ): if insidedge( currcross, c, d ): swapc := ccrossprod( a-b, f-currcross); swapc := ccrossprod(swapc,f-currcross); color somepo; numeric fract; (b-a)*fract = somepo-a; cdotprod(swapc,somepo) =cdotprod(swapc,f); if (fract>0) and (fract<1): crosspoin[counter] := somepo; counter := incr(counter); fi; fi; fi; fi; fi; endfor; if maycrossviewplanf( a, b, k ): currcross := crossingpointf( a, b, k ); if insidethisface( currcross, k ): if insidedge( currcross, a, b ): crosspoin[counter] := currcross; counter := incr(counter); fi; fi; fi; fi; endfor; crosspoin[counter] := b; sortangle[1] := 0; for k=2 upto counter: sortangle[k] := conorm(crosspoin[k]-a); endfor; forever: flag := true; for k=2 upto counter: if sortangle[k] < sortangle[k-1]: swapn := sortangle[k-1]; sortangle[k-1] := sortangle[k]; sortangle[k] := swapn; swapc := crosspoin[k-1]; crosspoin[k-1] := crosspoin[k]; crosspoin[k] := swapc; flag := false; fi; endfor; exitif flag; endfor; for k=2 upto counter: swapc := resolvec(crosspoin[k-1],crosspoin[k]); flag := themarkisinview( crosspoin[k-1]+swapc, i ); if flag and themarkisinview( crosspoin[k]-swapc, i ): draw rp(crosspoin[k-1])--rp(crosspoin[k]); fi; endfor; % if LabelCrossPoints: % for k=1 upto counter: % infolabel:=decimal(i)&","&decimal(j)&","&decimal(k); % infolabel := "0"; % dotlabelrand(infolabel,rp(crosspoin[k])); % endfor; % fi; endfor; endfor endgroup enddef; % Draw three-dimensional lines checking visibility. def lineraytrace(expr Press, Col) = begingroup save i, j, a, b; numeric i, j; color a, b; for i=1 upto NL: % scan all lines for j=1 upto npl[i]-1: a := L[i]p[j]; b := L[i]p[j+1]; coarse_line( a, b, 0, Press, Col); endfor; endfor endgroup enddef; % Draw only the faces, rigorously projecting the edges. def faceraytrace(expr Press, Col) = begingroup save i, j, a, b; numeric i, j; color a, b; for i=1 upto NF: % scan all faces for j=1 upto npf[i]: a := F[i]p[j]; if j <> npf[i]: b := F[i]p[j+1]; else: b := F[i]p1; fi; coarse_line( a, b, i, Press, Col); endfor; endfor endgroup enddef; % Fast test for your three-dimensional object def draw_all_test( expr AlsoDrawLines ) = begingroup save i, j, a, b; numeric i, j; color a, b; if ShadowOn: for i=1 upto NF: fill faceshadowpath( i ); endfor; if AlsoDrawLines: for i=1 upto NL: % scan all lines for j=1 upto npl[i]-1: a := L[i]p[j]; b := L[i]p[j+1]; drawsegment( cb(a), cb(b) ); endfor; endfor; fi; fi; for i=1 upto NF: % scan all faces for j=1 upto npf[i]: a := F[i]p[j]; if j <> npf[i]: b := F[i]p[j+1]; else: b := F[i]p1; fi; drawsegment( a, b ); endfor; endfor; if AlsoDrawLines: for i=1 upto NL: % scan all lines for j=1 upto npl[i]-1: a := L[i]p[j]; b := L[i]p[j+1]; drawsegment( a, b ); endfor; endfor; fi endgroup enddef; % Don't use SphericalDistortion here. def fill_faces( text LineAtribs ) = begingroup save i; numeric i; if ShadowOn: for i=1 upto NF: fill faceshadowpath( i ); endfor; fi; for i=1 upto NF: face_invisible( i, LineAtribs ); endfor endgroup enddef; def doitnow = begingroup save i, j, farone; numeric i, j, farone[]; for i=1 upto Nobjects: farone[i] := i; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shell's Method of sorting %%%%%%%% save inc, v, vv; numeric inc, v, vv; inc = 1; forever: inc := 3*inc+1; exitif inc > Nobjects; endfor; forever: inc := round(inc/3); for i=inc+1 upto Nobjects: v := RefDist[i]; vv:= farone[i]; j := i; forever: exitunless RefDist[j-inc] > v; RefDist[j] := RefDist[j-inc]; farone[j] := farone[j-inc]; j := j-inc; exitif j <= inc; endfor; RefDist[j] := v; farone[j] := vv; endfor; exitunless inc > 1; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=Nobjects downto 1: j := farone[i]; scantokens ostr[j]; endfor endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Nematic Liquid Crystal wise: def generatedirline(expr Lin, Phi, Theta, Long, Currpos ) = begingroup save longvec; color longvec; npl[Lin] := 2; longvec := Long*( cosd(Phi)*cosd(Theta), sind(Phi)*cosd(Theta), sind(Theta) ); L[Lin]p1 := Currpos-0.5*longvec; L[Lin]p2 := Currpos+0.5*longvec endgroup enddef; def generatedirface(expr Fen, Phi, Theta, Long, Base, Currpos ) = begingroup save basevec, longvec; color basevec, longvec; npf[Fen] := 3; longvec := Long*( cosd(Phi)*cosd(Theta), sind(Phi)*cosd(Theta), sind(Theta) ); basevec := Base*ncrossprod( Currpos-f, longvec ); F[Fen]p1 := Currpos-0.5*(longvec+basevec); F[Fen]p2 := Currpos+0.5*longvec; F[Fen]p3 := Currpos-0.5*(longvec-basevec) endgroup enddef; def generateonebiax(expr Lin, Phi, Theta, Long, SndDirAngl, Base, Currpos ) = begingroup save basevec, longvec, u, v; color basevec, longvec, u, v; npl[Lin] := 4; longvec := Long*( cosd(Phi)*cosd(Theta), sind(Phi)*cosd(Theta), sind(Theta) ); v = (-sind(Phi), cosd(Phi), 0); u = ( cosd(Phi)*cosd(Theta+90), sind(Phi)*cosd(Theta+90), sind(Theta+90) ); basevec := Base*( v*cosd(SndDirAngl)+u*sind(SndDirAngl) ); L[Lin]p1 := Currpos-0.5*longvec; L[Lin]p2 := Currpos+0.5*basevec; L[Lin]p3 := Currpos+0.5*longvec; L[Lin]p4 := Currpos-0.5*basevec endgroup enddef; def director_invisible( expr SortEmAll, ThickenFactor, CyclicLines ) = begingroup save i, j, k, farone, thisfar; save outerr, innerr, direc, ounum; numeric i, j, k, farone[], dist[], thisfar, ounum; pen actualpen, outerr, innerr; path direc; actualpen = currentpen; if SortEmAll: for i=1 upto NL: % scan all lines dist[i] := conorm( masscenter( npl[i], L[i]p ) - f ); farone[i] := i; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shell's Method of sorting %%%%%%%% save inc, v, vv; numeric inc, v, vv; inc = 1; forever: inc := 3*inc+1; exitif inc > NL; endfor; forever: inc := round(inc/3); for i=inc+1 upto NL: v := dist[i]; vv:= farone[i]; j := i; forever: exitunless dist[j-inc] > v; dist[j] := dist[j-inc]; farone[j] := farone[j-inc]; j := j-inc; exitif j <= inc; endfor; dist[j] := v; farone[j] := vv; endfor; exitunless inc > 1; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% else: for i=1 upto NL: farone[i] := i; endfor; fi; for i=NL downto 1: % draw all pathes j := farone[i]; direc := rp( L[j]p1 ) for k=2 upto npl[j]: --rp( L[j]p[k] ) endfor; if CyclicLines: direc := direc--cycle; fi; ounum := Spread*ps( masscenter(npl[j],L[j]p), ThickenFactor ); outerr := pencircle scaled ounum; innerr := pencircle scaled (0.8*ounum); %% DANGER %% pickup outerr; draw direc withcolor black; pickup innerr; draw direc withcolor background; endfor; pickup actualpen 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: % Define and draw surfaces with a triangular mesh. % On a hexagonal or triangular area. Without sorting (no need). def hexagonaltrimesh( expr BeHexa,theN,SideSize )( text SurFunc ) = begingroup save i, j, posx, posy, posz, higx, higy, counter, stepx, stepy, poi, lowx, lowy, newn, bola, bolb, bolc; numeric i, j, posx, posy, posz, higx, higy, counter, stepx, stepy, lowx, lowy, newn; color poi[][]; boolean bola, bolb, bolc; npf0 := 3; FCD0 := true; % this is used in the calls to fillfacewithlight ActuC := incr( ActuC ); if ActuC > TableColors: ActuC := 1; fi; FC0 := ActuC; %% counter = 0; stepy = SideSize/theN; stepx = 0.5*stepy*sqrt(3); lowy = -0.5*SideSize; lowx = -sqrt(3)*SideSize/6; higy = -lowy; higx = sqrt(3)*SideSize/3; for i=0 upto theN: for j=0 upto theN-i: posx := lowx + i*stepx; posy := lowy + i*stepx/sqrt(3) + j*stepy; posz := SurFunc( posx, posy ); poi[i][j] := ( posx, posy, posz ); endfor; endfor; if BeHexa: newn = round((theN+1)/3)+1; else: newn = 1; fi; for j=newn upto theN-newn+1: F0p1 := poi[0][j-1]; F0p2 := poi[0][j]; F0p3 := poi[1][j-1]; fillfacewithlight( 0 ); % see below endfor; for i=1 upto theN-1: for j=1 upto theN-i: bola := ( i < newn ) and ( j < newn-i ); bolb := ( i < newn ) and ( j > theN-newn+1 ); bolc := ( i > theN-newn ); if not ( bola or bolb or bolc ): F0p1 := poi[i-1][j]; F0p2 := poi[i][j-1]; F0p3 := poi[i][j]; fillfacewithlight( 0 ); F0p1 := poi[i+1][j-1]; fillfacewithlight( 0 ); fi; endfor; endfor; i := theN-newn+1; for j=1 upto newn-1: F0p1 := poi[i-1][j]; F0p2 := poi[i][j-1]; F0p3 := poi[i][j]; fillfacewithlight( 0 ); endfor; endgroup enddef; def fillfacewithlight( expr FaceN ) = begingroup save perpvec, reflectio, viewvec, inciden, refpos, projincid; save shiftv, fcol, lcol, theangle, ghost, pa, pb, pc, j, lowcolor; color perpvec, reflectio, viewvec, inciden, refpos, projincid; color shiftv, fcol, lcol, pa, pb, pc, lowcolor; numeric theangle, j; path ghost; ghost := rp( F[FaceN]p1 ) for j=2 upto npf[FaceN]: --rp( F[FaceN]p[j] ) endfor --cycle; if OverRidePolyhedricColor: unfill ghost; else: refpos = masscenter( npf[FaceN], F[FaceN]p ); pa = F[FaceN]p1; pb = F[FaceN]p2; pc = F[FaceN]p3; inciden = LightSource - refpos; viewvec = f - refpos; perpvec = ncrossprod( pa-pb, pb-pc ); if cdotprod( perpvec, blue ) < 0: perpvec := -perpvec; fi; projincid = perpvec*cdotprod( perpvec, inciden ); shiftv = inciden - projincid; reflectio = projincid - shiftv; theangle = getangle( reflectio, viewvec ); if FCD[FaceN]: lowcolor = TableC[FC[FaceN]]; else: lowcolor = TableC0; fi; lcol = (cosd( theangle ))[lowcolor,HigColor]; if cdotprod( viewvec, perpvec ) < 0: fcol = lcol - SubColor; else: fcol = lcol; fi; fill ghost withcolor fcol; fi; draw ghost; endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Part III (parametric plots and another renderer): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% Kindly contributed by Jens Schwaiger %%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % dmin_ is the minimal distance of all faces to f; % dmax_ is the maximal distance of all faces to f; % (both values are determined in "draw_invisible") % Facen is the number of the face to be filled; % ColAtrib is the color used for filling; % ColAtribone the color used for drawing; % Colordensity depends on distance of the face from f def face_drawfill( expr Facen, dmin_, dmax_ ,ColAtrib, ColAtribone ) = begingroup save j, ptmp, colfac_, coltmp_; path ghost; numeric j, colfac_; color ptmp, coltmp_; ghost := rp( F[Facen]p1 ); for j=2 upto npf[Facen]: ghost := ghost--rp( F[Facen]p[j] ); endfor; ghost := ghost--cycle; ptmp:= masscenter( npf[Facen], F[Facen]p ) - f; % 0<=colfac_<=1 colfac_ := ( conorm(ptmp)-dmin_ )/( dmax_ - dmin_ ); % color should be brighter, if distance to f is smaller colfac_ := 1 - colfac_; % color should be not to dark, i.e., >=0.1 (and <=1) colfac_ := 0.9colfac_ + 0.1; % now filling and drawing with appropriate color; fill ghost withcolor colfac_*ColAtrib; draw ghost withcolor colfac_*ColAtribone; endgroup enddef; % Now a much faster faces-only-ray-tracer based upon the unfill % command and the constraint of non-intersecting faces of similar % sizes. Faces are sorted by distance from nearest vertex or % masscenter to the point of view. This routine may be used to % draw graphs of 3D surfaces (use masscenters in this case). % % Option=true: test all vertices % Option=false: test masscenters of faces % DoJS=true: use face_drawfill by J. Schwaiger % DoJS=false: use fillfacewithlight by L. Nobre G. % ColAtrib=color for filling faces % ColAtribone=color for drawing edges 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_; color a, b, ptmp; for i=1 upto NF: % scan all faces if Option: % for distances of dist[i] = conorm( F[i]p1 - f ); % nearest vertices if i=1: distmin_ := dist1; % initialisation of distmax_ := dist1; % dmin_ and dmax_ fi; distmin_ := min( distmin_, dist[i] ); distmax_ := max( distmax_, dist[i] ); for j=2 upto npf[i]: thisfar := conorm( F[i]p[j] - f ); distmin_ := min( distmin_, thisfar ); distmax_ := max( distmax_, thisfar ); if thisfar < dist[i]: dist[i] := thisfar; fi; endfor; else: % for distances of centers of mass dist[i] := conorm( masscenter( npf[i], F[i]p ) - f ); if i=1: distmin_ := dist1; % initialisation of distmax_ := dist1; % dmin_ and dmax_ in this case fi; distmin_ := min( distmin_, dist[i] ); distmax_ := max( distmax_, dist[i] ); fi; farone[i] = i; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Shell's Method of sorting %%%%%%%% save inc, v, vv; numeric inc, v, vv; inc = 1; forever: inc := 3*inc+1; exitif inc > NF; endfor; forever: inc := round(inc/3); for i=inc+1 upto NF: v := dist[i]; vv:= farone[i]; j := i; forever: exitunless dist[j-inc] > v; dist[j] := dist[j-inc]; farone[j] := farone[j-inc]; j := j-inc; exitif j <= inc; endfor; dist[j] := v; farone[j] := vv; endfor; exitunless inc > 1; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=NF downto 1: % draw and fill all pathes j := farone[i]; if DoJS: face_drawfill( j, distmin_, distmax_, ColAtrib, ColAtribone ); else: fillfacewithlight( j ); fi; endfor; endgroup enddef; % Move to the good range (-1,1). def bracket( expr low, poi, hig ) = begingroup save zout; numeric zout; zout = (2*poi-hig-low)/(hig-low); if zout > 1: zout := 1; fi; if zout < -1: zout := -1; fi; ( zout ) endgroup enddef; % 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 ) = begingroup save i, j, k, posx, posy, posz; save counter, stept, steps, poss, post, tmpaux; save veca, vecb, vecc, vecd; numeric i, j, k, posx, posy, posz, counter, stept, steps; color poi[][], tmpaux, veca, vecb, vecc, vecd; counter := NF; % <-- NF must be initialized! ActuC := incr( ActuC ); if ActuC > TableColors: ActuC := 1; fi; steps = ( higs - lows )/ns; stept = ( higt - lowt )/nt; for i=0 upto ns: for j=0 upto nt: poss := lows + i*steps; post := lowt + j*stept; tmpaux := parSurFunc( poss, post ); posz := Z(tmpaux); posx := X(tmpaux); posy := Y(tmpaux); posx := bracket(lowx,posx,higx); posy := bracket(lowy,posy,higy); posz := bracket(lowz,posz,higz)/facz; poi[i][j] := ( posx, posy, posz ); endfor; endfor; for i=1 upto ns: for j=1 step 1 until nt: veca := poi[i][j]-poi[i-1][j]; vecb := poi[i][j]-poi[i-1][j-1]; vecc := poi[i][j]-poi[i][j-1]; if abs(cdotprod(ccrossprod(veca,vecb),vecc))<0.00005: %DANGER! counter := incr(counter); npf[counter] := 4; F[counter]p1 := poi[i-1][j-1]; F[counter]p2 := poi[i-1][j]; F[counter]p3 := poi[i][j]; F[counter]p4 := poi[i][j-1]; FC[counter] := ActuC; FCD[counter] := true; else: tmpaux:= 0.25*(poi[i-1][j-1]+poi[i-1][j]+poi[i][j]+poi[i][j-1]); veca := poi[i-1][j-1]-tmpaux; vecb := poi[i-1][j]-tmpaux; vecc := poi[i][j]-tmpaux; vecd := poi[i][j-1]-tmpaux; if getangle(vecb,vecd)>getangle(veca,vecc): for k=-1 upto 0: counter := incr(counter); npf[counter] := 3; F[counter]p1 := poi[i-1][j-1]; F[counter]p2 := poi[i+k][j-1-k]; F[counter]p3 := poi[i][j]; FC[counter] := ActuC; FCD[counter] := true; endfor; else: for k=-1 upto 0: counter := incr(counter); npf[counter] := 3; F[counter]p1 := poi[i+k][j+k]; F[counter]p2 := poi[i][j-1]; F[counter]p3 := poi[i-1][j]; FC[counter] := ActuC; FCD[counter] := true; endfor; fi; fi; endfor; endfor; NF := counter; endgroup enddef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Part IV (automatic perspective tuning and minimization): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% def randomfear = begingroup save a, b, c, i, h; numeric a, b, c, i, h; i = uniformdeviate( 360 ); h = uniformdeviate( 180 ); a = cosd( h )*cosd( i ); b = sind( h )*cosd( i ); c = sind( i ); ( (a,b,c) ) endgroup enddef; def renormalizevc( expr inF, inVC ) = begingroup save a; color a; cdotprod( inF, a ) = 0; inF - a = whatever*( inF - inVC ); ( a ) endgroup enddef; def calculatecost( expr TryF, TryVc, TrySp, TrySh ) = begingroup save sumsquares, i, difpair, xx, yy; numeric sumsquares, i, xx, yy; pair difpair; sumsquares = 0; f := TryF; viewcentr := TryVc; Spread := TrySp; ShiftV := TrySh; for i=1 upto PhotoMarks: difpair := 0.05*(rp(PhotoPoint[i]) - PhotoPair[i]); % DANGER xx := ((xpart difpair)**2)/PhotoMarks; yy := ((ypart difpair)**2)/PhotoMarks; sumsquares := sumsquares + xx + yy; endfor; ( sumsquares ) endgroup enddef; def forcepointinsidefear( text A ) = begingroup if abs( X(A) ) > 0.5*MaxFearLimit : A := 0.5*A*MaxFearLimit/abs( X(A) ); fi; if abs( Y(A) ) > 0.5*MaxFearLimit : A := 0.5*A*MaxFearLimit/abs( Y(A) ); fi; if abs( Z(A) ) > 0.5*MaxFearLimit : A := 0.5*A*MaxFearLimit/abs( Z(A) ); fi; ( A ) endgroup enddef; def forcepairinsidepage( text A ) = begingroup if xpart A > 2*(xpart OriginProjPagePos): A := ( 2*(xpart OriginProjPagePos), ypart A ); fi; if ypart A > 2*(ypart OriginProjPagePos): A := ( xpart A, 2*(ypart OriginProjPagePos) ); fi; if xpart A < 0: A := ( 0, ypart A ); fi; if ypart A < 0: A := ( xpart A, 0 ); fi; ( A ) endgroup enddef; def calculatejump( expr AverCost, PrevCost, RandCost, JumpLimit ) = begingroup save funfact, numer, denom; numeric funfact, numer, denom; if RandCost+PrevCost > 2*AverCost: numer := 3*PrevCost - 4*AverCost + RandCost; denom := 4*( RandCost + PrevCost - 2*AverCost ); if abs( denom )*JumpLimit < abs( numer ): funfact := 0; else: funfact := numer/denom; fi; else: funfact := 0; fi; ( funfact ) endgroup enddef; def photoreverse( expr IterNum, ExpTao, JumpFact ) = begingroup save j, auxvc, actfact, auxfac; save prevf, randf, averf, prevvc, randvc, avervc; save prevsp, randsp, aversp, prevsh, randsh, aversh; save prevcost, randcost, avercost, spreadlimit, expfact; numeric j, prevsp, randsp, aversp, actfact; numeric prevcost, randcost, avercost; numeric spreadlimit, expfact; color prevf, randf, averf, prevvc, randvc, avervc, auxvc; pair prevsh, randsh, aversh; spreadlimit = 20; % DANGER prevf = f; prevvc= viewcentr; prevsp= Spread; prevsh= ShiftV; prevcost = calculatecost( prevf, prevvc, prevsp, prevsh ); show prevcost; auxfac = -1280.0/ExpTao/IterNum; for j=0 upto IterNum: expfact := mexp(auxfac*j); randf := prevf + expfact*randomfear; randcost := calculatecost( randf, prevvc, prevsp, prevsh ); averf := 0.5[prevf,randf]; avercost := calculatecost( averf, prevvc, prevsp, prevsh ); actfact := calculatejump(avercost,prevcost,randcost,JumpFact); prevf := actfact[prevf,randf]; auxvc := prevvc + expfact*randomfear; randvc:= renormalizevc( randf, auxvc ); randcost := calculatecost( prevf, randvc, prevsp, prevsh ); auxvc := 0.5[prevvc,randvc]; avervc:= renormalizevc( averf, auxvc ); avercost := calculatecost( prevf, avervc, prevsp, prevsh ); actfact := calculatejump(avercost,prevcost,randcost,JumpFact); auxvc := actfact[prevvc,randvc]; prevvc:= renormalizevc( prevf, auxvc ); randsp:= prevsp+expfact*spreadlimit*(uniformdeviate( 1 ) - 0.5); randcost := calculatecost( prevf, prevvc, randsp, prevsh ); aversp:= 0.5[prevsp,randsp]; avercost := calculatecost( prevf, prevvc, aversp, prevsh ); actfact := calculatejump(avercost,prevcost,randcost,JumpFact); prevsp:= actfact[prevsp,randsp]; randsh:= prevsh + expfact*dir( uniformdeviate( 360 ) ); % DANGER randcost := calculatecost( prevf, prevvc, prevsp, randsh ); aversh:= 0.5[prevsh,randsh]; avercost := calculatecost( prevf, prevvc, prevsp, aversh ); actfact := calculatejump(avercost,prevcost,randcost,JumpFact); prevsh:= actfact[prevsh,randsh]; prevcost := calculatecost( prevf, prevvc, prevsp, prevsh ); %show (prevcost,expfact); endfor; show prevcost; endgroup enddef; % Minimization routine for scalar functions like y=f(x) where an initial % triplet (x1,x2,x3) with x1ya) 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 ((xaRay: 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Part VI (strictly two-dimensional): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Verify if a path is cyclic (written by Scott Pakin) def is_cyclic expr cpath = (point 0 of cpath = point (length cpath) of cpath) enddef; % Produce the schematics of a spring. def springpath( expr begp, endp, piturnum, piturnproj, spgfrac ) = begingroup boolean leftside; numeric counter, springwidth; pair stepdir, leftdir, rightdir, auxil, lastp, endvec; path thespring; leftside = true; stepdir = spgfrac*(endp-begp)/piturnum; endvec = (1-spgfrac)*(endp-begp)/2; springwidth = piturnproj +-+ abs( stepdir ); auxil = ( -ypart stepdir, xpart stepdir ); leftdir = springwidth*unitvector( auxil ); auxil := ( ypart stepdir, -xpart stepdir ); rightdir = springwidth*unitvector( auxil ); leftdir := (stepdir + leftdir)/2; rightdir := (stepdir + rightdir)/2; lastp = begp+endvec; thespring := begp--lastp; for counter=1 upto piturnum: if leftside: auxil := lastp + leftdir; else: auxil := lastp + rightdir; fi; lastp := lastp + stepdir; thespring := thespring--auxil--lastp; leftside := not leftside; endfor; thespring := thespring--endp; ( thespring ) endgroup enddef; % Summarize a great length in a zig-zag frontier line def zigzagfrontier( expr begp, endp, nzigs, dev, zthick, tthick, fthick, excol, incol ) = begingroup interim linecap := squared; interim linejoin := mitered; boolean leftside; numeric counter, tmpval; pair stepdir, leftdir, rightdir, auxil, lastp; path thefrontier; leftside = true; stepdir = (endp-begp)/nzigs; leftdir = unitvector( ( -ypart stepdir, xpart stepdir ) ); rightdir = unitvector( ( ypart stepdir, -xpart stepdir ) ); thefrontier = begp; lastp = begp; for counter=1 upto nzigs-1: lastp := lastp + stepdir; tmpval := zthick+normaldeviate*dev; if leftside: auxil := lastp + leftdir*tmpval; else: auxil := lastp + rightdir*tmpval; fi; thefrontier := thefrontier--auxil; leftside := not leftside; endfor; thefrontier := thefrontier--endp; draw thefrontier withcolor excol withpen pencircle scaled tthick; draw thefrontier withcolor incol withpen pencircle scaled fthick; ( thefrontier ) endgroup enddef; % The name says it all. def randomcirc( expr radi, stddev, numpois ) = begingroup numeric i, astep; path ranc; astep = 360/numpois; ranc = (radi+stddev*normaldeviate)*dir(-180); for i= -180+astep step astep until 180: ranc := ranc--((radi+stddev*normaldeviate)*dir(i)); endfor; ranc := ranc--cycle; ( ranc ) endgroup enddef; % Just label some 2D-place in a way similar to anglinen. def labeln(expr S, Pos, RelPos) = begingroup if RelPos = 0: label.rt( S, Pos ); elseif RelPos =1: label.urt( S, Pos ); elseif RelPos =2: label.top( S, Pos ); elseif RelPos =3: label.ulft( S, Pos ); elseif RelPos =4: label.lft( S, Pos ); elseif RelPos =5: label.llft( S, Pos ); elseif RelPos =6: label.bot( S, Pos ); elseif RelPos =7: label.lrt( S, Pos ); else: label( S, Pos ); fi endgroup enddef; % There must be a sort of planar cross-product. The z coordinate. def paircrossprod(expr A, B) = ( (xpart A)*(ypart B) - (xpart B)*(ypart A) ) enddef; % Just dot-label some 2D-place in a random way. def dotlabelrand(expr S, Pos ) = begingroup save RelPos; numeric RelPos; RelPos = floor( uniformdeviate(9) ); if RelPos = 0: dotlabel.rt( S, Pos ); elseif RelPos =1: dotlabel.urt( S, Pos ); elseif RelPos =2: dotlabel.top( S, Pos ); elseif RelPos =3: dotlabel.ulft( S, Pos ); elseif RelPos =4: dotlabel.lft( S, Pos ); elseif RelPos =5: dotlabel.llft( S, Pos ); elseif RelPos =6: dotlabel.bot( S, Pos ); elseif RelPos =7: dotlabel.lrt( S, Pos ); else: dotlabel( S, Pos ); fi endgroup enddef; 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: 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; 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: hide( va := unitvector( point i of Path - Point ); 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. 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+0.5*Beam*dir(aang); %%% DANGER %%%%%0.5%%%%%%% %draw b withpen pencircle scaled 3pt withcolor green; anew := bpoi+Beam*dir(90+bang); c := anew-0.5*Beam*dir(bang); %%% DANGER %%0.5%%%%%%% %draw c withpen pencircle scaled 3pt withcolor blue; 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; 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 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; % EOF