summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/metapost
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2019-03-08 22:22:36 +0000
committerKarl Berry <karl@freefriends.org>2019-03-08 22:22:36 +0000
commit56f9f18f9b7c498d77c5b2408b63ad4702ed50ba (patch)
tree85604697de75d43d73a78a6eabb060829c1694af /Master/texmf-dist/metapost
parent7bb37e17accc83d0fda43c70e5b7eb1368c76f65 (diff)
fiziko (8mar19)
git-svn-id: svn://tug.org/texlive/trunk@50293 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/metapost')
-rw-r--r--Master/texmf-dist/metapost/fiziko/fiziko.mp2163
1 files changed, 2163 insertions, 0 deletions
diff --git a/Master/texmf-dist/metapost/fiziko/fiziko.mp b/Master/texmf-dist/metapost/fiziko/fiziko.mp
new file mode 100644
index 00000000000..cf134b0fed0
--- /dev/null
+++ b/Master/texmf-dist/metapost/fiziko/fiziko.mp
@@ -0,0 +1,2163 @@
+% fiziko 0.1.3
+% MetaPost library for physics textbook illustrations
+% Copyright 2019 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 <http://www.gnu.org/licenses/>.
+
+% 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;
+ 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));
+ bp := ypart((ulcorner(p)--llcorner(p)) firstIntersectionTimes p);
+ invertedmask := (subpath (bp, length(p) + bp) 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;
+
+%
+% 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 := 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;
+
+ % here are some pens
+ pen thinpen, thickpen, fatpen;
+
+ thinpen := pencircle scaled minStrokeWidth;
+ thickpen := pencircle scaled 3minStrokeWidth;
+ fatpen := pencircle scaled 6minStrokeWidth;
+enddef;
+
+defineMinStrokeWidth(1/5pt);
+
+% here we set global light direction
+
+def defineLightDirection (expr ldx, ldy) =
+ pair lightDirection, lightDirectionVector;
+ lightDirection := (ldx, ldy);
+ lightDirectionVector := (sin(xpart(lightDirection)), sin(ypart(lightDirection)));
+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<length(p):
+ sp := subpath(i, i + subdivideStep) of p;
+ returnPath := returnPath .. controls (postcontrol 0 of sp) and (precontrol 1 of sp) .. pathSubdivideBase (p, subdivideStep, i + subdivideStep);
+ fi;
+ if (i = 0) and (cycle p):
+ (subpath(0, length(returnPath)-1) of returnPath) .. controls (postcontrol length(returnPath)-1 of returnPath) and (precontrol length(returnPath) of returnPath) .. cycle
+ else:
+ returnPath
+ fi
+enddef;
+
+vardef offsetPathSubdivide (expr p) =
+ pathSubdivideBase(p, 1/offsetPathSteps, 0)
+enddef;
+
+vardef pathSubdivide (expr p, n) =
+ pathSubdivideBase(p, 1/n, 0)
+enddef;
+
+%
+% This macro creates a template offset path to a straight line, so we can correct angles
+% It might appear as if we need to calculate derivative of the function somehow, instead of mocking it
+% but this function might be anything, function of coordinates of distance to some point etc.,
+% so consider this a lazy way to do the right thing.
+%
+% either offsetPathTime or offsetPathLength are intended to be used as arguments. offsetPathTime is for time and offsetPathLength is for distance
+%
+
+vardef offsetPathTemplate (expr p, i) (text offsetFunction) =
+ save returnPath, offsetPathTime, offsetPathLength, instantDirection, nextDirection;
+ numeric offsetPathTime, offsetPathLength, currentAngle;
+ pair instantDirection, nextDirection;
+ path returnPath;
+ if (i <= length(p)):
+ offsetPathTime := i;
+ else:
+ offsetPathTime := length(p);
+ fi;
+ if (arclength(p) > 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, bp, bt, t;
+ numeric w[], t[];
+ path bp[], 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):
+ bp1 := offsetPathGenerate (p, q yscaled 1/2, 0);
+ bp2 := offsetPathGenerate (p, q yscaled -1/2, 0);
+ fill bp1 -- reverse(bp2) -- 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
+%
+
+def thinBrushGenerate (expr p, q, i) =
+begingroup
+ save w, bp, bt, t, h, linecap;
+ numeric w[], t[];
+ path bp[], bt;
+ bt := q;
+ w0 := (ypart(urcorner(bt)));
+ w1 := (ypart(lrcorner(bt)));
+ w2 := floor((1/2(w0 + w1))/dashStrokeWidthStep)*dashStrokeWidthStep;
+ t := cutPathTime(bt, w2);
+ bp1 := subpath (0, t) of p;
+ bp2 := subpath (t, length(p)) of p;
+ if (((w0 - w1) > dashStrokeWidthStep) and (i < 15))
+ and ((arclength(bp1) > minDashStrokeLength)
+ or (arclength(bp2) > minDashStrokeLength)):
+ thinBrushGenerate (bp1, subpath (0, t) of q, i + 1);
+ thinBrushGenerate (bp2, subpath (t, length(q)) of q, i + 1);
+ else:
+ linecap := butt;
+ if (w2 > minStrokeWidth):
+ w2 := minStrokeWidth;
+ fi;
+ if (w2 >= minDashStrokeWidth) and (arclength(p) > 0):
+ draw p withpen thinpen dashed thinBrushPattern(w2, arclength(p));
+ 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
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle(
+ (angleRad(direction offsetPathTime of p), 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;
+
+%
+% macro that actually draws line of variable width
+%
+
+vardef brush (expr p) (text offsetFunction) =
+ image(
+ brushGenerate (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
+%
+
+def tubeGenerate (expr p, q, d, i) =
+begingroup
+ save w, bw, k, t, tubeWidth, sp, currentPath, currentTubePath, currentDepth;
+ numeric w[], bw[], t, currentDepth;
+ path tubeWidth, sp, currentPath, currentTubePath;
+ tubeWidth := d yscaled 2;
+ w0 := (ypart(urcorner(tubeWidth))) - 1/1000;
+ w1 := (ypart(lrcorner(tubeWidth))) + 1/1000;
+ w2 := ceiling(log(w0/shadingDensity, 2));
+ w3 := ceiling(log(w1/shadingDensity, 2));
+ if ((w2 > w3) and (i<20)):
+ t := cutPathTime(tubeWidth, shadingDensity*(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 := arccos(((k/bw1)*2) - 1);
+ currentDepth := -abs((1-sin(angleOnTube))*w0);
+ if shadowsEnabled:
+ currentPath := shadowCut(currentPath, currentDepth);
+ fi;
+ brushGenerate (currentPath,
+ offsetPathTemplate(currentPath, 0)(
+ maxShadingStrokeWidth
+ if odd (k): * (abs(ypart(point offsetPathTime of tubeWidth)/bw1) - 1/2shadingDensity) fi
+ %* orderFade(xpart(unitvector(direction offsetPathTime of tubeWidth yscaled 1/2cos(angleOnTube))), k) % why was it even here?
+ * angleToLightness(
+ tubeAngleToAbsoulteAngle((
+ angleOnTube,
+ angleRad(direction offsetPathTime of currentTubePath),
+ angleRad(direction offsetPathTime of tubeWidth yscaled 1/2)
+ )), currentDepth, point offsetPathTime of currentPath)
+ ), 0);
+ endfor;
+ fi;
+ fi;
+endgroup
+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)
+ * angleToLightness(
+ tubeAngleToAbsoulteAngle((
+ arccos(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 will eventually appear as a result
+% of some other macro.
+%
+
+vardef tubeAngleToAbsoulteAngle (expr p) =
+ save a;
+ numeric a[];
+ a1 := bluepart(p) + 1/2pi;
+ a2 := arccos(cos(redpart(p))*sin(a1));
+ a3 := greenpart(p) + 1/2pi;
+ a4 := arccos((cos(a1) * cos(a3) - cos(a2) * sin(a3))*(99/100));
+ a5 := arccos((cos(a1) * sin(a3) + cos(a2) * cos(a3))*(99/100));
+ (a5, a4)
+enddef;
+
+%
+% frontends to simplify tube drawing. tubeOutline variable changes on every call
+% of any tube frontend function and can be used afterwards.
+%
+
+path tubeOutline;
+boolean drawTubeEnds;
+drawTubeEnds := true;
+
+vardef tube.l (expr p)(text offsetFunction)=
+ save q;
+ path q[];
+ 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;
+ image(
+ tubeGenerate (q2, q3, q1, 0);
+ if (cycle p) or (not drawTubeEnds):
+ draw q2 withpen thinpen;
+ draw q3 withpen thinpen;
+ else:
+ draw q2--reverse(q3)--cycle withpen thinpen;
+ fi;
+ )
+enddef;
+
+vardef tube.t (expr p)(text offsetFunction)=
+ save q;
+ path q[];
+ 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;
+ image(
+ tubeGenerateAlt (q2, q3, q1);
+ if (cycle p) or (not drawTubeEnds):
+ draw q2 withpen thinpen;
+ draw q3 withpen thinpen;
+ else:
+ draw q2--reverse(q3)--cycle withpen thinpen;
+ fi;
+ )
+enddef;
+
+vardef tube.e (expr p)(text offsetFunction)=
+ save q;
+ path q[];
+ 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 not drawTubeEnds:
+ image(
+ draw q2 withpen thinpen;
+ draw q3 withpen thinpen;
+ )
+ else:
+ tubeOutline := q3--reverse(q2)--cycle;
+ tubeOutline
+ fi
+enddef;
+
+%
+% Sphere can be used as a cap for a tube, so it has same 2**n lines.
+%
+
+vardef sphere.c (expr d) =
+ save currentCircle, origCircle, currentRadius, currentDepth, order, circleThickness;
+ path currentCircle, origCircle;
+ numeric currentRadius, currentDepth, order;
+ origCircle := fullcircle;
+ order := 2**ceiling(log((1/2d)/shadingDensity, 2));
+ image(
+ draw fullcircle scaled d withpen thinpen;
+ for i := 1 upto order:
+ currentRadius := i/order;
+ currentCircle := origCircle scaled (currentRadius*d) rotated uniformdeviate (1/4pi);
+ if odd(i):
+ circleThickness := maxShadingStrokeWidth * ((abs(d - (shadingDensity*order)))/order);
+ else:
+ circleThickness := maxShadingStrokeWidth;
+ fi;
+ currentDepth:= -(1-sqrt(1-currentRadius**2))*(1/2d);
+ if shadowsEnabled:
+ currentCircle := shadowCut(currentCircle, currentDepth);
+ fi;
+ brushGenerate (currentCircle,
+ offsetPathTemplate (currentCircle, 0) (
+ circleThickness
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle(
+ (angleRad(direction offsetPathTime of currentCircle), currentRadius)
+ ), currentDepth, point offsetPathTime of currentCircle
+ )
+ ), 0);
+ endfor;
+ )
+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 sphere.l (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)
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle((
+ (
+ if (abs(point offsetPathTime of currentPath) > 0):
+ angleRad(point offsetPathTime of currentPath)
+ else:
+ 0
+ fi
+ + 1/2pi), 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 absolute angle
+%
+
+vardef sphereAngleToAbsoulteAngle (expr p) =
+ save a;
+ numeric a[];
+ a1 := xpart(p) - 1/2pi;
+ a2 := arcsin(ypart(p));
+ a3 := arccos(sin(a2)*cos(a1));
+ a4 := pi - arccos(sin(a2)*sin(a1));
+ (a3, a4)
+enddef;
+
+%
+% Once we get two angles at some point of some surface, we can compute light intensity there.
+%
+
+vardef angleToLightness (expr p, d, q) =
+ save returnValue, shiftedShadowPath;
+ path shiftedShadowPath;
+ if shadowsEnabled:
+ for i := 0 step 1 until numberOfShadows:
+ shiftedShadowPath := shadowPath[i] shifted (lightDirectionVector scaled (d-shadowDepth[i]));
+ if q isInside shiftedShadowPath:
+ returnValue := 1;
+ fi;
+ endfor;
+ fi;
+ if not known returnValue:
+ returnValue := (cos(xpart(p) + xpart(lightDirection))++cos(ypart(p) - ypart(lightDirection)));
+ returnValue := angleToLightnessPP(returnValue);
+ fi;
+ if returnValue > 1:
+ 1
+ else:
+ returnValue
+ fi
+enddef;
+
+vardef angleToLightnessPP (expr v) =
+ v**3
+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 (lightDirectionVector scaled (currentDepth - shadowDepth[j]));
+ 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 angleToLightnessPP;
+ vardef angleToLightnessPP (expr v) =
+ 1/2(v**3)
+ enddef;
+ p1 := image(draw sphere.l(2s, lat));
+ vardef angleToLightnessPP (expr v) =
+ if (abs(cos(sphlat)) > 7/8 + uniformdeviate (1/20)):
+ 1/4(v**2)
+ else:
+ 1/3(v**2) + 2/3
+ fi
+ enddef;
+ p2 := image(draw sphere.l(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
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle(
+ (angleRad(point offsetPathTime of lm) + 1/4pi, 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;
+ 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):
+ j := round(i*d1/(3/2shadingDensity));
+ if (j mod 4) = 0:
+ stripeskind := dashpattern (on 8shadingDensity off 4shadingDensity);
+ fi;
+ if ((j mod 4) = 1) or ((j mod 4) = 3):
+ stripeskind := dashpattern (off 1shadingDensity on 6shadingDensity off 5shadingDensity);
+ fi;
+ if (j mod 4) = 2:
+ stripeskind := dashpattern (on 0 off 12shadingDensity);
+ fi;
+ fi;
+ draw ((dir(a) scaled 1/2d2) -- (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);
+ thinBrushGenerate (q5,
+ offsetPathTemplate (q5, 0) (
+ minDashStrokeWidth + minStrokeWidth
+ * angleToLightness(
+ (arccos(1 - offsetPathLength*2), 1/2pi)
+ , 0, point offsetPathTime of q5
+ )
+ ), 0);
+ )
+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
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle(
+ (angleRad(direction offsetPathTime of fullcircle), i/4d)
+ ), 0, point offsetPathTime of fullcircle scaled i
+ )
+ ), 0);
+ else:
+ thinBrushGenerate (fullcircle scaled i,
+ offsetPathTemplate (fullcircle scaled i, 0) (
+ 1/4minStrokeWidth + minStrokeWidth
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle(
+ (angleRad(direction offsetPathTime of fullcircle) + pi, 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
+ * angleToLightness(
+ sphereAngleToAbsoulteAngle(
+ (angleRad(direction offsetPathTime of fullcircle) + pi, 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;
+
+%
+% 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;
+ valueExists := 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 omitDuplicates:
+ for i=scantokens(listName@#):
+ if (i = valueToAdd):
+ valueExists := true;
+ fi;
+ endfor;
+ 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 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;
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-26 02:06:32 +0000