% fiziko 0.2.0
% MetaPost library for physics textbook illustrations
% Copyright 2022 Sergey Slyusarev
%
% This program 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 program 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.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see .
% https://github.com/jemmybutton/fiziko
%
% Here we define some things of general interest
%
pi := 3.1415926;
radian := 180/pi;
vardef sin primary x = (sind(x*radian)) enddef;
vardef cos primary x = (cosd(x*radian)) enddef;
vardef log (expr n, b) =
save rv;
numeric rv;
if n > 0:
rv := (mlog(n)/mlog(b));
else:
rv := 0;
fi;
rv
enddef;
vardef arcsind primary x = angle((1+-+x,x)) enddef;
vardef arccosd primary x = angle((x,1+-+x)) enddef;
vardef arcsin primary x = ((arcsind(x))/radian) enddef;
vardef arccos primary x = ((arccosd(x))/radian) enddef;
vardef angleRad primary x = angle(x)/radian enddef;
vardef dirRad primary x = dir(x*radian) enddef;
% used here and there.
vardef sign (expr x)=
if x > 0: 1 fi
if x < 0: -1 fi
if x = 0: 1 fi
enddef;
% This is inverted `clip`
primarydef i maskedWith p =
begingroup
save q, invertedmask, resultimage, breakpoint;
pair q[];
path invertedmask;
picture resultimage;
resultimage := i;
q1 := ulcorner(i) shifted (-1, 1);
q3 := lrcorner(i) shifted (1, -1);
q2 := (xpart(q3), ypart(q1));
q4 := (xpart(q1), ypart(q3));
breakpoint := ypart((ulcorner(p)--llcorner(p)) firstIntersectionTimes p);
invertedmask := (subpath (breakpoint, length(p) + breakpoint) of p) -- q1 -- q2 -- q3 -- q4 -- q1 -- cycle;
clip resultimage to invertedmask;
resultimage
endgroup
enddef;
%
% Since metapost is somewhat unpredictable in determining where paths intersect, here's macro
% that returns first intersection times with first path (ray) priority.
% Actually, it is so in most cases, but sometimes second path can take precedence,
% so the macro just checks whether reversing 'q' changes something
%
primarydef p firstIntersectionTimes q =
begingroup
save t;
pair t[];
t1 := p intersectiontimes q;
t2 := p intersectiontimes reverse(q);
if xpart(t1) < xpart(t2):
t3 := t1;
else:
t3 := (xpart(t2), length(q) - ypart(t2));
fi;
if xpart(t1) < 0: t3 := t2; fi;
t3
endgroup
enddef;
% This checks if point a is inside of closed path p
primarydef a isInside p =
begingroup
save ang, v, i, rv, pp;
boolean rv;
pair pp[];
ang := 0;
for i := 0 step 1/4 until (length(p)):
pp1 := (point i of p) - a;
pp2 := (point i + 1/4 of p) - a;
if (pp1 <> (0, 0)) and (pp2 <> (0, 0)):
v := angle(pp1) - angle(pp2);
if v > 180: v := v - 360; fi; if v < -180: v := v + 360; fi;
ang := ang + v;
fi;
endfor;
if abs(ang) > 355:
rv := true;
else:
rv := false;
fi;
rv
endgroup
enddef;
% rotation in radians
primarydef somethingToRotate radRotated radAngle =
somethingToRotate rotated ((radAngle/pi)*180)
enddef;
%
% some 3D stuff
%
% this one's from byrne.mp
primarydef colorone dotprodXYZ colortwo =
begingroup
save xp, yp, zp;
numeric xp[], yp[], zp[];
xp1 := (redpart colorone);
yp1 := (greenpart colorone);
zp1 := (bluepart colorone);
xp2 := (redpart colortwo);
yp2 := (greenpart colortwo);
zp2 := (bluepart colortwo);
xp1*xp2 + yp1*yp2 + zp1*zp2
endgroup
enddef;
%
% sometimes it's useful to put some arrows along the path. this macro puts them
% in the middles of the segments that have length no less than midArrowLimit;
%
midArrowLimit := 1cm;
def drawmidarrow (expr p) text t =
begingroup
save i, j, q;
path q;
j := 0;
for i := 1 upto length(p):
if arclength(subpath(i-1, i) of p) >= midArrowLimit:
q := subpath(j, i - 1/2) of p;
j := i - 1/2;
draw q t;
filldraw arrowhead q t;
fi;
endfor;
draw subpath(j, length(p)) of p t;
endgroup
enddef;
% This macro marks angles, unsurprisingly
def markAngle (expr a, o, b) (text t) =
begingroup
save p, an, d;
numeric an[], d[];
pair p;
an1 := angle(a-o);
an2 := angle(b-o) - an1;
if (an2 < 0): an2 := an2 + 360; fi;
an3 := an1 + 1/2an2;
p := center(t);
d1 := abs(ulcorner(t)-lrcorner(t));
if (an2 < 90) and (an2 > 0):
d2 := max(1/3cm, (d1/(abs(sind(an2))*1/3cm))*1/3cm);
else:
d2 := 1/3cm;
fi;
draw subpath (0, 8an2/360) of fullcircle scaled 2d2 rotated an1 shifted o withpen thinpen;
draw (t) shifted -p shifted o shifted (dir(an3)*(d2 + d1));
endgroup
enddef;
%
% Here we define some auxilary global variables
%
% Offset path algorithm can subdivide original path in order to be more precise
offsetPathSteps := 4;
% The following macro sets all the values related to minimal stroke width at once.
% It can be used to easily redefine all of them.
def defineMinStrokeWidth (expr msw) =
% We don't want to display strokes that are too thin to print. Default value
% is subject to change when needed.
minStrokeWidth := msw;
maxShadingStrokeWidth := 3/2minStrokeWidth;
% At some point it's useless to display even dashes
minDashStrokeWidth := 1/3minStrokeWidth;
% this value corresponds to particular dashing algorithm and is subject to change whenever this algorithm changes
minDashStrokeLength := 3minStrokeWidth;
dashStrokeWidthStep := 1/5minDashStrokeWidth;
% all the shading algorithms need to know how close lines should be packed
shadingDensity := 3maxShadingStrokeWidth;
stippleSize := 3/2minStrokeWidth;
minStippleStep := 1/2stippleSize;
stippleShadingDensity := 3minStippleStep;
minStippleStrokeWidth := 1/20stippleSize;
% here are some pens
pen thinpen, thickpen, fatpen, stipplepen;
thinpen := pencircle scaled minStrokeWidth;
thickpen := pencircle scaled 3minStrokeWidth;
fatpen := pencircle scaled 6minStrokeWidth;
stipplepen := pencircle scaled stippleSize;
enddef;
defineMinStrokeWidth(1/5pt);
% here we set global light direction
def defineLightDirection (expr ldx, ldy) =
pair lightDirection;
color lightDirectionVectorXYZ;
lightDirection := (ldx, ldy);
lightDirectionVectorXYZ := (0, 0, 1);
lightDirectionVectorXYZ := rotateXYZaround.x(lightDirectionVectorXYZ, ldy);
lightDirectionVectorXYZ := rotateXYZaround.y(lightDirectionVectorXYZ, ldx);
enddef;
vardef rotateXYZaround@# (expr p, a) =
save partProj, rv;
pair partProj;
color rv;
if str @# = "x":
partProj := (greenpart(p), bluepart(p)) radRotated -a;
rv := (redpart(p), xpart(partProj), ypart(partProj));
elseif str @# = "y":
partProj := (redpart(p), bluepart(p)) radRotated -a;
rv := (xpart(partProj), greenpart(p), ypart(partProj));
elseif str @# = "z":
partProj := (redpart(p), greenpart(p)) radRotated -a;
rv := (xpart(partProj), ypart(partProj), bluepart(p));
else:
errmessage("What axis is " & str @# & "?");
fi;
rv
enddef;
defineLightDirection(-1/8pi, 1/8pi);
boolean shadowsEnabled;
shadowsEnabled := false;
%
% To simplify further calculations we need subdivided original path
%
vardef pathSubdivideBase (expr p, subdivideStep, i) =
save returnPath, sp;
path returnPath, sp;
returnPath := point i of p;
if i 0):
offsetPathLength := arclength(subpath (0, i) of p)/arclength(p);
else:
offsetPathLength := 0;
fi;
returnPath := (arclength(subpath (0, i) of p), offsetFunction);
if (i < length(p)):
% this thing is glitchy, but should be more accurate
%if (arclength(subpath (0, i) of p) < arclength(subpath (0, i + 1/4) of p)):
% offsetPathTime := i + 1/4;
% offsetPathLength := arclength(subpath (0, i + 1/4) of p)/arclength(p);
% instantDirection := unitvector((arclength(subpath (0, i + 1/4) of p), offsetFunction) - point 0 of returnPath);
% offsetPathTime := i + 1;
% offsetPathLength := arclength(subpath (0, i + 1) of p)/arclength(p);
% nextDirection := (arclength(subpath (0, i + 1) of p), offsetFunction);
% offsetPathTime := i + 3/4;
% offsetPathLength := arclength(subpath (0, i + 3/4) of p)/arclength(p);
% nextDirection := unitvector(nextDirection - (arclength(subpath (0, i + 3/4) of p), offsetFunction));
% returnPath := returnPath{instantDirection} .. {nextDirection}offsetPathTemplate(p, i + 1)(offsetFunction);
% returnPath := returnPath -- offsetPathTemplate(p, i + 1)(offsetFunction);
%else:
returnPath := returnPath -- offsetPathTemplate(p, i + 1)(offsetFunction);
%fi;
fi;
returnPath
enddef;
%
% This macro creates offset path p based on previously built template q, instead of function itself
% It is loosely based on something called Tiller-Hanson heuristic as described here:
% http://math.stackexchange.com/questions/465782/control-points-of-offset-bezier-curve
%
vardef offsetPathGenerate (expr p, q, i) =
save returnPath, c, d, pl, ps;
path returnPath, pl[];
pair c[], d[];
c1 := precontrol i of p;
c2 := point i of p;
c3 := postcontrol i of p;
if abs(c1-c2) = 0:
c1 := c2 shifted (c2-c3);
fi;
if abs(c3-c2) = 0:
c3 := c2 shifted (c2-c1);
fi;
if (abs(c1-c2) > 0) and (abs(c2-c3) > 0):
d1 := unitvector(c1-c2) rotated -90;
d2 := unitvector(c2-c3) rotated -90;
pl1 := (unitvector(c2-c1)--unitvector(c1-c2))
scaled arclength(subpath (i - 1/2, i + 1/2) of p)
shifted (point i of p shifted (d1 scaled ypart(point i of q)));
pl2 := (unitvector(c2-c3)--unitvector(c3-c2))
scaled arclength(subpath (i - 1/2, i + 1/2) of p)
shifted (point i of p shifted (d2 scaled ypart(point i of q)));
if (abs(angle(d1) - angle(d2)) > 2) and (xpart(pl1 intersectiontimes pl2) > 0):
c4 := pl1 intersectionpoint pl2;
else:
c4 := c2 shifted (d1 scaled ypart(point i of q));
fi;
returnPath := c4;
else:
returnPath := c2 shifted (unitvector( (point i-1 of p) - (point i+1 of p) rotated -90) scaled ypart (point i of q));
fi;
if i < length(p):
path ps;
ps := subpath (i, i + 1) of p;
c1 := point 0 of ps;
c2 := postcontrol 0 of ps;
c3 := precontrol 1 of ps;
c4 := point 1 of ps;
c5 := point 0 of returnPath;
if (abs(c3-c4)>0)
and (abs(c1-c2)>0)
and (abs(c1-c4)>0)
and (abs(direction i of q) > 0):
c6 := c4 shifted (unitvector(c4 - c3) rotated 90 scaled ypart(point i + 1 of q));
c7 := (c2 - c1) scaled (abs(c5-c6)/abs(c1-c4)) rotated angle(direction i of q) shifted c5;
c8 := (c3 - c4) scaled (abs(c5-c6)/abs(c1-c4)) rotated angle(direction i + 1 of q) shifted c6;
returnPath := returnPath .. controls c7 and c8 .. offsetPathGenerate (p, q, i + 1);
else:
returnPath := returnPath -- offsetPathGenerate (p, q, i + 1);
fi;
fi;
returnPath
enddef;
%
% Frontend for offsetPathGenerate and offsetPathTemplate
%
vardef offsetPath (expr p)(text offsetFunction) =
offsetPathGenerate (p, offsetPathTemplate(p, 0)(offsetFunction), 0)
enddef;
%
% Brush macro. It draws line with brush of variable width.
% For parts thicker than minStrokeWidth it uses offsetPath functions'
% results, for thiner parts it draws dashed lines of fixed width
%
def brushGenerate (expr p, q, i) =
begingroup
save w, brushPath, bt, t;
numeric w[], t[];
path brushPath[], bt;
bt := q;
w0 := (ypart(urcorner(bt)));
w1 := (ypart(lrcorner(bt)));
t := cutPathTime(bt, minStrokeWidth);
if ((w0 > minStrokeWidth)
and (w1 < minStrokeWidth)
and (t > 0)
and (t < length(p))
and (arclength(p) > minDashStrokeLength)
and (i < 10)):
brushGenerate (subpath (0, t) of p, subpath (0, t) of q, i + 1);
brushGenerate (subpath (t, length(p)) of p, subpath (t, length(q)) of q, i + 1);
elseif (arclength(p) > 0):
if (w0 > 99/100minStrokeWidth)
and (w1 > 99/100minStrokeWidth):
brushPath1 := offsetPathGenerate (p, q yscaled 1/2, 0);
brushPath2 := offsetPathGenerate (p, q yscaled -1/2, 0);
fill brushPath1 -- reverse(brushPath2) -- cycle;
elseif (w0 < 101/100minStrokeWidth) and (w1 < 101/100minStrokeWidth):
thinBrushGenerate (p, q, 0)
fi;
fi;
endgroup
enddef;
%
% macro for thin lines which are actually dashed
%
vardef thinBrushGenerate@#(expr p, q, i) =
begingroup
save w, brushPath, bt, t, h, minLength, minWidth, dashPatternImage;
numeric w[], t[];
path brushPath[], bt;
picture dashPatternImage;
bt := q;
w0 := (ypart(urcorner(bt)));
w1 := (ypart(lrcorner(bt)));
w2 := floor((1/2(w0 + w1))/dashStrokeWidthStep)*dashStrokeWidthStep;
t := cutPathTime(bt, w2);
brushPath1 := subpath (0, t) of p;
brushPath2 := subpath (t, length(p)) of p;
if (str @# = "") or (str @# = "hatches"):
minLength := minDashStrokeLength;
minWidth := minDashStrokeWidth;
elseif str @# = "stipples":
minLength := minStippleStep;
minWidth := minStippleStrokeWidth;
fi;
if (((w0 - w1) > dashStrokeWidthStep) and (i < 15))
and ((arclength(brushPath1) > minLength)
or (arclength(brushPath2) > minLength)):
thinBrushGenerate@#(brushPath1, subpath (0, t) of q, i + 1);
thinBrushGenerate@#(brushPath2, subpath (t, length(q)) of q, i + 1);
else:
if (str @# = "") or (str @# = "hatches"):
if (w2 > minStrokeWidth):
w2 := minStrokeWidth;
fi;
if (w2 >= minWidth) and (arclength(p) > 0):
draw p withpen thinpen dashed thinBrushPattern(w2, arclength(p));
fi;
elseif str @# = "stipples":
begingroup
interim linecap := rounded;
save stippleSizeVar;
stippleSizeVar := stippleSize;
save stippleSize;
w2 := 1/3w2;
if (w2 >= minWidth) and (arclength(p) > 0):
stippleSize := stippleSizeVar * (0.9 + uniformdeviate(0.3));
dashPatternImage := stipplesBrushPattern(w2, arclength(p));
if urcorner(dashPatternImage) <> (0,0):
brushPath1 := offsetPathGenerate (p, (q yscaled 0) shifted (0, 1/3stippleShadingDensity), 0);
draw brushPath1 withpen (pencircle scaled stippleSize) dashed dashPatternImage;
fi;
stippleSize := stippleSizeVar * (0.9 + uniformdeviate(0.3));
dashPatternImage := stipplesBrushPattern(w2, arclength(p));
if urcorner(dashPatternImage) <> (0,0):
brushPath2 := offsetPathGenerate (p, (q yscaled 0) shifted (0, -1/3stippleShadingDensity), 0);
draw brushPath2 withpen (pencircle scaled stippleSize) dashed dashPatternImage;
fi;
stippleSize := stippleSizeVar * (0.9 + uniformdeviate(0.3));
dashPatternImage := stipplesBrushPattern(w2, arclength(p));
if urcorner(dashPatternImage) <> (0,0):
draw p withpen (pencircle scaled stippleSize) dashed dashPatternImage;
fi;
fi;
endgroup
fi;
fi;
endgroup
enddef;
%
% this macro returns path as a shaded edge
%
vardef shadedEdge (expr p) =
image(
brushGenerate (p,
offsetPathTemplate (p, 0) (
1/2minStrokeWidth + 2*minStrokeWidth
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of p), arcsin(1/2))
), 0, point offsetPathTime of p
)
), 0);
)
enddef;
%
% Whenever we have brush thinner than minStrokeWidth we call this dash pattern macro
%
vardef thinBrushPattern (expr w, l) =
save d;
numeric d[];
d0 := w;
if d0 > minStrokeWidth: d0 := minStrokeWidth; fi;
% d1 is a result of some arbitrary function of line width
% we do not use simple linear function because minimal dash length
% also shouldn't be less than minStrokeWidth.
% After we get d1 other measurements are calculated,
% so filled area per unit length remains adequate and dashes are aligned
% with segments
d1 := (1/2minDashStrokeLength) + (((d0/minStrokeWidth)**2)*1/2minDashStrokeLength);
d1 := d1 + 1/2uniformdeviate(d1);
d2 := (minStrokeWidth - d0)*(d1/d0);
d3 := round(l/(d2 + d1));
if (d3 < 1): d3 := 1; fi;
d4 := (l/d3)/(d2 + d1);
d1 := d1*d4;
d2 := d2*d4;
if (uniformdeviate(2) > 1):
dashpattern (on d1 off 2d2)
else:
dashpattern (off 2d2 on d1)
fi
enddef;
%
% Stipples are also dashes
%
vardef stipplesBrushPattern (expr w, l) =
save d, n, rn, rv, ss;
ss := 1/1000;
numeric d[];
picture rv;
%if w > stippleSize:
% d0 := minStippleStep;
%else:
n := (w*l)/(stippleSize**2);
rn := floor(n);
if rn > 0:
d0 := l/rn;
fi;
%fi;
if rn > 0:
d1 := uniformdeviate(d0);
d2 := d0-d1;
if rn >=3:
d3 := uniformdeviate(d0);
d4 := d0-d3;
%if uniformdeviate(2) > 1:
% d5 := uniformdeviate(d0);
% d6 := d0-d5;
% rv := dashpattern (off d1 on ss off (d2+d5)-ss on ss off (d4+d6)-ss on ss off d3-ss);
%else:
rv := dashpattern (off d1 on ss off (d2+d3)-ss on ss off d4-ss);
%fi;
else:
rv := dashpattern (off d1 on ss off d2-ss);
fi;
else:
if uniformdeviate(1) < n:
rv := dashpattern (off uniformdeviate(l-ss) on ss off l);
else:
rv := image();
fi;
fi;
rv
enddef;
%
% macro that actually draws line of variable width
%
vardef brush (expr p) (text offsetFunction) =
image(
brushGenerate (p, offsetPathTemplate(p, 0)(offsetFunction), 0);
)
enddef;
%
% same, but only for thin brushes
%
vardef thinBrush@#(expr p) (text offsetFunction) =
image(
thinBrushGenerate@#(p, offsetPathTemplate(p, 0)(offsetFunction), 0);
)
enddef;
%
% This macro generates tube between paths p and q, of variable width d
% Tube is subdivided into segments in such a way that within every segment
% we need 2**n lines to generate even fill
%
vardef tubeGenerate@#(expr p, q, d, i) =
save w, bw, k, t, tubeWidth, sp, currentPath, currentTubePath, currentDepth, lineDensity;
numeric w[], bw[], t, currentDepth;
path tubeWidth, sp, currentPath, currentTubePath;
tubeWidth := d yscaled 2;
if (str @# = "") or (str @# = "hatches"):
lineDensity := shadingDensity;
elseif (str @# = "stipples"):
lineDensity := stippleShadingDensity;
fi;
w0 := (ypart(urcorner(tubeWidth))) - 1/1000;
w1 := (ypart(lrcorner(tubeWidth))) + 1/1000;
w2 := ceiling(log(w0/lineDensity, 2));
w3 := ceiling(log(w1/lineDensity, 2));
if ((w2 > w3) and (i<20)):
t := cutPathTime(tubeWidth, lineDensity*(2**(w2-1)));
tubeGenerate@#(subpath (0, t) of p, subpath (0, t) of q, subpath (0, t) of d, i + 1);
tubeGenerate@#(subpath (t, length(p)) of p, subpath (t, length(q)) of q, subpath (t, length(d)) of d, i + 1);
else:
if (arclength(p) > 0) and (arclength(q) > 0):
bw1 := 2**w2;
currentTubePath := interpath (1/2, q, p);
for k := 0 upto bw1:
currentPath := interpath (k/bw1, q, p);
angleOnTube := arcsin(((k/bw1)*2) - 1);
currentDepth := -abs((1-sin(angleOnTube + 1/2pi))*w0);
if shadowsEnabled:
currentPath := shadowCut(currentPath, currentDepth);
fi;
if (str @# = "") or (str @# = "hatches"):
brushGenerate (currentPath,
offsetPathTemplate(currentPath, 0)(
maxShadingStrokeWidth
if odd (k): * (abs(ypart(point offsetPathTime of tubeWidth)/bw1) - 1/2lineDensity) fi
* normalVectorToLightness(
tubeAnglesToNormalVector((
angleOnTube,
angleRad(direction offsetPathTime of currentTubePath),
angleRad(direction offsetPathTime of (tubeWidth yscaled 1/2))
)), currentDepth, point offsetPathTime of currentPath)
), 0);
elseif (str @# = "stipples"):
begingroup
save stippleShadingDensity;
if w2 > 0:
stippleShadingDensity := 2w0/(2**w2); % When the distance between the lines changes, wtipples should spread further apart
else:
stippleShadingDensity := w0;
fi;
thinBrushGenerate.@#(currentPath,
offsetPathTemplate(currentPath, 0)(
stippleSize
if odd (k): * (abs(ypart(point offsetPathTime of tubeWidth)/bw1) - 1/2lineDensity) fi
* normalVectorToLightness(
tubeAnglesToNormalVector((
angleOnTube,
angleRad(direction offsetPathTime of currentTubePath),
angleRad(direction offsetPathTime of (tubeWidth yscaled 1/2))
)), currentDepth, point offsetPathTime of currentPath)
), 0);
endgroup
fi;
endfor;
fi;
fi;
enddef;
%
% This macro is analogous to tubeGenerate, but draws transverse strokes
% result is somewhat suboptimal for now, but in simple cases it works ok
%
def tubeGenerateAlt (expr p, q, d) =
begingroup
save spth, lpth, currentPath, pos, t, pthdir, corr, o, l, i, j, k, tubeAngle, pathAngle, scorr, dt;
numeric l[];
path spth, lpth, currentPath;
pos := 0;
j := 0;
forever:
dt := (xpart(point pos of d) + 1/2shadingDensity);
scorr := cosd(angle(direction xpart(d intersectiontimes ((dt, ypart(lrcorner(d))) -- (dt, ypart(urcorner(d))))) of d));
t1 := arctime ((arclength(subpath(0, pos) of p)) + shadingDensity/scorr) of p;
t2 := arctime ((arclength(subpath(0, pos) of q)) + shadingDensity/scorr) of q;
if (arclength(subpath(pos, t1) of p) < arclength(subpath(pos, t1) of q)):
pthdir := -1;
t3 := t1;
else:
pthdir := 1;
t3 := t2;
fi;
corr := round(arclength(subpath(pos, t3) of if pthdir = 1: p else: q fi)/(shadingDensity/scorr));
if (corr < 1): corr := 1; fi;
corr := (arclength(subpath(pos, t3) of if pthdir = 1: p else: q fi) - (corr*(shadingDensity/scorr)))/corr;
t3 := arctime (arclength(subpath(0, t3) of if pthdir = 1: q else: p fi) - 1/3corr) of if pthdir = 1: q else: p fi;
spth := subpath(pos, t3) of if pthdir = 1: q else: p fi;
lpth := subpath(pos, t3) of if pthdir = 1: p else: q fi;
tubeAngle := angleRad(direction 1/2[pos, t3] of d);
pathAngle := angleRad(direction 1/2 of interpath (1/2, spth, lpth));
pos := t3;
l1 := round(arclength(lpth)/(shadingDensity/scorr));
if (l1 < 1): l1 := 1; fi;
l2 := arclength(lpth)/(l1*(shadingDensity/scorr));
for i := 0 upto l1 - 1:
j := j + 1;
k := i*(arclength(lpth)/l1);
currentPath := point (arctime k of lpth) of spth -- point (arctime k of lpth) of lpth;
currentPath := offsetPathSubdivide(currentPath);
brushGenerate (
currentPath,
offsetPathTemplate(currentPath, 0)(
maxShadingStrokeWidth
* orderFade(offsetPathLength[1/l1, l2], j)
* normalVectorToLightness(
tubeAnglesToNormalVector((
arcsin(pthdir*((offsetPathLength*2)-1)),
pathAngle,
tubeAngle)
), -2(1/2arclength(currentPath))+sqrt(1 - (2offsetPathLength - 1)**2)*(1/2arclength(currentPath)), point offsetPathTime of currentPath)
)
, 0);
endfor;
exitif pos >= length(p);
endfor;
endgroup
enddef;
%
% This macro converts some measurements of point on tube to absolute angle.
% Since there are three such measurements, macro gets them as as a single
% argument of "color" type, in case it would eventually appear as a result
% of some other macro.
%
% redpart is the angle on the tube's circumference
% greenpart is the angle of the tube path
% bluepart is the angle of the tube's outline
%
vardef tubeAnglesToNormalVector (expr p) =
save normalVector;
color normalVector;
normalVector := (0, 0, 1);
normalVector := rotateXYZaround.y(normalVector, -bluepart(p));
normalVector := rotateXYZaround.x(normalVector, redpart(p));
normalVector := rotateXYZaround.z(normalVector, -greenpart(p));
normalVector
enddef;
%
% frontend to simplify tube drawing. tubeOutline variable changes on every call
% of the function and can be used afterwards.
%
path tubeOutline;
boolean drawTubeEnds;
drawTubeEnds := true;
vardef tube@#(expr p)(text offsetFunction)=
save q, respic;
path q[];
picture respic;
q0 := offsetPathSubdivide(p);
q1 := offsetPathTemplate(q0, 0)(offsetFunction);
q2 := offsetPathGenerate (q0, q1, 0);
q3 := offsetPathGenerate (q0, q1 yscaled -1, 0);
tubeOutline := q3--reverse(q2)--cycle;
if str @# = "e":
if not drawTubeEnds:
image(
draw q2 withpen thinpen;
draw q3 withpen thinpen;
)
else:
tubeOutline
fi
else:
image(
if str @# = "l":
respic := image(tubeGenerate (q2, q3, q1, 0););
elseif str @# = "s":
respic := image(tubeGenerate.stipples(q2, q3, q1, 0);)
elseif str @# = "t":
respic := image(tubeGenerateAlt (q2, q3, q1););
fi;
if (cycle p) or (not drawTubeEnds):
draw q2 withpen thinpen;
draw q3 withpen thinpen;
else:
draw q2--reverse(q3)--cycle withpen thinpen;
fi;
clip respic to (q2--reverse(q3)--cycle);
draw respic;
)
fi
enddef;
%
% Sphere can be used as a cap for a tube, so it has same 2**n lines.
%
vardef sphere@#(expr d) =
save currentCircle, origCircle, currentRadius, currentDepth, order, circleThickness, lineDensity, shadingPicture;
path currentCircle, origCircle;
numeric currentRadius, currentDepth, order;
picture shadingPicture;
if (str @# = "") or (str @# = "c"):
lineDensity := shadingDensity;
elseif (str @# = "s"):
lineDensity := stippleShadingDensity;
fi;
origCircle := fullcircle;
order := 2**ceiling(log((1/2d)/lineDensity, 2));
image(
draw fullcircle scaled d withpen thinpen;
shadingPicture := image(
for i := 1 upto order:
currentRadius := i/order;
currentCircle := origCircle scaled (currentRadius*d) rotated uniformdeviate (1/4pi);
if odd(i):
circleThickness := maxShadingStrokeWidth * ((abs(d - (lineDensity*order)))/order);
else:
circleThickness := maxShadingStrokeWidth;
fi;
currentDepth:= -(1-sqrt(1-currentRadius**2))*(1/2d);
if shadowsEnabled:
currentCircle := shadowCut(currentCircle, currentDepth);
fi;
if (str @# = "") or (str @# = "c"):
brushGenerate (currentCircle,
offsetPathTemplate (currentCircle, 0) (
circleThickness
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of currentCircle), arcsin(currentRadius))
), currentDepth, point offsetPathTime of currentCircle
)
), 0);
elseif (str @# = "s"):
begingroup
save stippleShadingDensity;
if order > 0:
stippleShadingDensity := d/order; % When the distance between the lines changes, wtipples should spread further ap
else:
stippleShadingDensity := 1/2d;
fi;
thinBrushGenerate.stipples(currentCircle,
offsetPathTemplate (currentCircle, 0) (
circleThickness
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of currentCircle), arcsin(currentRadius))
), currentDepth, point offsetPathTime of currentCircle
)
), 0);
endgroup
fi;
endfor;
);
clip shadingPicture to (fullcircle scaled d);
draw shadingPicture;
)
enddef;
%
% Alternative sphere macro. It's all about latitudinal strokes.
% The idea is: when we have a sphere with evenly distributed parallel strokes
% we know how their density rises towards edge in a projection,
% so all we need to do is to fade lines correspondingly
%
vardef sphereLat (expr d, lat) =
save p, a, x, y, sphlat, latrad, n, c, currentPath, nline, tlat;
path p[], currentPath, currentArc;
sphlat := 0;
nline := 0;
latrad := (2pi*lat/360);
n := ceiling((pi*1/2d)/shadingDensity);
if (cosd(lat) <> 0): tlat := (sind(lat)/cosd(lat)); fi;
image(
draw fullcircle scaled d withpen thinpen;
p0 := fullcircle rotated 90;
for nline := 1 upto n-1:
sphlat := nline*(pi/n);
if (sphlat + latrad < pi) and (sphlat + latrad > 0):
if (cosd(lat) <> 0):
if (sin(sphlat) <> 0):
x := tlat*(cos(sphlat)/sin(sphlat));
else:
x := 0;
fi;
else:
if ((sphlat > 1/2pi) and (lat > 0)) or ((sphlat < 1/2pi) and (lat < 0)):
x := -2;
else:
x := 2;
fi;
fi;
if (abs(x) <= 1):
y := arcsin(x);
p1 := subpath(6 + 8y/2pi, 2 - 8y/2pi) of p0;
else:
p1 := p0;
fi;
if (x > -1) and (arclength(p1) > 0):
currentPath := (p1 scaled (d*sin(sphlat)) yscaled sind(lat)) shifted (0, 1/2d*cos(sphlat)*cosd(lat));
currentPath := offsetPathSubdivide(currentPath);
brushGenerate(currentPath,
offsetPathTemplate(currentPath, 0)(
maxShadingStrokeWidth * orderFade(
sqrt(1 -
abs(
ypart(point offsetPathTime of currentPath)/(1/2d),
1 - abs(
1 - abs(
xpart(point offsetPathTime of currentPath)
/(1/2d)
)
)**abs(sind(lat))
)**2)
, nline)
* normalVectorToLightness(
sphereAnglesToNormalVector((
(
if (abs(point offsetPathTime of currentPath) > 0):
angleRad(point offsetPathTime of currentPath)
else:
0
fi
), arcsin(2abs(point offsetPathTime of currentPath)/(d+1)))
), 0, point offsetPathTime of currentPath)
), 0);
fi;
fi;
endfor;
)
enddef;
vardef orderFade (expr v, n)=
save o;
if (v > 1/256):
o := 2**ceiling(log(1/v, 2));
if ((n mod 1/2o) = 0):
if ((n mod o) = 0):
1
else:
(v*o) - 1
fi
else:
0
fi
else:
0
fi
enddef;
%
% This one converts point location on sphere to normal vector
%
vardef sphereAnglesToNormalVector (expr p) =
save normalVector;
color normalVector;
normalVector := (0, 0, 1);
normalVector := rotateXYZaround.y(normalVector, ypart(p));
normalVector := rotateXYZaround.z(normalVector, -xpart(p));
normalVector
enddef;
%
% Once we get two angles at some point of some surface, we can compute light intensity there.
%
vardef normalVectorToLightness (expr normalVector, d, q) =
save returnValue, shiftedShadowPath;
path shiftedShadowPath;
if shadowsEnabled:
for i := 0 step 1 until numberOfShadows:
shiftedShadowPath := shadowPath[i] shifted
((redpart(lightDirectionVectorXYZ), greenpart(lightDirectionVectorXYZ))
scaled ((d-shadowDepth[i])*bluepart(lightDirectionVectorXYZ)));
if q isInside shiftedShadowPath:
returnValue := 1;
fi;
endfor;
fi;
if not known returnValue:
returnValue := 1 - (normalVector dotprodXYZ lightDirectionVectorXYZ);
returnValue := lightnessPP(returnValue);
fi;
if returnValue > 1:
1
else:
returnValue
fi
enddef;
vardef lightnessPP (expr v) =
v
enddef;
% Shadows are global
path shadowPath[];
numeric shadowDepth[];
% Shadows either require high path resolution, or some points
% on a path in just the right place for shadows.
% This macro adds such points.
vardef shadowCut (expr pathToCut, currentDepth)=
save shiftedShadowPath, pathShadowIntersection, pathShadowCut, currentPath;
path shiftedShadowPath, currentPath;
pair pathShadowIntersection;
numeric pathShadowCut;
currentPath := pathToCut;
for j := 0 step 1 until numberOfShadows:
shiftedShadowPath := shadowPath[j] shifted
((redpart(lightDirectionVectorXYZ), greenpart(lightDirectionVectorXYZ))
scaled ((currentDepth - shadowDepth[j])*bluepart(lightDirectionVectorXYZ)));
forever:
pathShadowIntersection := shiftedShadowPath firstIntersectionTimes currentPath;
pathShadowCut := ypart(pathShadowIntersection);
if (pathShadowCut > 1/10) and (pathShadowCut < length(currentPath) - 1/10):
currentPath := (subpath (0, pathShadowCut - 1/20) of currentPath) .. (subpath (pathShadowCut + 1/20, length(currentPath)) of currentPath);
fi;
shiftedShadowPath := subpath (xpart(pathShadowIntersection) + 1/5, length(shiftedShadowPath)) of shiftedShadowPath;
exitif (pathShadowCut = -1);
endfor;
endfor;
currentPath
enddef;
%
% Several macros rely on cutting offset path template at given height.
% Taking cutting points closer to the middle gives better results, and that's
% just what this macro tries to do.
%
vardef cutPathTime (expr p, h) =
save cutTime, d;
numeric cutTime[], d[];
d1 := xpart(urcorner(p));
d2 := xpart(ulcorner(p));
if (d2 < d1):
d3 := 1/2(d1 + d2);
cutTime1 := ypart(((d3, h) -- (d1, h)) firstIntersectionTimes p);
cutTime2 := ypart(((d3, h) -- (d2, h)) firstIntersectionTimes p);
d4 := xpart (point cutTime1 of p);
d5 := xpart (point cutTime2 of p);
if abs(d4-d3) < abs(d5-d3):
cutTime3 := cutTime1
else:
cutTime3 := cutTime2
fi;
else:
cutTime3 := -1;
fi;
cutTime3
enddef;
%
% This macro calculates ray angle after refraction. It takes raw angles (one of ray — p and one of surface — q)
% and refraction indices ratio. Whether ray comes from opticaly denser material is determined by direction of q
% relative to that of p
%
vardef refractionAngle (expr p, q, n) =
save a;
numeric a[];
a0 := p - q;
if (sin(a0) < 0):
a1 := cos(a0 + pi) * n;
a2 := pi;
else:
a1 := cos(a0) / n;
a2 := 0;
fi;
if abs(a1) <= 1:
a3 := arccos(a1) + q + a2;
else:
a3 := -1000;
fi;
a3
enddef;
%
% Same thing for reflection angle, just in case
%
vardef reflectionAngle (expr p, q) =
(2pi - p + 2q)
enddef;
%
% This macro returns path of ray 'sa' (which can actually be any path, but only ray from next to last to last point
% will count) refracted with coef. n through some shape p; if ray can't be refracted and, therefore, totally reflected,
% it will contunue as reflected from that point. i is total number of refractions to compute;
%
vardef refractionPathR (expr sa, p, n, i, mn) =
save ray, resultRay, d, s, a, iT;
path ray, resultRay;
pair s, iT;
numeric d[], a;
s := point (length(sa) - 1) of sa;
a := angleRad((point (length(sa)) of sa)-(point (length(sa) - 1) of sa));
ray := (s shifted (-dirRad(a) scaled 2)) -- s -- (s shifted (dirRad(a) scaled (abs(llcorner(p)-(urcorner(p))) + abs(s-(center(p))))));
if (i > 0): ray := subpath (1 + 1/1000, 2) of ray; fi;
iT := ray firstIntersectionTimes p;
d1 := xpart(iT);
d2 := ypart(iT);
d3 := a;
d4 := angleRad(direction d2 of p);
if (n > 0):
d5 := refractionAngle(d3, d4, n);
if (d5 < -100) and (d2 >= 0):
d5 := reflectionAngle(d3, d4);
fi;
else:
d5 := reflectionAngle(d3, d4);
fi;
if (d1 >= 0) and (i < mn) and (d5 > -100):
resultRay := (subpath (0, length(sa) - 1) of sa) -- refractionPathR(point d2 of p -- (point d2 of p shifted dirRad(d5)), p, n, i + 1, mn);
else:
if (d5 > -100) or (d1 < 0):
resultRay := subpath (0, 1/2) of ray;
else:
resultRay := subpath (0, d1) of ray;
fi;
fi;
resultRay
enddef;
vardef refractionPath (expr sa, p, n) =
refractionPathR(sa, p, n, 0, 10)
enddef;
%
% These macros are for isolines. cLine draws continuous line and is called by isoLines.
% For now they are only used to draw wood texture, but can be used elsewhere
%
%
% isoLines goes through i by j matrix of nodes (xy), looking for square, that has some of it's
% angles below zero and some - above, when found, it calls cLine, that tries to build segment of
% isoline, that happen to go through abovementioned square. Thickness of line is
% controlled by values in v array.
% All squares with lines already drawn through are ignored.
%
vardef isoLines (suffix xy)(expr cs, l, s) =
save xxyy, i, j, c, v, sqB, iL, lvl;
numeric xxyy[][], c[], v[], sqB, sqbM;
lvl := l;
path iL;
image(
for i := 0 step 1 until xpart(cs) - 1:
for j := 0 step 1 until ypart(cs) - 1:
if (unknown xxyy[i][j]):
c1 := xy[i][j]+lvl;
c2 := xy[i][j+1]+lvl;
c3 := xy[i+1][j]+lvl;
c4 := xy[i+1][j+1]+lvl;
sqB := 0;
sqBm := 0;
if (abs(sign(c1)+sign(c2)+sign(c3)+sign(c4)) < 4):
iL := cLine (xy)((i, j), (0, 0), 0, cs) scaled s;
brushGenerate (reverse(iL),
offsetPathTemplate(iL, 0)(
1/16minStrokeWidth
/(1/64 + 2(
if (offsetPathTime < length(iL) - 1):
(offsetPathTime - floor(offsetPathTime))
[v[floor(sqB + offsetPathTime)],
v[ceiling(sqB + offsetPathTime)]]
else:
1
fi
))
)
, 0);
fi;
fi;
endfor;
endfor;
draw (0,0);
)
enddef;
% cLine tries to generate continouos segment of an isoline
vardef cLine (suffix xy)(expr ij, dr, st, cs) =
save p, d, dd, n, i, j, k, outputPath, sqS, cp, dp, nd;
pair p[], d[], dd[], cp[], dp[];
path outputPath;
i := xpart(ij);
j := ypart(ij);
sqS := 0;
if (i >= 0) and (i <= xpart(cs)-1) and (j >= 0) and (j <= ypart(cs)-1):
n := 0;
c1 := xy[i][j] + lvl; d1:= (i, j);
c2 := xy[i][j+1] + lvl; d2:= (i, j+1);
c3 := xy[i+1][j] + lvl; d3:= (i+1, j);
c4 := xy[i+1][j+1] + lvl; d4:= (i+1, j+1);
cp1 := (1, 2); dp1 := (-1, 0);
cp2 := (1, 3); dp2 := (0, -1);
cp3 := (2, 4); dp3 := (0, 1);
cp4 := (3, 4); dp4 := (1, 0);
for k := 1 upto 4:
c5 := c[xpart(cp[k])]; d5 := d[xpart(cp[k])];
c6 := c[ypart(cp[k])]; d6 := d[ypart(cp[k])];
if (sign(c5)) <> (sign(c6)):
n := n + 1;
p[n] := (abs(-c5/(c6-c5)))[d5, d6];
dd[n] := dp[k];
fi;
endfor;
sqS := max(c1, c2, c3, c4) - min(c1, c2, c3, c4);
if (unknown xxyy[i][j]):
xxyy[i][j] := 1;
if (dr = (0, 0)):
outputPath := cLine (xy)(ij shifted dd2, dd2, st + 1,cs) -- p1 -- p2 -- cLine (xy)(ij shifted dd1, dd1, st - 1,cs);
else:
nd := 0;
if (unknown(xxyy[i + xpart(dd1)][j + ypart(dd1)])):
nd := 1;
elseif (unknown(xxyy[i + xpart(dd2)][j + ypart(dd2)])):
nd := 2;
fi;
if nd > 0:
outputPath := cLine (xy)(ij shifted dd[nd], dd[nd], st + sign(st),cs);
p3 := p[nd];
if (st > 0):
outputPath := outputPath -- p3;
else:
outputPath := p3 -- outputPath;
fi;
else:
outputPath := 1/2[p1, p2];
fi;
fi;
else:
outputPath := 1/2[p1, p2];
fi;
else:
outputPath := (i, j);
xxyy[i][j] := 1;
fi;
if (st < sqB): sqB := st; fi;
if (st > sqBm): sqBm := st; fi;
v[st] := sqS;
outputPath
enddef;
%
%
%
% In following section are gathered all the things for some commonly used
% real life objects
%
%
%
%
% Though there's a decent library to deal with geography for mp already
% (http://melusine.eu.org/lab/bmp/a_mp-geo), here's some basic globe-drawing
% routine for simple cases, note that latitude starts from the pole,
% not from the equator (for convenience sake)
% Below some landmasses are defined
%
path landmass[];
landmass1 := (206, 122.33)--(211.07, 116)--(213.3, 109.94)--(218.57, 106.03)--(218.38, 97.36)--(220.28, 91.28)--(229.75, 78.07)--(221.41, 78.29)--(220.78, 76.52)--(218.07, 74.48)--(213.8, 66.08)--(213.38, 62.04)--(222.31, 77.1)--(233.88, 72.27)--(237.79, 68.59)--(234.88, 64.69)--(229.83, 65.57)--(228.98, 64.73)--(227.37, 59.82)--(250.57, 68.12)--(254.63, 80.83)--(257.07, 80.93)--(257.38, 80.52)--(258.64, 75.5)--(266.4, 68.48)--(269.56, 67.49)--(271.88, 70.43)--(272.67, 74.49)--(275.36, 72.94)--(276.87, 78.6)--(276.68, 79.04)--(276.11, 79.28)--(276.3, 80.22)--(276.75, 79.96)--(276.56, 82.38)--(277.05, 82.04)--(280.5, 86.44)--(277.25, 85.56)--(276.55, 88.03)--(279.47, 92.77)--(283.29, 92.25)--(282.68, 90.91)--(283.74, 90.4)--(282.53, 89.58)--(283.03, 88.6)--(278.44, 80.08)--(279.15, 76.64)--(281.08, 78.25)--(282.29, 80.21)--(285.35, 79.72)--(288, 77.83)--(284.21, 71.22)--(287.94, 68.57)--(288, 68.6)--(288.74, 69.82)--(300.09, 61.89)--(300.86, 59.94)--(299.36, 59.63)--(297.64, 55.13)--(301.24, 52.55)--(296.1, 51.5)--(300.45, 49.51)--(299.83, 50.75)--(299.84, 50.82)--(299.44, 51.42)--(303.59, 50.57)--(302.72, 51.9)--(302.96, 52.12)--(304.97, 52.87)--(304.12, 55.13)--(307.89, 53.38)--(306.37, 50.11)--(308.65, 47.92)--(315.01, 45.12)--(319.69, 40.31)--(320.43, 44.25)--(321.66, 44.31)--(323.19, 41.66)--(320.37, 35.59)--(318.47, 37.21)--(315.99, 36.32)--(313.68, 35.16)--(320.43, 31.11)--(332.73, 30.38)--(338.5, 28.24)--(340.91, 28.61)--(334.92, 32.27)--(335, 39.2)--(340.58, 35.32)--(341.69, 32.15)--(340.43, 31.93)--(344.49, 29.68)--(352.49, 28.33)--(355.9, 25.35)--(358.67, 24.01)--(366.1, 25.61)--(368.78, 23.99)--(319.11, 17.34)--(309.82, 19)--(308.23, 18.4)--(307.69, 16.74)--(297.49, 16.63)--(290.61, 13.26)--(285.38, 13.37)--(284.06, 12.79)--(258.59, 16.61)--(260.79, 18.13)--(254.13, 18.01)--(253.53, 17.04)--(252.25, 17.02)--(252.44, 18.56)--(253.69, 19.64)--(251.71, 20.89)--(249.66, 16.97)--(245.54, 19.39)--(236.64, 19.73)--(239.08, 21)--(237.57, 21.46)--(232.4, 21.62)--(232.29, 21.34)--(225.16, 22.27)--(221.46, 21.23)--(218.52, 25.09)--(216.81, 24.69)--(214.76, 24.93)--(214.95, 25.52)--(213.66, 25.25)--(211.67, 23.33)--(215.44, 23.49)--(217.75, 21.65)--(200.52, 19.57)--(194.37, 21.1)--(186.19, 26.3)--(183.33, 30.4)--(187.61, 31.42)--(191.44, 33.88)--(194.61, 33.54)--(197.17, 30.74)--(196.08, 28.46)--(196.04, 27.66)--(203.54, 24.58)--(203.45, 24.88)--(200.38, 27.18)--(200.91, 27.54)--(200.05, 29.69)--(199.62, 29.91)--(201.03, 30.25)--(207.36, 29.93)--(205.2, 31.05)--(199.88, 30.97)--(199.94, 31.44)--(200.26, 32.2)--(202.19, 31.76)--(202.85, 32.2)--(199.62, 32.83)--(199.15, 34.61)--(189.46, 35.87)--(189.93, 35.46)--(191.12, 35.08)--(190.83, 34.56)--(188.1, 33.45)--(186.87, 34.75)--(187.11, 36.02)--(176.39, 40.19)--(176.65, 41.24)--(173.41, 42.02)--(176.82, 43.77)--(169.68, 46.56)--(169.15, 53.05)--(171.1, 53.62)--(173.12, 53.39)--(178.7, 51.26)--(183.17, 46.73)--(186.38, 46.75)--(192.72, 49.52)--(191.46, 52.64)--(193.74, 52.83)--(196.74, 50.32)--(190.71, 44.65)--(191.74, 44.4)--(198.11, 50.06)--(198.89, 52.03)--(200.95, 53.75)--(202.49, 51.99)--(201.15, 49.3)--(204.15, 49.28)--(206.54, 53.44)--(214.39, 53.25)--(211.18, 58.75)--(198.36, 57.7)--(197.88, 59.47)--(188.5, 55.87)--(189.63, 53.18)--(189.49, 52.79)--(173.31, 54.75)--(168.56, 58.21)--(161.34, 69.75)--(160.58, 75.18)--(161.25, 77.58)--(162.08, 79.09)--(163.71, 80.23)--(165.04, 82.46)--(168.88, 84.86)--(182.72, 83.77)--(184.88, 85.79)--(187.22, 85.99)--(186.79, 90.34)--(190.56, 95.97)--(190.23, 105.47)--(193.05, 115.74)--(196.18, 121.46)--(196.92, 124.65)--(206, 122.33)--(206, 122.33);
landmass2 := (111.44, 45.06)--(113.41, 44.75)--(111.77, 46)--(111.77, 46.07)--(118.69, 43.98)--(118.13, 42.88)--(116.49, 43.6)--(114.48, 42.7)--(114.1, 43.65)--(114.04, 41.9)--(113.28, 42.04)--(108.57, 42.18)--(114.57, 39.81)--(120.91, 38.84)--(119.04, 41.38)--(119.53, 41.59)--(122.9, 42.64)--(121.94, 42.77)--(121.82, 43.31)--(124.3, 42.48)--(125.29, 42.37)--(125.59, 41.84)--(125.41, 41.3)--(124.75, 41.38)--(122.58, 40.38)--(122.97, 39.95)--(121.71, 40.06)--(123.02, 38.38)--(121.65, 38.53)--(120.78, 35.69)--(116.8, 33.48)--(114.09, 29.59)--(110.59, 31.42)--(108.79, 28.81)--(104.18, 27.54)--(99.54, 29.42)--(100.85, 30.15)--(100.61, 31.83)--(100.51, 34.6)--(98.7, 35.56)--(99.11, 36.37)--(99.33, 38.41)--(96.3, 36.78)--(85.98, 32.83)--(83.79, 29.73)--(86.02, 28)--(87.63, 26.53)--(92.02, 23.6)--(94.53, 23.9)--(94.57, 23.59)--(95.6, 23.75)--(96.08, 21.63)--(96.05, 20.47)--(99.61, 19.73)--(103.5, 21.33)--(105.9, 22.76)--(103.75, 23.87)--(104.54, 24.37)--(101.78, 25.83)--(112.52, 27.74)--(113.4, 27.33)--(113.39, 27.23)--(110.6, 24.25)--(110.62, 24.18)--(111.27, 23.6)--(116.34, 23.81)--(117.17, 23.3)--(110.5, 21.34)--(111.78, 20.87)--(110.82, 20.4)--(111.3, 20.21)--(109.74, 19.84)--(110.11, 19.43)--(102.09, 17.29)--(97.38, 16.64)--(92.88, 17.67)--(92.21, 18.01)--(91.83, 17.28)--(89.16, 18.82)--(92.76, 20.05)--(93.22, 20.96)--(92.06, 21.98)--(91.69, 22.39)--(88.11, 21.59)--(87.79, 20.91)--(86.11, 20.22)--(86.94, 19.77)--(84.28, 18.1)--(84.81, 17.32)--(85.92, 16.37)--(82.62, 16.82)--(82.26, 18.46)--(82.22, 20.63)--(80.29, 21.42)--(73.81, 21.72)--(73.19, 21.56)--(73.02, 21.15)--(75.97, 21.21)--(75.92, 20.34)--(77.6, 20.64)--(73.05, 17.2)--(70.79, 18.11)--(67.81, 17.27)--(64.31, 17.23)--(54.27, 15.93)--(52.31, 18.03)--(60.06, 17.29)--(60.43, 18.73)--(59.79, 19.06)--(65.44, 20.05)--(71.86, 20.7)--(72.16, 20.95)--(69.43, 21.73)--(70.4, 21.96)--(70.06, 22.49)--(69.48, 22.4)--(63.3, 21.97)--(48.4, 19.77)--(41.64, 20.99)--(22.72, 18.59)--(11.06, 21.67)--(16.58, 23.75)--(14.57, 23.75)--(10.29, 24.54)--(17.36, 25.3)--(17.42, 26.18)--(12.61, 29.54)--(15.96, 29.89)--(15.52, 31.43)--(20.94, 31.31)--(13.31, 35.43)--(16.05, 35.66)--(19.1, 35.11)--(18.12, 34.75)--(27.83, 28.87)--(28.89, 30.18)--(33.87, 29.92)--(33.36, 30.35)--(38.31, 30.44)--(42.27, 33.16)--(42.87, 32.94)--(47.84, 35.37)--(50.85, 39.13)--(53.79, 42.3)--(54.39, 50.26)--(64.16, 61.59)--(63.13, 61.47)--(62.78, 61.95)--(64.93, 63.33)--(66.24, 65.62)--(67.78, 65.49)--(65.67, 61.75)--(65.09, 58.99)--(66.42, 61.44)--(68.81, 63.42)--(72.94, 69.36)--(81.5, 74.4)--(83.61, 73.92)--(85.99, 75.53)--(90.7, 76.75)--(93.55, 80.26)--(95.61, 82.43)--(96.7, 82.3)--(98.47, 82.54)--(101.05, 86.23)--(98.22, 93.14)--(100.55, 101.04)--(102.97, 105.27)--(107.45, 108.16)--(104.04, 132.27)--(104.08, 133.7)--(103.49, 134.62)--(102.5, 136.82)--(104.04, 136.98)--(103.44, 141.15)--(104.17, 143.63)--(108.08, 144.9)--(110.27, 145.35)--(111.35, 145.04)--(113.34, 144.63)--(109.31, 140.79)--(112.69, 137.21)--(111.47, 135.36)--(113.17, 133.67)--(113.34, 130.92)--(116.56, 129.1)--(119.97, 124.37)--(128.41, 119.87)--(133.19, 114.17)--(136.37, 113.08)--(138.56, 109.76)--(141.63, 100.78)--(139.94, 93.61)--(133.56, 91.17)--(129.89, 90.92)--(128.12, 89.29)--(123.75, 83.93)--(119.97, 83.01)--(117.24, 81.38)--(117.85, 80.79)--(116.9, 80.06)--(117.65, 79.02)--(114.24, 79.23)--(110.12, 78.93)--(108.35, 77.69)--(102.27, 80.15)--(102.6, 80.45)--(101.42, 81.66)--(96.3, 80.94)--(95.6, 75.16)--(91.88, 73.8)--(89.69, 73.89)--(91.19, 69.61)--(91.08, 68.3)--(87.73, 69.96)--(81.32, 69.32)--(81.63, 61.96)--(88.41, 60.66)--(88.78, 61.23)--(92.41, 60.31)--(95.93, 65.23)--(96.7, 65.47)--(97.12, 65.14)--(97.12, 59.81)--(100.34, 56.62)--(103.22, 54.85)--(104.21, 50.99)--(106.28, 48.84)--(108.51, 48.83)--(108.23, 48.05)--(111.68, 45.41);
landmass3 := (309.58, 101.89)--(307.85, 104.82)--(304.23, 104.36)--(301.98, 106.89)--(301.81, 107.2)--(301.39, 106.32)--(297.9, 109.91)--(292.61, 112.27)--(292.42, 116.24)--(291.76, 116.34)--(293.22, 123.3)--(295.54, 125.57)--(300.79, 123.84)--(313.31, 124.74)--(314.35, 125.13)--(314.72, 125.92)--(316.53, 125.74)--(318.45, 127.71)--(324.06, 128.78)--(326.88, 127.94)--(328.82, 125.64)--(331.64, 120.19)--(331.26, 115.24)--(328.34, 111.83)--(327.46, 109.78)--(327.32, 109.75)--(323.91, 106.24)--(320.57, 100.41)--(318.09, 106.52)--(315.29, 105.19)--(314.65, 104.27)--(314.4, 103.64)--(314.38, 101.88)--(314.02, 100.8)--(310.93, 100.72)--(310, 101.21)--(309.58, 101.89);
landmass4 := (360, 173.94)--(347.31, 173.02)--(337.83, 170.31)--(340.03, 168.66)--(345.63, 168.61)--(346.01, 167.92)--(341.09, 166.89)--(341.61, 165.5)--(343.88, 164.64)--(343.24, 163.51)--(349.37, 161.91)--(322.58, 157.73)--(323.11, 157.14)--(263.4, 156.39)--(263.26, 156.59)--(263.82, 157.03)--(245.18, 163.09)--(245.85, 157.68)--(234.93, 156.8)--(226.52, 156.65)--(194.69, 160.02)--(194.23, 160.12)--(182.27, 160.22)--(175.88, 161.21)--(175.09, 160.92)--(174.87, 161.24)--(171.99, 160.65)--(165.8, 161.33)--(163.91, 162.6)--(161.68, 163.73)--(141.54, 169.39)--(118.58, 172.13)--(116.78, 171.69)--(102.46, 169.67)--(94.78, 168.49)--(97.74, 167.88)--(101.28, 166.75)--(117.86, 164.33)--(118.52, 163.02)--(117.24, 161.46)--(117.41, 160.81)--(114.93, 158.64)--(113.38, 158.53)--(113.67, 158.17)--(112.94, 157.45)--(115.93, 156.83)--(115.81, 156.34)--(116.18, 155.88)--(116.84, 155.67)--(119.7, 154.63)--(119.77, 154.11)--(121.16, 153.99)--(121.76, 153.53)--(116.74, 154)--(113.84, 154.79)--(110.65, 157.28)--(111.18, 157.79)--(111.19, 159.12)--(104.73, 163.13)--(104.3, 163)--(90.27, 162.69)--(86.46, 162.59)--(74.77, 162.74)--(78.14, 164.63)--(64.48, 164.84)--(65.31, 164.46)--(50.04, 164.22)--(50.29, 164.61)--(32.44, 165.76)--(29.52, 166.6)--(27.21, 166.74)--(27.17, 166.97)--(22.41, 166.97)--(23.06, 168.33)--(21.43, 168.63)--(29.78, 169.97)--(27.58, 170.36)--(29.1, 170.81)--(23.15, 171.94)--(24.93, 172.46)--(17.69, 173.06)--(13.37, 172.69)--(3.71, 172.63)--(13.68, 173.98)--(0, 174.21);
landmass5 := (124.62, 19.08)--(123.98, 19.56)--(126.36, 20.09)--(127.24, 20.59)--(124.98, 23.19)--(130.77, 29.28)--(135.8, 28.99)--(137.84, 25.04)--(141.86, 24.34)--(146.13, 22.18)--(157, 19.51)--(155.91, 17.98)--(155.55, 16.8)--(159.25, 15.48)--(158.44, 14.93)--(159.9, 14.05)--(157.97, 12.46)--(159.52, 10.77)--(159.19, 10.3)--(167.06, 8.42)--(159.95, 8.2)--(156.76, 8.73)--(153.02, 7.87)--(131.44, 7.07)--(130.82, 7.37)--(127.12, 7.44)--(116.94, 8.32)--(111.03, 9.65)--(105.32, 11.75)--(107.03, 12.86)--(108.98, 13.43)--(112.23, 14.12)--(119.83, 14.46)--(121.54, 15.67)--(120.54, 15.91)--(122.91, 16.69)--(121.82, 17.23)--(124.62, 19.08);
landmass6 := (307.49, 56.47)--(307.06, 57.11)--(308.22, 57.5)--(310.39, 56.79)--(310.75, 57.43)--(312.59, 56.96)--(313.38, 56.17)--(316.84, 55.24)--(317.14, 55.51)--(319.46, 54.34)--(320.21, 51.88)--(320.32, 48.77)--(319.72, 48.29)--(319.35, 47.51)--(322.01, 46.93)--(323.56, 45.58)--(319.79, 44.82)--(319.05, 44.6)--(318.2, 46.06)--(317.67, 48.01)--(318.7, 48.63)--(315.73, 52.34)--(311.83, 54.71)--(307.49, 56.47);
landmass7 := (172.44, 33.8)--(171.76, 34.44)--(174.73, 35.22)--(173.52, 37.33)--(172.83, 38.17)--(172.56, 39.98)--(179.32, 38.52)--(178.63, 37.06)--(175.74, 33.85)--(176.19, 30.59)--(172.04, 32.2)--(171.52, 33.36)--(172.44, 33.8);
landmass8 := (222.04, 111.1)--(224.53, 115.35)--(228.6, 106.45)--(226.81, 103.35)--(222.04, 111.1);
%
% This macro draws contures for a landmass on globe centered on lon,lat
%
vardef drawLandMass (expr p, lon, lat) =
save i, j, k, l, horizon, currentPoint, horizonTimes, outHorizon, inHorizon, visibleContours, pathNumber, horizonArc, arcTimes, resultPath;
path resultPath, visibleContours[], horizon, horizonArc;
pair currentPoint, horizonTimes[];
numeric pathNumber, outHorizon, arcTimes[];
horizon := fullcircle scaled 2;
pathNumber := 0;
outHorizon := -1;
inHorizon := -1;
% In the following loop visible segments of landmass and points of
% horizon-crossing are calculated
% visibleContours are just what they are called and horizonTimes are
% times on horizon circle where visible segment should cross it
for i := 0 upto length(p):
currentPoint := pointOnGlobe (point i of p, lon, lat);
if (horizonOnGlobe (point i of p, lon, lat) < 0):
if (unknown visibleContours[pathNumber]):
visibleContours[pathNumber] := currentPoint;
if (i > 0):
outHorizon := xpart(horizon intersectiontimes ((0,0) -- (findHorizonPoint (subpath(i-1, i) of p, lon, lat, 0) scaled 5)));
if (outHorizon < inHorizon): outHorizon := outHorizon + 8; fi;
horizonTimes[pathNumber - 1] := (inHorizon, outHorizon);
fi;
else:
visibleContours[pathNumber] := visibleContours[pathNumber] -- currentPoint;
fi;
else:
if (known visibleContours[pathNumber]):
pathNumber := pathNumber + 1;
inHorizon := xpart(horizon intersectiontimes ((0,0) -- (findHorizonPoint (subpath(i, i-1) of p, lon, lat, 0) scaled 5)));
fi;
fi;
endfor;
if (unknown visibleContours0):
resultPath := (1,0);
else:
if (unknown visibleContours[pathNumber]): pathNumber := pathNumber - 1; fi;
if (unknown horizonTimes[-1]):
if (pathNumber > 0):
visibleContours0 := visibleContours[pathNumber] -- visibleContours0;
fi;
pathNumber := pathNumber - 1;
else:
if (ypart(horizonTimes[-1]) < inHorizon):
horizonTimes[pathNumber] := (inHorizon, ypart(horizonTimes[-1]) + 8);
else:
horizonTimes[pathNumber] := (inHorizon, ypart(horizonTimes[-1]));
fi;
fi;
% In these loops horizon arcs directions should be handled
% The idea is that when we have path with no self-intersections arcs should
% not cross one another, these conflicts are resolved here
% It's important to note, that horizon-time detecting algorithm is
% not absolutely precise for now, so at times in will work incorrect.
for i := 0 upto pathNumber - 1:
for j := i + 1 upto pathNumber:
l := 0;
for k := -8, 0, 8:
if (xpart(horizonTimes[j]) > xpart(horizonTimes[i]) + k) and (xpart(horizonTimes[j]) < ypart(horizonTimes[i]) + k): l := l + 1; fi;
if (xpart(horizonTimes[i]) > xpart(horizonTimes[j]) + k) and (xpart(horizonTimes[i]) < ypart(horizonTimes[j]) + k): l := l + 1; fi;
endfor;
if (l > 1): horizonTimes[j] := (xpart(horizonTimes[j]) + 8, ypart(horizonTimes[j])); fi;
endfor;
endfor;
% In the following loop previously calculated segments of a landmass and
% arcs of the horizon are sewed together in order
resultPath := visibleContours0
for i := 1 upto pathNumber + 1:
if (known horizonTimes[i-1]): -- (subpath horizonTimes[i-1] of horizon) fi
if (i <= pathNumber): -- visibleContours[i] fi
endfor
fi;
resultPath -- cycle
enddef;
% This thing just converts coordinates on globe rotated by lon, lat to screen coordinates
vardef pointOnGlobe (expr p, lon, lat) =
(cosd(lon + xpart(p)) * sind(ypart(p)),
cosd(ypart(p)) * cosd(lat)
+ sind(lon + xpart(p)) * sind(ypart(p)) * sind(lat))
enddef;
% This one is needed to check if point on globe is in view
vardef horizonOnGlobe (expr p, lon, lat) =
(sind(lon + xpart(p)) * cosd(lat) * sind(ypart(p))) - (cosd(ypart(p)) * sind(lat))
enddef;
% This macro calculates horizon crossing point with given precision (recursion depth).
% Likely, this could be done analytically, though.
vardef findHorizonPoint (expr p, lon, lat, i) =
save selecthalf, returnpoint;
pair selecthalf, returnpoint;
if (horizonOnGlobe (point 1/2 of p, lon, lat) < 0):
selecthalf := (0, 1/2);
else:
selecthalf := (1/2, 1);
fi;
if (i < 5):
returnpoint := findHorizonPoint (subpath selecthalf of p, lon, lat, i + 1)
else:
returnpoint := pointOnGlobe (point 1/2 of p, lon, lat);
fi;
returnpoint
enddef;
vardef globe (expr s, lon, lat) =
save i, p, lm;
picture p[];
path lm;
begingroup
save lightnessPP;
vardef lightnessPP (expr v) =
1/2v
enddef;
p1 := image(draw sphereLat(2s, lat));
vardef lightnessPP (expr v) =
if (abs(cos(sphlat)) > 7/8 + uniformdeviate (1/20)):
1/4v
else:
1/3v + 2/3
fi
enddef;
p2 := image(draw sphereLat(2s, lat));
endgroup;
image(
draw fullcircle scaled 2s withpen thinpen;
for i := 1 upto 8:
lm := drawLandMass(landmass[i], lon + 90, lat) scaled s;
p3 := p2;
clip p3 to lm;
draw p3;
thinBrushGenerate (lm,
offsetPathTemplate (lm, 0) (
2/3minStrokeWidth + 1/3minStrokeWidth
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of lm), arcsin(abs(point offsetPathTime of lm)/2s))
), 0, point offsetPathTime of lm
)
), 0);
p1 := p1 maskedWith lm;
endfor;
draw p1;
)
enddef;
%
% This macro draws an eye pointed in the direction a (in degrees)
% Eye is opened at random angle and pupil is scaled randomly by design
% Scaling below some level, dependent on minStrokeWidth, simplifies image
%
eyescale := 1/2cm;
vardef eye (expr a) =
save s, eyelids, pupil, eyeball, eyelash, loopstep, p, o;
path eyelids[], pupil[], eyelash;
pair p[];
picture eyeball;
numeric s, loopstep;
o := 10 + (15/(1 + 2**(3/2normaldeviate)));
s := eyescale;
p1 := (-3/4s, 0);
pupil1 := ((subpath (-2, 2) of fullcircle xscaled 3/5) .. (subpath (3, 5) of fullcircle xscaled 2/5) .. cycle) scaled 3/5s;
pupil2 := fullcircle scaled (1/3s + uniformdeviate(1/5s)) xscaled 1/3;
p2 := ((p1 -- ((1/2s, 0) rotatedabout (p1, o - 5))) intersectionpoint (subpath (0, 4) of pupil1));
eyelids1 := (p1 shifted (0, -1/16s)){dir(1/3o)} .. {dir(0)}p2;
eyelids2 := p1 {dir(-1/3o)} .. {dir(0)}((1/6s, 0) rotatedabout (p1, -2/3o - 5));
eyelids2 := subpath(xpart(eyelids2 intersectiontimes eyelids1), length(eyelids2)) of eyelids2;
eyelids3 := (p1){dir(2/3o)} .. tension 3/2 .. {dir(o-1/3o)}p2 rotatedabout (p1, 1/3o + 7);
eyelids3 := subpath (1/8, length(eyelids3)) of eyelids3;
eyelids4 := (p1 shifted (0, -1/16s)) .. {dir(1/4o - 2/3o)}((1/6s, -1/6s) rotatedabout (p1, -2/3o - 5));
eyelids4 := subpath (1/2, length(eyelids4)) of eyelids4;
loopstep := length(eyelids1)/20;
if (arclength(subpath(0, loopstep) of eyelids1) < 5minStrokeWidth):
loopstep := arctime (5minStrokeWidth) of eyelids1;
fi;
eyeball := image(
if (5loopstep <= 1):
draw pupil1 withpen thinpen;
fill pupil2;
for i := 0 step 5loopstep until length pupil1:
draw brush ((point i of pupil1) -- (point i + 6 of pupil2 scaled 5/6 yscaled 1/2))(cos(offsetPathLength*1/2pi)*2minStrokeWidth);
endfor;
else:
fill pupil1;
fi;
);
clip eyeball to (eyelids1{(1, 0)} .. (s, 0) .. {-1, 0}reverse(eyelids2) -- cycle);
eyeball := eyeball maskedWith ((fullcircle scaled (1/4s) xscaled 1/2 rotated 2 shifted (1/12s, 0) rotatedabout (p1, 1/3o + 2)));
image(
draw brush (eyelids1 .. tension 2.5 .. reverse (eyelids3))((1-offsetPathLength)*2minStrokeWidth);
draw brush (eyelids2)((offsetPathLength)*2minStrokeWidth);
draw brush (eyelids4)(sin(offsetPathLength*pi)*minStrokeWidth);
draw eyeball;
for i := length(eyelids1) step -loopstep until 0:
eyelash := (point i of eyelids1) {dir(angle(direction i of eyelids1) - 60 + 50*(i/length(eyelids1)))}
.. (point i of eyelids1) shifted (1/16s + (i/length(eyelids1))*1/4s, (i/length(eyelids1))*1/5s);
if (arclength(eyelash) > 2/3minDashStrokeLength):
draw brush (eyelash shifted ((-1/32s, uniformdeviate(1/12s)) scaled ((length(eyelids1)-i)/length(eyelids1))) )(minStrokeWidth + (1-offsetPathLength)*minStrokeWidth);
fi;
endfor;
for i := length(eyelids2) step -3/2loopstep until 0:
eyelash := (point i of eyelids2) {dir(angle(direction i of eyelids2) + 20 - 40*(i/length(eyelids2)))}
.. (point i of eyelids2) shifted (1/16s + (i/length(eyelids2))**2*1/7s, 1/16s - (i/length(eyelids2))*1/5s);
if (arclength(eyelash) > minDashStrokeLength):
draw brush (eyelash shifted (-1/32s, -uniformdeviate(1/24s)))(minStrokeWidth + (1-offsetPathLength)*minStrokeWidth);
fi;
endfor;
) if cosd(a) < 0: yscaled -1 rotated (a) else: rotated a fi
enddef;
%
% This macro draws solid surface
%
vardef solidSurface (expr p) =
save q, s, d, stripes;
path q, s;
picture stripes;
q := offsetPath(p) (-1/4cm);
s := p -- reverse(q) -- cycle;
image(
draw solid (s, 45, 0);
draw p withpen thinpen;
)
enddef;
vardef solid (expr p, a, t) =
save stripes, stripeskind, d, i, j, c, strokeVariation;
picture stripes, stripeskind;
pair c;
stripes := image(
d1 := abs(ulcorner(p rotated (90 - a)) - urcorner(p rotated (90 - a)));
d2 := abs(ulcorner(p rotated (90 - a)) - llcorner(p rotated (90 - a)));
stripeskind := dashpattern (on 1mm);
c := 1/2[ulcorner(p rotated (90 - a)), lrcorner(p rotated (90 - a))] rotated (a - 90);
for i:= 0 step (3/2shadingDensity)/d1 until 1:
if (t = 1):
strokeVariation := uniformdeviate(1)-1/2;
j := round(i*d1/(3/2shadingDensity));
if (j mod 4) = 0:
stripeskind := dashpattern (on (8-strokeVariation)*shadingDensity off (4+strokeVariation)*shadingDensity);
fi;
if ((j mod 4) = 1) or ((j mod 4) = 3):
stripeskind := dashpattern (off 1shadingDensity on (6-strokeVariation)*shadingDensity off (5+strokeVariation)*shadingDensity);
fi;
if (j mod 4) = 2:
stripeskind := dashpattern (on 0 off 12shadingDensity);
fi;
fi;
draw ((dir(a) scaled 1/2(d2-uniformdeviate(3shadingDensity))) -- (dir(a + 180) scaled 1/2d2)) shifted c shifted i[dir(a + 90) scaled 1/2d1, dir(a -90) scaled 1/2d1] withpen thinpen
dashed stripeskind;
endfor;
);
clip stripes to p;
stripes
enddef;
%
% Returns a picture of shaded gradient inside shape p at an angle a
%
vardef gradientShade (expr p, a) =
save stripes, stripeshd, d, i, j, s;
picture stripes;
path s;
stripes := image(
d1 := abs(ulcorner(p rotated (90 - a)) - urcorner(p rotated (90 - a)));
d2 := abs(ulcorner(p rotated (90 - a)) - llcorner(p rotated (90 - a)));
for i:= 0 step (shadingDensity)/d1 until 1:
s := ((dir(a) scaled 1/2d2) -- (dir(a + 180) scaled 1/2d2)) shifted 1/2[ulcorner(p), lrcorner(p)] shifted i[dir(a + 90) scaled 1/2d1, dir(a -90) scaled 1/2d1];
stripeshd := 1/4 + 3/4i;
draw brush(s)(minStrokeWidth*stripeshd);
endfor;
);
clip stripes to p;
stripes
enddef;
%
% This one draws a spring between points p and q with n steps
% if stretched more than possible, displayed as a straight line
% if compressed too much, displayed as having less steps
%
springwidth := 1/8cm;
vardef spring (expr p, q, n) =
save sp, ss, t, x, springstep, springcoef, springsegment;
transform t;
path ss[];
pair sp;
picture springsegment;
springstep := (arclength(p--q) - 2springwidth)/(n+1);
if (springstep < 6minStrokeWidth): springstep := arclength(p--q)/round(arclength(p--q)/6minStrokeWidth); fi;
if (springstep < (springwidth*pi)):
springcoef := (1-(springstep/(springwidth*pi)));
else:
springcoef := 0;
fi;
image(
for i := 0 step 30 until 360:
sp := ((cosd(i - 90), sind(i - 90)) scaled springwidth xscaled 1/4 yscaled springcoef) shifted (springstep*(i/360) - 1/2springstep, 0);
if (i = 0):
ss1 := sp;
else:
ss1 := ss1 -- sp;
fi;
endfor;
ss2 := subpath (0, 6) of ss1;
ss3 := subpath (6, 12) of ss1;
ss4 := subpath (0, ypart(ss2 shifted (-3minStrokeWidth, 0) intersectiontimes ss3)) of ss3;
x := ypart(ss2 shifted (3minStrokeWidth, 0) intersectiontimes ss3);
if (x > 0):
ss5 := subpath (x, length(ss3)) of ss3;
else:
ss5 := point length(ss3) of ss3;
fi;
if (xpart(llcorner(ss4)) - 3minStrokeWidth < xpart(urcorner(ss2 shifted (-springstep, 0)))):
x := ypart((subpath (3, 6) of ss2 shifted (-springstep + 3minStrokeWidth, 0)) intersectiontimes ss4);
if (x > 0):
ss6 := subpath (0, x) of ss4;
else:
ss6 := point 0 of ss4;
fi;
x := ypart((subpath (0, 3) of ss2 shifted (-springstep + 3minStrokeWidth, 0)) intersectiontimes ss4);
if (x > 0):
ss7 := subpath (x, length(ss4)) of ss4;
else:
ss7 := point length(ss4) of ss4;
fi;
fi;
springsegment := image(
draw brush (ss2)(minStrokeWidth + sin(offsetPathLength*pi)*minStrokeWidth);
if (unknown(ss6)):
if (arclength(ss4) > minStrokeWidth): draw ss4 withpen thinpen; fi;
else:
if (arclength(ss6) > minStrokeWidth): draw ss6 withpen thinpen; fi
if (arclength(ss7) > minStrokeWidth): draw ss7 withpen thinpen; fi
fi;
if (arclength(ss5) > minStrokeWidth): draw ss5 withpen thinpen; fi;
);
ss8 := (ss2 shifted (-2minStrokeWidth, 0)) rotated angle(q-p) shifted ((springwidth + 1/2springstep)/arclength(p--q))[p, q];
for i := springwidth + 1/2springstep step springstep until arclength(p--q) - springwidth - 1/2springstep + 1:
t := identity rotated angle(q-p) shifted (i/arclength(p--q))[p, q];
draw springsegment transformed t;
if (i <= springwidth + 1/2springstep + 2/3springwidth):
ss9 := ss3 transformed t;
ss10 := subpath (
xpart(ss9 intersectiontimes (subpath (0, 3) of ss8)),
xpart(ss9 intersectiontimes (subpath (3, 6) of ss8))
) of ss9;
if (arclength(ss10) > minStrokeWidth): draw ss10 withpen thinpen; fi;
fi;
endfor;
draw brush (((springwidth)/arclength(p--q))[p, q] shifted (dir(angle(p-q) + 90)*springwidth * springcoef){(p-q)} .. {(p-q)}p)(minStrokeWidth);
draw brush (((springwidth)/arclength(p--q))[q, p] shifted (dir(angle(p-q) + 90)*springwidth * springcoef){(q-p)} .. {(q-p)}q)(minStrokeWidth);
)
enddef;
%
% This macro draws some kind of weight. Not very nice one at the moment
%
vardef weight.h (expr h) =
save auricle, q, r;
path q[];
auricle.d := 2mm;
auricle.t := 2shadingDensity;
r := 2/5(h-auricle.d);
image(
q0 := offsetPathSubdivide((0, -h) -- (0, -auricle.d - 1/6h));
q1 := offsetPathTemplate(q0, 0)(r-(offsetPathLength*(1/8r)));
q2 := offsetPathGenerate (q0, q1, 0);
q3 := offsetPathGenerate (q0, q1 yscaled -1, 0);
tubeGenerate (q2, q3, q1, 0);
draw reverse(q2) -- q3 withpen thinpen;
q0 := offsetPathSubdivide((0, -auricle.d - 1/6h) -- (0, -auricle.d));
q1 := offsetPathTemplate(q0, 0)(2/8r + 5/8r*(sqrt(1-offsetPathLength**2)));
q2 := offsetPathGenerate (q0, q1, 0);
q3 := offsetPathGenerate (q0, q1 yscaled -1, 0);
tubeGenerateAlt (q2, q3, q1);
draw q3 -- reverse(q2) withpen thinpen;
q5 := offsetPathSubdivide(point 0 of q3 -- point 0 of q2);
q6 := tube.e((0,0) -- (0, -3/2auricle.t))(1/2auricle.t);
draw image(
q4 := (((0, -1/2) {dir(90)} .. (1/2, 1/2) .. (0, 1) .. {dir(-90)}(-1/2, 1/4)) shifted (0, -1)) scaled 2/3auricle.d;
draw shadedEdge(tube.e(q4) (1/4auricle.t)) shifted (0, -1/2auricle.t);
) maskedWith (q3 -- reverse(q2) -- (q2 yscaled 0 shifted (0, -h)) -- (reverse(q3) yscaled 0 shifted (0, -h)) -- cycle)
maskedWith q6;
draw shadedEdge(q6);
)
enddef;
vardef weight.s (expr h) =
save q,r;
path q[];
r := 2/5h;
image(
q0 := offsetPathSubdivide((0, 0) -- (0, h-2/3r));
q1 := offsetPathTemplate(q0, 0)(r);
q2 := offsetPathGenerate (q0, q1, 0);
q3 := offsetPathGenerate (q0, q1 yscaled -1, 0);
tubeGenerate (q2, q3, q1, 0);
draw reverse(q2)--q3 withpen thinpen;
q0 := offsetPathSubdivide((0, h-2/3r) -- (0, h - 1/8r));
q1 := offsetPathTemplate(q0, 0)(r-sqrt(1- (1- offsetPathLength*2)**2)*1/3r - 1/6r*offsetPathLength);
q2 := offsetPathGenerate (q0, q1, 0);
q3 := offsetPathGenerate (q0, q1 yscaled -1, 0);
tubeGenerateAlt (q2, q3, q1);
draw q2 withpen thinpen;
draw q3 withpen thinpen;
q0 := offsetPathSubdivide((0, h-1/8r) -- (0, h));
q1 := offsetPathTemplate(q0, 0)((r-1/6r)+sqrt(1- (1- offsetPathLength*2)**2)*1/16r);
q2 := offsetPathGenerate (q0, q1, 0);
q3 := offsetPathGenerate (q0, q1 yscaled -1, 0);
tubeGenerateAlt (q2, q3, q1);
draw q2 --reverse(q3) withpen thinpen;
)
enddef;
%
% This macro makes a lens-shaped clockwise path with radii
% r = (left radius, right radius), thickness t and diameter d
%
vardef lens (expr r, t, d, s) =
save p, q, m, c;
pair c[];
path p[], q[];
if (xpart(r) = infinity):
p1 := (0, d) -- (0, -d);
else:
p1 := subpath (2, 6) of fullcircle scaled 2xpart(r);
fi;
if (ypart(r) = infinity):
p2 := (0, d) -- (0, -d);
else:
p2 := subpath (-2, 2) of fullcircle scaled 2ypart(r);
fi;
q1 := (min(xpart(ulcorner(p1)), xpart(ulcorner(p2))) - 1, -1/2d)--(max(xpart(urcorner(p1)), xpart(urcorner(p2))) + 1,-1/2d);
q2 := q1 shifted (0, d);
q3 := q1 shifted (0, 1/2d);
c1 := p1 intersectiontimes q1;
c2 := p1 intersectiontimes q2;
c3 := p2 intersectiontimes q2;
c4 := p2 intersectiontimes q1;
if (xpart(c1) > 0):
p1 := subpath (xpart(c1), xpart(c2)) of p1;
fi;
if (xpart(c3) > 0):
p2 := subpath (xpart(c3), xpart(c4)) of p2;
fi;
p1 := p1 shifted (-xpart(point xpart(p1 intersectiontimes q3) of p1), 0);
p2 := p2 shifted (-xpart(point xpart(p2 intersectiontimes q3) of p2) + t, 0);
reverse(p1--p2--cycle) scaled s
enddef;
%
% This one returns a picture of a pulley with diameter d and it's support rotated
% at angle a. pulleyOutline path changes every time
%
path pulleyOutline;
numeric pulleySupportSize;
pulleySupportSize := 2/3;
vardef pulley (expr d, a) =
save pw, p, r;
picture pw;
path p[];
r1 := 3/5d;
r2 := 1/6d;
image(
p0 := fullcircle scaled d;
p1 := (subpath (3, 9) of fullcircle) scaled r1;
p0 := subpath (
xpart(p0 intersectiontimes ((point 0 of p1) -- (point 0 of p1 shifted (0, d)))),
8 + xpart(p0 intersectiontimes ((point 0 of reverse(p1)) -- (point 0 of reverse(p1) shifted (0, d)))))
of p0;
p0 := ((xpart(point 0 of reverse(p0)), pulleySupportSize*d) -- (xpart(point 0 of p0), pulleySupportSize*d) -- p0 -- cycle) rotated a;
pulleyOutline := p0;
p1 := (p1 -- (xpart(point 0 of reverse(p1)), pulleySupportSize*d) -- (xpart(point 0 of p1), pulleySupportSize*d) -- cycle) rotated a;
pw := pulleyWheel(d) maskedWith p1;
draw pw;
draw p1 withpen thinpen;
draw shadedEdge(reverse(fullcircle) scaled r2);
)
enddef;
vardef pulleyWheel (expr d) =
save pw, r, i;
picture pw;
r1 := 7/9d;
r2 := 8/9d;
r3 := 3/5d;
r4 := 1/8d;
if (r2-r1) > shadingDensity:
pw := image(
for i := r3 step 2shadingDensity until r2:
if (i <= r1) or (i >= r2):
thinBrushGenerate (fullcircle scaled i,
offsetPathTemplate (fullcircle scaled i, 0) (
2/3minStrokeWidth + minStrokeWidth
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of fullcircle), arcsin(i/4d))
), 0, point offsetPathTime of fullcircle scaled i
)
), 0);
else:
thinBrushGenerate (fullcircle scaled i,
offsetPathTemplate (fullcircle scaled i, 0) (
1/4minStrokeWidth + minStrokeWidth
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of fullcircle) + pi, arcsin(1/2))
), 0, point offsetPathTime of fullcircle scaled i
)
), 0);
fi;
endfor;
draw brush (fullcircle scaled r1)(minStrokeWidth);
draw brush (fullcircle scaled r2)(minStrokeWidth);
draw fullcircle scaled d withpen thinpen;
draw fullcircle scaled r3 withpen thinpen;
draw shadedEdge(reverse(fullcircle) scaled r4);
);
else:
pw := image(
draw shadedEdge (reverse(fullcircle) scaled r1);
draw fullcircle scaled d withpen thinpen;
draw fullcircle scaled r4 withpen thickpen;
)
fi;
pw
enddef;
%
% This macro draws a wheel
%
vardef wheel (expr d, a) =
save pc, p, r, i;
picture pc[];
path p[];
r1 := 2/8d;
r2 := d-6minStrokeWidth;
r3 := 7/8d-2shadingDensity;
if r1 > 3shadingDensity:
pc1 := image(
for i := r1 step 2shadingDensity until r3:
thinBrushGenerate (fullcircle scaled i,
offsetPathTemplate (fullcircle scaled i, 0) (
2/3minStrokeWidth + minStrokeWidth
* normalVectorToLightness(
sphereAnglesToNormalVector(
(angleRad(point offsetPathTime of fullcircle) + pi, arcsin(i/4d))
), 0, point offsetPathTime of fullcircle scaled i
)
), 0);
endfor;
draw sphere.c(r1);
draw shadedEdge (reverse(fullcircle) scaled r3);
);
else:
pc1 := image(
draw shadedEdge (fullcircle scaled d);
draw shadedEdge (fullcircle scaled r1);
draw shadedEdge (reverse(fullcircle) scaled r2);
);
fi;
pc2 := image(fill fullcircle scaled d;) maskedWith (fullcircle scaled r2);
pc3 := image(
p1 := reverse((1/3d, 1/4d) -- (-1/3d, 1/4d) -- (-1/2d, 1/2d) -- (1/2d, 1/2d) -- cycle) rotated a;
draw p1 withpen thinpen;
draw gradientShade(p1, 180);
);
pc3 := pc3 maskedWith (fullcircle scaled d);
image(
draw pc1;
draw pc2;
draw pc3;
)
enddef;
%
% This one roughly finds the center of mass for a closed shape
%
vardef shapeCenterOfMass (expr p) =
save i, a, aTotal, q;
numeric i, a, aTotal;
pair rv, q[];
q0 := point 0 of p;
aTotal := 0;
rv := (0, 0);
for i := 1 step 1 until (length(p)-2):
q1 := point i of p;
q2 := point (i+1) of p;
if (xpart(q1-q0) = 0):
a := (abs(ypart(q1-q0)/cm)*abs(xpart(q2-q0)/cm))/2;
elseif (ypart(q1-q0) = 0):
a := (abs(xpart(q1-q0)/cm)*abs(ypart(q2-q0)/cm))/2;
elseif (abs(xpart(q1-q0)) > abs(ypart(q1-q0))):
begingroup
save qA;
pair qA;
qA = whatever[q0, q0 + (0,1)];
qA = whatever[q2, q2 + (q1-q0)];
a := (abs(ypart(qA-q0)/cm)*abs(xpart(q1-q0)/cm))/2;
endgroup
else:
begingroup
save qA;
pair qA;
qA = whatever[q0, q0 + (1,0)];
qA = whatever[q2, q2 + (q1-q0)];
a := (abs(xpart(qA-q0)/cm)*abs(ypart(q1-q0)/cm))/2;
endgroup
fi;
%a := 2;
aTotal := aTotal + a;
rv := rv + ((q0 + q1 + q2) scaled (a/3));
endfor;
if (aTotal > 0):
rv := rv scaled (1/aTotal);
else:
rv := (0, 0);
fi;
rv
enddef;
%
% This macro is for drawing shaded flat surfaces
%
vardef flatSurface@#(expr surfPath, normalVector, hatchAngle) =
save p, aHatch, hatchImage, surfLight, totalHeight, totalWidth, lineDensity, distFromEdge, hatchLength;
path p, aHatch;
picture hatchImage;
numeric distFromEdge, hatchLength;
surfLight := normalVectorToLightness(normalVector, 0, (0, 0));
p := surfPath rotated -hatchAngle;
if (str @# = "") or (str @# = "hatches"):
lineDensity := shadingDensity;
elseif (str @# = "stipples"):
lineDensity := stippleShadingDensity;
fi;
hatchImage := image(
totalHeight := abs(ypart(llcorner(p)) - ypart(urcorner(p)));
totalWidth := abs(xpart(llcorner(p)) - xpart(urcorner(p)));
%fill surfPath withcolor (surfLight, surfLight, surfLight);
for i := ypart(llcorner(p)) step (totalHeight/round(totalHeight/lineDensity)) until ypart(urcorner(p)):
distFromEdge := abs((i - ypart(llcorner(p)))/(ypart(urcorner(p))-ypart(llcorner(p))));
aHatch := (xpart(llcorner(p)),i) -- (xpart(lrcorner(p)), i);
aHatch := (point xpart(aHatch firstIntersectionTimes p) of aHatch) --
(point xpart(reverse(aHatch) firstIntersectionTimes p) of reverse(aHatch));
hatchLength := arclength(aHatch);
if arclength(aHatch) > 0:
aHatch := aHatch rotated hatchAngle;
aHatch := pathSubdivide(aHatch,3+round(uniformdeviate(2)));
if (str @# = "") or (str @# = "hatches"):
draw brush(aHatch)(
(2minStrokeWidth*(1/15+surfLight))*
(1-6/7sqrt(
sin(offsetPathLength*pi)*(hatchLength/totalWidth)*
sin(distFromEdge*pi)
))
);
elseif (str @# = "stipples"):
draw thinBrush.stipples(aHatch)(
(stippleSize*(1/15+surfLight))*
(1-5/6sqrt(
sin(offsetPathLength*pi)*(hatchLength/totalWidth)*
sin(distFromEdge*pi)
))
);
fi;
fi;
endfor;
);
clip hatchImage to surfPath;
image(
draw hatchImage;
)
enddef;
%
% These macros are for drawing wood texture. A bunch of wood-related global
% variables are also here.
%
woodBlockRes := 1/3mm;
woodBlockDetail := 1/5;
woodBlockYRdensity := 1/24;
woodKnotAngle := 1/3pi;
woodKnotRadius := 1/8cm;
woodKnotDensity := 2cm;
% wField returns a value on a scalar field, that is surface of year ring
vardef wField (suffix wk)(expr i, j, cs, nwk)=
save x, y, csx, csy, ba, a, r, outputValue, k, bd;
csx := xpart(cs);
csy := ypart(cs);
x := i*woodBlockDetail;
y := j*woodBlockDetail;
bd := 1/2(woodKnotRadius/woodBlockRes)*woodBlockDetail;
outputValue := 0;
for k := 0 upto nwk:
r := (abs(((x, y) shifted -wk[k]) yscaled (sin(woodKnotAngle))));
a := angleRad((x, y) shifted -wk[k]);
if (r >= 2bd) and (1/2r-bd < 8):
outputValue := outputValue + ((sqrt(1/2r-bd)/(3**(1/2r-bd)))*(sin(woodKnotAngle)*(1/2sin(-a)) + 1)*(1/3 + 2/3abs(cos(a))))**2;
elseif (r < 2bd):
outputValue := outputValue - 10sqrt(1-(r/2bd)**4);
fi;
outputValue := outputValue + 1/64cos(2pi*1/30r);
endfor;
outputValue := outputValue - (i*woodBlockYRdensity)/5 + 1/32sin(pi*i/xpart(cs) + 4pi*j/ypart(cs));
outputValue
enddef;
% woodBlock generates coordinates of knots, calls wField to
% generate matrix of heights (one for all years, that's simplification, of course)
% and then calls isoLines for each year ring (shifting values in matrix by woodBlockYRdensity)
vardef woodBlock (expr w, l) =
save wood, wKnot, nwKnot, p, q, i, j, k, cl, tr, wS, lS;
numeric wood[][];
pair wKnot[];
path q;
picture p;
if (l > w):
wS := round(w/woodBlockRes);
lS := round(l/woodBlockRes);
else:
wS := round(l/woodBlockRes);
lS := round(w/woodBlockRes);
fi;
image(
p := image(
i := -1;
j := 0;
forever:
wKnot[-1] := (uniformdeviate(wS*woodBlockDetail + woodKnotRadius*woodBlockDetail/woodBlockRes), uniformdeviate(lS*woodBlockDetail + woodKnotRadius*woodBlockDetail/woodBlockRes)) shifted (-1/2woodKnotRadius*woodBlockDetail/woodBlockRes, -1/2woodKnotRadius*woodBlockDetail/woodBlockRes);
cl := 0;
if (i > -1):
for k := 0 step 1 until i:
if (woodBlockRes*abs(wKnot[k]-wKnot[-1])/woodBlockDetail) < woodKnotDensity:
cl := 1;
fi;
endfor;
fi;
if cl > 0:
j := j + 1;
else:
i := i + 1;
wKnot[i] := wKnot[-1];
for tr := 2woodKnotRadius step -8minStrokeWidth until 2minStrokeWidth:
draw brush(
((fullcircle scaled tr yscaled (1/sin(woodKnotAngle))) shifted (woodBlockRes*wKnot[i]/woodBlockDetail))
)(
minStrokeWidth*(1/2+abs(sin(woodKnotAngle))*abs(sin(angleRad(point offsetPathTime of fullcircle))))
);
endfor;
fi;
exitif (j >= 10) or (i >= 1/2((w/cm)*(l*cm))/(woodKnotDensity/cm));
endfor;
nwKnot := i;
for i := 0 step 1 until wS:
k := uniformdeviate(1/8woodBlockYRdensity);
for j := 0 step 1 until lS:
wood[i][j] := wField (wKnot)(i, j, (wS, lS), nwKnot) + k;
endfor;
endfor;
for i := -wS/(80woodBlockYRdensity) - 2 upto wS/(80woodBlockYRdensity) + 2:
draw isoLines(wood)((wS, lS), woodBlockYRdensity*i + uniformdeviate(1/8woodBlockYRdensity), woodBlockRes);
endfor;
);
q := (0, 0) -- (w, 0) -- (w, l) -- (0, l) -- cycle;
if (w > l): q := q xscaled -1 rotated -90; fi;
clip p to q;
draw p;
draw q withpen thinpen;
) if (l < w): xscaled -1 rotated -90 fi
enddef;
% woodenThing fits a woodBlock into thingOutline at a given woodAngle
vardef woodenThing (expr thingOutline, woodAngle) =
save shiftedThing, thingOrigin, thingWoodBlock;
path shiftedThing;
pair thingOrigin;
picture thingWoodBlock;
shiftedThing := thingOutline rotated -woodAngle;
thingOrigin := llcorner(shiftedThing);
shiftedThing := shiftedThing shifted -thingOrigin;
thingWoodBlock := woodBlock(xpart(urcorner(shiftedThing)), ypart(urcorner(shiftedThing)));
clip thingWoodBlock to shiftedThing;
thingWoodBlock := (thingWoodBlock shifted thingOrigin) rotated woodAngle;
image(
draw thingOutline withpen thinpen;
draw thingWoodBlock;
)
enddef;
%
% This part is related to knots
%
% lists are only used in knots so far.
% The following two macros are taken from byrne.mp
vardef appendList@#(suffix listName)(expr valueToAdd, whereToAdd, omitDuplicates) =
save v, valueExists;
string v;
boolean valueExists;
if str @# = "":
if not string listName:
string listName;
fi;
else:
if not string listName0:
string listName[];
fi;
fi;
if unknown listName@#:
listName@# := "";
fi;
valueExists := false;
if omitDuplicates:
valueExists := isInList@#(valueToAdd, listName)
fi;
if not valueExists:
if string valueToAdd:
v := valueToAdd;
else:
v := decimal(valueToAdd);
fi;
if length(listName@#) = 0:
listName@# := v;
else:
if (whereToAdd = 1):
listName@# := listName@# & ", " & v;
else:
listName@# := v & ", " & listName@#;
fi;
fi;
fi;
enddef;
vardef isInList@#(expr valueToLookFor)(suffix listName) =
save rv, i;
boolean rv;
rv := false;
if str @# = "":
if not string listName:
string listName;
fi;
else:
if not string listName0:
string listName[];
fi;
fi;
if unknown listName@#:
listName@# := "";
fi;
if string valueToLookFor:
forsuffixes i=scantokens(listName@#):
if (str i = valueToLookFor):
rv := true;
fi;
endfor;
else:
for i=scantokens(listName@#):
if (i = valueToLookFor):
rv := true;
fi;
endfor;
fi;
rv
enddef;
vardef sortList (expr listToSort, ascending) =
save nPre, nPost, pivot, isSorted, lastValue, preList, postList, rv;
numeric nPre, nPost, pivot;
boolean isSorted;
string preList, postList, rv;
nPre := 0;
nPost := 0;
isSorted := true;
if ascending:
lastValue := -infinity;
else:
lastValue := infinity;
fi;
for i=scantokens(listToSort):
if (unknown pivot):
pivot := i;
fi;
if ((i < pivot) and ascending) or ((i > pivot) and not ascending):
appendList (preList, i, -1, false);
nPre := nPre + 1;
else:
appendList (postList, i, -1, false);
nPost := nPost + 1;
fi;
if ((lastValue > i) and ascending) or ((lastValue < i) and not ascending):
isSorted := false;
fi;
lastValue := i;
endfor;
if isSorted:
rv := listToSort;
else:
if nPre > 1:
preList := sortList(preList, ascending);
fi;
if nPre > 0:
preList := preList & ", ";
else:
preList := "";
fi;
if nPost > 1:
postList := sortList(postList, ascending);
fi;
rv := preList & postList;
fi;
rv
enddef;
% When looking for intersections, knot is browsed with knotStepSize step.
% Affects nothing interesting.
numeric knotStepSize;
knotStepSize := 1/8;
vardef knotFromStrands (suffix knotName) =
save slidingSegment, timeToAdd, pathSegments, intTimes, tSegWidth, tSegStyle, numberOfSegments, segmentWidth, totalNumberOfSegments, intersections, intersectionsList, layerContents, layersList, b, e, n, layerContents, segmentPicture;
save shadowsEnabled, allShadowPaths, totalNumberOfShadows, numberOfShadows, shadowPath, tmpShadows;
boolean shadowsEnabled;
path shadowPath[], allShadowPaths[];
numeric timeToAdd, numberOfShadows, totalNumberOfShadows;
totalNumberOfShadows := -1;
tmpShadows := -1;
path inspectedPath, slidingSegment, pathSegments[];
pair intTimes[];
numeric intersections[], numberOfSegments[], segmentWidth[], totalNumberOfSegments, tSegWidth, b, e, n;
picture layerPicture[];
string layerContents[], segmentStyle[], tSegStyle;
totalNumberOfSegments := 0;
for sNa := 1 step 1 until knotName.nStrands:
inspectedPath := knotName.strandPath[sNa];
tSegWidth := knotName.strandWidth[sNa];
tSegStyle := knotName.strandStyle[sNa];
for sNb := 1 step 1 until knotName.nStrands:
for i := -knotStepSize step knotStepSize until length(knotName.strandPath[sNb]):
slidingSegment := subpath (i, i + 2knotStepSize) of knotName.strandPath[sNb];
intTimes0 := (inspectedPath firstIntersectionTimes slidingSegment);
intTimes1 := (reverse(inspectedPath) firstIntersectionTimes slidingSegment);
if (sNb = sNa):
if (xpart(intTimes0) >= i)
and (xpart(intTimes0) <= i + 2knotStepSize):
intTimes0 := (-1, -1);
fi;
if ((length(inspectedPath) - xpart(intTimes1)) >= i)
and ((length(inspectedPath) - xpart(intTimes1)) <= i + 2knotStepSize):
intTimes1 := (-1, -1);
fi;
if (i+knotStepSize >= length(inspectedPath)):
if (xpart(intTimes0) <= knotStepSize)
or (xpart(intTimes0) = length(inspectedPath)):
intTimes0 := (-1, -1);
fi;
if ((length(inspectedPath) - xpart(intTimes1)) <= knotStepSize)
or (xpart(intTimes1) = 0):
intTimes1 := (-1, -1);
fi;
fi;
if (i-knotStepSize <= length(inspectedPath)):
if (xpart(intTimes0) >= (length(inspectedPath) - knotStepSize))
or (xpart(intTimes0) = 0):
intTimes0 := (-1, -1);
fi;
if ((length(inspectedPath) - xpart(intTimes1)) >= (length(inspectedPath) - knotStepSize))
or (xpart(intTimes1) = length(inspectedPath)):
intTimes1 := (-1, -1);
fi;
fi;
fi;
timeToAdd := -1;
if ((ypart(intTimes0) > 0) and (ypart(intTimes0) < length(slidingSegment)))
and ((xpart(intTimes0) >= 0) and (xpart(intTimes0) <= length(inspectedPath)))
and ((sNb <> sNa) or ((xpart(intTimes0) < i) or (xpart(intTimes0) > i + 1))):
timeToAdd := xpart(intTimes0);
elseif (sNb = sNa):
if ((ypart(intTimes1) > 0) and (ypart(intTimes1) < length(slidingSegment)))
and ((xpart(intTimes1) >= 0) and (xpart(intTimes1) <= length(inspectedPath)))
and ((length(inspectedPath) - xpart(intTimes1) < i) or (length(inspectedPath) - xpart(intTimes1) > i + 1)):
timeToAdd := length(inspectedPath) - xpart(intTimes1);
fi;
fi;
if (timeToAdd >= 0):
if (timeToAdd = length(inspectedPath)):
timeToAdd := 0;
fi;
appendList[sNa](intersectionsList, round(timeToAdd*10)/10, 1, true);
fi;
endfor;
endfor;
if known intersectionsList[sNa]:
numberOfSegments[sNa] := 0;
for i := scantokens(sortList(intersectionsList[sNa], true)):
numberOfSegments[sNa] := numberOfSegments[sNa] + 1;
intersections[numberOfSegments[sNa]] := i;
endfor;
intersections[0] := 0;
if (not cycle knotName.strandPath[sNa]):
intersections[numberOfSegments[sNa] + 1] := length(knotName.strandPath[sNa]);
fi;
for i := 1 step 1 until numberOfSegments[sNa]:
totalNumberOfSegments := totalNumberOfSegments + 1;
if (i > 1):
b := 1/2[intersections[i-1], intersections[i]];
else:
if (not cycle knotName.strandPath[sNa]):
b := 0;
else:
b := 1/2[intersections[numberOfSegments[sNa]] - length(knotName.strandPath[sNa]), intersections[i]];
fi;
fi;
if (i < numberOfSegments[sNa]):
e := 1/2[intersections[i], intersections[i+1]];
else:
if (not cycle knotName.strandPath[sNa]):
e := length(inspectedPath);
else:
e := 1/2[intersections[i], intersections[1] + length(knotName.strandPath[sNa])];
fi;
fi;
pathSegments[totalNumberOfSegments] := subpath (b, e) of inspectedPath;
if (length(pathSegments[totalNumberOfSegments])<=2):
pathSegments[totalNumberOfSegments] := pathSubdivide(pathSegments[totalNumberOfSegments], 2);
fi;
segmentWidth[totalNumberOfSegments] := tSegWidth;
segmentStyle[totalNumberOfSegments] := tSegStyle;
endfor;
else:
totalNumberOfSegments := totalNumberOfSegments + 1;
numberOfSegments[sNa] := 1;
pathSegments[totalNumberOfSegments] := inspectedPath;
segmentWidth[totalNumberOfSegments] := tSegWidth;
segmentStyle[totalNumberOfSegments] := tSegStyle;
fi;
n := 0;
for i := scantokens(knotName.intLayers[sNa]):
n := n + 1;
if (n <= numberOfSegments[sNa]):
appendList[i](layerContents, totalNumberOfSegments - numberOfSegments[sNa] + n, 1, true);
appendList(layersList, i, 1, true);
fi;
endfor;
if n > 0:
if n < numberOfSegments[sNa]:
for i := n + 1 step 1 until numberOfSegments[sNa]:
appendList0(layerContents, totalNumberOfSegments - numberOfSegments[sNa] + i, 1, true);
endfor;
appendList(layersList, 0, 1, true);
fi;
else:
for i := 1 step 1 until numberOfSegments[sNa]:
appendList0(layerContents, totalNumberOfSegments - numberOfSegments[sNa] + i, 1, true);
endfor;
appendList(layersList, 0, 1, true);
fi;
endfor;
image(
for i := scantokens(sortList(layersList, false)):
layerPicture[i] := image(
for j := scantokens(layerContents[i]):
numberOfShadows := -1;
shadowsEnabled := false;
for k := 0 step 1 until totalNumberOfShadows:
if xpart((subpath (1/10, length(pathSegments[j]) - 1/10) of pathSegments[j])
intersectiontimes allShadowPaths[k]) > 0:
shadowsEnabled := true;
numberOfShadows := numberOfShadows + 1;
shadowPath[numberOfShadows] := allShadowPaths[k];
shadowDepth[numberOfShadows] := 3/2segmentWidth[j];
fi;
endfor;
save drawTubeEnds;
boolean drawTubeEnds;
drawTubeEnds := false;
draw tube.scantokens(segmentStyle[j])(pathSegments[j])(segmentWidth[j]) if segmentStyle[j] = "e": withpen thickpen fi;
tmpShadows := tmpShadows + 1;
allShadowPaths[tmpShadows] := tubeOutline;
endfor;
);
for j := 0 step 1 until totalNumberOfShadows:
layerPicture[i] := layerPicture[i] maskedWith allShadowPaths[j];
endfor;
totalNumberOfShadows := tmpShadows;
endfor;
for i := scantokens(sortList(layersList, true)):
draw layerPicture[i];
endfor;
)
enddef;
vardef addStrandToKnot (suffix knotName) (expr p, w, s, intersectionLayers) =
save n;
if not known knotName.nStrands:
numeric knotName.nStrands;
knotName.nStrands := 0;
fi;
if not path knotName.strandPath0:
path knotName.strandPath[];
fi;
if not numeric knotName.strandWidth0:
numeric knotName.strandWidth[];
fi;
if not string knotName.strandStyle0:
string knotName.strandStyle[];
fi;
if not string knotName.intLayers0:
string knotName.intLayers[];
fi;
knotName.nStrands := knotName.nStrands + 1;
n := knotName.nStrands;
knotName.strandPath[n] := p;
knotName.strandWidth[n] := w;
knotName.strandStyle[n] := s;
knotName.intLayers[n] := intersectionLayers;
enddef;