% featpost3Dplus2D.mp
% L. Nobre G.,
[email protected],
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<lena) and (actang>=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<a:
sp = 1 +-+ cp;
else:
sp = -( 1 +-+ cp );
fi;
bdm = (k*cm,sm)*Ray;
bdp = (k*cp,sp)*Ray;
bdh = 0.5[bdp,bdm];
tmpc := Ray*( 1 +-+ ((xpart bdh)/a) );
tmpb := tmpc +-+ (ypart bdh);
bdv = bdm-bdp;
wax = 0.5*( xax*( xpart (bdv) ) + yax*( ypart (bdv) ) );
cac = CentrPoi+ xax*( xpart (bdh) ) + yax*( ypart (bdh) );
fill ellipticshadowpath( cac, wax, zax*tmpb );
endgroup
enddef;
def spheroid( expr CentrPoi, NorthPoleVec, Ray ) =
begingroup
save a, k, fx, fy, tmpa, tmpb, tmpc, ep;
save 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;
if ShadowOn:
spheroidshadow( CentrPoi, NorthPoleVec, Ray );
fi;
vax = f-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( (-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<a:
sp = 1 +-+ cp;
else:
sp = -( 1 +-+ cp );
fi;
bdm = (k*cm,sm)*Ray;
bdp = (k*cp,sp)*Ray;
bdh = 0.5[bdp,bdm];
tmpc := Ray*( 1 +-+ ((xpart bdh)/a) );
tmpb := tmpc +-+ (ypart bdh);
bdv = bdm-bdp;
wax = 0.5*( xax*( xpart (bdv) ) + yax*( ypart (bdv) ) );
cac = CentrPoi+ xax*( xpart (bdh) ) + yax*( ypart (bdh) );
ep = ellipticpath( cac, wax, zax*tmpb );
unfill ep;
draw ep;
draw spatialhalfcircle( CentrPoi, NorthPoleVec, Ray, true );
endgroup
enddef;
% Another brute-force algorythm.
% It's advisable to use three orthogonal axes.
def ellipsoid( expr Centr, AxOne, AxTwo, AxThr ) =
begingroup
save count, i, j, axx, axy, cyc, cy, di, leng;
numeric count, i, j, leng;
color axx, axy, di[];
path cy, cyc[];
di1 = AxOne;
di2 = AxTwo;
di3 = AxThr;
count = 0;
for i=1 upto 3:
if i=1:
axx := di2;
axy := di3;
elseif i=2:
axx := di1;
axy := di3;
else:
axx := di1;
axy := di2;
fi;
for j=5 step 10 until 175:
cyc[incr(count)] = ellipticpath( Centr, di[i],
Centr + planarrotation( axx, axy, j ) );
endfor;
endfor;
cy = cyc1;
for i=2 upto count-1:
cy := twocyclestogether( cy, cyc[i] );
endfor;
leng = (length cy) - 1;
( point 0 of cy for i = 1 upto leng: ..point i of cy endfor ..cycle )
endgroup
enddef;
def ellipsoidshadow( expr Centr, AxOne, AxTwo, AxThr ) =
begingroup
save count, i, j, axx, axy, cyc, cy, di, leng;
numeric count, i, j, leng;
color axx, axy, di[];
path cy, cyc[];
di1 = AxOne;
di2 = AxTwo;
di3 = AxThr;
count = 0;
for i=1 upto 3:
if i=1:
axx := di2;
axy := di3;
elseif i=2:
axx := di1;
axy := di3;
else:
axx := di1;
axy := di2;
fi;
for j=5 step 10 until 175:
cyc[incr(count)] = ellipticshadowpath( Centr, di[i],
Centr + planarrotation( axx, axy, j ) );
endfor;
endfor;
cy = cyc1;
for i=2 upto count-1:
cy := twocyclestogether( cy, cyc[i] );
endfor;
leng = (length cy) - 1;
( point 0 of cy for i = 1 upto leng: ..point i of cy endfor ..cycle )
endgroup
enddef;
def revolparab( expr BaseCenter, ParabTip, BaseRay ) =
begingroup
save bcpt, conetip, fakex, fakey, fakez, tipview, ellicenter, coneview;
save tanefe, majorvec, minorvec, cutvec, auxpoi, crux, auxx, baseview;
save xzero, yzero, xdelta, fakea, xpos, ypos, xneg, yneg, l, apertur;
save auxy, maxy, ellmaxang, ymin, xefe, auxray, auxcos, auxsin, ste, a;
save auxpath, tippath;
save conda, condb, condc;
color bcpt, conetip, fakex, fakey, fakez, tipview, ellicenter, coneview;
color tanefe, majorvec, minorvec, cutvec, auxpoi, crux, auxx, baseview;
numeric xzero, yzero, xdelta, fakea, xpos, ypos, xneg, yneg, l, apertur;
numeric auxy, maxy, ellmaxang, ymin, xefe, auxray, auxcos, auxsin, ste;
numeric a;
path auxpath, tippath;
boolean conda, condb, condc;
bcpt = BaseCenter-ParabTip;
conetip = BaseCenter-2*bcpt;
maxy = conorm(bcpt);
apertur = angle(2*maxy,BaseRay);
fakey = N(bcpt);
coneview = f-conetip;
tipview = f-ParabTip;
baseview = f-BaseCenter;
a = getangle(coneview,bcpt);
conda = a<=apertur;
condb = a>=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 yneg<ypos:
ymin=yneg;
xefe=xneg;
else:
ymin=ypos;
xefe=xpos;
fi;
tanefe = ParabTip +fakex*xefe+fakey*ymin;
minorvec = xdelta*fakez;
majorvec = tanefe-ellicenter;
if (yneg<maxy) and (ypos<maxy):
draw ellipticpath( ellicenter, majorvec, minorvec );
else:
auxpoi = lineintersectplan(conetip,conetip-f,BaseCenter,bcpt);
auxray = conorm( BaseCenter-auxpoi );
auxcos = BaseRay/auxray;
auxsin = 1 +-+ auxcos;
if cdotprod(fakex,auxpoi-BaseCenter)>0:
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 (dist<Bray+dray) and (dist>Bray-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 x1<x2<x3 is given as a parabolic squeleton that
% provides a way to search for the smallest value of y (if iterated)
def minimizestep( expr Abcisscolor )( text PlainFunc ) =
begingroup
save xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
save colout;
numeric xa, xb, xc, xd, ya, yb, yc, yd, aux, coeb, coec, den;
color colout;
xa = X( Abcisscolor );
xb = Y( Abcisscolor );
xc = Z( Abcisscolor );
ya = PlainFunc(xa);
yb = PlainFunc(xb);
yc = PlainFunc(xc);
if ya = yb:
colout = (-0.125[xa,xb],xb,xc);
elseif yb = yc:
colout = (xa,xb,1.125[xb,xc]);
else:
if (yb>ya) or (yb>yc):
show Abcisscolor;
message " Unable to minimizestep!";
fi;
den = (xb-xc)*((xa**2)-(xb**2))-(xa-xb)*((xb**2)-(xc**2));
if abs(den) < 0.0005:
show den;
message " Unable to minimizestep!";
fi;
coeb = ((yb-yc)*((xa**2)-(xb**2))-(ya-yb)*((xb**2)-(xc**2)))/den;
coec = ((xb-xc)*(ya-yb)-(xa-xb)*(yb-yc))/den;
xd = -0.5*coeb/coec;
yd = PlainFunc( xd );
if ((xa<xd) and (xd<xb)):
if (yd<yb):
colout = (xa,xd,xb);
else:
colout = (xd,xb,xc);
fi;
elseif ((xb<xd) and (xd<xc)):
if (yd<yb):
colout = (xb,xd,xc);
else:
colout = (xa,xb,xd);
fi;
else:
aux := 0.125[xb,xc]-0.125[xb,xa];
colout = (xa,0.125[xb,xa]+uniformdeviate(aux),xc);
fi;
fi;
( colout )
endgroup
enddef;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Part V:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Intersections:
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Check lineintersectplan above %%%%%%%%%%%%%%%
def calculatecostver(expr VerA,DisA,VerB,DisB,VerC,DisC,TryV) =
begingroup
save a, b, c;
numeric a, b, c;
a = conorm( TryV - VerA ) - DisA;
b = conorm( TryV - VerB ) - DisB;
c = conorm( TryV - VerC ) - DisC;
( ( a ++ b ++ c )**2 )
endgroup
enddef;
% Approximation of the intersection of three spheres.
% Be aware of the next three danger parameters.
def improvertex( expr VerA, DisA, VerB, DisB, VerC, DisC, IniV ) =
begingroup
save j, iternum, prevcost, factry, actpoi, actfact, auxa, auxb;
save jumplim, trypoi, halfpo, expfact, randcost, avercost, auxc;
numeric j, iternum, prevcost, factry, randcost, avercost, actfact;
numeric jumplim, expfact;
color actpoi, trypoi, halfpo, auxa, auxb, auxc;
trypoi = IniV;
iternum = 24; % DANGER!
factry = 0.25; % DANGER!
jumplim = 2.5; % DANGER!
prevcost = calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
message "improvertex initial cost: " & decimal(prevcost)
& " at point " & cstr(trypoi);
for j=1 upto 5:
auxa := VerA-trypoi;
auxa := N(auxa)*(conorm(auxa)-DisA);
auxb := VerB-trypoi;
auxb := N(auxb)*(conorm(auxb)-DisB);
auxc := VerC-trypoi;
auxc := N(auxc)*(conorm(auxc)-DisC);
trypoi := trypoi+factry*(auxa+auxb+auxc);
endfor;
prevcost := calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
message "phase one cost: " & decimal(prevcost)
& " at point " & cstr(trypoi);
for j=0 upto iternum:
expfact := mexp(-j*200/iternum); % DANGER!
actpoi := trypoi+factry*expfact*randomfear;
randcost :=
calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,actpoi);
halfpo := 0.5[trypoi,actpoi];
avercost :=
calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,halfpo);
actfact := calculatejump(avercost,prevcost,randcost,jumplim);
trypoi := actfact[trypoi,actpoi];
prevcost :=
calculatecostver(VerA,DisA,VerB,DisB,VerC,DisC,trypoi);
% show prevcost;
endfor;
message "final cost: " & decimal(prevcost)
& " at point " & cstr(trypoi);
( trypoi )
endgroup;
enddef;
% The approximation of the intersection of a plane, an infinite cylinder (tube)
% and a spheroid (an oblate spheroid, a prolate spheroid or a perfect sphere)
def ultraimprovertex( expr PlanPoi, PlanDir, BaseCenter, Radius, LenVec,
CentrPoi, NorthPoleVec, Ray, IniV ) =
begingroup
save trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, factry, focn, focs;
save norpoldi, major, majortmp, nordi, tmpdi;
color trypoi, auxa, auxb, auxc, auxd, plandi, cyldi, focn, focs;
color norpoldi, nordi, tmpdi;
numeric factry, major, majortmp;
trypoi = IniV;
factry = 0.25;
plandi = N(PlanDir);
cyldi = N(LenVec);
norpoldi = N(NorthPoleVec);
major = conorm(NorthPoleVec);
if major>Ray:
focn := CentrPoi+(major+-+Ray)*norpoldi;
focs := CentrPoi-(major+-+Ray)*norpoldi;
fi;
for j=1 upto 50:
if major<Ray:
tmpdi := trypoi-CentrPoi;
nordi := N(tmpdi-cdotprod(tmpdi,norpoldi)*norpoldi);
focn := CentrPoi+(Ray+-+major)*nordi;
focs := CentrPoi-(Ray+-+major)*nordi;
elseif major=Ray:
focn := CentrPoi;
focs := CentrPoi;
fi;
auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
auxb := cyldi*cdotprod(trypoi-BaseCenter,cyldi)-trypoi+BaseCenter;
auxb := N(auxb)*(conorm(auxb)-Radius);
auxc := focn-trypoi;
auxd := focs-trypoi;
if major>Ray:
auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*major);
else:
auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*Ray);
fi;
trypoi := trypoi+factry*(auxa+auxb+auxc);
endfor;
( trypoi )
endgroup
enddef;
% The approximation of the intersection of a plane and two prolate spheroids
def necplusimprovertex( expr PlanPoi, PlanDir,
CentrPoiA, NorthPoleVecA, RayA,
CentrPoiB, NorthPoleVecB, RayB, IniV ) =
begingroup
save trypoi, auxa, auxb, auxc, auxd, plandi, factry, focni, focsi;
save norpoldi, norpoldj, maior, major, focnj, focsj, auxe;
color trypoi, auxa, auxb, auxc, auxd, plandi, focni, focsi, focnj, focsj;
color norpoldi, norpoldj, auxe;
numeric factry, maior, major;
trypoi = IniV;
factry = 0.25;
plandi = N(PlanDir);
norpoldi = N(NorthPoleVecA);
maior = conorm(NorthPoleVecA);
norpoldj = N(NorthPoleVecB);
major = conorm(NorthPoleVecB);
focni = CentrPoiA+(maior+-+RayA)*norpoldi;
focsi = CentrPoiA-(maior+-+RayA)*norpoldi;
focnj = CentrPoiB+(major+-+RayB)*norpoldj;
focsj = CentrPoiB-(major+-+RayB)*norpoldj;
for j=1 upto 50:
auxa := plandi*cdotprod(PlanPoi-trypoi,plandi);
auxb := focni-trypoi;
auxe := focsi-trypoi;
auxc := focnj-trypoi;
auxd := focsj-trypoi;
auxb := (N(auxb)+N(auxe))*(conorm(auxb)+conorm(auxe)-2*maior);
auxc := (N(auxc)+N(auxd))*(conorm(auxc)+conorm(auxd)-2*major);
trypoi := trypoi+factry*(auxa+auxb+auxc);
endfor;
( trypoi )
endgroup
enddef;
% The approximation of the intersection of a straight line and
% one prolate spheroid
def intersectprolatespheroid( expr CentrPoi, NorthPoleVec, Ray,
LinePoi, LineDir ) =
begingroup
save trypoi, factry, linedi, norpol, focn, focs, maior, j;
save auxa, auxb, auxc, auxd;
color trypoi, linedi, norpol, focn, focs;
color auxa, auxb, auxc, auxd;
numeric factry, maior, j;
trypoi = LinePoi;
factry = 0.25;
linedi = N(LineDir);
norpol = N(NorthPoleVec);
maior = conorm(NorthPoleVec);
focn = CentrPoi+(maior+-+Ray)*norpol;
focs = CentrPoi-(maior+-+Ray)*norpol;
for j=1 upto 50:
auxb := focn-trypoi;
auxd := focs-trypoi;
auxa := (N(auxb)+N(auxd))*(conorm(auxb)+conorm(auxd)-2*maior);
auxc := linedi*cdotprod(auxa,linedi);
trypoi := trypoi+factry*auxc;
endfor;
( trypoi )
endgroup
enddef;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 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
