asymptote 2.33

git-svn-id: svn://tug.org/texlive/trunk@37365 c570f23f-e606-0410-a88d-b1316a301751

12 files changed, 1619 insertions, 57 deletions

diff --git a/Master/texmf-dist/asymptote/GUI/xasyVersion.py b/Master/texmf-dist/asymptote/GUI/xasyVersion.py
index 0dbd3d3614d..05ca2c3a8c1 100755
--- a/Master/texmf-dist/asymptote/GUI/xasyVersion.py
+++ b/Master/texmf-dist/asymptote/GUI/xasyVersion.py
@@ -1,2 +1,2 @@
 #!/usr/bin/env python
-xasyVersion = "2.31"
+xasyVersion = "2.33"
diff --git a/Master/texmf-dist/asymptote/asy-keywords.el b/Master/texmf-dist/asymptote/asy-keywords.el
index f9ef7269cc6..ea6435749f5 100644
--- a/Master/texmf-dist/asymptote/asy-keywords.el
+++ b/Master/texmf-dist/asymptote/asy-keywords.el
@@ -2,16 +2,16 @@
 ;; This file is automatically generated by asy-list.pl.
 ;; Changes will be overwritten.
 ;;
-(defvar asy-keywords-version "2.31")
+(defvar asy-keywords-version "2.33")
 
 (defvar asy-keyword-name '(
 and controls tension atleast curl if else while for do return break continue struct typedef new access import unravel from include quote static public private restricted this explicit true false null cycle newframe operator ))
 
 (defvar asy-type-name '(
-Braid FitResult Label Legend Solution TreeNode abscissa arc arrowhead binarytree binarytreeNode block bool bool3 bounds bqe circle conic coord coordsys cputime ellipse file filltype frame grid3 guide horner hsv hyperbola indexedTransform int inversion key light line linefit marginT marker mass object pair parabola path path3 pen picture point position projection real revolution scaleT scientific segment side slice solution splitface string surface tensionSpecifier ticklocate ticksgridT tickvalues transform transformation tree triangle trilinear triple vector vertex void ))
+Braid FitResult Label Legend Solution TreeNode abscissa arc arrowhead binarytree binarytreeNode block bool bool3 bounds bqe circle conic coord coordsys cputime ellipse evaluatedpoint file filltype frame grid3 guide horner hsv hyperbola indexedTransform int inversion key light line linefit marginT marker mass object pair parabola patch path path3 pen picture point position positionedvector projection real revolution scaleT scientific segment side slice solution splitface string surface tensionSpecifier ticklocate ticksgridT tickvalues transform transformation tree triangle trilinear triple vector vertex void ))
 
 (defvar asy-function-name '(
-AND Arc ArcArrow ArcArrows Arrow Arrows Automatic AvantGarde BBox BWRainbow BWRainbow2 Bar Bars BeginArcArrow BeginArrow BeginBar BeginDotMargin BeginMargin BeginPenMargin Blank Bookman Bottom BottomTop Bounds Break Broken BrokenLog CLZ CTZ Ceil Circle CircleBarIntervalMarker Cos Courier CrossIntervalMarker DOSendl DOSnewl DefaultFormat DefaultLogFormat Degrees Dir DotMargin DotMargins Dotted Draw Drawline Embed EndArcArrow EndArrow EndBar EndDotMargin EndMargin EndPenMargin Fill FillDraw Floor Format Full Gaussian Gaussrand Gaussrandpair Gradient Grayscale Helvetica Hermite HookHead InOutTicks InTicks Jn Label Landscape Left LeftRight LeftTicks Legend Linear Log LogFormat Margin Margins Mark MidArcArrow MidArrow NOT NewCenturySchoolBook NoBox NoMargin NoModifier NoTicks NoTicks3 NoZero NoZeroFormat None OR OmitFormat OmitTick OmitTickInterval OmitTickIntervals OutTicks Ox Oy Palatino PaletteTicks Pen PenMargin PenMargins Pentype Portrait RadialShade RadialShadeDraw Rainbow Range Relative Right RightTicks Rotate Round SQR Scale ScaleX ScaleY ScaleZ Seascape Shift Sin Slant Spline StickIntervalMarker Straight Symbol Tan TeXify Ticks Ticks3 TildeIntervalMarker TimesRoman Top TrueMargin UnFill UpsideDown Wheel X XEquals XOR XY XYEquals XYZero XYgrid XZEquals XZZero XZero XZgrid Y YEquals YXgrid YZ YZEquals YZZero YZero YZgrid Yn Z ZX ZXgrid ZYgrid ZapfChancery ZapfDingbats _begingroup3 _cputime _draw _eval _image _labelpath _projection _strokepath _texpath aCos aSin aTan abort abs accel acos acosh acot acsc activatequote add addArrow addMargins addSaveFunction addpenarc addpenline addseg adjust alias align all altitude angabscissa angle angpoint animate annotate anticomplementary antipedal apply approximate arc arcarrowsize arccircle arcdir arcfromcenter arcfromfocus arclength arcnodesnumber arcpoint arcsubtended arcsubtendedcenter arctime arctopath array arrow arrow2 arrowbase arrowbasepoints arrowsize ascii asec asin asinh ask assert asy asycode asydir asyfigure asyfilecode asyinclude asywrite atan atan2 atanh atbreakpoint atexit attach attract atupdate autoformat autoscale autoscale3 axes axes3 axialshade axis axiscoverage azimuth babel background bangles bar barmarksize barsize basealign baseline bbox beep begin beginclip begingroup beginpoint between bevel bezier bezierP bezierPP bezierPPP bezulate bibliography bibliographystyle binarytree binarytreeNode binomial bins bisector bisectorpoint bispline blend blockconnector box bqe brace breakpoint breakpoints brick buildRestoreDefaults buildRestoreThunk buildcycle bulletcolor byte calculateScaling canonical canonicalcartesiansystem cartesiansystem case1 case2 case3 cbrt cd ceil center centerToFocus centroid cevian change2 changecoordsys checkSegment checkconditionlength checker checkincreasing checklengths checkposition checktriangle choose circle circlebarframe circlemarkradius circlenodesnumber circumcenter circumcircle clamped clear clip clipdraw close cmyk code colatitude collect collinear color colorless colors colorspace comma compassmark complement complementary concat concurrent cone conic conicnodesnumber conictype conj connect containmentTree contains contour contour3 controlSpecifier convert coordinates coordsys copy copyPairOrTriple cos cosh cot countIntersections cputime crop cropcode cross crossframe crosshatch crossmarksize csc cubicroots curabscissa curlSpecifier curpoint currentarrow currentexitfunction currentmomarrow currentpolarconicroutine curve cut cutafter cutbefore cyclic cylinder deactivatequote debugger deconstruct defaultdir defaultformat defaultpen defined degenerate degrees delete deletepreamble determinant diagonal diamond diffdiv dir dirSpecifier dirtime display distance divisors do_overpaint dot dotframe dotsize downcase draw drawAll drawDoubleLine drawFermion drawGhost drawGluon drawMomArrow drawPRCcylinder drawPRCdisk drawPRCsphere drawPRCtube drawPhoton drawScalar drawVertex drawVertexBox drawVertexBoxO drawVertexBoxX drawVertexO drawVertexOX drawVertexTriangle drawVertexTriangleO drawVertexX drawarrow drawarrow2 drawline drawpixel drawtick duplicate elle ellipse ellipsenodesnumber embed embed3 embedplayer empty enclose end endScript endclip endgroup endgroup3 endl endpoint endpoints eof eol equation equations erase erasestep erf erfc error errorbar errorbars eval excenter excircle exit exitXasyMode exitfunction exp expfactors expi expm1 exradius extend extension extouch fabs factorial fermat fft fhorner figure file filecode fill filldraw filloutside fillrule filltype find finite finiteDifferenceJacobian firstcut firstframe fit fit2 fixedscaling floor flush fmdefaults fmod focusToCenter font fontcommand fontsize foot format frac frequency fromCenter fromFocus fspline functionshade gamma generate_random_backtrace generateticks gergonne getc getint getpair getreal getstring gettriple gluon gouraudshade graph graphic gray grestore grid grid3 gsave halfbox hatch hdiffdiv hermite hex histogram history hline hprojection hsv hyperbola hyperbolanodesnumber hyperlink hypot identity image incenter incentral incircle increasing incrementposition indexedTransform indexedfigure initXasyMode initdefaults input inradius insert inside integrate interactive interior interp interpolate intersect intersection intersectionpoint intersectionpoints intersections intouch inverse inversion invisible is3D isDuplicate isnan isogonal isogonalconjugate isotomic isotomicconjugate isparabola italic item jobname key kurtosis kurtosisexcess label labelaxis labelmargin labelpath labels labeltick labelx labelx3 labely labely3 labelz labelz3 lastcut latex latitude latticeshade layer layout ldexp leastsquares legend legenditem length lexorder lift light limits line linear linecap lineinversion linejoin linemargin lineskip linetype linewidth link list lm_enorm lm_evaluate_default lm_lmdif lm_lmpar lm_minimize lm_print_default lm_print_quiet lm_qrfac lm_qrsolv locale locate locatefile location log log10 log1p logaxiscoverage longitude lookup makeNode makedraw makepen map margin markangle markangleradius markanglespace markarc marker markinterval marknodes markrightangle markthin markuniform mass masscenter massformat math max max3 maxAfterTransform maxbezier maxbound maxcoords maxlength maxratio maxtimes mean medial median midpoint min min3 minAfterTransform minbezier minbound minipage minratio mintimes miterlimit mktemp momArrowPath momarrowsize monotonic multifigure nativeformat natural needshipout newl newpage newslide newton newtree nextframe nextnormal nextpage nib nodabscissa none norm normalvideo notaknot nowarn numberpage nurb object offset onpath opacity opposite orientation origin orthic orthocentercenter outformat outline outname outprefix output overloadedMessage overwrite pack pad pairs palette parabola parabolanodesnumber parallel parallelogram partialsum path path3 pattern pause pdf pedal periodic perp perpendicular perpendicularmark phantom phi1 phi2 phi3 photon piecewisestraight point polar polarconicroutine polargraph polygon postcontrol postscript pow10 ppoint prc prc0 prconly precision precontrol prepend printBytecode print_random_addresses project projection purge pwhermite quadrant quadraticroots quantize quarticroots quotient radialshade radians radicalcenter radicalline radius rand randompath rd readline realmult realquarticroots rectangle rectangular rectify reflect relabscissa relative relativedistance reldir relpoint reltime remainder remark removeDuplicates rename replace report resetdefaultpen restore restoredefaults reverse reversevideo rf rfind rgb rgba rgbint rms rotate rotateO rotation round roundbox roundedpath roundrectangle samecoordsys sameside sample save savedefaults saveline scale scale3 scaleO scaleT scaleless scientific search searchtree sec secondaryX secondaryY seconds section sector seek seekeof segment segmentlimits sequence setpens sgn sgnd sharpangle sharpdegrees shift shiftless shipout shipout3 show side simeq simpson sin sinh size size3 skewness skip slant sleep slope slopefield solve solveBVP sort sourceline sphere split sqrt square srand standardizecoordsys startScript stdev step stickframe stickmarksize stickmarkspace stop straight straightness string stripdirectory stripextension stripfile stripsuffix strokepath subdivide subitem subpath substr sum surface symmedial symmedian system tab tableau tan tangent tangential tangents tanh tell tensionSpecifier tensorshade tex texcolor texify texpath texpreamble texreset texshipout texsize textpath thick thin tick tickMax tickMax3 tickMin tickMin3 ticklabelshift ticklocate tildeframe tildemarksize tile tiling time times title titlepage topbox transform transformation transpose trembleFuzz triangle triangleAbc triangleabc triangulate tricoef tridiagonal trilinear trim truepoint tube uncycle unfill uniform unique unit unitrand unitsize unityroot unstraighten upcase updatefunction uperiodic upscale uptodate usepackage usersetting usetypescript usleep value variance variancebiased vbox vector vectorfield verbatim view vline vperiodic vprojection warn warning windingnumber write xaxis xaxis3 xaxis3At xaxisAt xequals xlimits xpart xscale xscaleO xtick xtick3 xtrans yaxis yaxis3 yaxis3At yaxisAt yequals ylimits ypart yscale yscaleO ytick ytick3 ytrans zaxis3 zaxis3At zero zero3 zlimits zpart ztick ztick3 ztrans ))
+AND Arc ArcArrow ArcArrows Arrow Arrows AtA Automatic AvantGarde B03 B13 B23 B33 BBox BWRainbow BWRainbow2 Bar Bars BeginArcArrow BeginArrow BeginBar BeginDotMargin BeginMargin BeginPenMargin Blank Bookman Bottom BottomTop Bounds Break Broken BrokenLog CLZ CTZ Ceil Circle CircleBarIntervalMarker Cos Courier CrossIntervalMarker DOSendl DOSnewl DefaultFormat DefaultLogFormat Degrees Dir DotMargin DotMargins Dotted Draw Drawline Embed EndArcArrow EndArrow EndBar EndDotMargin EndMargin EndPenMargin Fill FillDraw Floor Format Full Gaussian Gaussrand Gaussrandpair Gradient Grayscale Helvetica Hermite HookHead InOutTicks InTicks Jn Label Landscape Left LeftRight LeftTicks Legend Linear Log LogFormat Margin Margins Mark MidArcArrow MidArrow NOT NewCenturySchoolBook NoBox NoMargin NoModifier NoTicks NoTicks3 NoZero NoZeroFormat None OR OmitFormat OmitTick OmitTickInterval OmitTickIntervals OutTicks Ox Oy Palatino PaletteTicks Pen PenMargin PenMargins Pentype Portrait RadialShade RadialShadeDraw Rainbow Range Relative Right RightTicks Rotate Round SQR Scale ScaleX ScaleY ScaleZ Seascape Shift Sin Slant Spline StickIntervalMarker Straight Symbol Tan TeXify Ticks Ticks3 TildeIntervalMarker TimesRoman Top TrueMargin UnFill UpsideDown Wheel X XEquals XOR XY XYEquals XYZero XYgrid XZEquals XZZero XZero XZgrid Y YEquals YXgrid YZ YZEquals YZZero YZero YZgrid Yn Z ZX ZXgrid ZYgrid ZapfChancery ZapfDingbats _begingroup3 _cputime _draw _eval _image _labelpath _projection _strokepath _texpath aCos aSin aTan abort abs accel acos acosh acot acsc activatequote add addArrow addMargins addSaveFunction addpenarc addpenline addseg adjust alias align all altitude angabscissa angle angledegrees angpoint animate annotate anticomplementary antipedal apply applytranspose approximate arc arcarrowsize arccircle arcdir arcfromcenter arcfromfocus arclength arcnodesnumber arcpoint arcsubtended arcsubtendedcenter arctime arctopath array arrow arrow2 arrowbase arrowbasepoints arrowsize ascii asec asin asinh ask assert asy asycode asydir asyfigure asyfilecode asyinclude asywrite atan atan2 atanh atbreakpoint atexit attach attract atupdate autoformat autoscale autoscale3 axes axes3 axialshade axis axiscoverage azimuth babel background bangles bar barmarksize barsize basealign baseline bbox beep begin beginclip begingroup beginpoint between bevel bezier bezierP bezierPP bezierPPP bezulate bibliography bibliographystyle binarytree binarytreeNode binomial bins bisector bisectorpoint bispline blend blockconnector box bqe brace breakpoint breakpoints brick buildRestoreDefaults buildRestoreThunk buildcycle bulletcolor byte calculateScaling canonical canonicalcartesiansystem cartesiansystem case1 case2 case3 cbrt cd ceil center centerToFocus centroid cevian change2 changecoordsys checkSegment check_fpt_zero checkconditionlength checker checkincreasing checklengths checkposition checkpt checkptincube checktriangle choose circle circlebarframe circlemarkradius circlenodesnumber circumcenter circumcircle clamped clear clip clipdraw close cmyk code colatitude collect collinear color colorless colors colorspace comma compassmark complement complementary concat concurrent cone conic conicnodesnumber conictype conj connect containmentTree contains contour contour3 controlSpecifier convert coordinates coordsys copy copyPairOrTriple cos cosh cot countIntersections cputime crop cropcode cross crossframe crosshatch crossmarksize csc cubicroots curabscissa curlSpecifier curpoint currentarrow currentexitfunction currentmomarrow currentpolarconicroutine curve cut cutafter cutbefore cyclic cylinder deactivatequote debugger deconstruct defaultdir defaultformat defaultpen defined degenerate degrees delete deletepreamble determinant diagonal diamond diffdiv dir dirSpecifier dirtime display distance divisors do_overpaint dot dotframe dotsize downcase draw drawAll drawDoubleLine drawFermion drawGhost drawGluon drawMomArrow drawPRCcylinder drawPRCdisk drawPRCsphere drawPRCtube drawPhoton drawScalar drawVertex drawVertexBox drawVertexBoxO drawVertexBoxX drawVertexO drawVertexOX drawVertexTriangle drawVertexTriangleO drawVertexX drawarrow drawarrow2 drawline drawpixel drawtick duplicate elle ellipse ellipsenodesnumber embed embed3 embedplayer empty enclose end endScript endclip endgroup endgroup3 endl endpoint endpoints eof eol equation equations erase erasestep erf erfc error errorbar errorbars eval excenter excircle exit exitXasyMode exitfunction exp expfactors expi expm1 exradius extend extension extouch fabs factorial fermat fft fhorner figure file filecode fill filldraw filloutside fillrule filltype find findroot finite finiteDifferenceJacobian firstcut firstframe fit fit2 fixedscaling floor flush fmdefaults fmod focusToCenter font fontcommand fontsize foot format frac frequency fromCenter fromFocus fspline functionshade gamma generate_random_backtrace generateticks gergonne getc getint getpair getreal getstring gettriple gluon gouraudshade graph graphic gray grestore grid grid3 gsave halfbox hatch hdiffdiv hermite hex histogram history hline hprojection hsv hyperbola hyperbolanodesnumber hyperlink hypot identity image implicitsurface incenter incentral incircle increasing incrementposition indexedTransform indexedfigure initXasyMode initdefaults input inradius insert inside insphere integrate interactive interior interp interpolate intersect intersection intersectionpoint intersectionpoints intersections intouch inverse inversion invisible is3D isDuplicate isnan isogonal isogonalconjugate isotomic isotomicconjugate isparabola italic item jobname key kurtosis kurtosisexcess label labelaxis labelmargin labelpath labels labeltick labelx labelx3 labely labely3 labelz labelz3 lastcut latex latitude latticeshade layer layout ldexp leastsquares legend legenditem length lexorder lift light limits line linear linecap lineinversion linejoin linemargin lineskip linetype linewidth link list lm_enorm lm_evaluate_default lm_lmdif lm_lmpar lm_minimize lm_print_default lm_print_quiet lm_qrfac lm_qrsolv locale locate locatefile location log log10 log1p logaxiscoverage longitude lookup make3dgrid makeNode makecircle makedraw makepen map margin markangle markangleradius markanglespace markarc marker markinterval marknodes markrightangle markthin markuniform mass masscenter massformat math max max3 maxAfterTransform maxbezier maxbound maxcoords maxlength maxratio maxtimes mean medial median midpoint min min3 minAfterTransform minbezier minbound minipage minratio mintimes miterlimit mktemp momArrowPath momarrowsize monotonic multifigure nGrad nativeformat natural needshipout newl newpage newslide newton newtree nextframe nextnormal nextpage nib nodabscissa none norm normalout normalvideo notaknot nowarn numberpage nurb object offset onpath opacity opposite orient orientation origin orthic orthocentercenter outformat outline outname outprefix output overloadedMessage overwrite pack pad pairs palette parabola parabolanodesnumber parallel parallelogram partialsum patchwithnormals path path3 pathbetween pathinface pattern pause pdf pedal periodic perp perpendicular perpendicularmark phantom phi1 phi2 phi3 photon piecewisestraight point polar polarconicroutine polargraph polygon postcontrol postscript pow10 ppoint prc prc0 prconly precision precontrol prepend printBytecode print_random_addresses project projection projecttospan projecttospan_findcoeffs purge pwhermite quadpatches quadrant quadraticroots quantize quarticroots quotient radialshade radians radicalcenter radicalline radius rand randompath rd readline realmult realquarticroots rectangle rectangular rectify reflect relabscissa relative relativedistance reldir relpoint reltime remainder remark removeDuplicates rename replace report resetdefaultpen restore restoredefaults reverse reversevideo rf rfind rgb rgba rgbint rms rotate rotateO rotation round roundbox roundedpath roundrectangle samecoordsys sameside sample save savedefaults saveline scale scale3 scaleO scaleT scaleless scientific search searchtree sec secondaryX secondaryY seconds section sector seek seekeof segment segmentlimits sequence setpens sgn sgnd sharpangle sharpdegrees shift shiftless shipout shipout3 show simeq simpson sin sinh size size3 skewness skip slant sleep slice slope slopefield solve solveBVP sort sourceline sphere split sqrt square srand standardizecoordsys startScript stdev step stickframe stickmarksize stickmarkspace stop straight straightness string stripdirectory stripextension stripfile stripsuffix strokepath subdivide subitem subpath substr sum surface symmedial symmedian system tab tableau tan tangent tangential tangents tanh tell tensionSpecifier tensorshade tex texcolor texify texpath texpreamble texreset texshipout texsize textpath thick thin tick tickMax tickMax3 tickMin tickMin3 ticklabelshift ticklocate tildeframe tildemarksize tile tiling time times title titlepage topbox transform transformation transpose trembleFuzz triangle triangleAbc triangleabc triangletoquads triangulate tricoef tridiagonal trilinear trim truepoint tube uncycle unfill uniform unique unit unitrand unitsize unityroot unstraighten upcase updatefunction uperiodic upscale uptodate usepackage usersetting usetypescript usleep value variance variancebiased vbox vector vectorfield verbatim view vline vperiodic vprojection warn warning windingnumber write xaxis xaxis3 xaxis3At xaxisAt xequals xlimits xpart xscale xscaleO xtick xtick3 xtrans yaxis yaxis3 yaxis3At yaxisAt yequals ylimits ypart yscale yscaleO ytick ytick3 ytrans zaxis3 zaxis3At zero zero3 zlimits zpart ztick ztick3 ztrans ))
 
 (defvar asy-variable-name '(
-AliceBlue Align Allow AntiqueWhite Apricot Aqua Aquamarine Aspect Azure BeginPoint Beige Bisque Bittersweet Black BlanchedAlmond Blue BlueGreen BlueViolet Both Break BrickRed Brown BurlyWood BurntOrange CCW CW CadetBlue CarnationPink Center Centered Cerulean Chartreuse Chocolate Coeff Coral CornflowerBlue Cornsilk Crimson Crop Cyan Dandelion DarkBlue DarkCyan DarkGoldenrod DarkGray DarkGreen DarkKhaki DarkMagenta DarkOliveGreen DarkOrange DarkOrchid DarkRed DarkSalmon DarkSeaGreen DarkSlateBlue DarkSlateGray DarkTurquoise DarkViolet DeepPink DeepSkyBlue DefaultHead DimGray DodgerBlue Dotted Down Draw E ENE EPS ESE E_Euler E_PC E_RK2 E_RK3BS Emerald EndPoint Euler Fill FillDraw FireBrick FloralWhite ForestGreen Fuchsia Gainsboro GhostWhite Gold Goldenrod Gray Green GreenYellow Honeydew HookHead Horizontal HotPink I IgnoreAspect IndianRed Indigo Ivory JOIN_IN JOIN_OUT JungleGreen Khaki LM_DWARF LM_MACHEP LM_SQRT_DWARF LM_SQRT_GIANT LM_USERTOL Label Lavender LavenderBlush LawnGreen Left LeftJustified LeftSide LemonChiffon LightBlue LightCoral LightCyan LightGoldenrodYellow LightGreen LightGrey LightPink LightSalmon LightSeaGreen LightSkyBlue LightSlateGray LightSteelBlue LightYellow Lime LimeGreen Linear Linen Log Logarithmic Magenta Mahogany Mark MarkFill MarkPath Maroon Max MediumAquamarine MediumBlue MediumOrchid MediumPurple MediumSeaGreen MediumSlateBlue MediumSpringGreen MediumTurquoise MediumVioletRed Melon MidPoint MidnightBlue Min MintCream MistyRose Moccasin Move MoveQuiet Mulberry N NE NNE NNW NW NavajoWhite Navy NavyBlue NoAlign NoCrop NoFill NoSide OldLace Olive OliveDrab OliveGreen Orange OrangeRed Orchid Ox Oy PC PaleGoldenrod PaleGreen PaleTurquoise PaleVioletRed PapayaWhip Peach PeachPuff Periwinkle Peru PineGreen Pink Plum PowderBlue ProcessBlue Purple RK2 RK3 RK3BS RK4 RK5 RK5DP RK5F RawSienna Red RedOrange RedViolet Rhodamine Right RightJustified RightSide RosyBrown RoyalBlue RoyalPurple RubineRed S SE SSE SSW SW SaddleBrown Salmon SandyBrown SeaGreen Seashell Sepia Sienna Silver SimpleHead SkyBlue SlateBlue SlateGray Snow SpringGreen SteelBlue Suppress SuppressQuiet Tan TeXHead Teal TealBlue Thistle Ticksize Tomato Turquoise UnFill Up VERSION Value Vertical Violet VioletRed W WNW WSW Wheat White WhiteSmoke WildStrawberry XYAlign YAlign Yellow YellowGreen YellowOrange addpenarc addpenline align allowstepping angularsystem animationdelay appendsuffix arcarrowangle arcarrowfactor arrow2sizelimit arrowangle arrowbarb arrowdir arrowfactor arrowhookfactor arrowlength arrowsizelimit arrowtexfactor authorpen axis axiscoverage axislabelfactor background backgroundcolor backgroundpen barfactor barmarksizefactor basealign baselinetemplate beveljoin bigvertexpen bigvertexsize black blue bm bottom bp bracedefaultratio braceinnerangle bracemidangle braceouterangle brown bullet byfoci byvertices camerafactor chartreuse circlemarkradiusfactor circlenodesnumberfactor circleprecision circlescale cm codefile codepen codeskip colorPen coloredNodes coloredSegments conditionlength conicnodesfactor count cputimeformat crossmarksizefactor currentcoordsys currentlight currentpatterns currentpen currentpicture currentposition currentprojection curvilinearsystem cuttings cyan darkblue darkbrown darkcyan darkgray darkgreen darkgrey darkmagenta darkolive darkred dashdotted dashed datepen dateskip debuggerlines debugging deepblue deepcyan deepgray deepgreen deepgrey deepmagenta deepred default defaultControl defaultS defaultbackpen defaultcoordsys defaultexcursion defaultfilename defaultformat defaultmassformat defaultpen defaultseparator diagnostics differentlengths dot dotfactor dotframe dotted doublelinepen doublelinespacing down duplicateFuzz ellipsenodesnumberfactor eps epsgeo epsilon evenodd expansionfactor extendcap fermionpen figureborder figuremattpen file3 firstnode firststep foregroundcolor fuchsia fuzz gapfactor ghostpen gluonamplitude gluonpen gluonratio gray green grey hatchepsilon havepagenumber heavyblue heavycyan heavygray heavygreen heavygrey heavymagenta heavyred hline hwratio hyperbolanodesnumberfactor identity4 ignore inXasyMode inch inches includegraphicscommand inf infinity institutionpen intMax intMin invert invisible itempen itemskip itemstep labelmargin landscape lastnode left legendhskip legendlinelength legendmargin legendmarkersize legendmaxrelativewidth legendvskip lightblue lightcyan lightgray lightgreen lightgrey lightmagenta lightolive lightred lightyellow linemargin lm_infmsg lm_shortmsg longdashdotted longdashed magenta magneticRadius mantissaBits markangleradius markangleradiusfactor markanglespace markanglespacefactor maxrefinements mediumblue mediumcyan mediumgray mediumgreen mediumgrey mediummagenta mediumred mediumyellow middle minDistDefault minblockheight minblockwidth mincirclediameter minipagemargin minipagewidth minvertexangle miterjoin mm momarrowfactor momarrowlength momarrowmargin momarrowoffset momarrowpen monoPen morepoints nCircle newbulletcolor ngraph nil nmesh nobasealign nodeMarginDefault nodesystem nomarker nopoint noprimary nullpath nullpen numarray ocgindex oldbulletcolor olive orange origin overpaint page pageheight pagemargin pagenumberalign pagenumberpen pagenumberposition pagewidth paleblue palecyan palegray palegreen palegrey palemagenta palered paleyellow parabolanodesnumberfactor perpfactor phi photonamplitude photonpen photonratio pi pink plain plain_bounds plain_scaling plus preamblenodes pt purple r3 r4a r4b randMax realDigits realEpsilon realMax realMin red relativesystem reverse right roundcap roundjoin royalblue salmon saveFunctions scalarpen sequencereal settings shipped signedtrailingzero solid springgreen sqrtEpsilon squarecap squarepen startposition stdin stdout stepfactor stepfraction steppagenumberpen stepping stickframe stickmarksizefactor stickmarkspacefactor swap textpen ticksize tildeframe tildemarksizefactor tinv titlealign titlepagepen titlepageposition titlepen titleskip top trailingzero treeLevelStep treeMinNodeWidth treeNodeStep trembleAngle trembleFrequency trembleRandom undefined unitcircle unitsquare up urlpen urlskip version vertexpen vertexsize viewportmargin viewportsize vline white wye xformStack yellow ylabelwidth zerotickfuzz zerowinding ))
+AliceBlue Align Allow AntiqueWhite Apricot Aqua Aquamarine Aspect Azure BeginPoint Beige Bisque Bittersweet Black BlanchedAlmond Blue BlueGreen BlueViolet Both Break BrickRed Brown BurlyWood BurntOrange CCW CW CadetBlue CarnationPink Center Centered Cerulean Chartreuse Chocolate Coeff Coral CornflowerBlue Cornsilk Crimson Crop Cyan Dandelion DarkBlue DarkCyan DarkGoldenrod DarkGray DarkGreen DarkKhaki DarkMagenta DarkOliveGreen DarkOrange DarkOrchid DarkRed DarkSalmon DarkSeaGreen DarkSlateBlue DarkSlateGray DarkTurquoise DarkViolet DeepPink DeepSkyBlue DefaultHead DimGray DodgerBlue Dotted Down Draw E ENE EPS ESE E_Euler E_PC E_RK2 E_RK3BS Emerald EndPoint Euler Fill FillDraw FireBrick FloralWhite ForestGreen Fuchsia Gainsboro GhostWhite Gold Goldenrod Gray Green GreenYellow Honeydew HookHead Horizontal HotPink I IgnoreAspect IndianRed Indigo Ivory JOIN_IN JOIN_OUT JungleGreen Khaki LM_DWARF LM_MACHEP LM_SQRT_DWARF LM_SQRT_GIANT LM_USERTOL Label Lavender LavenderBlush LawnGreen Left LeftJustified LeftSide LemonChiffon LightBlue LightCoral LightCyan LightGoldenrodYellow LightGreen LightGrey LightPink LightSalmon LightSeaGreen LightSkyBlue LightSlateGray LightSteelBlue LightYellow Lime LimeGreen Linear Linen Log Logarithmic Magenta Mahogany Mark MarkFill MarkPath Maroon Max MediumAquamarine MediumBlue MediumOrchid MediumPurple MediumSeaGreen MediumSlateBlue MediumSpringGreen MediumTurquoise MediumVioletRed Melon MidPoint MidnightBlue Min MintCream MistyRose Moccasin Move MoveQuiet Mulberry N NE NNE NNW NULL_VERTEX NW NavajoWhite Navy NavyBlue NoAlign NoCrop NoFill NoSide OldLace Olive OliveDrab OliveGreen Orange OrangeRed Orchid Ox Oy PC PaleGoldenrod PaleGreen PaleTurquoise PaleVioletRed PapayaWhip Peach PeachPuff Periwinkle Peru PineGreen Pink Plum PowderBlue ProcessBlue Purple RK2 RK3 RK3BS RK4 RK5 RK5DP RK5F RawSienna Red RedOrange RedViolet Rhodamine Right RightJustified RightSide RosyBrown RoyalBlue RoyalPurple RubineRed S SE SSE SSW SW SaddleBrown Salmon SandyBrown SeaGreen Seashell Sepia Sienna Silver SimpleHead SkyBlue SlateBlue SlateGray Snow SpringGreen SteelBlue Suppress SuppressQuiet Tan TeXHead Teal TealBlue Thistle Ticksize Tomato Turquoise UnFill Up VERSION Value Vertical Violet VioletRed W WNW WSW Wheat White WhiteSmoke WildStrawberry XHIGH XLOW XYAlign YAlign YHIGH YLOW Yellow YellowGreen YellowOrange ZHIGH ZLOW addpenarc addpenline align allowstepping angularsystem animationdelay appendsuffix arcarrowangle arcarrowfactor arrow2sizelimit arrowangle arrowbarb arrowdir arrowfactor arrowhookfactor arrowlength arrowsizelimit arrowtexfactor authorpen axis axiscoverage axislabelfactor background backgroundcolor backgroundpen barfactor barmarksizefactor basealign baselinetemplate bernstein beveljoin bigvertexpen bigvertexsize black blue bm bottom bp bracedefaultratio braceinnerangle bracemidangle braceouterangle brown bullet byfoci byvertices camerafactor chartreuse circlemarkradiusfactor circlenodesnumberfactor circleprecision circlescale cm codefile codepen codeskip colorPen coloredNodes coloredSegments conditionlength conicnodesfactor count cputimeformat crossmarksizefactor currentcoordsys currentlight currentpatterns currentpen currentpicture currentposition currentprojection curvilinearsystem cuttings cyan darkblue darkbrown darkcyan darkgray darkgreen darkgrey darkmagenta darkolive darkred dashdotted dashed datepen dateskip debuggerlines debugging deepblue deepcyan deepgray deepgreen deepgrey deepmagenta deepred default defaultControl defaultS defaultbackpen defaultcoordsys defaultexcursion defaultfilename defaultformat defaultmassformat defaultpen defaultseparator diagnostics differentlengths dot dotfactor dotframe dotted doublelinepen doublelinespacing down duplicateFuzz ellipsenodesnumberfactor eps epsgeo epsilon evenodd expansionfactor extendcap fermionpen figureborder figuremattpen file3 firstnode firststep foregroundcolor fuchsia fuzz gapfactor ghostpen gluonamplitude gluonpen gluonratio gray green grey hatchepsilon havepagenumber heavyblue heavycyan heavygray heavygreen heavygrey heavymagenta heavyred hline hwratio hyperbolanodesnumberfactor identity4 ignore inXasyMode inch inches includegraphicscommand inf infinity institutionpen intMax intMin invert invisible itempen itemskip itemstep labelmargin landscape lastnode left legendhskip legendlinelength legendmargin legendmarkersize legendmaxrelativewidth legendvskip lightblue lightcyan lightgray lightgreen lightgrey lightmagenta lightolive lightred lightyellow linemargin lm_infmsg lm_shortmsg longdashdotted longdashed magenta magneticRadius mantissaBits markangleradius markangleradiusfactor markanglespace markanglespacefactor maxrefinements mediumblue mediumcyan mediumgray mediumgreen mediumgrey mediummagenta mediumred mediumyellow middle minDistDefault minblockheight minblockwidth mincirclediameter minipagemargin minipagewidth minvertexangle miterjoin mm momarrowfactor momarrowlength momarrowmargin momarrowoffset momarrowpen monoPen morepoints nCircle newbulletcolor ngraph nil nmesh nobasealign nodeMarginDefault nodesystem nomarker nopoint noprimary nullpath nullpen numarray ocgindex oldbulletcolor olive orange origin overpaint page pageheight pagemargin pagenumberalign pagenumberpen pagenumberposition pagewidth paleblue palecyan palegray palegreen palegrey palemagenta palered paleyellow parabolanodesnumberfactor perpfactor phi photonamplitude photonpen photonratio pi pink plain plain_bounds plain_scaling plus preamblenodes pt purple r3 r4a r4b randMax realDigits realEpsilon realMax realMin red relativesystem reverse right roundcap roundjoin royalblue salmon saveFunctions scalarpen sequencereal settings shipped signedtrailingzero solid springgreen sqrtEpsilon squarecap squarepen startposition stdin stdout stepfactor stepfraction steppagenumberpen stepping stickframe stickmarksizefactor stickmarkspacefactor swap textpen ticksize tildeframe tildemarksizefactor tinv titlealign titlepagepen titlepageposition titlepen titleskip top trailingzero treeLevelStep treeMinNodeWidth treeNodeStep trembleAngle trembleFrequency trembleRandom undefined unitcircle unitsquare up urlpen urlskip version vertexpen vertexsize viewportmargin viewportsize vline white wye xformStack yellow ylabelwidth zerotickfuzz zerowinding ))
diff --git a/Master/texmf-dist/asymptote/geometry.asy b/Master/texmf-dist/asymptote/geometry.asy
index 6dc10674bd9..580c49e2e69 100644
--- a/Master/texmf-dist/asymptote/geometry.asy
+++ b/Master/texmf-dist/asymptote/geometry.asy
@@ -7182,7 +7182,9 @@ void perpendicular(picture pic = currentpicture, pair z, pair align, path g,
 // If r < 0, return the corresponding exterior arc of radius |r|.
 path arc(explicit pair B, explicit pair A, explicit pair C, real r)
 {
-  return arc(A, r, degrees(B - A), degrees(C - A));
+  real BA = degrees(B - A);
+  real CA = degrees(C - A);
+  return arc(A, abs(r), BA, CA, (r < 0) ^ ((BA-CA) % 360 < 180) ? CW : CCW);
 }
 
 // *.......End of compatibility routines........*
diff --git a/Master/texmf-dist/asymptote/ode.asy b/Master/texmf-dist/asymptote/ode.asy
index 0b88060e081..712196a3018 100644
--- a/Master/texmf-dist/asymptote/ode.asy
+++ b/Master/texmf-dist/asymptote/ode.asy
@@ -258,7 +258,7 @@ void write(solution S)
 
 // Integrate dy/dt+cy=f(t,y) from a to b using initial conditions y,
 // specifying either the step size h or the number of steps n.
-solution integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
+solution integrate(real y, real c=0, real f(real t, real y), real a, real b=a,
                    real h=0, int n=0, bool dynamic=false, real tolmin=0,
                    real tolmax=0, real dtmin=0, real dtmax=realMax,
                    RKTableau tableau, bool verbose=false)
@@ -273,8 +273,8 @@ solution integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
     else h=(b-a)/n;
   }
 
