pgfplots (8may21)

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

11 files changed, 1116 insertions, 59 deletions

diff --git a/Master/texmf-dist/tex/generic/pgfplots/libs/tikzlibrarypgfplots.contourlua.code.tex b/Master/texmf-dist/tex/generic/pgfplots/libs/tikzlibrarypgfplots.contourlua.code.tex
new file mode 100644
index 00000000000..f7821013eb5
--- /dev/null
+++ b/Master/texmf-dist/tex/generic/pgfplots/libs/tikzlibrarypgfplots.contourlua.code.tex
@@ -0,0 +1,91 @@
+%--------------------------------------------
+%
+% Package pgfplots
+%
+% Provides a user-friendly interface to create function plots (normal
+% plots, semi-logplots and double-logplots).
+% 
+% It is based on Till Tantau's PGF package.
+%
+% Copyright 2007-2021 by Christian Feuersänger.
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+% 
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+% 
+% You should have received a copy of the GNU General Public License
+% along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+\ifpgfplots@LUA@supported
+	\pgfplotsutil@directlua{require("prepcontour");}
+	\pgfplotsset{
+		contour lua/.style={
+			contour external={%
+				% we decided to NOT precompute contour/levels:
+				% @before cmd=\pgfplotscontourpopulateallkeys,
+				scanline marks=required,
+				script type=lua,
+				%
+				% this library consumes FOUR coordinates...
+				script input format=x y z meta,
+				% ... and it also generates FOUR. Always read the last one
+				% as it is never subject to rounding issues within the
+				% same contour line:
+				output point meta=explicit,
+				%
+				script={
+					copylines = 0;
+					nblocks   = \if 0\ordering \pgfkeysvalueof{/pgfplots/mesh/rows}\else \pgfkeysvalueof{/pgfplots/mesh/cols}\fi;
+					nlines    = \if 0\ordering \pgfkeysvalueof{/pgfplots/mesh/cols}\else \pgfkeysvalueof{/pgfplots/mesh/rows}\fi;
+					if copylines == nil or nblocks == nil or nlines == nil then
+						io.stderr:write("failed to parse arguments");
+						os.exit(2);
+					end;
+					yvaries = \if 0\ordering 'x'\else 'y'\fi;
+					yvaries = (yvaries == 'y');
+					infile = io.open(\"\infile\", \"r\");
+					outfile = io.open(\"\outfile\", \"w\");
+					mesh = PrepcMesh.new(yvaries, nblocks, nlines, copylines, infile, outfile);
+					\ifpgfplotsplothandlercontour@corners
+					corners = true;
+					\else
+					corners = false;
+					\fi
+					\ifx\thecontourlevels\pgfutil@empty
+					mesh:autocontour(\thecontournumber, \pgfplotsmetamin, \pgfplotsmetamax, corners);
+					\else
+					levels = {\thecontourlevels};
+					for n = 1, \#levels do
+						level = tonumber(levels[n]);
+						if level then
+							mesh:contour(level, corners);
+						end;
+					end;
+					\fi
+					io.close(infile);
+					mesh:printcontours();
+					io.close(outfile);
+				},%
+				#1,%
+			},
+		}
+	}
+\else
+	% do not throw the exception at load time ... this library is
+	% included as default
+	\pgfplotsset{
+		contour lua/.code={
+			\pgfplots@error{Library 'contour lua' requires you to run in LUA mode (i.e. lualatex instead of pdflatex)}%
+		}
+	}
+\fi
+
+
+
diff --git a/Master/texmf-dist/tex/generic/pgfplots/lua/prepcontour.lua b/Master/texmf-dist/tex/generic/pgfplots/lua/prepcontour.lua
new file mode 100755
index 00000000000..cfed282c248
--- /dev/null
+++ b/Master/texmf-dist/tex/generic/pgfplots/lua/prepcontour.lua
@@ -0,0 +1,708 @@
+#!/usr/bin/luatex
+
+--[[
+
+prepcontour [Lua variant] - prepare contour lines (for pgfplots)
+
+Version:  1.4 (2021-02-22)
+
+Copyright (C) 2020-2021  Francesco Poli <invernomuto@paranoici.org>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+--]]
+
+-- class creation function
+function FpolyCreateClass()
+    local newclass = {}   -- this table is the created class (which will act
+                          -- as metatable for the instantiated objects)
+
+    newclass.__index = newclass   -- set the __index metamethod (so that
+                                  -- members missing from objects are searched
+                                  -- for in the class)
+
+    newclass.new = function (...)   -- member function that will be used to
+                                    -- instantiate a new object of this class
+        local object = setmetatable({}, newclass)   -- a table with this class
+                                                    -- set as metatable
+
+        if object.init then
+            object:init(...)   -- run the init method, if present
+        end
+
+        return object
+    end
+
+    return newclass
+end
+
+-- class for the coordinates (and meta value) of a point in 3D space
+Coord = FpolyCreateClass()
+
+function Coord:init()
+    self.x = { nil, nil, nil }  -- point x,y,z coordinates
+    self.meta = nil             -- point meta value
+end
+
+-- class for a 2D mesh of points in 3D space
+PrepcMesh = FpolyCreateClass()
+
+function PrepcMesh:init(yvaries, nblocks, nlines, copylines,
+                        input, output)
+    if yvaries then
+        self.ni = nblocks - 1  -- number of cells in the i direction
+        self.nj = nlines  - 1  -- number of cells in the j direction
+    else
+        self.ni = nlines  - 1  -- number of cells in the i direction
+        self.nj = nblocks - 1  -- number of cells in the j direction
+    end
+    self.st = self.nj + 1  -- stride (number of nodes in the j direction)
+    self.si = self.nj + 1  -- stride for done_i
+    self.sj = self.nj      -- stride for done_j
+    local cl = tonumber(copylines) or 0
+    self.is =  input or io.stdin
+    self.os = output or io.stdout
+
+    self.debug = false
+
+    local st = self.st    -- shorter name
+
+    self.coords = {}  -- mesh of nodes
+                      -- use single-index table for efficiency
+                      -- (i,j) will be mapped to [i*st+j]
+    --[[
+        The mesh of nodes represented by this table is:
+
+            nj    o  o  o  o  o  o  o  o  o  o  o  o
+
+              ^   o  o  o  o  o  o  o  o  o  o  o  o
+              |
+             2|   o  o  o  o  o  o  o  o  o  o  o  o
+              |
+             1|   o  o  o  o  o  o  o  o  o  o  o  o
+              |
+           j=0+   o  o  o  o  o  o  o  o  o  o  o  o
+
+                  +----------------------------->
+                i=0  1  2                          ni
+    --]]
+
+    -- copy the first "cl" text lines from self.is to self.os
+    for l = 1, cl do
+        self.os:write(self.is:read('L'))
+    end
+
+    -- read all the nodes
+    for bl = 0, nblocks - 1 do
+        if bl > 0 then
+            self.is:read('L')  -- this text line should be empty
+        end
+        for ln = 0, nlines - 1 do
+            if yvaries then
+                i, j = bl, ln
+            else
+                i, j = ln, bl
+            end
+            n = Coord.new()
+            n.x[1], n.x[2], n.x[3], n.meta = self.is:read('n', 'n', 'n', 'n')
+            self.coords[i*st+j] = n
+        end
+    end
+end
+
+
+--[[ begin core of the program logic ]]--
+
+-- build N contour lines between meta_min and meta_max
+function PrepcMesh:autocontour(N, meta_min, meta_max, corners, tolerance)
+    -- subdivide the meta_min÷meta_max interval into N equal sub-intervals
+    -- and pick the midpoints of those sub-intervals
+    meta_min = tonumber(meta_min)
+    meta_max = tonumber(meta_max)
+    local step     = (meta_max - meta_min)/N
+    local meta_mid = (meta_max + meta_min)/2
+    local n_mid    = (N + 1)/2
+    for n = 1, N do
+        local isoval = meta_mid + (n - n_mid)*step
+        self:contour(isoval, corners, tolerance)
+    end
+end
+
+-- build contour lines for meta==isoval
+function PrepcMesh:contour(isoval, corners, tolerance)
+    -- set relative tolerance
+    local tol = tolerance or 1e-3
+
+    -- short names
+    local st = self.st
+    local si = self.si
+    local sj = self.sj
+
+    -- local variables
+    local i, j
+
+    if self.done_i == nil then
+        self.done_i = {}  -- table of markers for i-sides
+                          -- again single-index table
+                          -- (i,j) will be mapped to [i*si+j]
+        --[[
+            The corresponding i-sides are:
+
+                nj    o--o--o--o--o--o--o--o--o--o--o--
+
+                  ^   o--o--o--o--o--o--o--o--o--o--o--
+                  |
+                 2|   o--o--o--o--o--o--o--o--o--o--o--
+                  |
+                 1|   o--o--o--o--o--o--o--o--o--o--o--
+                  |
+               j=0+   o--o--o--o--o--o--o--o--o--o--o--
+
+                      +-------------------------->
+                    i=0  1  2                      ni-1
+        --]]
+    end
+
+    if self.done_j == nil then
+        self.done_j = {}  -- table of markers for j-sides
+                          -- again single-index table
+                          -- (i,j) will be mapped to [i*sj+j]
+        --[[
+            The corresponding j-sides are:
+
+                      |  |  |  |  |  |  |  |  |  |  |  |
+              nj-1    o  o  o  o  o  o  o  o  o  o  o  o
+                      |  |  |  |  |  |  |  |  |  |  |  |
+                  ^   o  o  o  o  o  o  o  o  o  o  o  o
+                  |   |  |  |  |  |  |  |  |  |  |  |  |
+                 1|   o  o  o  o  o  o  o  o  o  o  o  o
+                  |   |  |  |  |  |  |  |  |  |  |  |  |
+               j=0+   o  o  o  o  o  o  o  o  o  o  o  o
+
+                      +----------------------------->
+                    i=0  1  2                          ni
+        --]]
+    end
+
+    if self.newcoords == nil then
+        self.newcoords = {}  -- nodes for contour lines will be placed here
+    end
+
+    -- scan all the i-sides, searching for
+    -- intersections with the level surface
+    for i = 0, self.ni - 1 do
+        for j = 0, self.nj do
+            self.done_i[i*si+j] = ((self.coords[    i*st+j].meta > isoval) ==
+                                   (self.coords[(i+1)*st+j].meta > isoval))
+            --[[
+                This is equivalent to testing the two nodes of the i-side:
+
+                    j o-----o
+                      i    i+1
+
+                If meta > isoval in both nodes, or in none of them, then
+                set done_i to true (never look again at this i-side).
+                Otherwise, set done_i to false (we have not finished with
+                this i-side, since it contains one point of a contour line).
+           --]]
+        end
+    end
+
+    -- scan all the j-sides, searching for
+    -- intersections with the level surface
+    for i = 0, self.ni do
+        for j = 0, self.nj - 1 do
+            self.done_j[i*sj+j] = ((self.coords[i*st+j  ].meta > isoval) ==
+                                   (self.coords[i*st+j+1].meta > isoval))
+            --[[
+                This is equivalent to testing the two nodes of the j-side:
+
+                  j+1 o
+                      |
+                      |
+                    j o
+                      i
+
+                Similarly, set done_j accordingly...
+           --]]
+        end
+    end
+
+    if self.debug then
+        self:show_sides(isoval)
+    end
+
+    --[[
+        We now have all the sides containing points of the meta==isoval
+        contour lines. We need to connect the points into distinct lines.
+
+        First of all, scan all the boundary sides, searching for sides
+        with done==false. For each boundary side with done==false, set
+        its done to true and build one contour line, starting from
+        the side itself and the real cell adjacent to it.
+
+        How can we specify a cell and one of its sides?
+        We can specify cell (i,j) and a side (South, West, North, East)
+
+                       N
+                 j o-------o
+                   |       |
+                 W | (i,j) | E
+                   |       |
+               j-1 o-------o
+                  i-1  S   i
+    --]]
+
+    j = 0  -- i-sides with j=0
+    for i = 0, self.ni - 1 do
+        if not self.done_i[i*si+j] then
+            self.done_i[i*si+j] = true
+            self:build_line(isoval, tol, i+1, j+1, 'S', corners)
+            if self.debug then
+                self:show_sides(isoval)
+            end
+        end
+    end
+
+    i = 0  -- j-sides with i=0
+    for j = 0, self.nj - 1 do
+        if not self.done_j[i*sj+j] then
+            self.done_j[i*sj+j] = true
+            self:build_line(isoval, tol, i+1, j+1, 'W', corners)
+            if self.debug then
+                self:show_sides(isoval)
+            end
+        end
+    end
+
+    j = self.nj  -- i-sides with j=nj
+    for i = 0, self.ni - 1 do
+        if not self.done_i[i*si+j] then
+            self.done_i[i*si+j] = true
+            self:build_line(isoval, tol, i+1, j  , 'N', corners)
+            if self.debug then
+                self:show_sides(isoval)
+            end
+        end
+    end
+
+    i = self.ni  -- j-sides with i=ni
+    for j = 0, self.nj - 1 do
+        if not self.done_j[i*sj+j] then
+            self.done_j[i*sj+j] = true
+            self:build_line(isoval, tol, i  , j+1, 'E', corners)
+            if self.debug then
+                self:show_sides(isoval)
+            end
+        end
+    end
+
+    --[[
+        Finally, scan all the internal sides, again searching for sides
+        with done==false. For each internal side with done==false, do _not_
+        alter its done value and build one contour line, starting from
+        the side itself and one of the two real cells adjacent to it.
+    --]]
+
+    -- i-sides with j=1,...,nj-1  (i-sides with j=0 are already done)
+    -- j-sides with i=1,...,ni-1  (j-sides with i=0 are already done)
+    for j = 0, self.nj - 1 do
+        for i = 0, self.ni - 1 do
+            if not self.done_i[i*si+j] then
+                self:build_line(isoval, tol, i+1, j+1, 'S', corners)
+                if self.debug then
+                    self:show_sides(isoval)
+                end
+            end
+            if not self.done_j[i*sj+j] then
+                self:build_line(isoval, tol, i  , j+1, 'E', corners)
+                if self.debug then
+                    self:show_sides(isoval)
+                end
+            end
+        end
+    end
+end
+
+-- build a single contour line
+function PrepcMesh:build_line(isoval, tol, ic, jc, side, corners)
+    -- short names
+    local st = self.st
+    local si = self.si
+    local sj = self.sj
+
+    -- local variables
+    local ia, ja, na, wa
+    local ib, jb, nb, wb
+    local ie, je, ne
+    local id, jd, nd
+    local pt, k, count, next_side
+    local cc, xi_k, eta_k, ck, above, sign_ck
+    local xi_v, eta_v, xi_0, xi_1, eta_0, eta_1
+    local next_done, next_done_idx
+
+    while true do
+        --[[
+            Start from cell (ic,jc) and its South, West, North, or East side
+
+              e         d    b         d    a    N    b    d         b
+               o-------o      o-------o      o-*-----o      o-------o
+               |       |      |       |      |       |      |       *
+               |   c   |    W |   c   |      |   c   |      |   c   | E
+               |       |      *       |      |       |      |       |
+               o----*--o      o-------o      o-------o      o-------o
+              a    S    b    a         e    e         d    e         a
+        --]]
+
+        -- find nodes a and b
+        if side == 'E' then ia = ic     else ia = ic - 1 end
+        if side == 'N' then ja = jc     else ja = jc - 1 end
+        if side == 'W' then ib = ic - 1 else ib = ic     end
+        if side == 'S' then jb = jc - 1 else jb = jc     end
+        if ia < 0 or ib > self.ni or ja < 0 or jb > self.nj then
+            self.newcoords[#self.newcoords + 1] = Coord.new()
+            break  -- out of range: abort
+        end
+        na = self.coords[ia*st+ja]                   -- node a
+        nb = self.coords[ib*st+jb]                   -- node b
+
+        -- compute the intersection point between the side
+        -- and the meta==isoval level surface
+        wa = (isoval - nb.meta)/(na.meta - nb.meta)  -- weights
+        wb = 1 - wa
+
+        pt = Coord.new()
+        for k = 1, 3 do
+            pt.x[k] = wa*na.x[k] + wb*nb.x[k]        -- intersection point
+        end
+        pt.meta = isoval
+        -- add the intersection point to the contour line
+        self.newcoords[#self.newcoords + 1] = pt
+
+        -- check the cell: if it is a phantom cell, then stop
+        if ic < 1 or ic > self.ni or jc < 1 or jc > self.nj then
+            self.newcoords[#self.newcoords + 1] = Coord.new()
+            break  -- the contour line is finished
+        end
+
+        -- find nodes e and d
+        if side == 'W' then ie = ic     else ie = ic - 1 end
+        if side == 'S' then je = jc     else je = jc - 1 end
+        if side == 'E' then id = ic - 1 else id = ic     end
+        if side == 'N' then jd = jc - 1 else jd = jc     end
+        ne = self.coords[ie*st+je]                   -- node e
+        nd = self.coords[id*st+jd]                   -- node d
+
+        -- look at the other three sides of the cell:
+        -- how many of them have done==false ?
+        count = 0
+        if side ~= 'S' and not self.done_i[(ic-1)*si+(jc-1)] then
+            count = count + 1
+            next_side = 'N'  -- next side could be the North side of the
+                             -- adjacent south cell
+        end
+        if side ~= 'W' and not self.done_j[(ic-1)*sj+(jc-1)] then
+            count = count + 1
+            next_side = 'E'  -- next side could be the East side of the
+                             -- adjacent west cell
+        end
+        if side ~= 'N' and not self.done_i[(ic-1)*si+ jc   ] then
+            count = count + 1
+            next_side = 'S'  -- next side could be the South side of the
+                             -- adjacent north cell
+        end
+        if side ~= 'E' and not self.done_j[ ic   *sj+(jc-1)] then
+            count = count + 1
+            next_side = 'W'  -- next side could be the West side of the
+                             -- adjacent east cell
+        end
+
+        --[[
+            if zero other sides have done==false, then stop
+
+                 o     o
+
+
+                 o---*-o
+        --]]
+        if count == 0 then
+            self.newcoords[#self.newcoords + 1] = Coord.new()
+            break  -- the contour line is finished
+        end
+
+        --[[
+            if exactly one other side has done==false, then
+            next_side already specifies the next side to start from
+
+                 o-*---o
+                   |
+                    \
+                 o---*-o
+
+            nothing to be done for the time being...
+        --]]
+
+        --[[
+            the surface inside the cell is thought as the bilinear
+            interpolation of the four nodes and hence mapped onto
+            a unit square in the xi,eta plane:
+
+               eta ^
+                   |               pt.x[k] = (1-xi)*(1-eta)*na.x[k] +
+                  1+------+                     xi *(1-eta)*nb.x[k] +
+                   |e    d|                  (1-xi)*   eta *ne.x[k] +
+                   |      |                     xi *   eta *nd.x[k]
+                   |a    b|
+                  0+------+------>
+                   0      1      xi
+
+            the equation of the contour line is therefore:
+            cc*xi*eta - eta_k*xi - xi_k*eta == ck
+        --]]
+        cc    = na.meta - nb.meta + nd.meta - ne.meta
+        eta_k = na.meta - nb.meta
+        xi_k  = na.meta - ne.meta
+        ck    = isoval  - na.meta
+
+        if math.abs(cc) <=
+           math.max(math.abs(tol*eta_k), math.abs(tol*xi_k)) then
+            --[[
+                cc is negligible and the contour line is
+                basically a straight line:
+                eta_k*xi + xi_k*eta + ck == 0
+            --]]
+            cc = 0
+        else
+            --[[
+                cc is non-zero (and non-negligible) and we can divide
+                both sides of the equation by cc:
+                xi*eta - (eta_k/cc)*xi - (xi_k/cc)*eta == ck/cc
+            --]]
+            xi_k  = xi_k /cc        -- let's rename the coefficients
+            eta_k = eta_k/cc
+            ck    = xi_k*eta_k + ck/cc
+            --[[
+                the contour line is an equilateral hyperbola (in the
+                xi,eta plane) with asymptotes xi == xi_k and eta == eta_k
+                and constant product ck:
+                (xi - xi_k)*(eta - eta_k) == ck
+            --]]
+            if math.abs(ck) <= math.abs(tol/2) then
+                --[[
+                    ck is negligible and the hyperbola basically
+                    degenerates into its asymptotes
+                --]]
+                ck = 0
+            end
+        end
+
+        --[[
+            if two or three other sides have done==false, then
+            choose the next side
+
+                   ?
+                 o-*---o
+                 |     * ?
+               ? *     |
+                 o---*-o
+
+        --]]
+        if count >= 2 then
+            -- how can we choose the next side?
+            --[[
+                if there's more than one side to choose from, cc is
+                necessarily non-zero (a straight line could not cross
+                more than two sides of the unit square) and the
+                center k of the hyperbola is necessarily inside the
+                cell (otherwise one hyperbola branch would be completely
+                outside the cell and the other branch could not cross
+                more than two sides)
+
+                hence, choose the next side by checking the quadrant
+                where the hyperbola branch crossing side a-b lies
+
+                   eta ^
+                       |  :
+                      1+--:---+
+                       |e :  d|
+                    ······+k······
+                       |a :  b|
+                      0+--:---+------>
+                       0  :   1      xi
+            --]]
+            if ck == 0 then      -- math.abs(ck) <= math.abs(tol/2)
+                -- choose side e-d
+                next_side = side
+            elseif ck < 0 then   -- ck < -tol/2
+                -- choose side b-d
+                if side == 'N' or side == 'S' then
+                    next_side = 'W'
+                else
+                    next_side = 'S'
+                end
+            else                 -- ck > +tol/2
+                -- choose side a-e
+                if side == 'N' or side == 'S' then
+                    next_side = 'E'
+                else
+                    next_side = 'N'
+                end
+            end
+        end
+
+        if corners then   -- enhanced corner algorithm
+
+            --[[
+                if, within the cell, the contour line is not straight
+                we can improve its representation by computing an
+                additional point: the vertex v of the hyperbola branch,
+                as long as it lies inside the cell
+            --]]
+            if cc ~= 0 then
+                -- depending on where the center k is, consider the
+                -- vertex above or below k
+                if eta_k < 0 then above = 1 else above = -1 end
+                -- also check the sign of ck
+                if ck < 0 then sign_ck = -1 else sign_ck = 1 end
+                -- compute the vertex v
+                xi_v  =  xi_k + sign_ck*above*math.sqrt(math.abs(ck))
+                eta_v = eta_k +         above*math.sqrt(math.abs(ck))
+
+                --[[
+                    the vertex will be considered inside the cell, as long as
+                    xi_0 <= xi_v <= xi_1   and    eta_0 <= eta_v <= eta_1
+                --]]
+                xi_0  =   - math.abs(tol/10)
+                xi_1  = 1 + math.abs(tol/10)
+                eta_0 =     math.abs(tol)
+                eta_1 = 1 + math.abs(tol/10)
+                -- eta_0 is meant to reject a vertex too close to side a-b
+                -- (which already has a contour line point!)
+                -- we should also reject a vertex too close to the next side
+                if next_side == side then
+                  -- next side is e-d
+                  eta_1 = 1 - math.abs(tol)
+                elseif next_side == 'W' or next_side == 'S' then
+                  -- next side is b-d
+                  xi_1  = 1 - math.abs(tol)
+                else  -- next_side == 'E' or next_side == 'N'
+                  -- next side is a-e
+                  xi_0  =     math.abs(tol)
+                end
+
+                if xi_0  <= xi_v  and xi_v  <= xi_1  and
+                   eta_0 <= eta_v and eta_v <= eta_1 then
+                    -- v is inside the cell: compute its actual coordinates
+                    pt = Coord.new()
+                    for k = 1, 3 do
+                        pt.x[k] = (1-xi_v)*(1-eta_v)*na.x[k] +
+                                     xi_v *(1-eta_v)*nb.x[k] +
+                                  (1-xi_v)*   eta_v *ne.x[k] +
+                                     xi_v *   eta_v *nd.x[k]
+                    end
+                    pt.meta = isoval
+                    -- add the vertex to the contour line
+                    self.newcoords[#self.newcoords + 1] = pt
+                end
+            end
+
+        end
+
+        -- point to done value for the next side and
+        -- determine indexes for the next cell
+        -- (the other cell adjacent to the next side)
+        if next_side == 'W' then ic = ic + 1 end
+        if next_side == 'S' then jc = jc + 1 end
+        if next_side == 'N' or next_side == 'S' then
+            next_done = self.done_i
+            next_done_idx = (ic-1)*si+(jc-1)
+        else
+            next_done = self.done_j
+            next_done_idx = (ic-1)*sj+(jc-1)
+        end
+        if next_side == 'E' then ic = ic - 1 end
+        if next_side == 'N' then jc = jc - 1 end
+
+        -- if next side is already done, then stop
+        if next_done[next_done_idx] then
+            self.newcoords[#self.newcoords + 1] = Coord.new()
+            break  -- the contour line is finished
+        end
+
+        -- set next side done to true and
+        -- iterate starting from this new side and the next cell
+        next_done[next_done_idx] = true
+        side = next_side
+    end
+end
+
+--[[ end core of the program logic ]]--
+
+
+-- print out all the contour lines
+function PrepcMesh:printcontours()
+    for n = 1, #self.newcoords do
+        if self.newcoords[n].x[1] then
+            self.os:write(string.format("%14.6g%14.6g%14.6g%14.6g\n",
+                          self.newcoords[n].x[1],
+                          self.newcoords[n].x[2],
+                          self.newcoords[n].x[3],
+                          self.newcoords[n].meta))
+        else
+            self.os:write('\n')
+        end
+    end
+end
+
+-- print a debug representation of mesh and not-done sides
+function PrepcMesh:show_sides(isoval)
+    -- short names
+    local st = self.st
+    local si = self.si
+    local sj = self.sj
+
+    self.os:write(string.format("# contour level = %.6g \n#\n", isoval))
+    for j = self.nj, 0, -1 do
+        if j <= self.nj - 1 then
+            self.os:write("#    ")
+            for i = 0, self.ni do
+                if self.done_j[i*sj+j] then
+                    self.os:write("  ")
+                else
+                    self.os:write(" *")
+                end
+            end
+            self.os:write("\n")
+        end
+        self.os:write("#     ")
+        for i = 0, self.ni do
+            if self.coords[i*st+j].meta > isoval then
+                self.os:write("•")
+            else
+                self.os:write("o")
+            end
+            if i <= self.ni - 1 then
+                if self.done_i[i*si+j] then
+                    self.os:write(" ")
+                else
+                    self.os:write("*")
+                end
+            end
+        end
+        self.os:write("\n")
+    end
+    self.os:write("#\n")
+end
+
diff --git a/Master/texmf-dist/tex/generic/pgfplots/lua/prepcontour_cli.lua b/Master/texmf-dist/tex/generic/pgfplots/lua/prepcontour_cli.lua
new file mode 100755
index 00000000000..7584228d95c
--- /dev/null
+++ b/Master/texmf-dist/tex/generic/pgfplots/lua/prepcontour_cli.lua
@@ -0,0 +1,97 @@
+#!/usr/bin/luatex
+
+--[[
+
+prepcontour_cli [Lua variant] - prepare contour lines (for pgfplots)
+
+Version:  1.4 (2021-02-22)
+
+Copyright (C) 2020-2021  Francesco Poli <invernomuto@paranoici.org>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+--]]
+
+-- directory separator ('\' for Windows, '/' for all other OSes)
+dsep = package.config:sub(1, 1)
+-- directory where the script lives
+mydir = arg[0]:match("(.*" .. dsep .. ")")
+-- add mydir to the search path for Lua loaders
+package.path = package.path .. ";" .. mydir .. "?.lua"
+-- load the module
+require "prepcontour"
+
+-- main program
+usage = "Usage: " .. arg[0] .. [[
+ CL NB NL {x|y} [VAL]... < IN > OUT
+
+where:
+CL  number of text lines to be copied from the input beginning to the output
+NB  number of text blocks to be read from the input
+NL  number of text lines included in each input block
+x|y ordering ('x' if x varies, 'y' if y varies)
+VAL level for contour lines or
+    min:n:max for n levels between min and max or
+    corners=true  to enable enhanced corner algorithm or
+    corners=false to disable it (this is the default initial state)
+IN  input file
+OUT output file
+]]
+    
+-- simplistic command-line argument parsing
+-- TODO: use a better approach (argparse? some other library?)
+if arg[4] == nil then
+    io.stderr:write(arg[0] .. ": too few arguments\n" .. usage)
+    os.exit(1)
+end
+copylines = tonumber(arg[1])
+nblocks   = tonumber(arg[2])
+nlines    = tonumber(arg[3])
+if copylines == nil or nblocks == nil or nlines == nil then
+    io.stderr:write(arg[0] .. ": failed to parse arguments\n" .. usage)
+    os.exit(2)
+end
+yvaries = arg[4]:sub(1, 1)
+if yvaries ~= 'x' and yvaries ~= 'y' then
+    io.stderr:write(arg[0] .. ": failed to parse argument '"
+                 .. arg[4] .. "'!\n" .. usage)
+    os.exit(3)
+end
+yvaries = (yvaries == 'y')
+
+mesh = PrepcMesh.new(yvaries, nblocks, nlines, copylines)
+corners = false
+for n = 5, #arg do
+    if arg[n] == 'corners=true' then
+        corners = true
+    elseif arg[n] == 'corners=false' then
+        corners = false
+    else
+        level = tonumber(arg[n])
+        if level then
+            mesh:contour(level, corners)
+        else
+            lmin, nlevels, lmax = arg[n]:match('([^%:]+)%:([^%:]+)%:([^%:]+)')
+            nlevels = math.tointeger(tonumber(nlevels))
+            lmin = tonumber(lmin)
+            lmax = tonumber(lmax)
+            if nlevels and lmin and lmax then
+                mesh:autocontour(nlevels, lmin, lmax, corners)
+            end
+        end
+    end
+end
+
+mesh:printcontours()
diff --git a/Master/texmf-dist/tex/generic/pgfplots/numtable/pgfplotstable.code.tex b/Master/texmf-dist/tex/generic/pgfplots/numtable/pgfplotstable.code.tex
index c169f730d42..f8f28883059 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/numtable/pgfplotstable.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/numtable/pgfplotstable.code.tex
@@ -675,6 +675,10 @@
 	/pgfplots/table/create col/linear regression/variance/.initial=,%
 	/pgfplots/table/create col/linear regression/variance list/.initial=,%
 	/pgfplots/table/create col/linear regression/variance src/.initial=,%
+	/pgfplots/table/create col/linear regression/variance format/.is choice,%
+	/pgfplots/table/create col/linear regression/variance format/linear/.code={\def\pgfplotstable@regression@variance@format{l}},%
+	/pgfplots/table/create col/linear regression/variance format/log/.code={\def\pgfplotstable@regression@variance@format{m}},%
+	/pgfplots/table/create col/linear regression/variance format/log,%
 	/pgfplots/table/create col/linear regression/xmode/.initial=,% auto
 	/pgfplots/table/create col/linear regression/ymode/.initial=,% auto
 	/pgfplots/table/create col/quotient/.style={%
@@ -2789,48 +2793,83 @@
 		\pgfplotslistpopfront\pgfplotstable@X\to\pgfplotstable@x
 		\pgfplotslistpopfront\pgfplotstable@Y\to\pgfplotstable@y
 		%
-		\pgfplotstableparsex{\pgfplotstable@x}%
-		\let\pgfplotstable@x=\pgfmathresult
-		\expandafter\pgfplotslistpushback\pgfmathresult\to\pgfplotstable@Xparsed
-		\pgfplotstableparsey{\pgfplotstable@y}%
-		\let\pgfplotstable@y=\pgfmathresult
-		%
-		\pgfplotslistcheckempty\pgfplotstable@VARIANCE
-		\ifpgfplotslistempty
-			\pgfmathfloatcreate{1}{1.0}{0}%
-			\let\pgfplotstable@invsqr=\pgfmathresult
+		% if either the x or y entry is empty, skip this point
+		\ifx\pgfplotstable@x\pgfutil@empty
+		\else\ifx\pgfplotstable@y\pgfutil@empty
 		\else
-			\pgfplotslistpopfront\pgfplotstable@VARIANCE\to\pgfplotstable@variance
-			\pgfmathfloatparsenumber{\pgfplotstable@variance}%
-			\let\pgfplotstable@variance=\pgfmathresult
+			\pgfplotstableparsex{\pgfplotstable@x}%
+			\let\pgfplotstable@x=\pgfmathresult
+			\expandafter\pgfplotslistpushback\pgfmathresult\to\pgfplotstable@Xparsed
+			\pgfplotstableparsey{\pgfplotstable@y}%
+			\let\pgfplotstable@y=\pgfmathresult
+			%
+			\pgfplotslistcheckempty\pgfplotstable@VARIANCE
+			\ifpgfplotslistempty
+				\pgfmathfloatcreate{1}{1.0}{0}%
+				\let\pgfplotstable@variance=\pgfmathresult
+			\else
+				\pgfplotslistpopfront\pgfplotstable@VARIANCE\to\pgfplotstable@variance
+				\pgfmathfloatparsenumber{\pgfplotstable@variance}%
+				\let\pgfplotstable@variance=\pgfmathresult
+			\fi
+			%
+			% if we're in y log mode, then instead of y +- dy, we have log_b(y +- dy)
+			% which Taylor expands to log_b(y) +- dy/(y ln(b)) + O(dy^2)
+			% However, this might introduce a division-by-zero in the exceptional case
+			% where b=1 or log_b(0) = -infinity, so instead of computing 1/(dy/(y ln b))^2
+			% we will compute (y ln b)^2*1/dy^2.  Note that we don't factor out the ^2
+			% because floating point operations don't commute and the historical behavior
+			% was to compute 1/dy^2.
+			%
 			\pgfmathfloatmultiply@{\pgfplotstable@variance}{\pgfplotstable@variance}%
 			\let\pgfplotstable@sqr=\pgfmathresult
 			\pgfmathfloatreciprocal@{\pgfplotstable@sqr}%
 			\let\pgfplotstable@invsqr=\pgfmathresult
-		\fi
-		%
-		\pgfmathfloatadd@{\pgfplotstable@S}{\pgfplotstable@invsqr}%
-		\let\pgfplotstable@S=\pgfmathresult
-		%
-		\pgfmathfloatmultiply@{\pgfplotstable@x}{\pgfplotstable@invsqr}%
-		\let\pgfplots@table@accum=\pgfmathresult
-		\pgfmathfloatadd@{\pgfplotstable@Sx}{\pgfplots@table@accum}%
-		\let\pgfplotstable@Sx=\pgfmathresult
-		%
-		\pgfmathfloatmultiply@{\pgfplotstable@x}{\pgfplots@table@accum}%
-		\let\pgfplots@table@accum=\pgfmathresult
-		\pgfmathfloatadd@{\pgfplotstable@Sxx}{\pgfplots@table@accum}%
-		\let\pgfplotstable@Sxx=\pgfmathresult
-		%
-		\pgfmathfloatmultiply@{\pgfplotstable@y}{\pgfplotstable@invsqr}%
-		\let\pgfplots@table@accum=\pgfmathresult
-		\pgfmathfloatadd@{\pgfplotstable@Sy}{\pgfplots@table@accum}%
-		\let\pgfplotstable@Sy=\pgfmathresult
-		%
-		\pgfmathfloatmultiply@{\pgfplotstable@x}{\pgfplots@table@accum}%
-		\let\pgfplots@table@accum=\pgfmathresult
-		\pgfmathfloatadd@{\pgfplotstable@Sxy}{\pgfplots@table@accum}%
-		\let\pgfplotstable@Sxy=\pgfmathresult
+			\ifx\pgfplotstableparseylogbase\pgfutil@empty
+			\else
+				\if m\pgfplotstable@regression@variance@format
+					% variance format=log
+				\else
+					% variance format=linear
+					% N.B. at this point, \pgfplotstable@invsqr, \pgfplotstable@sqr, and \pgfplotstable@variance are
+					% all missing log-scaling.  We update \pgfplotstable@invsqr now, but leave \pgfplotstable@sqr
+					% and \pgfplotstable@variance alone because they are not used below this point.
+					\pgfplotstableparseyinv@{\pgfplotstable@y}% y
+					\let\pgfplotstable@tmp=\pgfmathresult
+					\pgfmathfloatmultiply@{\pgfplotstable@tmp}{\pgfplotstableparseylogofbasefloat}% y * ln(b)
+					\let\pgfplotstable@ylnb=\pgfmathresult
+					\pgfmathfloatmultiply@{\pgfplotstable@ylnb}{\pgfplotstable@ylnb}% (y * ln(b))^2
+					\let\pgfplotstable@ylnbsqr=\pgfmathresult
+					% since \pgfplotstable@invsqr is used, but \pgfplotstable@sqr and \pgfplotstable@variance are not
+					% we only need to update \pgfplotstable@invsqr
+					\pgfmathfloatmultiply@{\pgfplotstable@invsqr}{\pgfplotstable@ylnbsqr}% (y * ln(b))^2 / dy^2
+					\let\pgfplotstable@invsqr=\pgfmathresult
+				\fi
+			\fi
+			%
+			\pgfmathfloatadd@{\pgfplotstable@S}{\pgfplotstable@invsqr}%
+			\let\pgfplotstable@S=\pgfmathresult
+			%
+			\pgfmathfloatmultiply@{\pgfplotstable@x}{\pgfplotstable@invsqr}%
+			\let\pgfplots@table@accum=\pgfmathresult
+			\pgfmathfloatadd@{\pgfplotstable@Sx}{\pgfplots@table@accum}%
+			\let\pgfplotstable@Sx=\pgfmathresult
+			%
+			\pgfmathfloatmultiply@{\pgfplotstable@x}{\pgfplots@table@accum}%
+			\let\pgfplots@table@accum=\pgfmathresult
+			\pgfmathfloatadd@{\pgfplotstable@Sxx}{\pgfplots@table@accum}%
+			\let\pgfplotstable@Sxx=\pgfmathresult
+			%
+			\pgfmathfloatmultiply@{\pgfplotstable@y}{\pgfplotstable@invsqr}%
+			\let\pgfplots@table@accum=\pgfmathresult
+			\pgfmathfloatadd@{\pgfplotstable@Sy}{\pgfplots@table@accum}%
+			\let\pgfplotstable@Sy=\pgfmathresult
+			%
+			\pgfmathfloatmultiply@{\pgfplotstable@x}{\pgfplots@table@accum}%
+			\let\pgfplots@table@accum=\pgfmathresult
+			\pgfmathfloatadd@{\pgfplotstable@Sxy}{\pgfplots@table@accum}%
+			\let\pgfplotstable@Sxy=\pgfmathresult
+		\fi\fi
 	\pgfutil@repeat
 	%
 	\pgfmathparse{\pgfplotstable@S * \pgfplotstable@Sxx - \pgfplotstable@Sx *\pgfplotstable@Sx}%
diff --git a/Master/texmf-dist/tex/generic/pgfplots/pgfcontrib/tikzlibrarydecorations.softclip.code.tex b/Master/texmf-dist/tex/generic/pgfplots/pgfcontrib/tikzlibrarydecorations.softclip.code.tex
index 7544d66d313..3b604760ca6 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/pgfcontrib/tikzlibrarydecorations.softclip.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/pgfcontrib/tikzlibrarydecorations.softclip.code.tex
@@ -433,6 +433,8 @@
 }
 
 \def\tikzlibsoftclip@setkey@@{%
+	\edef\tikzlibsoftclip@setkey@@reset{\noexpand\tikz@expandcount=\the\tikz@expandcount\relax}%
+	\tikz@resetexpandcount
 	%
 	\tikz@scan@one@point\tikzlibsoftclip@setkey@bb@scan@a
 }%
@@ -455,6 +457,7 @@
 \def\tikzlibsoftclip@setkey@bb@scan@b#1{%
 	\def\tikzlibsoftclip@setkey@bb@b{#1}%
 	\pgfutil@ifnextchar \pgf@stop{%
+		\tikzlibsoftclip@setkey@@reset
 		\tikzlibsoftclip@setkey@@activate
 		\tikz@gobble@until@stop
 	}{%
diff --git a/Master/texmf-dist/tex/generic/pgfplots/pgfplots.code.tex b/Master/texmf-dist/tex/generic/pgfplots/pgfplots.code.tex
index 6b5b8b808da..b02206937e7 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/pgfplots.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/pgfplots.code.tex
@@ -45,6 +45,7 @@
 
 
 \usetikzlibrary{decorations,decorations.pathmorphing,decorations.pathreplacing}
+\usepgfplotslibrary{contourlua}%
 
 % FIXME: reduce number of variables!
 
@@ -4318,8 +4319,10 @@
 	/pgfplots/compat/anchors/1.15/.style= 	{/pgfplots/compat/anchors/1.13},%
 	/pgfplots/compat/anchors/1.16/.style= 	{/pgfplots/compat/anchors/1.13},%
 	/pgfplots/compat/anchors/1.17/.style= 	{/pgfplots/compat/anchors/1.13},%
+	/pgfplots/compat/anchors/1.18/.style= 	{/pgfplots/compat/anchors/1.13},%
 	/pgfplots/compat/anchors/default/.style={/pgfplots/compat/anchors/1.3},%
 	%
+	%
 	/pgfplots/compat/empty line/.is choice,
 	/pgfplots/compat/empty line/pre 1.3/.code={\pgfplots@emptyline@compattrue},% FIXME: WAS \global
 	/pgfplots/compat/empty line/1.3/.style= {/pgfplots/compat/empty line/pre 1.3},%
@@ -4338,6 +4341,7 @@
 	/pgfplots/compat/empty line/1.15/.style= 	{/pgfplots/compat/empty line/1.4},%
 	/pgfplots/compat/empty line/1.16/.style= 	{/pgfplots/compat/empty line/1.4},%
 	/pgfplots/compat/empty line/1.17/.style= 	{/pgfplots/compat/empty line/1.4},%
+	/pgfplots/compat/empty line/1.18/.style= 	{/pgfplots/compat/empty line/1.4},%
 	/pgfplots/compat/empty line/default/.style=	{/pgfplots/compat/empty line/1.4},%
 	%
 	/pgfplots/compat/path replacement/.is choice,
@@ -4358,6 +4362,7 @@
 	/pgfplots/compat/path replacement/1.15/.style=   	{/pgfplots/compat/path replacement/1.5.1},
 	/pgfplots/compat/path replacement/1.16/.style=   	{/pgfplots/compat/path replacement/1.5.1},
 	/pgfplots/compat/path replacement/1.17/.style=   	{/pgfplots/compat/path replacement/1.5.1},
+	/pgfplots/compat/path replacement/1.18/.style=   	{/pgfplots/compat/path replacement/1.5.1},
 	/pgfplots/compat/path replacement/default/.style=	{/pgfplots/compat/path replacement/pre 1.3},%
 	%
 	/pgfplots/compat/pgfpoint substitution/.is choice,
@@ -4378,6 +4383,7 @@
 	/pgfplots/compat/pgfpoint substitution/1.15/.style=   	{/pgfplots/compat/pgfpoint substitution/1.11},
 	/pgfplots/compat/pgfpoint substitution/1.16/.style=   	{/pgfplots/compat/pgfpoint substitution/1.11},
 	/pgfplots/compat/pgfpoint substitution/1.17/.style=   	{/pgfplots/compat/pgfpoint substitution/1.11},
+	/pgfplots/compat/pgfpoint substitution/1.18/.style=   	{/pgfplots/compat/pgfpoint substitution/1.11},
 	/pgfplots/compat/pgfpoint substitution/default/.style=	{/pgfplots/compat/pgfpoint substitution/pre 1.3},%
 	%ellipse/.is if=pgfplots@path@replace@ellipse,
 	%
@@ -4408,6 +4414,7 @@
 	/pgfplots/compat/labels/1.15/.style=		{/pgfplots/compat/labels/1.8},%
 	/pgfplots/compat/labels/1.16/.style=		{/pgfplots/compat/labels/1.8},%
 	/pgfplots/compat/labels/1.17/.style=		{/pgfplots/compat/labels/1.8},%
+	/pgfplots/compat/labels/1.18/.style=		{/pgfplots/compat/labels/1.8},%
 	/pgfplots/compat/labels/default/.style=	{/pgfplots/compat/labels/pre 1.3},% maintain backwards compatibility
 	%
 	/pgfplots/compat/bar nodes/.is choice,
@@ -4463,6 +4470,7 @@
 	/pgfplots/compat/bar nodes/1.15/.style=		{/pgfplots/compat/bar nodes/1.13},%
 	/pgfplots/compat/bar nodes/1.16/.style=		{/pgfplots/compat/bar nodes/1.13},%
 	/pgfplots/compat/bar nodes/1.17/.style=		{/pgfplots/compat/bar nodes/1.13},%
+	/pgfplots/compat/bar nodes/1.18/.style=		{/pgfplots/compat/bar nodes/1.13},%
 	/pgfplots/compat/bar nodes/default/.style=	{/pgfplots/compat/bar nodes/pre 1.3},% maintain backwards compatibility
 	/pgfplots/compat/bar nodes/default,
 	%
@@ -4498,6 +4506,7 @@
 	/pgfplots/compat/scaling/1.15/.style={/pgfplots/compat/scaling/1.6},
 	/pgfplots/compat/scaling/1.16/.style={/pgfplots/compat/scaling/1.6},
 	/pgfplots/compat/scaling/1.17/.style={/pgfplots/compat/scaling/1.6},
+	/pgfplots/compat/scaling/1.18/.style={/pgfplots/compat/scaling/1.6},
 	/pgfplots/compat/scaling/default/.style={/pgfplots/compat/scaling/1.4},
 	%
 	/pgfplots/compat/scale mode/.is choice,
@@ -4518,6 +4527,7 @@
 	/pgfplots/compat/scale mode/1.15/.style={/pgfplots/compat/scale mode/1.6},
 	/pgfplots/compat/scale mode/1.16/.style={/pgfplots/compat/scale mode/1.6},
 	/pgfplots/compat/scale mode/1.17/.style={/pgfplots/compat/scale mode/1.6},
+	/pgfplots/compat/scale mode/1.18/.style={/pgfplots/compat/scale mode/1.6},
 	/pgfplots/compat/scale mode/default/.style={/pgfplots/compat/scale mode/1.5},
 	%
 	%
@@ -4539,6 +4549,7 @@
 	/pgfplots/compat/plot3graphics/1.15/.style=	{/pgfplots/compat/plot3graphics/1.6},%
 	/pgfplots/compat/plot3graphics/1.16/.style=	{/pgfplots/compat/plot3graphics/1.6},%
 	/pgfplots/compat/plot3graphics/1.17/.style=	{/pgfplots/compat/plot3graphics/1.6},%
+	/pgfplots/compat/plot3graphics/1.18/.style=	{/pgfplots/compat/plot3graphics/1.6},%
 	/pgfplots/compat/plot3graphics/default/.style=	{/pgfplots/compat/plot3graphics/1.5},%
 	%
 	% ATTENTION: there is a compatibility issue which slipped through
@@ -4586,6 +4597,7 @@
 	/pgfplots/compat/bar width by units/1.15/.style=		{/pgfplots/compat/bar width by units/1.7},
 	/pgfplots/compat/bar width by units/1.16/.style=		{/pgfplots/compat/bar width by units/1.7},
 	/pgfplots/compat/bar width by units/1.17/.style=		{/pgfplots/compat/bar width by units/1.7},
+	/pgfplots/compat/bar width by units/1.18/.style=		{/pgfplots/compat/bar width by units/1.7},
 	/pgfplots/compat/bar width by units/default/.style={/pgfplots/compat/bar width by units/pre 1.3},
 	%
 	/pgfplots/compat/BB/.is choice,
@@ -4606,6 +4618,7 @@
 	/pgfplots/compat/BB/1.15/.style=	{/pgfplots/compat/BB/1.8},%
 	/pgfplots/compat/BB/1.16/.style=	{/pgfplots/compat/BB/1.8},%
 	/pgfplots/compat/BB/1.17/.style=	{/pgfplots/compat/BB/1.8},%
+	/pgfplots/compat/BB/1.18/.style=	{/pgfplots/compat/BB/1.8},%
 	/pgfplots/compat/BB/default/.style={/pgfplots/compat/BB/pre 1.3},
 	%
 	/pgfplots/compat/general/.is choice,
@@ -4670,6 +4683,7 @@
 		/pgfplots/compat/general/1.15,
 		%/pgfplots/table/percent is letter=true,
 	},%
+	/pgfplots/compat/general/1.18/.style=	{/pgfplots/compat/general/1.17},%
 	/pgfplots/compat/general/default/.style={/pgfplots/compat/general/1.11},
 	%
 	/pgfplots/enable tick line clipping/.is if=pgfplots@enable@tick@line@clipping,
@@ -4687,7 +4701,7 @@
 	/pgfplots/compat/show suggested version/.is if=pgfplots@show@suggested@version,
 	/pgfplots/compat/show suggested version=true,
 	/pgfplots/compat/current/.initial=,
-	/pgfplots/compat/mostrecent/.initial=1.17,
+	/pgfplots/compat/mostrecent/.initial=1.18,
 	/utils/exec={%
 		\pgfplotsutilforeachcommasep{%
 			/pgfplots/compat/current,%
@@ -7925,13 +7939,13 @@
 			{\csname pgfplots@environment@#2axis@\endcsname[]}%
 	}%
 	\expandafter\long\expandafter\def\csname pgfplots@environment@#2axis@\endcsname[##1]{%
-		\begin{axis}[%
+		\axis[%
 			##1,%
 			data cs=#1,%  assume that coordinates provided by \addplot are in the correct coordsystem
 			axis type=#2% active the correct axis type
 		]%
 	}%
-	\expandafter\def\csname endpgfplots@environment@#2axis\endcsname{\end{axis}}%
+	\expandafter\def\csname endpgfplots@environment@#2axis\endcsname{\endaxis}%
 	\edef\pgfplots@glob@TMPc{\expandafter\noexpand\csname #2axis\endcsname\expandafter\noexpand\csname pgfplots@environment@#2axis\endcsname}%
 	\expandafter
 	\tikzaddtikzonlycommandshortcutlet\pgfplots@glob@TMPc
@@ -11571,7 +11585,7 @@
 			\pgfplotsutil@directlua{%
 				pgfplots.gca = pgfplots.Axis.new();
 				pgfplots.gca.is3d = \pgfplots@boolval{pgfplots@threedim};
-				pgfplots.gca.clipLimits = { 
+				pgfplots.gca.clipLimits = {
 					\pgfplots@boolval{pgfplots@clip@limits@x},
 					\pgfplots@boolval{pgfplots@clip@limits@y},
 					\pgfplots@boolval{pgfplots@clip@limits@z} };
diff --git a/Master/texmf-dist/tex/generic/pgfplots/pgfplots.revision.tex b/Master/texmf-dist/tex/generic/pgfplots/pgfplots.revision.tex
index 2bf436fe808..be95dfac86a 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/pgfplots.revision.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/pgfplots.revision.tex
@@ -5,10 +5,10 @@
 \catcode`\:=12
 \catcode`\+=12
 \catcode`\-=12
-\gdef\pgfplotsrevision{1.17}
-\gdef\pgfplotsversion{1.17}
-\gdef\pgfplotsversiondatetime{2020-02-29 09:43:25 +0100}
-\gdef\pgfplotsrevisiondatetime{2020-02-29 09:43:25 +0100}
+\gdef\pgfplotsrevision{1.18}
+\gdef\pgfplotsversion{1.18}
+\gdef\pgfplotsversiondatetime{2021-05-08 12:48:32 +0200}
+\gdef\pgfplotsrevisiondatetime{2021-05-08 12:48:32 +0200}
 \gdef\pgfplots@glob@TMPa#1-#2-#3 #4\relax{#1/#2/#3}
 \xdef\pgfplotsversiondate{\expandafter\pgfplots@glob@TMPa\pgfplotsversiondatetime\relax}
 \xdef\pgfplotsrevisiondate{\expandafter\pgfplots@glob@TMPa\pgfplotsrevisiondatetime\relax}
diff --git a/Master/texmf-dist/tex/generic/pgfplots/pgfplotsmeshplotimage.code.tex b/Master/texmf-dist/tex/generic/pgfplots/pgfplotsmeshplotimage.code.tex
index b5a51cd3b2a..3b5115bf9a2 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/pgfplotsmeshplotimage.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/pgfplotsmeshplotimage.code.tex
@@ -106,7 +106,10 @@
 		\else
 			\ifx\pgfplotsplothandlermesh@image@lastA\pgfutil@empty
 			\else
-				\pgfplotscoordmath{x}{parse}{\pgfplotsplothandlermesh@image@lastA - 0.5*(\pgfplotsplothandlermesh@cur - \pgfplotsplothandlermesh@image@lastA)}%
+				\pgfplotscoordmath{\pgfplotsplothandlermesh@image@A@dir}{parse}{\pgfplotsplothandlermesh@image@lastA - 0.5*(\pgfplotsplothandlermesh@cur - \pgfplotsplothandlermesh@image@lastA)}%
+				\pgfplots@if{pgfplots@\pgfplotsplothandlermesh@image@A@dir islinear}{}{%
+					\pgfplotscoordmath{\pgfplotsplothandlermesh@image@A@dir}{parsenumber}{\pgfmathresult}%
+				}
 				\let\pgfplots@current@point@xh=\pgfmathresult
 				\pgfplotsplothandlermesh@image@updatelimits@AB\pgfplots@current@point@xh\pgfplotsplothandlermesh@cur@B\pgfplots@current@point@z
 				%
@@ -120,12 +123,18 @@
 		\else
 			\ifx\pgfplotsplothandlermesh@image@lastB\pgfutil@empty
 			\else
-				\pgfplotscoordmath{y}{parse}{\pgfplotsplothandlermesh@image@lastB - 0.5*(\pgfplotsplothandlermesh@cur@B - \pgfplotsplothandlermesh@image@lastB)}%
+				\pgfplotscoordmath{\pgfplotsplothandlermesh@image@B@dir}{parse}{\pgfplotsplothandlermesh@image@lastB - 0.5*(\pgfplotsplothandlermesh@cur@B - \pgfplotsplothandlermesh@image@lastB)}%
+				\pgfplots@if{pgfplots@\pgfplotsplothandlermesh@image@B@dir islinear}{}{%
+					\pgfplotscoordmath{\pgfplotsplothandlermesh@image@B@dir}{parsenumber}{\pgfmathresult}%
+				}
 				\let\pgfplots@current@point@yh=\pgfmathresult
 				\pgfplotsplothandlermesh@image@updatelimits@AB\pgfplotsplothandlermesh@cur\pgfplots@current@point@yh\pgfplots@current@point@z
 				%
 				%
-				\pgfplotscoordmath{y}{parse}{\pgfplotsplothandlermesh@cur@B + 0.5*(\pgfplotsplothandlermesh@cur@B - \pgfplotsplothandlermesh@image@lastB)}%
+				\pgfplotscoordmath{\pgfplotsplothandlermesh@image@B@dir}{parse}{\pgfplotsplothandlermesh@cur@B + 0.5*(\pgfplotsplothandlermesh@cur@B - \pgfplotsplothandlermesh@image@lastB)}%
+				\pgfplots@if{pgfplots@\pgfplotsplothandlermesh@image@B@dir islinear}{}{%
+					\pgfplotscoordmath{\pgfplotsplothandlermesh@image@B@dir}{parsenumber}{\pgfmathresult}%
+				}
 				\let\pgfplots@current@point@yh=\pgfmathresult
 				\pgfplotsplothandlermesh@image@updatelimits@AB\pgfplotsplothandlermesh@cur\pgfplots@current@point@yh\pgfplots@current@point@z
 				%
@@ -145,7 +154,10 @@
 		\ifx\pgfplotsplothandlermesh@image@lastA\pgfutil@empty
 		\else
 			\expandafter\let\expandafter\pgfplotsplothandlermesh@cur@B\csname pgfplots@current@point@\pgfplotsplothandlermesh@image@B@dir\endcsname
-			\pgfplotscoordmath{x}{parse}{\pgfplotsplothandlermesh@image@lastA + 0.5*(\pgfplotsplothandlermesh@image@lastA - \pgfplotsplothandlermesh@image@lastlastA)}%
+			\pgfplotscoordmath{\pgfplotsplothandlermesh@image@A@dir}{parse}{\pgfplotsplothandlermesh@image@lastA + 0.5*(\pgfplotsplothandlermesh@image@lastA - \pgfplotsplothandlermesh@image@lastlastA)}%
+			\pgfplots@if{pgfplots@\pgfplotsplothandlermesh@image@A@dir islinear}{}{%
+				\pgfplotscoordmath{\pgfplotsplothandlermesh@image@A@dir}{parsenumber}{\pgfmathresult}%
+			}
 			\let\pgfplots@current@point@xh=\pgfmathresult
 			\pgfplotsplothandlermesh@image@updatelimits@AB\pgfplots@current@point@xh\pgfplotsplothandlermesh@cur@B\pgfplots@current@point@z
 		\fi
diff --git a/Master/texmf-dist/tex/generic/pgfplots/pgfplotsplothandlers.code.tex b/Master/texmf-dist/tex/generic/pgfplots/pgfplotsplothandlers.code.tex
index bad6354b2f3..4f138c8077e 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/pgfplotsplothandlers.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/pgfplotsplothandlers.code.tex
@@ -1952,6 +1952,7 @@
 
 \newif\ifpgfplotsplothandlercontour@labels
 \newif\ifpgfplotsplothandlercontour@filled
+\newif\ifpgfplotsplothandlercontour@corners
 
 \def\pgfplotsplothandlercontour@axis@set@inverse#1#2#3#4{%
 	\pgfplotsutilifstringequal{#1}{#4}{%
@@ -2041,6 +2042,8 @@
 	contour/contour label style/.style={
 		/pgfplots/contour/every contour label/.append style={#1}},
 	%
+	contour/corners/.is if=pgfplotsplothandlercontour@corners,
+	contour/corners=false,
 	%
 	%  Styles to actually *compute* the contour.
 	%  These are mostly placeholders here: As long as the
@@ -2090,6 +2093,16 @@
 	contour external/script extension/.initial=script,
 	contour external/script/.initial=,% not yet initialised
 	contour external/cmd/.initial=,% not yet initialised
+	contour external/script type/.is choice=,%
+	contour external/script type/system call/.code={\def\pgfplotsplothandlercontourexternal@script@type{S}},%
+	contour external/script type/lua/.code={\def\pgfplotsplothandlercontourexternal@script@type{L}\pgfkeyssetvalue{/pgfplots/contour external/script extension}{lua}},%
+	contour external/script type/system call,%
+	contour external/script input format/.is choice,%
+	contour external/script input format/x y meta meta/.code={\def\pgfplotsplothandlercontourexternal@script@input@format{M}},%
+	contour external/script input format/x y z meta/.code={\def\pgfplotsplothandlercontourexternal@script@input@format{Z}},%
+	contour external/script input format=x y meta meta,%
+	contour external/@before cmd/.initial=,
+	contour external/@after cmd/.initial=,
 	contour external/.search also=/pgfplots/contour,
 	contour gnuplot/.style={
 		contour external={%
@@ -2684,6 +2697,7 @@
 	\else
 		\pgfplotsset{plot to file/scanline marks/if in input}% the choice 'always' is unaware of existing end-of-scanline marks
 	\fi
+	%
 	\pgfplotsplothandlertofile{\pgfplotsplothandlercontourexternal@file.dat}%
 	\pgfplotsplothandlersurveystart
 	\let\pgfplotsplothandlersurveypoint@tofile=\pgfplotsplothandlersurveypoint
@@ -2694,6 +2708,10 @@
 	\pgfplotsplothandlersurveypoint@contour@prepare@axes x%
 	\pgfplotsplothandlersurveypoint@contour@prepare@axes y%
 	\pgfplotsplothandlersurveypoint@contour@prepare@axes z%
+	%
+	\global\let\pgfplotsplothandlercontourexternal@metamin=\relax
+	\global\let\pgfplotsplothandlercontourexternal@metamax=\relax
+	%
 }%
 
 % #1: either x, y, or z
@@ -2745,22 +2763,39 @@
 }%
 
 \def\pgfplotsplothandlersurveypoint@contourexternal{%
-	% temporarily disable updates to point meta limits. They should be
-	% updated during the contour prepared processing.
-	\let\pgfplotsplothandlercontourexternal@orig@perpointlimits@limits=\pgfplots@set@perpointmeta@limits
-	\let\pgfplots@set@perpointmeta@limits=\relax
+	\ifx\pgfplotsplothandlercontourexternal@metamin\relax
+		% remember values. We want to restore them right before we parse the output of contour external:
+		\global\let\pgfplotsplothandlercontourexternal@metamin=\pgfplots@metamin
+		\global\let\pgfplotsplothandlercontourexternal@metamax=\pgfplots@metamax
+	\fi
+	\pgfplotsplothandlersurveypoint@contour@axes@std@to@reordered
 	%
 	\pgfplotsaxissurveysetpointmeta
-	\let\pgfplots@current@point@z=\pgfplots@current@point@meta
 	%
-	\let\pgfplots@set@perpointmeta@limits=\pgfplotsplothandlercontourexternal@orig@perpointlimits@limits
+	\if M\pgfplotsplothandlercontourexternal@script@input@format
+		\let\pgfplots@current@point@z=\pgfplots@current@point@meta
+	\fi
 	%
-	\pgfplotsplothandlersurveypoint@contour@axes@std@to@reordered
 	%
-%\message{contour gnuplot: collecting point (\pgfplots@current@point@x,\pgfplots@current@point@y,\pgfplots@current@point@z) [\pgfplots@current@point@meta]...}%
+%\message{contour external collecting point (\pgfplots@current@point@x,\pgfplots@current@point@y,\pgfplots@current@point@z) [\pgfplots@current@point@meta]...}%
 	\pgfplotsplothandlersurveypoint@tofile
 }%
 
+% Utility function which is NOT normally invoked. It can be invoked inside of 'contour external/@before cmd'. Should it become the default?
+\def\pgfplotscontourpopulateallkeys{%
+	\pgfkeysgetvalue{/pgfplots/contour/levels}\pgfplots@loc@TMPa
+	\ifx\pgfplots@loc@TMPa\pgfutil@empty
+		\ifx\pgfplots@metamax\pgfutil@empty
+		\else
+			\edef\pgfplots@loc@TMPa{\pgfplots@metamin:\pgfplots@metamax}%
+			\edef\thecontournumber{\pgfkeysvalueof{/pgfplots/contour/number}}%
+			\expandafter\pgfplots@domain@to@foreach\pgfplots@loc@TMPa\relax{\thecontournumber}%
+			\expandafter\pgfplots@foreach@to@list@macro\expandafter{\pgfplotsretval}%
+			\pgfkeyslet{/pgfplots/contour/levels}\pgfplotsretval
+		\fi
+	\fi
+}%
+
 {
 \catcode`\`=12
 \catcode`\'=12
@@ -2777,6 +2812,7 @@
 		\fi
 	\fi
 	%
+	\pgfkeysvalueof{/pgfplots/contour external/@before cmd}\relax%
 	\begingroup
 	\let\numcoords=\pgfplots@current@point@coordindex%
 	\pgfplotsautocompletemeshkeys
@@ -2790,16 +2826,46 @@
 	\edef\ordering{\pgfplots@plot@mesh@ordering}%
 	\edef\infile{\pgfplotsplothandlercontourexternal@file.dat}%
 	\edef\outfile{\pgfplotsplothandlercontourexternal@file.table}%
-	\edef\script{\pgfplotsplothandlercontourexternal@file.\pgfkeysvalueof{/pgfplots/contour external/script extension}}%
+	%
+	\pgfkeysgetvalue{/pgfplots/contour external/cmd}\cmd
+	\pgfkeysgetvalue{/pgfplots/contour external/script extension}\scriptext
+	\if S\pgfplotsplothandlercontourexternal@script@type
+	\else
+		% script type=lua
+		\ifx\cmd\pgfutil@empty
+			\pgfkeyssetvalue{/pgfplots/contour external/cmd}{require(\"\scriptbase\")}%
+			\pgfkeysgetvalue{/pgfplots/contour external/cmd}\cmd
+		\fi
+	\fi
+	%
+	\pgfplotscoordmath{meta}{tofixed}{\pgfplots@metamin}%
+	\let\pgfplotsmetamin=\pgfmathresult
+	\pgfplotscoordmath{meta}{tofixed}{\pgfplots@metamax}%
+	\let\pgfplotsmetamax=\pgfmathresult
+	%
+	\edef\script{\pgfplotsplothandlercontourexternal@file.\scriptext}%
+	\edef\scriptbase{\pgfplotsplothandlercontourexternal@file}%
 	\edef\thecontourlevels{\pgfkeysvalueof{/pgfplots/contour/levels}}%
 	\edef\thecontournumber{\pgfkeysvalueof{/pgfplots/contour/number}}%
 	\immediate\openout\w@pgf@writea=\script
 	\immediate\write\w@pgf@writea{\pgfkeysvalueof{/pgfplots/contour external/script}}%
 	\immediate\closeout\w@pgf@writea
 	%
-	\pgfplots@shellescape{\pgfkeysvalueof{/pgfplots/contour external/cmd}}%
+	\if S\pgfplotsplothandlercontourexternal@script@type
+		\pgfplots@shellescape{\cmd}%
+	\else
+		\ifpgfplots@LUA@supported
+			\pgfplotsutil@directlua{\pgfkeysvalueof{/pgfplots/contour external/cmd}}%
+		\else
+			\pgfplots@error{Calls to 'lua' are currently unsupported. Please invoke lualatex instead.}%
+		\fi
+	\fi
+	\pgfkeysvalueof{/pgfplots/contour external/@after cmd}\relax%
 	\endgroup
 	%
+	% RESET point meta computation here such that it can be collected from the 'contour prepared' stream (i.e. the output of contour external):
+	\global\let\pgfplots@metamin=\pgfplotsplothandlercontourexternal@metamin
+	\global\let\pgfplots@metamax=\pgfplotsplothandlercontourexternal@metamax
 	%
 	\pgfplotsplothandlercontourprepared
 	% the PREPARE steps in the coord stream start/end have already
diff --git a/Master/texmf-dist/tex/generic/pgfplots/pgfplotsstackedplots.code.tex b/Master/texmf-dist/tex/generic/pgfplots/pgfplotsstackedplots.code.tex
index 7f45bfb5bef..c9e1c15ac20 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/pgfplotsstackedplots.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/pgfplotsstackedplots.code.tex
@@ -686,3 +686,8 @@
 	\def\pgf@plotyzerolevelstreamend{%
 	}%
 }%
+
+% Initialize \pgfplots@stacked@diff to 0
+% Otherwise, if the first point in a plot is NaN \pgfplots@stacked@diff will be
+% undefined and an error is raised, see issue #219
+\def\pgfplots@stacked@diff{0}
diff --git a/Master/texmf-dist/tex/generic/pgfplots/util/pgfplotsutil.code.tex b/Master/texmf-dist/tex/generic/pgfplots/util/pgfplotsutil.code.tex
index 353a04ea945..5a014d44f23 100644
--- a/Master/texmf-dist/tex/generic/pgfplots/util/pgfplotsutil.code.tex
+++ b/Master/texmf-dist/tex/generic/pgfplots/util/pgfplotsutil.code.tex
@@ -412,6 +412,28 @@
 	\pgfplotslistcopy\pgfplots@glob@TMPa\to#2\relax
 }
 
+% converts a foreach list to a fully expanded comma-separated list and assigns the result to \pgfplotsretval
+\long\def\pgfplots@foreach@to@list@macro#1{%
+	\begingroup
+	\pgfplotsapplistXnewempty\pgfplots@tmpa
+	\pgfplots@loc@tmptrue
+	\pgfkeys{/pgf/fpu}%
+	%
+	\pgfplotsforeachungrouped \pgfplots@i in {#1} {%
+		\ifpgfplots@loc@tmp
+			\pgfplots@loc@tmpfalse
+		\else
+			\pgfplotsapplistXpushback,\to\pgfplots@tmpa
+		\fi
+		\pgfmathfloattosci\pgfplots@i
+		\expandafter\pgfplotsapplistXpushback\pgfmathresult\to\pgfplots@tmpa
+	}%
+	\pgfplotsapplistXlet\pgfplots@loc@TMPa=\pgfplots@tmpa
+	\global\let\pgfplots@glob@TMPa=\pgfplots@loc@TMPa
+	\endgroup
+	\let\pgfplotsretval=\pgfplots@glob@TMPa
+}
+
 
 % Removes duplicates in a comma separated list and creates a new list
 % into the macro \pgfplotsretval. The list doesn't need to be sorted,