-  real f(real t, real y)=(c == 0 || tableau.exponential) ? g :
-    new real(real t, real y) {return g(t,y)-c*y;};
+  real F(real t, real y)=(c == 0 || tableau.exponential) ? f :
+    new real(real t, real y) {return f(t,y)-c*y;};
 
   tableau.stepDependence(h,c,tableau.a);
       
@@ -283,7 +283,7 @@ solution integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
   if(tableau.a.lowOrderWeights.length == 0) dynamic=false;
   bool fsal=dynamic &&
     (tableau.a.lowOrderWeights.length > tableau.a.highOrderWeights.length);
-  if(fsal) f0=f(t,y);
+  if(fsal) f0=F(t,y);
 
   real dt=h;
   while(t < b) {
@@ -295,9 +295,9 @@ solution integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
       dt=h;
     }
  
-    real[] predictions={fsal ? f0 : f(t,y)};
+    real[] predictions={fsal ? f0 : F(t,y)};
     for(int i=0; i < tableau.a.steps.length; ++i)
-      predictions.push(f(t+h*tableau.a.steps[i],
+      predictions.push(F(t+h*tableau.a.steps[i],
                          tableau.a.factors[i]*y+h*dot(tableau.a.weights[i],
                                                       predictions)));
 
@@ -306,7 +306,7 @@ solution integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
     if(dynamic) {
       real f1;
       if(fsal) {
-        f1=f(t+h,y0+highOrder);
+        f1=F(t+h,y0+highOrder);
         predictions.push(f1);
       }
       real lowOrder=h*dot(tableau.a.lowOrderWeights,predictions);
diff --git a/Master/texmf-dist/asymptote/palette.asy b/Master/texmf-dist/asymptote/palette.asy
index 62ca83f1c9a..96bbf1550ec 100644
--- a/Master/texmf-dist/asymptote/palette.asy
+++ b/Master/texmf-dist/asymptote/palette.asy
@@ -118,11 +118,11 @@ bounds image(picture pic=currentpicture, real f(real, real),
   real ymax=pic.scale.y.T(final.y);
   real[][] data=new real[ny][nx];
   for(int j=0; j < ny; ++j) {
-    real y=pic.scale.y.Tinv(interp(ymin,ymax,(j+0.5)/nx));
+    real y=pic.scale.y.Tinv(interp(ymin,ymax,(j+0.5)/ny));
     scalefcn Tinv=pic.scale.x.Tinv;
     // Take center point of each bin
     data[j]=sequence(new real(int i) {
-        return f(Tinv(interp(xmin,xmax,(i+0.5)/ny)),y);
+        return f(Tinv(interp(xmin,xmax,(i+0.5)/nx)),y);
       },nx);
   }
   return image(pic,data,range,initial,final,palette,transpose=false,
diff --git a/Master/texmf-dist/asymptote/plain_arrows.asy b/Master/texmf-dist/asymptote/plain_arrows.asy
index d4220103ce0..d7d79217dd4 100644
--- a/Master/texmf-dist/asymptote/plain_arrows.asy
+++ b/Master/texmf-dist/asymptote/plain_arrows.asy
@@ -83,17 +83,9 @@ SimpleHead.head=new path(path g, position position=EndPoint, pen p=currentpen,
   path r=subpath(g,position,0);
   pair x=point(r,0);
   real t=arctime(r,size);
-  pair y=point(r,t);
-  path base=arrowbase(r,y,t,size);
   path left=rotate(-angle,x)*r;
   path right=rotate(angle,x)*r;
-  real[] T=arrowbasepoints(base,left,right,1);
-  pair denom=point(right,T[1])-y;
-  real factor=denom != 0 ? length((point(left,T[0])-y)/denom) : 1;
-  path left=rotate(-angle*factor,x)*r;
-  path right=rotate(angle*factor,x)*r;
-  real[] T=arrowbasepoints(base,left,right,1);
-  return subpath(left,T[0],0)--subpath(right,0,T[1]);
+  return subpath(left,t,0)--subpath(right,0,t);
 };
 
 arrowhead HookHead(real dir=arrowdir, real barb=arrowbarb)
diff --git a/Master/texmf-dist/asymptote/slide.asy b/Master/texmf-dist/asymptote/slide.asy
index 0a3e0f65fab..433426c2f37 100644
--- a/Master/texmf-dist/asymptote/slide.asy
+++ b/Master/texmf-dist/asymptote/slide.asy
@@ -432,10 +432,13 @@ void indexedfigure(string prefix, int first, int last,
                    pair align=S, pen p=itempen, pen figuremattpen=figuremattpen,
                    bool step=itemstep)
 {
+  bool Stepping=stepping;
+  stepping=true;
   string[] s;
   for(int i=first; i <= last; ++i)
     s.push(prefix+string(i));
   multifigure(s,options,caption,align,p,figuremattpen,step=step);
+  stepping=Stepping;
 }
 
 string[] codefile;
diff --git a/Master/texmf-dist/asymptote/smoothcontour3.asy b/Master/texmf-dist/asymptote/smoothcontour3.asy
new file mode 100644
index 00000000000..216f3269d92
--- /dev/null
+++ b/Master/texmf-dist/asymptote/smoothcontour3.asy
@@ -0,0 +1,1563 @@
+// Copyright 2015 Charles Staats III
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// 
+//     http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// smoothcontour3
+// An Asymptote module for drawing smooth implicitly defined surfaces
+// author: Charles Staats III
+// charles dot staats dot iii at gmail dot com
+
+import graph_settings;  // for nmesh
+import three;
+
+/***********************************************/
+/******** LINEAR ALGEBRA ROUTINES **************/
+/********     LEAST-SQUARES       **************/
+/***********************************************/
+
+// Apply a matrix to a vector.
+real[] apply(real[][] matrix, real[] v) {
+  real[] ans = new real[matrix.length];
+  for (int r = 0; r < matrix.length; ++r) {
+    ans[r] = 0;
+    for (int c = 0; c < v.length; ++c) {
+      ans[r] += matrix[r][c] * v[c];
+    }
+  }
+  return ans;
+}
+
+// Apply the transpose of a matrix to a vector,
+// without actually forming the transpose.
+real[] applytranspose(real[][] matrix, real[] v) {
+  real[] ans = new real[matrix[0].length];
+  for (int r = 0; r < ans.length; ++r) ans[r] = 0;
+  for (int c = 0; c < matrix.length; ++c) {
+    for (int r = 0; r < ans.length; ++r) {
+      ans[r] += matrix[c][r] * v[c];
+    }
+  }
+  return ans;
+}
+
+// For a matrix A, returns the matrix product
+// (A transposed) * A.
+// The transpose of A is never actually formed.
+real[][] AtA(real[][] matrix) {
+  real[][] toreturn = new real[matrix[0].length][matrix[0].length];
+  for (int i = 0; i < toreturn.length; ++i) {
+    for (int j = 0; j < toreturn.length; ++j) {
+      toreturn [i][j] = 0;
+    }
+  }
+  for (int k = 0; k < matrix.length; ++k) {
+    for (int i = 0; i < toreturn.length; ++i) {
+      for (int j = 0; j < toreturn.length; ++j) {
+	toreturn[i][j] += matrix[k][i] * matrix[k][j];
+      }
+    }
+  }
+  return toreturn;
+}
+
+// Assuming A is a matrix with independent columns, returns
+// the unique vector y minimizing |Ay - b|^2 (the L2 norm).
+// If the columns of A are not linearly independent,
+// throws an error (if warn == true) or returns an empty array
+// (if warn == false).
+real[] leastsquares(real[][] A, real[] b, bool warn = true) {
+  real[] solution = solve(AtA(A), applytranspose(A, b), warn=false);
+  if (solution.length == 0 && warn)
+    abort("Cannot compute least-squares approximation for " +
+	  "a matrix with linearly dependent columns.");
+  return solution;
+}
+
+/***********************************************/
+/******** CREATING BEZIER PATCHES **************/
+/******** WITH SPECIFIED NORMALS  **************/
+/***********************************************/
+
+// The weight given to minimizing the sum of squares of
+// the mixed partials at the corners of the bezier patch.
+// If this weight is zero, the result is undefined in
+// places and can be rather wild even where it is
+// defined.
+// The struct is used to as a namespace.
+struct pathwithnormals_settings {
+  static real wildnessweight = 1e-3;
+}
+private from pathwithnormals_settings unravel wildnessweight;
+
+// The Bernstein basis polynomials of degree 3:
+real B03(real t) { return (1-t)^3; }
+real B13(real t) { return 3*t*(1-t)^2; }
+real B23(real t) { return 3*t^2*(1-t); }
+real B33(real t) { return t^3; }
+
+private typedef real function(real);
+function[] bernstein = new function[] {B03, B13, B23, B33};
+
+// This function attempts to produce a Bezier patch
+// with the specified boundary path and normal directions.
+// For instance, the patch should be normal to
+// u0normals[0] at (0, 0.25),
+// normal to u0normals[1] at (0, 0.5), and
+// normal to u0normals[2] at (0, 0.75).
+// The actual normal (as computed by the patch.normal() function)
+// may be parallel to the specified normal, antiparallel, or
+// even zero.
+//
+// A small amount of deviation is allowed in order to stabilize
+// the algorithm (by keeping the mixed partials at the corners from
+// growing too large).
+//
+// Note that the specified normals are projected to be orthogonal to
+// the specified boundary path. However, the entries in the array
+// remain intact.
+patch patchwithnormals(path3 external, triple[] u0normals, triple[] u1normals,
+		       triple[] v0normals, triple[] v1normals)
+{
+  assert(cyclic(external));
+  assert(length(external) == 4);
+  assert(u0normals.length == 3);
+  assert(u1normals.length == 3);
+  assert(v0normals.length == 3);
+  assert(v1normals.length == 3);
+  
+  triple[][] controlpoints = new triple[4][4];
+  controlpoints[0][0] = point(external,0);
+  controlpoints[1][0] = postcontrol(external,0);
+  controlpoints[2][0] = precontrol(external,1);
+  controlpoints[3][0] = point(external,1);
+  controlpoints[3][1] = postcontrol(external,1);
+  controlpoints[3][2] = precontrol(external,2);
+  controlpoints[3][3] = point(external,2);
+  controlpoints[2][3] = postcontrol(external,2);
+  controlpoints[1][3] = precontrol(external,3);
+  controlpoints[0][3] = point(external,3);
+  controlpoints[0][2] = postcontrol(external,3);
+  controlpoints[0][1] = precontrol(external, 4);
+
+  real[][] matrix = new real[24][12];
+  for (int i = 0; i < matrix.length; ++i)
+    for (int j = 0; j < matrix[i].length; ++j)
+      matrix[i][j] = 0;
+  real[] rightvector = new real[24];
+  for (int i = 0; i < rightvector.length; ++i)
+    rightvector[i] = 0;
+
+  void addtocoeff(int i, int j, int count, triple coeffs) {
+    if (1 <= i && i <= 2 && 1 <= j && j <= 2) {
+      int position = 3 * (2 * (i-1) + (j-1));
+      matrix[count][position] += coeffs.x;
+      matrix[count][position+1] += coeffs.y;
+      matrix[count][position+2] += coeffs.z;
+    } else {
+      rightvector[count] -= dot(controlpoints[i][j], coeffs);
+    }
+  }
+
+  void addtocoeff(int i, int j, int count, real coeff) {
+    if (1 <= i && i <= 2 && 1 <= j && j <= 2) {
+      int position = 3 * (2 * (i-1) + (j-1));
+      matrix[count][position] += coeff;
+      matrix[count+1][position+1] += coeff;
+      matrix[count+2][position+2] += coeff;
+    } else {
+      rightvector[count] -= controlpoints[i][j].x * coeff;
+      rightvector[count+1] -= controlpoints[i][j].y * coeff;
+      rightvector[count+2] -= controlpoints[i][j].z * coeff;
+    }
+  }
+
+  int count = 0;
+
+  void apply_u0(int j, real a, triple n) {
+    real factor = 3 * bernstein[j](a);
+    addtocoeff(0,j,count,-factor*n);
+    addtocoeff(1,j,count,factor*n);
+  }
+  void apply_u0(real a, triple n) {
+    triple tangent = dir(external, 4-a);
+    n -= dot(n,tangent)*tangent;
+    n = unit(n);
+    for (int j = 0; j < 4; ++j) {
+      apply_u0(j,a,n);
+    }
+    ++count;
+  }
+  apply_u0(0.25, u0normals[0]);
+  apply_u0(0.5, u0normals[1]);
+  apply_u0(0.75, u0normals[2]);
+
+  void apply_u1(int j, real a, triple n) {
+    real factor = 3 * bernstein[j](a);
+    addtocoeff(3,j,count,factor*n);
+    addtocoeff(2,j,count,-factor*n);
+  }
+  void apply_u1(real a, triple n) {
+    triple tangent = dir(external, 1+a);
+    n -= dot(n,tangent)*tangent;
+    n = unit(n);
+    for (int j = 0; j < 4; ++j)
+      apply_u1(j,a,n);
+    ++count;
+  }
+  apply_u1(0.25, u1normals[0]);
+  apply_u1(0.5, u1normals[1]);
+  apply_u1(0.75, u1normals[2]);
+
+  void apply_v0(int i, real a, triple n) {
+    real factor = 3 * bernstein[i](a);
+    addtocoeff(i,0,count,-factor*n);
+    addtocoeff(i,1,count,factor*n);
+  }
+  void apply_v0(real a, triple n) {
+    triple tangent = dir(external, a);
+    n -= dot(n,tangent) * tangent;
+    n = unit(n);
+    for (int i = 0; i < 4; ++i)
+      apply_v0(i,a,n);
+    ++count;
+  }
+  apply_v0(0.25, v0normals[0]);
+  apply_v0(0.5, v0normals[1]);
+  apply_v0(0.75, v0normals[2]);
+
+  void apply_v1(int i, real a, triple n) {
+    real factor = 3 * bernstein[i](a);
+    addtocoeff(i,3,count,factor*n);
+    addtocoeff(i,2,count,-factor*n);
+  }
+  void apply_v1(real a, triple n) {
+    triple tangent = dir(external, 3-a);
+    n -= dot(n,tangent)*tangent;
+    n = unit(n);
+    for (int i = 0; i < 4; ++i)
+      apply_v1(i,a,n);
+    ++count;
+  }
+  apply_v1(0.25, v1normals[0]);
+  apply_v1(0.5, v1normals[1]);
+  apply_v1(0.75, v1normals[2]);
+
+  addtocoeff(0,0,count,9*wildnessweight);
+  addtocoeff(1,1,count,9*wildnessweight);
+  addtocoeff(0,1,count,-9*wildnessweight);
+  addtocoeff(1,0,count,-9*wildnessweight);
+  count+=3;
+  addtocoeff(3,3,count,9*wildnessweight);
+  addtocoeff(2,2,count,9*wildnessweight);
+  addtocoeff(3,2,count,-9*wildnessweight);
+  addtocoeff(2,3,count,-9*wildnessweight);
+  count+=3;
+  addtocoeff(0,3,count,9*wildnessweight);
+  addtocoeff(1,2,count,9*wildnessweight);
+  addtocoeff(1,3,count,-9*wildnessweight);
+  addtocoeff(0,2,count,-9*wildnessweight);
+  count += 3;
+  addtocoeff(3,0,count,9*wildnessweight);
+  addtocoeff(2,1,count,9*wildnessweight);
+  addtocoeff(3,1,count,-9*wildnessweight);
+  addtocoeff(2,0,count,-9*wildnessweight);
+  count += 3;
+
+  real[] solution = leastsquares(matrix, rightvector, warn=false);
+  if (solution.length == 0) { // if the matrix was singular
+    write("Warning: unable to solve matrix for specifying normals "
+	  + "on bezier patch. Using standard method.");
+    return patch(external);
+  }
+  
+  for (int i = 1; i <= 2; ++i) {
+    for (int j = 1; j <= 2; ++j) {
+      int position = 3 * (2 * (i-1) + (j-1));
+      controlpoints[i][j] = (solution[position],
+			     solution[position+1],
+			     solution[position+2]);
+    }
+  }
+
+  return patch(controlpoints);
+}
+
+// A wrapper for the previous function when the normal direction
+// is given as a function of direction. The wrapper can also
+// accommodate cyclic boundary paths of between one and four
+// segments, although the results are best by far when there
+// are four segments.
+patch patchwithnormals(path3 external, triple normalat(triple)) {
+  assert(cyclic(external));
+  assert(1 <= length(external) && length(external) <= 4);
+  while (length(external) < 4) external = external -- cycle;
+  triple[] u0normals = new triple[3];
+  triple[] u1normals = new triple[3];
+  triple[] v0normals = new triple[3];
+  triple[] v1normals = new triple[3];
+  for (int i = 1; i <= 3; ++i) {
+    v0normals[i-1] = unit(normalat(point(external, i/4)));
+    u1normals[i-1] = unit(normalat(point(external, 1 + i/4)));
+    v1normals[i-1] = unit(normalat(point(external, 3 - i/4)));
+    u0normals[i-1] = unit(normalat(point(external, 4 - i/4)));
+  }
+  return patchwithnormals(external, u0normals, u1normals, v0normals, v1normals);
+}
+
+/***********************************************/
+/*********** ROOT-FINDER UTILITY ***************/
+/***********************************************/
+
+// Namespace
+struct rootfinder_settings {
+  static real roottolerance = 1e-4;
+}
+
+// Find a root for the specified continuous (but not
+// necessarily differentiable) function. Whatever
+// value t is returned, it is guaranteed that either
+// t is within tolerance of a sign change, or 
+// abs(f(t)) <= 0.1 tolerance.
+// An error is thrown if fa and fb are both positive
+// or both negative.
+//
+// In the current implementation, binary search is interleaved
+// with a modified version of linear interpolation.
+real findroot(real f(real), real a, real b,
+	       real tolerance = rootfinder_settings.roottolerance,
+	       real fa = f(a), real fb = f(b))
+{
+  if (fa == 0) return a;
+  if (fb == 0) return b;
+  real g(real);
+  if (fa < 0) {
+    assert(fb > 0);
+    g = f;
+  } else {
+    assert(fb < 0);
+    fa = -fa;
+    fb = -fb;
+    g = new real(real t) { return -f(t); };
+  }
+
+  real t = a;
+  real ft = fa;
+
+  while (b - a > tolerance && abs(ft) > 0.1*tolerance) {
+    t = a + (b - a) / 2;
+    ft = g(t);
+    if (ft == 0) return t;
+    else if (ft > 0) {
+      b = t;
+      fb = ft;
+    } else if (ft < 0) {
+      a = t;
+      fa = ft;
+    }
+    
+    // linear interpolation
+    t = a - (b - a) / (fb - fa) * fa;
+
+    // If the interpolated value is close to one edge of
+    // the interval, move it farther away from the edge in
+    // an effort to catch the root in the middle.
+    if (t - a < (b-a)/8) t = a + 2*(t-a);
+    else if (b - t < (b-a)/8) t = b - 2*(b-t);
+
+    assert(t >= a && t <= b);
+    
+    ft = g(t);
+    if (ft == 0) return t;
+    else if (ft > 0) {
+      b = t;
+      fb = ft;
+    } else if (ft < 0) {
+      a = t;
+      fa = ft;
+    }
+
+  }
+  return a - (b - a) / (fb - fa) * fa;
+}
+
+/***********************************************/
+/********* DUAL CUBE GRAPH UTILITY *************/
+/***********************************************/
+
+// Suppose a plane intersects a (hollow) cube, and
+// does not intersect any vertices. Then its intersection
+// with cube forms a cycle. The goal of the code below
+// is to reconstruct the order of the cycle
+// given only an unordered list of which edges the plane
+// intersects.
+//
+// Basically, the question is this: If we know the points
+// in which a more-or-less planar surface intersects the
+// edges of cube, how do we connect those points?
+//
+// When I wrote the code, I was thinking in terms of the
+// dual graph of a cube, in which "vertices" are really
+// faces of the cube and "edges" connect those "vertices."
+
+// An enum for the different "vertices" (i.e. faces)
+// available. NULL_VERTEX is primarily intended as a
+// return value to indicate the absence of a desired
+// vertex.
+private int NULL_VERTEX = -1;
+private int XHIGH = 0;
+private int XLOW = 1;
+private int YHIGH = 2;
+private int YLOW = 3;
+private int ZHIGH = 4;
+private int ZLOW = 5;
+
+// An unordered set of nonnegative integers.
+// Since the intent is to use
+// only the six values from the enum above, no effort
+// was made to use scalable algorithms.
+struct intset {
+  private bool[] ints = new bool[0];
+  private int size = 0;
+
+  bool contains(int item) {
+    assert(item >= 0);
+    if (item >= ints.length) return false;
+    return ints[item];
+  }
+
+  // Returns true if the item was added (i.e., was
+  // not already present).
+  bool add(int item) {
+    assert(item >= 0);
+    while (item >= ints.length) ints.push(false);
+    if (ints[item]) return false;
+    ints[item] = true;
+    ++size;
+    return true;
+  }
+
+  int[] elements() {
+    int[] toreturn;
+    for (int i = 0; i < ints.length; ++i) {
+      if (ints[i]) toreturn.push(i);
+    }
+    return toreturn;
+  }
+
+  int size() { return size; }
+}
+
+// A map from integers to sets of integers. Again, no
+// attempt is made to use scalable data structures.
+struct int_to_intset {
+  int[] keys = new int[0];
+  intset[] values = new intset[0];
+
+  void add(int key, int value) {
+    for (int i = 0; i < keys.length; ++i) {
+      if (keys[i] == key) {
+	values[i].add(value);
+	return;
+      }
+    }
+    keys.push(key);
+    intset newset;
+    values.push(newset);
+    newset.add(value);
+  }
+
+  private int indexOf(int key) {
+    for (int i = 0; i < keys.length; ++i) {
+      if (keys[i] == key) return i;
+    }
+    return -1;
+  }
+
+  int[] get(int key) {
+    int i = indexOf(key);
+    if (i < 0) return new int[0];
+    else return values[i].elements();
+  }
+
+  int numvalues(int key) {
+    int i = indexOf(key);
+    if (i < 0) return 0;
+    else return values[i].size();
+  }
+
+  int numkeys() {
+    return keys.length;
+  }
+}
+
+// A struct intended to represent an undirected edge between
+// two "vertices."
+struct edge {
+  int start;
+  int end;
+  void operator init(int a, int b) {
+    start = a;
+    end = b;
+  }
+  bool bordersvertex(int v) { return start == v || end == v; }
+}
+
+string operator cast(edge e) {
+  int a, b;
+  if (e.start <= e.end) {a = e.start; b = e.end;}
+  else {a = e.end; b = e.start; }
+  return (string)a + " <-> " + (string)b;
+}
+
+bool operator == (edge a, edge b) {
+  if (a.start == b.start && a.end == b.end) return true;
+  if (a.start == b.end && a.end == b.start) return true;
+  return false;
+}
+
+string operator cast(edge[] edges) {
+  string toreturn = "{ ";
+  for (int i = 0; i < edges.length; ++i) {
+    toreturn += edges[i];
+    if (i < edges.length-1) toreturn += ", ";
+  }
+  return toreturn + " }";
+}
+
+// Finally, the function that strings together a list of edges
+// into a cycle. It makes assumptions that hold true if the
+// list of edges did in fact come from a plane intersection
+// containing no vertices of the cube. For instance, such a
+// plane can contain at most two noncollinear points of any
+// one face; consequently, no face can border more than two of
+// the selected edges.
+//
+// If the underlying assumptions prove to be false, the function
+// returns null.
+int[] makecircle(edge[] edges) {
+  if (edges.length == 0) return new int[0];
+  int_to_intset graph;
+  for (edge e : edges) {
+    graph.add(e.start, e.end);
+    graph.add(e.end, e.start);
+  }
+  int currentvertex = edges[0].start;
+  int startvertex = currentvertex;
+  int lastvertex = NULL_VERTEX;
+  int[] toreturn = new int[0];
+  do {
+    toreturn.push(currentvertex);
+    int[] adjacentvertices = graph.get(currentvertex);
+    if (adjacentvertices.length != 2) return null;
+    for (int v : adjacentvertices) {
+      if (v != lastvertex) {
+	lastvertex = currentvertex;
+	currentvertex = v;
+	break;
+      }
+    }
+  } while (currentvertex != startvertex);
+  if (toreturn.length != graph.numkeys()) return null;
+  toreturn.cyclic = true;
+  return toreturn;
+}
+
+/***********************************************/
+/********** PATHS BETWEEN POINTS ***************/
+/***********************************************/
+// Construct paths between two points with additional
+// constraints; for instance, the path must be orthogonal
+// to a certain vector at each of the endpoints, must
+// lie within a specified plane or a specified face
+// of a rectangular solid,....
+
+// A vector (typically a normal vector) at a specified position.
+struct positionedvector {
+  triple position;
+  triple direction;
+  void operator init(triple position, triple direction) {
+    this.position = position;
+    this.direction = direction;
+  }
+}
+
+string operator cast(positionedvector vv) {
+  return "position: " + (string)(vv.position) + " vector: " + (string)vv.direction;
+}
+
+// The angle, in degrees, between two vectors.
+real angledegrees(triple a, triple b) {
+  real lengthprod = abs(a) * abs(b);
+  if (lengthprod == 0) return 0;
+  return aCos(dot(a,b) / lengthprod);
+}
+
+// A path (single curved segment) between two points. At each point
+// is specified a vector orthogonal to the path.
+path3 pathbetween(positionedvector v1, positionedvector v2) {
+  triple n1 = unit(v1.direction);
+  triple n2 = unit(v2.direction);
+
+  triple p1 = v1.position;
+  triple p2 = v2.position;
+  triple delta = p2-p1;
+
+  triple dir1 = delta - dot(delta, n1)*n1;
+  triple dir2 = delta - dot(delta, n2)*n2;
+  return p1 {dir1} .. {dir2} p2;
+}
+
+// Assuming v1 and v2 are linearly independent, returns an array {a, b}
+// such that a v1 + b v2 is the orthogonal projection of toproject onto
+// the span of v1 and v2. If v1 and v2 are dependent, returns an empty array
+// (if warn==false) or throws an error (if warn==true).
+real[] projecttospan_findcoeffs(triple toproject, triple v1, triple v2,
+				   bool warn=false) {
+  real[][] matrix = {{v1.x, v2.x},
+		     {v1.y, v2.y},
+		     {v1.z, v2.z}};
+  real[] desiredanswer = {toproject.x, toproject.y, toproject.z};
+  return leastsquares(matrix, desiredanswer, warn=warn);
+}
+
+// Project the triple toproject into the span of a and b, but restrict
+// to the quarter-plane of linear combinations a v1 + b v2 such that
+// a >= mincoeff and b >= mincoeff. If v1 and v2 are linearly dependent,
+// return a random (positive) linear combination.
+triple projecttospan(triple toproject, triple v1, triple v2,
+		       real mincoeff = 0.05) {
+  real[] coeffs = projecttospan_findcoeffs(toproject, v1, v2, warn=false);
+  real a, b;
+  if (coeffs.length == 0) {
+    a = mincoeff + unitrand();
+    b = mincoeff + unitrand();
+  } else {
+    a = max(coeffs[0], mincoeff);
+    b = max(coeffs[1], mincoeff);
+  }
+  return a*v1 + b*v2;
+}
+
+// A path between two specified vertices of a cyclic path. The
+// path tangent at each endpoint is guaranteed to lie within the
+// quarter-plane spanned by positive linear combinations of the
+// tangents of the two outgoing paths at that endpoint.
+path3 pathbetween(path3 edgecycle, int vertex1, int vertex2) {
+  triple point1 = point(edgecycle, vertex1);
+  triple point2 = point(edgecycle, vertex2);
+  
+  triple v1 = -dir(edgecycle, vertex1, sign=-1);
+  triple v2 =  dir(edgecycle, vertex1, sign= 1);
+  triple direction1 = projecttospan(unit(point2-point1), v1, v2);
+
+  v1 = -dir(edgecycle, vertex2, sign=-1);
+  v2 =  dir(edgecycle, vertex2, sign= 1);
+  triple direction2 = projecttospan(unit(point1-point2), v1, v2);
+
+  return point1 {direction1} .. {-direction2} point2;
+}
+
+// This function applies a heuristic to choose two "opposite"
+// vertices (separated by three segments) of edgecycle, which
+// is required to be a cyclic path consisting of 5 or 6 segments.
+// The two chosen vertices are pushed to savevertices.
+//
+// The function returns a path between the two chosen vertices. The
+// path tangent at each endpoint is guaranteed to lie within the
+// quarter-plane spanned by positive linear combinations of the
+// tangents of the two outgoing paths at that endpoint.
+path3 bisector(path3 edgecycle, int[] savevertices) {
+  real mincoeff = 0.05;
+  assert(cyclic(edgecycle));
+  int n = length(edgecycle);
+  assert(n >= 5 && n <= 6);
+  triple[] forwarddirections = sequence(new triple(int i) {
+      return dir(edgecycle, i, sign=1);
+    }, n);
+  forwarddirections.cyclic = true;
+  triple[] backwarddirections = sequence(new triple(int i) {
+      return -dir(edgecycle, i, sign=-1);
+    }, n);
+  backwarddirections.cyclic = true;
+  real[] angles = sequence(new real(int i) {
+      return angledegrees(forwarddirections[i], backwarddirections[i]);
+    }, n);
+  angles.cyclic = true;
+  int lastindex = (n == 5 ? 4 : 2);
+  real maxgoodness = 0;
+  int chosenindex = -1;
+  triple directionout, directionin;
+  for (int i = 0; i <= lastindex; ++i) {
+    int opposite = i + 3;
+    triple vec = unit(point(edgecycle, opposite) - point(edgecycle, i));
+    real[] coeffsbegin = projecttospan_findcoeffs(vec, forwarddirections[i],
+						    backwarddirections[i]);
+    if (coeffsbegin.length == 0) continue;
+    coeffsbegin[0] = max(coeffsbegin[0], mincoeff);
+    coeffsbegin[1] = max(coeffsbegin[1], mincoeff);
+
+    real[] coeffsend = projecttospan_findcoeffs(-vec, forwarddirections[opposite],
+						    backwarddirections[opposite]);
+    if (coeffsend.length == 0) continue;
+    coeffsend[0] = max(coeffsend[0], mincoeff);
+    coeffsend[1] = max(coeffsend[1], mincoeff);
+
+    real goodness = angles[i] * angles[opposite] * coeffsbegin[0] * coeffsend[0]
+      * coeffsbegin[1] * coeffsend[1];
+    if (goodness > maxgoodness) {
+      maxgoodness = goodness;
+      directionout = coeffsbegin[0] * forwarddirections[i] +
+	coeffsbegin[1] * backwarddirections[i];
+      directionin = -(coeffsend[0] * forwarddirections[opposite] +
+		       coeffsend[1] * backwarddirections[opposite]);
+      chosenindex = i;
+    }
+  }
+  if (chosenindex == -1) {
+    savevertices.push(0);
+    savevertices.push(3);
+    return pathbetween(edgecycle, 0, 3);
+  } else {
+    savevertices.push(chosenindex);
+    savevertices.push(chosenindex+3);
+    return point(edgecycle, chosenindex) {directionout} ..
+      {directionin} point(edgecycle, chosenindex + 3);
+  }
+}
+
+// A path between two specified points (with specified normals) that lies
+// within a specified face of a rectangular solid.
+path3 pathinface(positionedvector v1, positionedvector v2,
+		    triple facenorm, triple edge1normout, triple edge2normout)
+{
+  triple dir1 = cross(v1.direction, facenorm);
+  real dotprod = dot(dir1, edge1normout);
+  if (dotprod > 0) dir1 = -dir1;
+  // Believe it or not, this "tiebreaker" is actually relevant at times,
+  // for instance, when graphing the cone x^2 + y^2 = z^2 over the region
+  // -1 <= x,y,z <= 1.
+  else if (dotprod == 0 && dot(dir1, v2.position - v1.position) < 0) dir1 = -dir1;
+
+  triple dir2 = cross(v2.direction, facenorm);
+  dotprod = dot(dir2, edge2normout);
+  if (dotprod < 0) dir2 = -dir2;
+  else if (dotprod == 0 && dot(dir2, v2.position - v1.position) < 0) dir2 = -dir2;
+
+  return v1.position {dir1} .. {dir2} v2.position;
+}
+
+triple normalout(int face) {
+  if (face == XHIGH) return X;
+  else if (face == YHIGH) return Y;
+  else if (face == ZHIGH) return Z;
+  else if (face == XLOW) return -X;
+  else if (face == YLOW) return -Y;
+  else if (face == ZLOW) return -Z;
+  else return O;
+}
+
+// A path between two specified points (with specified normals) that lies
+// within a specified face of a rectangular solid.
+path3 pathinface(positionedvector v1, positionedvector v2,
+		    int face, int edge1face, int edge2face) {
+  return pathinface(v1, v2, normalout(face), normalout(edge1face),
+		      normalout(edge2face));
+}
+
+/***********************************************/
+/******** DRAWING IMPLICIT SURFACES ************/
+/***********************************************/
+
+// Quadrilateralization:
+// Produce a surface (array of *nondegenerate* Bezier patches) with a
+// specified three-segment boundary. The surface should approximate the
+// zero locus of the specified f with its specified gradient.
+//
+// If it is not possible to produce the desired result without leaving the
+// specified rectangular region, returns a length-zero array.
+//
+// Dividing a triangle into smaller quadrilaterals this way is opposite
+// the usual trend in mathematics. However, the pathwithnormals algorithm
+// does a poor job of choosing a good surface when the boundary path does
+// not consist of four positive-length segments.
+patch[] triangletoquads(path3 external, real f(triple), triple grad(triple),
+			  triple a, triple b) {
+  static real epsilon = 1e-3;
+  assert(length(external) == 3);
+  assert(cyclic(external));
+
+  triple c0 = point(external, 0);
+  triple c1 = point(external, 1);
+  triple c2 = point(external, 2);
+
+  triple center = (c0 + c1 + c2) / 3;
+  triple n = unit(cross(c1-c0, c2-c0));
+
+  real g(real t) { return f(center + t*n); }
+
+  real tmin = -realMax, tmax = realMax;
+  void absorb(real t) {
+    if (t < 0) tmin = max(t,tmin);
+    else tmax = min(t,tmax);
+  }
+  if (n.x != 0) {
+    absorb((a.x - center.x) / n.x);
+    absorb((b.x - center.x) / n.x);
+  }
+  if (n.y != 0) {
+    absorb((a.y - center.y) / n.y);
+    absorb((b.y - center.y) / n.y);
+  }
+  if (n.z != 0) {
+    absorb((a.z - center.z) / n.z);
+    absorb((b.z - center.z) / n.z);
+  }
+
+  real fa = g(tmin);
+  real fb = g(tmax);
+  if ((fa > 0 && fb > 0) || (fa < 0 && fb < 0)) {
+    return new patch[0];
+  } else {
+    real t = findroot(g, tmin, tmax, fa=fa, fb=fb);
+    center += t * n;
+  }
+
+  n = unit(grad(center));
+
+  triple m0 = point(external, 0.5);
+  positionedvector m0 = positionedvector(m0, unit(grad(m0)));
+  triple m1 = point(external, 1.5);
+  positionedvector m1 = positionedvector(m1, unit(grad(m1)));
+  triple m2 = point(external, 2.5);
+  positionedvector m2 = positionedvector(m2, unit(grad(m2)));
+  positionedvector c = positionedvector(center, unit(grad(center)));
+
+  path3 pathto_m0 = pathbetween(c, m0);
+  path3 pathto_m1 = pathbetween(c, m1);
+  path3 pathto_m2 = pathbetween(c, m2);
+
+  path3 quad0 = subpath(external, 0, 0.5)
+    & reverse(pathto_m0)
+    & pathto_m2
+    & subpath(external, -0.5, 0)
+    & cycle;
+  path3 quad1 = subpath(external, 1, 1.5)
+    & reverse(pathto_m1)
+    & pathto_m0
+    & subpath(external, 0.5, 1)
+    & cycle;
+  path3 quad2 = subpath(external, 2, 2.5)
+    & reverse(pathto_m2)
+    & pathto_m1
+    & subpath(external, 1.5, 2)
+    & cycle;
+
+  return new patch[] {patchwithnormals(quad0, grad),
+      patchwithnormals(quad1, grad),
+      patchwithnormals(quad2, grad)};   
+}
+
+// Returns true if the point is "nonsingular" (in the sense that the magnitude
+// of the gradient is not too small) AND very close to the zero locus of f
+// (assuming f is locally linear).
+bool check_fpt_zero(triple testpoint, real f(triple), triple grad(triple)) {
+  real testval = f(testpoint);
+  real slope = abs(grad(testpoint));
+  static real tolerance = 2*rootfinder_settings.roottolerance;
+  return !(slope > tolerance && abs(testval) / slope > tolerance);
+}
+
+// Returns true if pt lies within the rectangular solid with
+// opposite corners at a and b.
+bool checkptincube(triple pt, triple a, triple b) {
+  real xmin = a.x;
+  real xmax = b.x;
+  real ymin = a.y;
+  real ymax = b.y;
+  real zmin = a.z;
+  real zmax = b.z;
+  if (xmin > xmax) { real t = xmax; xmax=xmin; xmin=t; }
+  if (ymin > ymax) { real t = ymax; ymax=ymin; ymin=t; }
+  if (zmin > zmax) { real t = zmax; zmax=zmin; zmin=t; }
+
+  return ((xmin <= pt.x) && (pt.x <= xmax) &&
+	  (ymin <= pt.y) && (pt.y <= ymax) &&
+	  (zmin <= pt.z) && (pt.z <= zmax));
+
+}
+
+// A convenience function for combining the previous two tests.
+bool checkpt(triple testpt, real f(triple), triple grad(triple),
+	      triple a, triple b) {
+  return checkptincube(testpt, a, b) &&
+    check_fpt_zero(testpt, f, grad);
+}
+
+// Attempts to fill in the boundary cycle with a collection of
+// patches to approximate smoothly the zero locus of f. If unable to
+// do so while satisfying certain checks, returns null.
+// This is distinct from returning an empty
+// array, which merely indicates that the boundary cycle is too small
+// to be worth filling in.
+patch[] quadpatches(path3 edgecycle, positionedvector[] corners,
+		     real f(triple), triple grad(triple),
+		     triple a, triple b) {
+  assert(corners.cyclic);
+
+  // The tolerance for considering two points "essentially identical."
+  static real tolerance = 2.5 * rootfinder_settings.roottolerance;
+  
+  // If there are two neighboring vertices that are essentially identical,
+  // unify them into one.
+  for (int i = 0; i < corners.length; ++i) {
+    if (abs(corners[i].position - corners[i+1].position) < tolerance) {
+      if (corners.length == 2) return new patch[0];
+      corners.delete(i);
+      edgecycle = subpath(edgecycle, 0, i)
+	& subpath(edgecycle, i+1, length(edgecycle))
+	& cycle;
+      --i;
+      assert(length(edgecycle) == corners.length);
+    }
+  }
+
+  static real areatolerance = tolerance^2;
+
+  assert(corners.length >= 2);
+  if (corners.length == 2) {
+    // If the area is too small, just ignore it; otherwise, subdivide.
+    real area0 = abs(cross(-dir(edgecycle, 0, sign=-1, normalize=false),
+			   dir(edgecycle, 0, sign=1, normalize=false)));
+    real area1 = abs(cross(-dir(edgecycle, 1, sign=-1, normalize=false),
+			   dir(edgecycle, 1, sign=1, normalize=false)));
+    if (area0 < areatolerance && area1 < areatolerance) return new patch[0];
+    else return null;
+  }
+  if (length(edgecycle) > 6) abort("too many edges: not possible.");
+
+  for (int i = 0; i < length(edgecycle); ++i) {
+    if (angledegrees(dir(edgecycle,i,sign=1),
+		     dir(edgecycle,i+1,sign=-1)) > 80) {
+      return null;
+    }
+  }
+  
+  if (length(edgecycle) == 3) {
+    patch[] toreturn = triangletoquads(edgecycle, f, grad, a, b);
+    if (toreturn.length == 0) return null;
+    else return toreturn;
+  }
+  if (length(edgecycle) == 4) {
+    return new patch[] {patchwithnormals(edgecycle, grad)};
+  }
+
+  int[] bisectorindices;
+  path3 middleguide = bisector(edgecycle, bisectorindices);
+  
+  triple testpoint = point(middleguide, 0.5);
+  if (!checkpt(testpoint, f, grad, a, b)) {
+    return null;
+  }
+  
+  patch[] toreturn = null;
+  path3 firstpatch = subpath(edgecycle, bisectorindices[0], bisectorindices[1])
+    & reverse(middleguide) & cycle;
+  if (length(edgecycle) == 5) {
+    path3 secondpatch = middleguide
+	& subpath(edgecycle, bisectorindices[1], 5+bisectorindices[0]) & cycle;
+    toreturn = triangletoquads(secondpatch, f, grad, a, b);
+    if (toreturn.length == 0) return null;
+    toreturn.push(patchwithnormals(firstpatch, grad));
+  } else {
+    // now length(edgecycle) == 6
+    path3 secondpatch = middleguide
+      & subpath(edgecycle, bisectorindices[1], 6+bisectorindices[0])
+      & cycle;
+    toreturn = new patch[] {patchwithnormals(firstpatch, grad),
+	patchwithnormals(secondpatch, grad)};
+  }
+  return toreturn;
+}
+
+// Numerical gradient of a function
+typedef triple vectorfunction(triple);
+vectorfunction nGrad(real f(triple)) {
+  static real epsilon = 1e-3;
+  return new triple(triple v) {
+    return ( (f(v + epsilon*X) - f(v - epsilon*X)) / (2 epsilon),
+			 (f(v + epsilon*Y) - f(v - epsilon*Y)) / (2 epsilon),
+			 (f(v + epsilon*Z) - f(v - epsilon*Z)) / (2 epsilon) );
+  };
+}
+
+// A point together with a value at that location.
+struct evaluatedpoint {
+  triple pt;
+  real value;
+  void operator init(triple pt, real value) {
+    this.pt = pt;
+    this.value = value;
+  }
+}
+
+triple operator cast(evaluatedpoint p) { return p.pt; }
+
+// Compute the values of a function at every vertex of an nx by ny by nz
+// array of rectangular solids.
+evaluatedpoint[][][] make3dgrid(triple a, triple b, int nx, int ny, int nz,
+				   real f(triple), bool allowzero = false)
+{
+  evaluatedpoint[][][] toreturn = new evaluatedpoint[nx+1][ny+1][nz+1];
+  for (int i = 0; i <= nx; ++i) {
+    for (int j = 0; j <= ny; ++j) {
+      for (int k = 0; k <= nz; ++k) {
+	triple pt = (interp(a.x, b.x, i/nx),
+		     interp(a.y, b.y, j/ny),
+		     interp(a.z, b.z, k/nz));
+	real value = f(pt);
+	if (value == 0 && !allowzero) value = 1e-5;
+	toreturn[i][j][k] = evaluatedpoint(pt, value);
+      }
+    }
+  }
+  return toreturn;
+}
+
+// The following utilities make, for instance, slice(A, i, j, k, l)
+// equivalent to what A[i:j][k:l] ought to mean for two- and three-
+// -dimensional arrays of evaluatedpoints and of positionedvectors.
+typedef evaluatedpoint T;
+T[][] slice(T[][] a, int start1, int end1, int start2, int end2) {
+  T[][] toreturn = new T[end1-start1][];
+  for (int i = start1; i < end1; ++i) {
+    toreturn[i-start1] = a[i][start2:end2];
+  }
+  return toreturn;
+}
+T[][][] slice(T[][][] a, int start1, int end1,
+	      int start2, int end2,
+	      int start3, int end3) {
+  T[][][] toreturn = new T[end1-start1][][];
+  for (int i = start1; i < end1; ++i) {
+    toreturn[i-start1] = slice(a[i], start2, end2, start3, end3);
+  }
+  return toreturn;
+}
+typedef positionedvector T;
+T[][] slice(T[][] a, int start1, int end1, int start2, int end2) {
+  T[][] toreturn = new T[end1-start1][];
+  for (int i = start1; i < end1; ++i) {
+    toreturn[i-start1] = a[i][start2:end2];
+  }
+  return toreturn;
+}
+T[][][] slice(T[][][] a, int start1, int end1,
+	      int start2, int end2,
+	      int start3, int end3) {
+  T[][][] toreturn = new T[end1-start1][][];
+  for (int i = start1; i < end1; ++i) {
+    toreturn[i-start1] = slice(a[i], start2, end2, start3, end3);
+  }
+  return toreturn;
+}
+
+// An object of class gridwithzeros stores the values of a function at each vertex
+// of a three-dimensional grid, together with zeros of the function along edges
+// of the grid and the gradient of the function at each such zero.
+struct gridwithzeros {
+  int nx, ny, nz;
+  evaluatedpoint[][][] corners;
+  positionedvector[][][] xdirzeros;
+  positionedvector[][][] ydirzeros;
+  positionedvector[][][] zdirzeros;
+  triple grad(triple);
+  real f(triple);
+  int maxdepth;
+
+  // Populate the edges with zeros that have a sign change and are not already
+  // populated.
+  void fillzeros() {
+    for (int j = 0; j < ny+1; ++j) {
+      for (int k = 0; k < nz+1; ++k) {
+	real y = corners[0][j][k].pt.y;
+	real z = corners[0][j][k].pt.z;
+	real f_along_x(real t) { return f((t, y, z)); }
+	for (int i = 0; i < nx; ++i) {
+	  if (xdirzeros[i][j][k] != null) continue;
+	  evaluatedpoint start = corners[i][j][k];
+	  evaluatedpoint end = corners[i+1][j][k];
+	  if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0))
+	    xdirzeros[i][j][k] = null;
+	  else {
+	    triple root = (0,y,z);
+	    root += X * findroot(f_along_x, start.pt.x, end.pt.x,
+				    fa=start.value, fb=end.value);
+	    triple normal = grad(root);
+	    xdirzeros[i][j][k] = positionedvector(root, normal);
+	  }
+	}
+      }
+    }
+
+    for (int i = 0; i < nx+1; ++i) {
+      for (int k = 0; k < nz+1; ++k) {
+	real x = corners[i][0][k].pt.x;
+	real z = corners[i][0][k].pt.z;
+	real f_along_y(real t) { return f((x, t, z)); }
+	for (int j = 0; j < ny; ++j) {
+	  if (ydirzeros[i][j][k] != null) continue;
+	  evaluatedpoint start = corners[i][j][k];
+	  evaluatedpoint end = corners[i][j+1][k];
+	  if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0))
+	    ydirzeros[i][j][k] = null;
+	  else {
+	    triple root = (x,0,z);
+	    root += Y * findroot(f_along_y, start.pt.y, end.pt.y,
+				    fa=start.value, fb=end.value);
+	    triple normal = grad(root);
+	    ydirzeros[i][j][k] = positionedvector(root, normal);
+	  }
+	}
+      }
+    }
+
+    for (int i = 0; i < nx+1; ++i) {
+      for (int j = 0; j < ny+1; ++j) {
+	real x = corners[i][j][0].pt.x;
+	real y = corners[i][j][0].pt.y;
+	real f_along_z(real t) { return f((x, y, t)); }
+	for (int k = 0; k < nz; ++k) {
+	  if (zdirzeros[i][j][k] != null) continue;
+	  evaluatedpoint start = corners[i][j][k];
+	  evaluatedpoint end = corners[i][j][k+1];
+	  if ((start.value > 0 && end.value > 0) || (start.value < 0 && end.value < 0))
+	    zdirzeros[i][j][k] = null;
+	  else {
+	    triple root = (x,y,0);
+	    root += Z * findroot(f_along_z, start.pt.z, end.pt.z,
+			       fa=start.value, fb=end.value);
+	    triple normal = grad(root);
+	    zdirzeros[i][j][k] = positionedvector(root, normal);
+	  }
+	}
+      }
+    }
+  }
+
+  // Fill in the grid vertices and the zeros along edges. Each cube starts at
+  // depth one and the depth increases each time it subdivides; maxdepth is the
+  // maximum subdivision depth. When a cube at maxdepth cannot be resolved to
+  // patches, it is left empty.
+  void operator init(int nx, int ny, int nz,
+		     real f(triple), triple a, triple b,
+		     int maxdepth = 6) {
+    this.nx = nx;
+    this.ny = ny;
+    this.nz = nz;
+    grad = nGrad(f);
+    this.f = f;
+    this.maxdepth = maxdepth;
+    corners = make3dgrid(a, b, nx, ny, nz, f);
+    xdirzeros = new positionedvector[nx][ny+1][nz+1];
+    ydirzeros = new positionedvector[nx+1][ny][nz+1];
+    zdirzeros = new positionedvector[nx+1][ny+1][nz];
+
+    for (int i = 0; i <= nx; ++i) {
+      for (int j = 0; j <= ny; ++j) {
+	for (int k = 0; k <= nz; ++k) {
+	  if (i < nx) xdirzeros[i][j][k] = null;
+	  if (j < ny) ydirzeros[i][j][k] = null;
+	  if (k < nz) zdirzeros[i][j][k] = null;
+	}
+      }
+    }
+
+    fillzeros();
+  }
+
+  // Doubles nx, ny, and nz by halving the sizes of the cubes along the x, y, and z
+  // directions (resulting in 8 times as many cubes). Already existing data about
+  // function values and zeros is copied; vertices and edges with no such pre-existing
+  // data are populated.
+  //
+  // Returns true if subdivide succeeded, false if it failed (because maxdepth
+  // was exceeded).
+  bool subdivide() {
+    if (maxdepth <= 1) {
+      return false;
+    }
+    --maxdepth;
+    triple a = corners[0][0][0];
+    triple b = corners[nx][ny][nz];
+    nx *= 2;
+    ny *= 2;
+    nz *= 2;
+    evaluatedpoint[][][] oldcorners = corners;
+    corners = new evaluatedpoint[nx+1][ny+1][nz+1];
+    for (int i = 0; i <= nx; ++i) {
+      for (int j = 0; j <= ny; ++j) {
+	for (int k = 0; k <= nz; ++k) {
+	  if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0) {
+	    corners[i][j][k] = oldcorners[quotient(i,2)][quotient(j,2)][quotient(k,2)];
+	  } else {
+	    triple pt = (interp(a.x, b.x, i/nx),
+			 interp(a.y, b.y, j/ny),
+			 interp(a.z, b.z, k/nz));
+	    real value = f(pt);
+	    if (value == 0) value = 1e-5;
+	    corners[i][j][k] = evaluatedpoint(pt, value);
+	  }
+	}
+      }
+    }
+
+    positionedvector[][][] oldxdir = xdirzeros;
+    xdirzeros = new positionedvector[nx][ny+1][nz+1];
+    for (int i = 0; i < nx; ++i) {
+      for (int j = 0; j < ny + 1; ++j) {
+	for (int k = 0; k < nz + 1; ++k) {
+	  if (j % 2 != 0 || k % 2 != 0) {
+	    xdirzeros[i][j][k] = null;
+	  } else {
+	    positionedvector zero = oldxdir[quotient(i,2)][quotient(j,2)][quotient(k,2)];
+	    if (zero == null) {
+	      xdirzeros[i][j][k] = null;
+	      continue;
+	    }
+	    real x = zero.position.x;
+	    if (x > interp(a.x, b.x, i/nx) && x < interp(a.x, b.x, (i+1)/nx)) {
+	      xdirzeros[i][j][k] = zero;
+	    } else {
+	      xdirzeros[i][j][k] = null;
+	    }
+	  }
+	}
+      }
+    }
+
+    positionedvector[][][] oldydir = ydirzeros;
+    ydirzeros = new positionedvector[nx+1][ny][nz+1];
+    for (int i = 0; i < nx+1; ++i) {
+      for (int j = 0; j < ny; ++j) {
+	for (int k = 0; k < nz + 1; ++k) {
+	  if (i % 2 != 0 || k % 2 != 0) {
+	    ydirzeros[i][j][k] = null;
+	  } else {
+	    positionedvector zero = oldydir[quotient(i,2)][quotient(j,2)][quotient(k,2)];
+	    if (zero == null) {
+	      ydirzeros[i][j][k] = null;
+	      continue;
+	    }
+	    real y = zero.position.y;
+	    if (y > interp(a.y, b.y, j/ny) && y < interp(a.y, b.y, (j+1)/ny)) {
+	      ydirzeros[i][j][k] = zero;
+	    } else {
+	      ydirzeros[i][j][k] = null;
+	    }
+	  }
+	}
+      }
+    }
+
+    positionedvector[][][] oldzdir = zdirzeros;
+    zdirzeros = new positionedvector[nx+1][ny+1][nz];
+    for (int i = 0; i < nx + 1; ++i) {
+      for (int j = 0; j < ny + 1; ++j) {
+	for (int k = 0; k < nz; ++k) {
+	  if (i % 2 != 0 || j % 2 != 0) {
+	    zdirzeros[i][j][k] = null;
+	  } else {
+	    positionedvector zero = oldzdir[quotient(i,2)][quotient(j,2)][quotient(k,2)];
+	    if (zero == null) {
+	      zdirzeros[i][j][k] = null;
+	      continue;
+	    }
+	    real z = zero.position.z;
+	    if (z > interp(a.z, b.z, k/nz) && z < interp(a.z, b.z, (k+1)/nz)) {
+	      zdirzeros[i][j][k] = zero;
+	    } else {
+	      zdirzeros[i][j][k] = null;
+	    }
+	  }
+	}
+      }
+    }
+
+    fillzeros();
+    return true;
+  }
+
+  // Forward declaration of the draw method, which will be called by drawcube().
+  patch[] draw(bool[] reportactive = null);
+
+  // Construct the patches, assuming that we are working
+  // with a single cube (nx = ny = nz = 1). This method will subdivide the
+  // cube if necessary. The parameter reportactive should be an array of
+  // length 6. Setting an entry to true indicates that the surface abuts the
+  // corresponding face (according to the earlier enum), and thus that the
+  // algorithm should be sure that something is drawn in the cube sharing
+  // that face--even if all the vertices of that cube have the same sign.
+  patch[] drawcube(bool[] reportactive = null) {
+    // First, determine which edges (if any) actually have zeros on them.
+    edge[] zeroedges = new edge[0];
+    positionedvector[] zeros = new positionedvector[0];
+
+    int currentface, nextface;
+
+    void pushifnonnull(positionedvector v) {
+      if (v != null) {
+	zeroedges.push(edge(currentface, nextface));
+	zeros.push(v);
+      }
+    }
+    positionedvector findzero(int face1, int face2) {
+      edge e = edge(face1, face2);
+      for (int i = 0; i < zeroedges.length; ++i) {
+	if (zeroedges[i] == e) return zeros[i];
+      }
+      return null;
+    }
+    
+    currentface = XLOW;
+    nextface = YHIGH;
+    pushifnonnull(zdirzeros[0][1][0]);
+    nextface = YLOW;
+    pushifnonnull(zdirzeros[0][0][0]);
+    nextface = ZHIGH;
+    pushifnonnull(ydirzeros[0][0][1]);
+    nextface = ZLOW;
+    pushifnonnull(ydirzeros[0][0][0]);
+
+    currentface = XHIGH;
+    nextface = YHIGH;
+    pushifnonnull(zdirzeros[1][1][0]);
+    nextface = YLOW;
+    pushifnonnull(zdirzeros[1][0][0]);
+    nextface = ZHIGH;
+    pushifnonnull(ydirzeros[1][0][1]);
+    nextface = ZLOW;
+    pushifnonnull(ydirzeros[1][0][0]);
+
+    currentface = YHIGH;
+    nextface = ZHIGH;
+    pushifnonnull(xdirzeros[0][1][1]);
+    currentface = ZHIGH;
+    nextface = YLOW;
+    pushifnonnull(xdirzeros[0][0][1]);
+    currentface = YLOW;
+    nextface = ZLOW;
+    pushifnonnull(xdirzeros[0][0][0]);
+    currentface = ZLOW;
+    nextface = YHIGH;
+    pushifnonnull(xdirzeros[0][1][0]);
+
+    //Now, string those edges together to make a circle.
+
+    patch[] subdividecube() {
+      if (!subdivide()) {
+	return new patch[0];
+      }
+      return draw(reportactive);
+    }
+    if (zeroedges.length < 3) {
+      return subdividecube();
+    }
+    int[] faceorder = makecircle(zeroedges);
+    if (alias(faceorder,null)) {
+      return subdividecube();
+    }
+    positionedvector[] patchcorners = new positionedvector[0];
+    for (int i = 0; i < faceorder.length; ++i) {
+      patchcorners.push(findzero(faceorder[i], faceorder[i+1]));
+    }
+    patchcorners.cyclic = true;
+
+    //Now, produce the cyclic path around the edges.
+    path3 edgecycle;
+    for (int i = 0; i < faceorder.length; ++i) {
+      path3 currentpath = pathinface(patchcorners[i], patchcorners[i+1],
+					faceorder[i+1], faceorder[i],
+					faceorder[i+2]);
+      triple testpoint = point(currentpath, 0.5);
+      if (!checkpt(testpoint, f, grad, corners[0][0][0], corners[1][1][1])) {
+	return subdividecube();
+      }
+       
+      edgecycle = edgecycle & currentpath;
+    }
+    edgecycle = edgecycle & cycle;
+
+    patch[] toreturn = quadpatches(edgecycle, patchcorners, f, grad,
+				     corners[0][0][0], corners[1][1][1]);
+    if (alias(toreturn, null)) return subdividecube();
+    return toreturn;
+  }
+
+  // Extracts the specified cube as a gridwithzeros object with
+  // nx = ny = nz = 1.
+  gridwithzeros getcube(int i, int j, int k) {
+    gridwithzeros cube = new gridwithzeros;
+    cube.grad = grad;
+    cube.f = f;
+    cube.nx = 1;
+    cube.ny = 1;
+    cube.nz = 1;
+    cube.maxdepth = maxdepth;
+    cube.corners = slice(corners,i,i+2,j,j+2,k,k+2);
+    cube.xdirzeros = slice(xdirzeros,i,i+1,j,j+2,k,k+2);
+    cube.ydirzeros = slice(ydirzeros,i,i+2,j,j+1,k,k+2);
+    cube.zdirzeros = slice(zdirzeros,i,i+2,j,j+2,k,k+1);
+    return cube;
+  }
+
+  // Returns an array of patches representing the surface.
+  // The parameter reportactive should be an array of
+  // length 6. Setting an entry to true indicates that the surface abuts the
+  // corresponding face of the cube that bounds the entire grid.
+  //
+  // If reportactive == null, it is assumed that this is a top-level call;
+  // a dot is printed to stdout for each cube drawn as a very rough
+  // progress indicator.
+  //
+  // If reportactive != null, then it is assumed that the caller had a strong
+  // reason to believe that this grid contains a part of the surface; the
+  // grid will subdivide all the way to maxdepth if necessary to find points
+  // on the surface.
+  draw = new patch[](bool[] reportactive = null) {
+    // A list of all the patches not already drawn but known
+    // to contain part of the surface. This "queue" is
+    // actually implemented as stack for simplicity, since
+    // it does not make any difference. In a multi-threaded
+    // version of the algorithm, a queue (shared across all threads)
+    // would make more sense than a stack.
+    triple[] queue = new triple[0];
+    bool[][][] enqueued = new bool[nx][ny][nz];
+    for (int i = 0; i < enqueued.length; ++i) {
+      for (int j = 0; j < enqueued[i].length; ++j) {
+	for (int k = 0; k < enqueued[i][j].length; ++k) {
+	  enqueued[i][j][k] = false;
+	}
+      }
+    }
+
+    void enqueue(int i, int j, int k) {
+      if (i >= 0 && i < nx
+	  && j >= 0 && j < ny
+	  && k >= 0 && k < nz
+	  && !enqueued[i][j][k]) {
+	queue.push((i,j,k));
+	enqueued[i][j][k] = true;
+      }
+      if (!alias(reportactive, null)) {
+	if (i < 0) reportactive[XLOW] = true;
+	if (i >= nx) reportactive[XHIGH] = true;
+	if (j < 0) reportactive[YLOW] = true;
+	if (j >= ny) reportactive[YHIGH] = true;
+	if (k < 0) reportactive[ZLOW] = true;
+	if (k >= nz) reportactive[ZHIGH] = true;
+      }
+    }
+    
+    for (int i = 0; i < nx+1; ++i) {
+      for (int j = 0; j < ny+1; ++j) {
+	for (int k = 0; k < nz+1; ++k) {
+	  if (i < nx && xdirzeros[i][j][k] != null) {
+	    for (int jj = j-1; jj <= j; ++jj)
+	      for (int kk = k-1; kk <= k; ++kk)
+		enqueue(i, jj, kk);
+	  }
+	  if (j < ny && ydirzeros[i][j][k] != null) {
+	    for (int ii = i-1; ii <= i; ++ii)
+	      for (int kk = k-1; kk <= k; ++kk)
+		enqueue(ii, j, kk);
+	  }
+	  if (k < nz && zdirzeros[i][j][k] != null) {
+	    for (int ii = i-1; ii <= i; ++ii)
+	      for (int jj = j-1; jj <= j; ++jj)
+		enqueue(ii, jj, k);
+	  }
+	}
+      }
+    }
+
+    if (!alias(reportactive, null) && queue.length == 0) {
+      if (subdivide()) return draw(reportactive);
+    }
+
+    patch[] surface = new patch[0];
+
+    while (queue.length > 0) {
+      triple coord = queue.pop();
+      int i = floor(coord.x);
+      int j = floor(coord.y);
+      int k = floor(coord.z);
+      bool[] reportface = array(6, false);
+      patch[] toappend = getcube(i,j,k).drawcube(reportface);
+      if (reportface[XLOW]) enqueue(i-1,j,k);
+      if (reportface[XHIGH]) enqueue(i+1,j,k);
+      if (reportface[YLOW]) enqueue(i,j-1,k);
+      if (reportface[YHIGH]) enqueue(i,j+1,k);
+      if (reportface[ZLOW]) enqueue(i,j,k-1);
+      if (reportface[ZHIGH]) enqueue(i,j,k+1);
+      surface.append(toappend);
+      if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '.');
+    }
+    if (settings.verbose > 1 && alias(reportactive, null)) write(stdout, '\n');
+    return surface;
+  };
+}
+
+// The external interface of this whole module. Accepts exactly one 
+// function (throws an error if two or zero functions are specified).
+// The function should be differentiable. (Whatever you do, do not
+// pass in an indicator function!) Ideally, the zero locus of the
+// function should be smooth; singularities will significantly slow
+// down the algorithm and potentially give bad results.
+//
+// Returns a plot of the zero locus of the function within the
+// rectangular solid with opposite corners at a and b.
+//
+// Additional parameters:
+// n - the number of initial segments in each of the x, y, z directions.
+// overlapedges - if true, the patches of the surface are slightly enlarged
+//     to compensate for an artifact in which the viewer can see through the
+//     boundary between patches. (Some of this may actually be a result of
+//     edges not lining up perfectly, but I'm fairly sure a lot of it arises
+//     purely as a rendering artifact.)
+// nx - override n in the x direction
+// ny - override n in the y direction
+// nz - override n in the z direction
+// maxdepth - the maximum depth to which the algorithm will subdivide in
+//     an effort to find patches that closely approximate the true surface.
+surface implicitsurface(real f(triple) = null, real ff(real,real,real) = null,
+			 triple a, triple b,
+			 int n = nmesh,
+			 bool keyword overlapedges = false,
+			 int keyword nx=n, int keyword ny=n,
+			 int keyword nz=n,
+			 int keyword maxdepth = 8) {
+  if (f == null && ff == null)
+    abort("implicitsurface called without specifying a function.");
+  if (f != null && ff != null)
+    abort("Only specify one function when calling implicitsurface.");
+  if (f == null) f = new real(triple w) { return ff(w.x, w.y, w.z); };
+  gridwithzeros grid = gridwithzeros(nx, ny, nz, f, a, b, maxdepth=maxdepth);
+  patch[] patches = grid.draw();
+  if (overlapedges) {
+    for (int i = 0; i < patches.length; ++i) {
+      triple center = patches[i].point(1/2,1/2);
+      patches[i] = shift(center) * scale3(1.01) * shift(-center) * patches[i];
+    }
+  }
+  return surface(...patches);
+}
diff --git a/Master/texmf-dist/asymptote/three_surface.asy b/Master/texmf-dist/asymptote/three_surface.asy
index bf070e2a4a9..faf56b546f4 100644
--- a/Master/texmf-dist/asymptote/three_surface.asy
+++ b/Master/texmf-dist/asymptote/three_surface.asy
@@ -1115,6 +1115,21 @@ patch subpatch(patch s, pair a, pair b)
   return patch(subpatch(s.P,a,b),s.straight,s.planar);
 }
 
+// return an array containing an intersection times of path p and surface s.
+real[] intersect(path3 p, patch s, real fuzz=-1)
+{
+  return intersect(p,s.P,fuzz);
+}
+
+// return an array containing an intersection times of path p and surface s.
+real[] intersect(path3 p, surface s, real fuzz=-1)
+{
+  for(int i=0; i < s.s.length; ++i) {
+    real[] T=intersect(p,s.s[i].P,fuzz);
+    if(T.length > 0) return T;
+  }
+  return new real[];
+}
 
 // return an array containing all intersection times of path p and patch s.
 real[][] intersections(path3 p, patch s, real fuzz=-1)
@@ -1168,8 +1183,8 @@ bool overlap(triple[][] p, triple[][] q, real fuzz=-1)
   triple qmin=minbound(q);
   triple qmax=maxbound(q);
 
-  static real Fuzz=1000*realEpsilon;
-  real fuzz=max(10*fuzz,Fuzz*max(abs(pmin),abs(pmax)));
+  if(fuzz == -1)
+    fuzz=1000*realEpsilon*max(abs(pmin),abs(pmax),abs(qmin),abs(qmax));
   
   return
     pmax.x+fuzz >= qmin.x &&
@@ -1363,7 +1378,7 @@ nullpens.cyclic=true;
 
 void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
           material[] surfacepen, pen[] meshpen=nullpens,
-          light light=currentlight, light meshlight=light, string name="",
+          light light=currentlight, light meshlight=nolight, string name="",
           render render=defaultrender, projection P=currentprojection)
 {
   bool is3D=is3D();
@@ -1450,7 +1465,7 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
 
 void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
           material surfacepen=currentpen, pen meshpen=nullpen,
-          light light=currentlight, light meshlight=light, string name="",
+          light light=currentlight, light meshlight=nolight, string name="",
           render render=defaultrender, projection P=currentprojection)
 {
   material[] surfacepen={surfacepen};
@@ -1462,7 +1477,7 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
 
 void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
           material[] surfacepen, pen[] meshpen=nullpens,
-          light light=currentlight, light meshlight=light, string name="",
+          light light=currentlight, light meshlight=nolight, string name="",
           render render=defaultrender)
 {
   if(s.empty()) return;
@@ -1507,7 +1522,7 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
 
 void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
           material surfacepen=currentpen, pen meshpen=nullpen,
-          light light=currentlight, light meshlight=light, string name="",
+          light light=currentlight, light meshlight=nolight, string name="",
           render render=defaultrender)
 {
   material[] surfacepen={surfacepen};
@@ -1519,7 +1534,7 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
 
 void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
           material[] surfacepen, pen meshpen,
-          light light=currentlight, light meshlight=light, string name="",
+          light light=currentlight, light meshlight=nolight, string name="",
           render render=defaultrender)
 {
   pen[] meshpen={meshpen};
diff --git a/Master/texmf-dist/asymptote/three_tube.asy b/Master/texmf-dist/asymptote/three_tube.asy
index 8b621e92dda..bbbd5a65633 100644
--- a/Master/texmf-dist/asymptote/three_tube.asy
+++ b/Master/texmf-dist/asymptote/three_tube.asy
@@ -83,8 +83,10 @@ private real[][][] bispline0(real[][] z, real[][] p, real[][] q, real[][] r,
     count=0;
     for(int i=0; i < n; ++i) {
       bool[] condi=cond[i];
+      bool[] condp=cond[i+1];
       for(int j=0; j < m; ++j)
-        if(condi[j]) ++count;
+        if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1])) 
+          ++count;
     }
   }
 
@@ -104,8 +106,9 @@ private real[][][] bispline0(real[][] z, real[][] p, real[][] q, real[][] r,
     real[] qi=q[i];
     real[] qp=q[ip];
     bool[] condi=all ? null : cond[i];
+    bool[] condp=all ? null : cond[i+1];
     for(int j=0; j < m; ++j) {
-      if(all || condi[j]) {
+      if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1])) {
         real yj=y[j];
         int jp=j+1;
         real yp=y[jp];
@@ -214,7 +217,8 @@ surface surface(triple f(pair z), real[] u, real[] v,
     bool[] activei=all ? null : active[i];
     for(int j=0; j <= nv; ++j) {
       pair z=(ui,v[j]);
-      triple f=(all || (activei[j]=cond(z))) ? f(z) : O;
+      if(!all) activei[j]=cond(z);
+      triple f=f(z);
       fxi[j]=f.x;
       fyi[j]=f.y;
       fzi[j]=f.z;
diff --git a/Master/texmf-dist/asymptote/trembling.asy b/Master/texmf-dist/asymptote/trembling.asy
index 6f6f5973332..f9f0d6a7474 100644
--- a/Master/texmf-dist/asymptote/trembling.asy
+++ b/Master/texmf-dist/asymptote/trembling.asy
@@ -29,10 +29,7 @@
 real magneticRadius=1; // unit is bp in postscript coordinates.
 real trembleFuzz(){return min(1e-3,magneticRadius/10);}
 
-/* <asyxml> <variable type="real" signature="trembleAngle"> <code> </asyxml> */
 real trembleAngle=4, trembleFrequency=0.5, trembleRandom=2;
-/* <asyxml> </code> <documentation> Variables used by the routine 'tremble'.
-   </documentation> </variable> </asyxml> */ 
 
 struct tremble
 {
@@ -58,13 +55,8 @@ struct tremble
     return op;
   }
 
-  /* <asyxml> <function type="pair" signature="attract(pair,path,real)"> <code> 
-     </asyxml> */
   real atime(pair m, path g, real fuzz=trembleFuzz())
-  {/* <asyxml> </code> <documentation> Return the time of the point on path g
-      nearest to 'm'.
-      'fuzz' is the argument 'fuzz' of 'intersect'. </documentation> </function> 
-      </asyxml> */
+  {// Return the time of the point on path g nearest to m, within fuzz.
     if(length(g) == 0) return 0.0;
     real[] t=intersect(m,g,fuzz);
     if(t.length > 0) return t[1];
@@ -91,11 +83,8 @@ struct tremble
     return ot;
   }
 
-  /* <asyxml> <function type="path" signature="addnode(path,real)"> <code>
-     </asyxml> */
   path addnode(path g, real t)
-  {/* <asyxml> </code> <documentation> Add a node to 'g' at point(g,t).
-      </documentation> </function> </asyxml> */
+  {// Add a node to 'g' at point(g,t).
     real l=length(g);
     real rt=t % 1;
     if(l == 0 || (t > l && !cyclic(g)) || rt == 0) return g;
@@ -129,11 +118,8 @@ struct tremble
     return og;
   }
 
-  /* <asyxml> <function type="path" signature="addnodes(path,int)"> <code>
-     </asyxml> */
   path addnodes(path g, int n)
-  {/* <asyxml> </code> <documentation> Add 'n' nodes between each node of 'g'.
-      </documentation> </function> </asyxml> */
+  {// Add 'n' nodes between each node of 'g'.
     real l=length(g);
     if(n == 0 || l == 0) return g;
     path og=g;
@@ -148,8 +134,6 @@ struct tremble
     return og;
   }
 
-  /* <asyxml> <function type="path"
-     signature="tremble(path,real,real,real,real,pair[] )"> <code> </asyxml> */
   void operator init(real angle=trembleAngle, real frequency=trembleFrequency,
                      real random=trembleRandom, real fuzz=trembleFuzz()) {
     this.angle=angle;
@@ -159,7 +143,7 @@ struct tremble
   }
   
   path deform(path g...pair[] magneticPoints) {
-    /* <asyxml> </code> <documentation> Return g as it was handwriting.
+    /* Return g as it was handwriting.
        The postcontrols and precontrols of the nodes of g will be rotated
        by an angle proportional to 'angle'(in degrees).
        If frequency < 1, floor(1/frequency) nodes will be added to g to
@@ -169,8 +153,7 @@ struct tremble
        'random' controls the randomized coefficient which will be multiplied
        by 'angle'.
        random is 0 means don't use randomized coefficient;
-       The higher 'random' is, the more the trembling is randomized.
-       </documentation> </function> </asyxml> */ 
+       The higher 'random' is, the more the trembling is randomized. */
     if(length(g) == 0) return g;
     g=addnodes(g,fuzz*abs(max(g)-min(g))...magneticPoints);
     path tg=g;
diff --git a/Master/texmf-dist/asymptote/version.asy b/Master/texmf-dist/asymptote/version.asy
index 455a0f4aa07..04d0fe337f2 100644
--- a/Master/texmf-dist/asymptote/version.asy
+++ b/Master/texmf-dist/asymptote/version.asy
@@ -1 +1 @@
-string VERSION="2.31";
+string VERSION="2.33";