% numodel-plot.dtx
%
% Docstrip source for numodel-plot.sty.
% Run
%   tex numodel-plot.ins
% to extract the derived file.  User-facing documentation lives
% in numodel-plot-manual.tex (a stand-alone LaTeX file).
%
% Copyright (C) 2026 Paul Zuurbier <mail@paulzuurbier.nl>
%
% This work may be distributed and/or modified under the conditions
% of the LaTeX Project Public License, either version 1.3c of this
% license or (at your option) any later version.  The latest version
% of this license is in https://www.latex-project.org/lppl.txt
%
% This work has the LPPL maintenance status 'maintained'.
% The Current Maintainer of this work is Paul Zuurbier.
%
% This work consists of the files numodel-plot.dtx and numodel-plot.ins,
% and the derived file numodel-plot.sty.
% \section{Implementation}
%
% All internal helper macros use the \texttt{\textbackslash nmp@\ldots}
% prefix (standard LaTeX internal convention).  Because the package
% may be loaded through |\input| as well as |\usepackage|, an
% explicit |\makeatletter| wraps the body so that |@| is a letter
% regardless of caller.
%
%    \begin{macrocode}
\NeedsTeXFormat{LaTeX2e}
\ProvidesPackage{numodel-plot}[2026/09/07 v0.9.1 Auto-sizing PGFPlots engine]

\makeatletter

\RequirePackage{expl3}
\RequirePackage{xparse}
\RequirePackage{l3keys2e}
% \qtyPlain drives siunitx through its code-level interface.  Every
% function it uses has been public since siunitx v3.0.0 (2021-05-17),
% but the axis labels also rely on `print-unity-mantissa = false',
% whose handling of a bare 1 was corrected in v3.3.8, and on the
% empty-exponent case of `prefix-mode = combine-exponent', fixed in
% v3.3.7.  Hence the v3.3.8 date: with anything older LaTeX warns
% rather than letting the labels come out subtly wrong.
\RequirePackage{siunitx}[2023-11-06]
\RequirePackage{pgfplots}
\RequirePackage{pdfrender}
\usepgfplotslibrary{fillbetween}

\pgfplotsset{
  compat=1.18,
  numodel/grid/.style={
    grid=both,
    grid style={
      line width=0.5pt,
      draw=black,
      line cap=round,
      dash pattern=on 0pt off 1mm
      },
    major grid style={
      line width=0.5pt,
      draw=black,
      line cap=round,
      dash pattern=on 0pt off 1mm}
      },
  numodel/ticks/.style={
    tick style={black}
  },
  % Semi-transparent white halo behind tick labels.  Only applied to an
  % axis that \calcplotdims draws *through the middle* of the plot, where
  % curves pass behind the numbers; edge-placed axes leave their labels on
  % the white margin and get no halo (see \calcplotdims).
  %
  % The halo follows the glyph outlines instead of boxing the label in a
  % rectangle: the label text is wrapped in \nmp@ticklabelhalo (defined
  % below), which emits it twice -- first stroke-only with a fat white
  % pen (the halo), then normally on top of it.  \axisdefaultticklabel
  % is pgfplots' stock label body, so number formatting and tick
  % scaling stay exactly as without the halo.
  %
  % The halo has to cover the curves as well as the grid, so the haloed
  % labels are lifted onto the `axis descriptions' layer, which sits
  % above `main' (the curves) while the grid stays on `axis grid' below
  % them.  This only takes effect when layered graphics is active;
  % \calcplotdims switches that on (`set layers') in the same branches
  % that apply these styles.
  numodel/xticklabel halo/.style={
    xticklabel={\nmp@ticklabelhalo{\axisdefaultticklabel}},
    xticklabel style={/pgfplots/on layer=axis descriptions}},
  numodel/yticklabel halo/.style={
    yticklabel={\nmp@ticklabelhalo{\axisdefaultticklabel}},
    yticklabel style={/pgfplots/on layer=axis descriptions}},
  numodel/axis/.style={
    xticklabel style={/pgf/number format/.cd, fixed, precision=6, fixed zerofill=false, use comma, 1000 sep={\,}, /tikz/.cd, text depth=0.3ex, inner sep=1pt},
    yticklabel style={/pgf/number format/.cd, fixed, precision=6, fixed zerofill=false, use comma, 1000 sep={\,}, /tikz/.cd, text depth=0.3ex, inner sep=1pt},
    hide obscured x ticks=false,
    hide obscured y ticks=false,
    axis lines=left,
    xlabel near ticks,
    ylabel near ticks,
    ylabel style={rotate=-90},
    axis line style={-}
  }
}

% Glyph-following halo used by the `numodel/...ticklabel halo' styles.
% The label is emitted twice: first with pdfrender's stroke-only text
% rendering mode -- a fat pen traces the glyph outlines, at 80%
% opacity so grid and curves still shimmer through in the halo -- and
% then the untouched original on top.  The \rlap keeps the halo copy
% out of the width, so tick-label spacing stays identical to the
% halo-less labels on edge-placed axes.
%
% The halo colour is the named colour `numodelhalo'.  It is governed
% by the `halo-color' setup key: an explicit colour pins it, the
% default `auto' re-resolves it at every \drawplot against the
% surrounding background.  LaTeX has no universal "what is behind me"
% query, so `auto' probes the places where common environments
% register their background, cascading from general to specific --
% each later probe overrides the earlier ones when it applies:
%
%   white         -- base fallback (LaTeX's plain \pagecolor leaves
%                    no readable trace);
%   \thepagecolor -- the `pagecolor' package's record of \pagecolor;
%   bg            -- beamer's background colour of the current
%                    colour theme;
%   tcbcolback    -- the current tcolorbox's background, but only
%                    when actually inside a box: tcolorbox keeps the
%                    (documented) counter `tcblayer' at the current
%                    box nesting depth, and tcbcolback itself is
%                    already defined at load time, so definedness
%                    alone would misfire outside boxes.
\colorlet{numodelhalo}{white}
\newif\ifnmp@haloauto
\nmp@haloautotrue
\newcommand{\nmp@resolvehalocolor}{%
  \ifnmp@haloauto
    \colorlet{numodelhalo}{white}%
    \ifdefined\thepagecolor
      \colorlet{numodelhalo}{\thepagecolor}%
    \fi
    \ifcsname\string\color@bg\endcsname
      \colorlet{numodelhalo}{bg}%
    \fi
    \ifcsname c@tcblayer\endcsname
      \ifnum\value{tcblayer}>0
        \colorlet{numodelhalo}{tcbcolback}%
      \fi
    \fi
  \fi
}
\newcommand{\nmp@ticklabelhalo}[1]{%
  \rlap{%
    \pgfsetstrokeopacity{0.8}%
    \textpdfrender{TextRenderingMode=Stroke, LineWidth=1.6pt,%
      LineJoinStyle=Round, LineCapStyle=Round, StrokeColor=numodelhalo}{#1}%
    \pgfsetstrokeopacity{1}%
  }%
  #1%
}

\ExplSyntaxOn

\tl_new:N \g__numodelplot_axislabel_tl
\tl_gset:Nn \g__numodelplot_axislabel_tl {ieee}

\providecommand{\xcmmax}{12}
\providecommand{\ycmmax}{10}

\keys_define:nn { numodel-plot }
  {
    axis-label-format .choice:,
    axis-label-format / iso       .code:n =
      \tl_gset:Nn \g__numodelplot_axislabel_tl {iso},
    axis-label-format / ieee      .code:n =
      \tl_gset:Nn \g__numodelplot_axislabel_tl {ieee},
    axis-label-format / brackets  .code:n =
      \tl_gset:Nn \g__numodelplot_axislabel_tl {brackets},
    axis-label-format / qty-only  .code:n =
      \tl_gset:Nn \g__numodelplot_axislabel_tl {qty-only},
    axis-label-format / unit-only .code:n =
      \tl_gset:Nn \g__numodelplot_axislabel_tl {unit-only},

    grid .choice:,
    grid / mm-dots .code:n =
      \pgfplotsset{ numodel/grid/.style={
        grid=both,
        grid~style={line~width=0.5pt, draw=black, line~cap=round,
          dash~pattern=on~0pt~off~1mm},
        major~grid~style={line~width=0.5pt, draw=black, line~cap=round,
          dash~pattern=on~0pt~off~1mm}
      } },
    grid / none .code:n =
      \pgfplotsset{ numodel/grid/.style={grid=none} },
    grid / unknown .code:n =
      \pgfplotsset{ numodel/grid/.style={#1} },

    xcmmax .code:n = \def\xcmmax{#1},
    ycmmax .code:n = \def\ycmmax{#1},

    % Colour of the tick-label halo.  `auto' (the default) matches the
    % surrounding background at each \drawplot (tcolorbox, beamer,
    % pagecolor -- see \nmp@resolvehalocolor); any other value is a
    % colour expression that pins the halo to that colour.
    halo-color .code:n =
      \str_if_eq:nnTF {#1} {auto}
        { \nmp@haloautotrue }
        { \nmp@haloautofalse
          \colorlet{numodelhalo}{#1} },

    % How a scaled axis presents its power of ten in the axis label:
    % `prefix' (the default) folds it into an SI prefix on the unit --
    % t (ks) instead of t (10^3 s) -- falling back to the exponent
    % form whenever no engineering prefix fits (see
    % \__numodelplot_prefix_check:nn); `exponent' always shows the
    % power of ten.
    scale-format .choice:,
    scale-format / prefix .code:n =
      \tl_gset:Nn \g__numodelplot_scaleformat_tl {prefix},
    scale-format / exponent .code:n =
      \tl_gset:Nn \g__numodelplot_scaleformat_tl {exponent},
  }

\NewDocumentCommand{\numodelplotsetup}{m}
  { \keys_set:nn { numodel-plot } {#1} }

\ExplSyntaxOff

\ExplSyntaxOn

% \qtyPlain[<options>]{<number>}{<unit>}
%
% Like siunitx's \qty, except that a mantissa of exactly 1 with no
% exponent left over prints as the bare unit ("m", not "1 m") -- the
% form an axis label wants when the scale factor has been absorbed
% into the unit's own prefix.
%
% Up to v0.8.0 this was done by overriding two siunitx internals
% (\__siunitx_number_output_integer:nnn and
% \__siunitx_quantity_print_unit:n) and reading four internal
% variables.  The investigation behind siunitx issue #864 established
% that siunitx's documented code-level interface already exposes the
% whole pipeline, so nothing below reaches into siunitx internals.
%
% The pipeline mirrors \__siunitx_quantity_parsed:nn: parse the
% number, run the stage selected by prefix-mode, then print.  All
% three prefix-mode values are handled, so \qtyPlain remains a
% drop-in replacement for \qty rather than being tied to
% extract-exponent.
\tl_new:N \l__nmp_qty_number_tl
\tl_new:N \l__nmp_qty_unit_tl
\tl_new:N \l__nmp_qty_input_tl
\tl_new:N \l__nmp_qty_rawunit_tl
\fp_new:N \l__nmp_qty_exponent_fp

% prefix-mode = extract-exponent: strip the SI prefixes off the unit
% and fold their net power of ten into the number, then process
% (rounding, exponent-mode, ...).  The order matters: prefixes come
% off before the exponent is adjusted, processing runs last.
\cs_new_protected:cpn { __nmp_qty_prefix_extract-exponent:n } #1
  {
    \siunitx_unit_format_extract_prefixes:nNN {#1}
      \l__nmp_qty_unit_tl \l__nmp_qty_exponent_fp
    \tl_set:Ne \l__nmp_qty_number_tl
      {
        \siunitx_number_adjust_exponent:Nn
          \l__nmp_qty_number_tl \l__nmp_qty_exponent_fp
      }
    \siunitx_number_process:NN \l__nmp_qty_number_tl \l__nmp_qty_number_tl
  }

% prefix-mode = combine-exponent: process first, then lift the
% exponent out of the number and hand it to the unit formatter, which
% turns it into an SI prefix.
\cs_new_protected:cpn { __nmp_qty_prefix_combine-exponent:n } #1
  {
    \siunitx_number_process:NN \l__nmp_qty_number_tl \l__nmp_qty_number_tl
    \exp_after:wN \__nmp_qty_drop_exponent:nnnnnnn \l__nmp_qty_number_tl
    \siunitx_unit_format_combine_exponent:nnN {#1}
      \l__nmp_qty_exponent_fp \l__nmp_qty_unit_tl
  }
\cs_new_protected:Npn \__nmp_qty_drop_exponent:nnnnnnn #1#2#3#4#5#6#7
  {
    \fp_set:Nn \l__nmp_qty_exponent_fp { #6 0#7 }
    \tl_set:Nn \l__nmp_qty_number_tl { {#1} {#2} {#3} {#4} {#5} { } { } }
  }

% prefix-mode = input: prefixes stay where the user wrote them.
\cs_new_protected:Npn \__nmp_qty_prefix_input:n #1
  {
    \siunitx_number_process:NN \l__nmp_qty_number_tl \l__nmp_qty_number_tl
    \siunitx_unit_format:nN {#1} \l__nmp_qty_unit_tl
  }

% The prefix-mode currently in force, as a string.
%
% siunitx v3.5.0 renamed the variable holding it from
% \l_siunitx_quantity_prefix_mode_tl -- a documented public variable
% up to and including v3.4.14 -- to
% \l_siunitx_quantity_prefix_mode_str, and left the old name
% declared but no longer assigned.  Reading the stale one therefore
% yields an *empty* mode and no error at all, and a \cs_if_exist
% guard does not catch it either: on siunitx 3.5 and newer, v0.9.0 of
% this package silently degraded every \qtyPlain to a plain \qty, so
% a scaled axis label read `t (1 ks)' instead of `t (ks)'.
%
% Prefer the str, fall back to the tl, and leave the result empty if
% neither carries anything -- an empty mode matches no
% \__nmp_qty_prefix_...:n and so takes the \qty branch below, which is
% what an unknown mode should do anyway.
\cs_new:Npn \__nmp_qty_mode:
  {
    \bool_lazy_and:nnTF
      { \cs_if_exist_p:N \l_siunitx_quantity_prefix_mode_str }
      { ! \tl_if_empty_p:N \l_siunitx_quantity_prefix_mode_str }
      { \tl_use:N \l_siunitx_quantity_prefix_mode_str }
      {
        \cs_if_exist:NT \l_siunitx_quantity_prefix_mode_tl
          { \tl_use:N \l_siunitx_quantity_prefix_mode_tl }
      }
  }

% #1-#7 are the seven groups \siunitx_number_parse:nN produces --
% comparator, sign, integer, decimal, uncertainty, exponent sign,
% exponent -- and #8 is the unit.  A number carrying an uncertainty
% is handed straight to \qty: separate-uncertainty-units has three
% layouts that are not worth reproducing here, and such a number
% never reduces to the bare mantissa 1 that \qtyPlain exists to
% suppress.  An unrecognised or unreadable prefix-mode falls back the
% same way.
\cs_new_protected:Npn \__nmp_qty_parsed:nnnnnnnn #1#2#3#4#5#6#7#8
  {
    \bool_lazy_and:nnTF
      { \tl_if_blank_p:n {#5} }
      { \cs_if_exist_p:c { __nmp_qty_prefix_ \__nmp_qty_mode: :n } }
      {
        \use:c { __nmp_qty_prefix_ \__nmp_qty_mode: :n } {#8}
        \exp_after:wN \__nmp_qty_output:nnnnnnn \l__nmp_qty_number_tl
      }
      { \__nmp_qty_fallback: }
  }

% The one thing \qtyPlain does differently from \qty: a number that
% has boiled down to a bare 1 prints as nothing at all, leaving the
% unit to stand alone -- and with it goes the quantity-product symbol
% that would otherwise separate the two.
\cs_new_protected:Npn \__nmp_qty_output:nnnnnnn #1#2#3#4#5#6#7
  {
    \bool_lazy_all:nTF
      {
        { \tl_if_blank_p:n {#1} }           % no comparator
        { \tl_if_blank_p:n {#2} }           % no sign
        { \str_if_eq_p:nn {#3} { 1 } }      % integer part is exactly 1
        { \tl_if_blank_p:n {#4} }           % no decimal part
        { \tl_if_blank_p:n {#5} }           % no uncertainty
        { \fp_compare_p:n { #6 0#7 = 0 } }  % exponent absent or zero
      }
      { \siunitx_print_unit:V \l__nmp_qty_unit_tl }
      {
        \siunitx_quantity_print:eV
          { \siunitx_number_output:N \l__nmp_qty_number_tl }
          \l__nmp_qty_unit_tl
      }
  }

% Everything \qtyPlain does not special-case goes to \qty unchanged.
% The options are already in force in the enclosing group, so they
% need not be passed again.
\cs_new_protected:Npn \__nmp_qty_fallback:
  { \exp_args:NVV \qty \l__nmp_qty_input_tl \l__nmp_qty_rawunit_tl }

\NewDocumentCommand{\qtyPlain}{O{} m m}
  {
    \group_begin:
      \tl_set:Nn \l__nmp_qty_input_tl {#2}
      \tl_set:Nn \l__nmp_qty_rawunit_tl {#3}
      \keys_set:nn { siunitx } {#1}
      \siunitx_unit_options_apply:n {#3}
      \bool_lazy_and:nnTF
        { \l_siunitx_number_parse_bool }
        { ! \tl_if_blank_p:n {#2} }
        {
          \siunitx_number_parse:nN {#2} \l__nmp_qty_number_tl
          \tl_if_empty:NTF \l__nmp_qty_number_tl
            { \__nmp_qty_fallback: }
            {
              \exp_after:wN \__nmp_qty_parsed:nnnnnnnn
                \l__nmp_qty_number_tl {#3}
            }
        }
        { \__nmp_qty_fallback: }
    \group_end:
  }
\ExplSyntaxOff

\ExplSyntaxOn

% \pzuIfUnitNonEngTF{<unit macro>}{<true>}{<false>}
%
% True when the unit carries a power of ten that is not a multiple of
% three -- cm, dm, daN, hPa and friends.  Such a unit already encodes
% its own order of magnitude, so the label builders leave the axis
% unscaled rather than stacking a second power of ten on top of it.
%
% The power of ten is the one siunitx itself extracts, i.e. the net
% exponent of the whole unit expression.  It therefore sees through
% abbreviations (\cm), through powers (\cm\squared is 10^{-4}, still
% non-engineering, while \cm\cubed is 10^{-6}, which is not) and
% through cancellation (\cm\per\cs is 10^{0}).  The hard-coded token
% list this replaces could do none of that, and roughly half of the
% control sequences it named (\cg, \dam, \hPa, ...) are not defined
% by siunitx at all.
%
% The argument is expanded once, as it always was, so call sites pass
% \xlabelunit rather than its contents.
\tl_new:N \l__nmp_noneng_unit_tl
\fp_new:N \l__nmp_noneng_exponent_fp
\prg_new_protected_conditional:Npnn \__nmp_unit_noneng:n #1 { TF }
  {
    \group_begin:
      \siunitx_unit_options_apply:n {#1}
      \siunitx_unit_format_extract_prefixes:nNN {#1}
        \l__nmp_noneng_unit_tl \l__nmp_noneng_exponent_fp
      \fp_compare:nNnTF
        {
          \l__nmp_noneng_exponent_fp
            - 3 * floor ( \l__nmp_noneng_exponent_fp / 3 )
        }
        = \c_zero_fp
        { \group_end: \prg_return_false: }
        { \group_end: \prg_return_true: }
  }
\prg_generate_conditional_variant:Nnn \__nmp_unit_noneng:n { o } { TF }
\NewDocumentCommand{\pzuIfUnitNonEngTF}{m m m}
  { \__nmp_unit_noneng:oTF {#1} {#2} {#3} }
\ExplSyntaxOff

\ExplSyntaxOn

\cs_new_protected:Npn \__numodelplot_xlabel_emit_unscaled:
  {
    \str_case:Vn \g__numodelplot_axislabel_tl
      {
        {iso}
          { \def \xlabel { $\nmp@savedxlabelqty \, /\, ( \si{\nmp@savedxlabelunit} )$ } }
        {ieee}
          { \def \xlabel { $\nmp@savedxlabelqty \,( \si{\nmp@savedxlabelunit} )$ } }
        {brackets}
          { \def \xlabel { $\nmp@savedxlabelqty \,[ \si{\nmp@savedxlabelunit} ]$ } }
        {qty-only}
          { \def \xlabel { $\nmp@savedxlabelqty$ } }
        {unit-only}
          { \def \xlabel { $\si{\nmp@savedxlabelunit}$ } }
      }
  }

\cs_new_protected:Npn \__numodelplot_xlabel_emit_scaled:
  {
    \str_case:Vn \g__numodelplot_axislabel_tl
      {
        {iso}
          {
            \edef\xlabel{\noexpand$\noexpand\nmp@savedxlabelqty
              \noexpand\,/\noexpand\,(\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@xlabelexpuse}{\noexpand\nmp@savedxlabelunit})\noexpand$}%
          }
        {ieee}
          {
            \edef\xlabel{\noexpand$\noexpand\nmp@savedxlabelqty
              \noexpand\,(\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@xlabelexpuse}{\noexpand\nmp@savedxlabelunit})\noexpand$}%
          }
        {brackets}
          {
            \edef\xlabel{\noexpand$\noexpand\nmp@savedxlabelqty
              \noexpand\,[\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@xlabelexpuse}{\noexpand\nmp@savedxlabelunit}]\noexpand$}%
          }
        {qty-only}
          { \def \xlabel { $\nmp@savedxlabelqty$ } }
        {unit-only}
          {
            \edef\xlabel{\noexpand$\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@xlabelexpuse}{\noexpand\nmp@savedxlabelunit}\noexpand$}%
          }
      }
  }

\cs_new_protected:Npn \__numodelplot_ylabel_emit_unscaled:
  {
    \str_case:Vn \g__numodelplot_axislabel_tl
      {
        {iso}
          { \def \ylabel { $\nmp@savedylabelqty \, /\, ( \si{\nmp@savedylabelunit} )$ } }
        {ieee}
          { \def \ylabel { $\nmp@savedylabelqty \,( \si{\nmp@savedylabelunit} )$ } }
        {brackets}
          { \def \ylabel { $\nmp@savedylabelqty \,[ \si{\nmp@savedylabelunit} ]$ } }
        {qty-only}
          { \def \ylabel { $\nmp@savedylabelqty$ } }
        {unit-only}
          { \def \ylabel { $\si{\nmp@savedylabelunit}$ } }
      }
  }

\cs_new_protected:Npn \__numodelplot_ylabel_emit_scaled:
  {
    \str_case:Vn \g__numodelplot_axislabel_tl
      {
        {iso}
          {
            \edef\ylabel{\noexpand$\noexpand\nmp@savedylabelqty
              \noexpand\,/\noexpand\,(\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@ylabelexpuse}{\noexpand\nmp@savedylabelunit})\noexpand$}%
          }
        {ieee}
          {
            \edef\ylabel{\noexpand$\noexpand\nmp@savedylabelqty
              \noexpand\,(\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@ylabelexpuse}{\noexpand\nmp@savedylabelunit})\noexpand$}%
          }
        {brackets}
          {
            \edef\ylabel{\noexpand$\noexpand\nmp@savedylabelqty
              \noexpand\,[\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@ylabelexpuse}{\noexpand\nmp@savedylabelunit}]\noexpand$}%
          }
        {qty-only}
          { \def \ylabel { $\nmp@savedylabelqty$ } }
        {unit-only}
          {
            \edef\ylabel{\noexpand$\noexpand\qtyPlain
              [evaluate-expression=false, round-mode=none,
               prefix-mode=\g__numodelplot_prefixmode_tl,
               extract-mass-in-kilograms=true,
               print-zero-exponent=false,
               print-unity-mantissa=false]%
              {1e\nmp@ylabelexpuse}{\noexpand\nmp@savedylabelunit}\noexpand$}%
          }
      }
  }

\ExplSyntaxOff

\newcommand{\pznmpAppendXScaleTicks}[1]{%
  \pgfplotsset{numodel/axis/.append style={scaled x ticks=base 10:#1}}%
}
\newcommand{\pznmpAppendYScaleTicks}[1]{%
  \pgfplotsset{numodel/axis/.append style={scaled y ticks=base 10:#1}}%
}
\newcommand{\pznmpSuppressXScaleLabel}{%
  \pgfplotsset{numodel/axis/.append style={xtick scale label code/.code={}}}%
}
\newcommand{\pznmpSuppressYScaleLabel}{%
  \pgfplotsset{numodel/axis/.append style={ytick scale label code/.code={}}}%
}

\ExplSyntaxOn

\cs_new_protected:Npn \__numodelplot_suppress_xscale_label:
  {
    \str_if_eq:VnF \g__numodelplot_axislabel_tl {qty-only}
      { \pznmpSuppressXScaleLabel }
  }
\cs_new_protected:Npn \__numodelplot_suppress_yscale_label:
  {
    \str_if_eq:VnF \g__numodelplot_axislabel_tl {qty-only}
      { \pznmpSuppressYScaleLabel }
  }

\ExplSyntaxOff

\ExplSyntaxOn

% How the power of ten of a scaled axis reaches the axis label.
% \g__numodelplot_scaleformat_tl holds the user's wish (`prefix' or
% `exponent', set through the scale-format key); the label builders
% translate that per axis into \g__numodelplot_prefixmode_tl, the
% actual siunitx prefix-mode value spliced into the \qtyPlain calls of
% the emit-scaled routines: `combine-exponent' folds the exponent into
% an SI prefix on the unit, `extract-exponent' keeps the 10^k form.
\tl_new:N \g__numodelplot_scaleformat_tl
\tl_gset:Nn \g__numodelplot_scaleformat_tl {prefix}
\tl_new:N \g__numodelplot_prefixmode_tl
\tl_gset:Nn \g__numodelplot_prefixmode_tl {extract-exponent}

% Can 10^{#1} be folded into an SI prefix on unit #2?  siunitx's
% combine-exponent divides the exponent by the power of the unit's
% first atom and adds any existing prefix; if the result is not in its
% prefix table it raises an *error*, and results like 10^1 (deca)
% would be legal SI but unwanted engineering-wise.  So this check
% mirrors that computation via the siunitx parser and additionally
% requires a multiple of 3: 10^6 m^2 -> km^2 passes, 10^3 m^2 (k=1.5)
% and 10^3 m^3 (deca) fall back to the exponent form.  Everything
% lives behind \cs_if_exist guards because it reaches into siunitx
% internals; if a future siunitx renames them the check just returns
% false and the label keeps the always-safe 10^k form.
%
% This is now the only place in the package that touches siunitx
% internals, and it is the one thing siunitx issue #864 asks for that
% is genuinely unreachable today: \siunitx_unit_format_extract_-
% prefixes:nNN reports the net power of ten of the whole unit but not
% *which* prefix sat on *which* atom, which is exactly what deciding
% the fold needs.  \l__siunitx_unit_parsed_prop holds it (prefix-n,
% unit-n, power-n, ...); if siunitx ever exposes that, this routine
% becomes internals-free like the rest.
\bool_new:N \g__numodelplot_prefix_ok_bool
\tl_new:N \l__numodelplot_power_tl
\tl_new:N \l__numodelplot_prefix_tl
\tl_new:N \l__numodelplot_target_tl
\cs_new_protected:Npn \__numodelplot_prefix_check:nn #1#2
  {
    \bool_gset_false:N \g__numodelplot_prefix_ok_bool
    \bool_lazy_all:nT
      {
        { \cs_if_exist_p:N \__siunitx_unit_parse:n }
        { \cs_if_exist_p:N \__siunitx_unit_if_symbolic:nTF }
        { \cs_if_exist_p:N \l_siunitx_unit_symbolic_seq }
        { \cs_if_exist_p:N \l__siunitx_unit_parsed_prop }
        { \cs_if_exist_p:N \l__siunitx_unit_prefixes_forward_prop }
        { \cs_if_exist_p:N \l__siunitx_unit_prefixes_reverse_prop }
      }
      { \__numodelplot_prefix_check_aux:nn {#1} {#2} }
  }
\cs_new_protected:Npn \__numodelplot_prefix_check_aux:nn #1#2
  {
        \group_begin:
          % unit macros only gain their parser meaning inside a
          % formatting group; replicate what siunitx's own
          % \__siunitx_unit_format:nNN sets up before parsing (the
          % symbolic test, too, only answers correctly after this)
          \seq_map_inline:Nn \l_siunitx_unit_symbolic_seq
            { \cs_set_eq:Nc ##1 { __siunitx_unit_ \token_to_str:N ##1 :w } }
          \__siunitx_unit_if_symbolic:nTF {#2}
            { \__siunitx_unit_parse:n {#2} }
            { \prop_clear:N \l__siunitx_unit_parsed_prop }
          \prop_if_empty:NF \l__siunitx_unit_parsed_prop
            {
              \prop_get:NnNF \l__siunitx_unit_parsed_prop { power-1 }
                \l__numodelplot_power_tl
                { \tl_set:Nn \l__numodelplot_power_tl { 1 } }
              \prop_get:NnNTF \l__siunitx_unit_parsed_prop { prefix-1 }
                \l__numodelplot_prefix_tl
                {
                  \prop_get:NVNF \l__siunitx_unit_prefixes_forward_prop
                    \l__numodelplot_prefix_tl \l__numodelplot_prefix_tl
                    { \tl_set:Nn \l__numodelplot_prefix_tl { 0 } }
                }
                { \tl_set:Nn \l__numodelplot_prefix_tl { 0 } }
              \tl_set:Nx \l__numodelplot_target_tl
                {
                  \fp_eval:n
                    { #1 / \l__numodelplot_power_tl + \l__numodelplot_prefix_tl }
                }
              \fp_compare:nNnTF { \l__numodelplot_target_tl } = \c_zero_fp
                { \bool_gset_true:N \g__numodelplot_prefix_ok_bool }
                {
                  % engineering only: target must be a multiple of 3
                  \fp_compare:nNnT
                    {
                      \l__numodelplot_target_tl
                      - 3 * trunc( \l__numodelplot_target_tl / 3 )
                    } = \c_zero_fp
                    {
                      \prop_get:NVNT \l__siunitx_unit_prefixes_reverse_prop
                        \l__numodelplot_target_tl \l__numodelplot_target_tl
                        { \bool_gset_true:N \g__numodelplot_prefix_ok_bool }
                    }
                }
            }
        \group_end:
  }

% Pick the prefix-mode for one axis: combine-exponent when the user
% asked for prefixes AND the fold is feasible, extract-exponent
% otherwise.
\cs_new_protected:Npn \__numodelplot_set_prefixmode:nn #1#2
  {
    \tl_gset:Nn \g__numodelplot_prefixmode_tl {extract-exponent}
    \str_if_eq:VnT \g__numodelplot_scaleformat_tl {prefix}
      {
        \__numodelplot_prefix_check:nn {#1} {#2}
        \bool_if:NT \g__numodelplot_prefix_ok_bool
          { \tl_gset:Nn \g__numodelplot_prefixmode_tl {combine-exponent} }
      }
  }

\NewDocumentCommand{\xlabelbuild}{}{
  \let\nmp@savedxlabelqty\xlabelqty
  \let\nmp@savedxlabelunit\xlabelunit
  \pzuIfUnitNonEngTF{\xlabelunit}{
    \__numodelplot_xlabel_emit_unscaled:
    \def\nmp@xlabelexp{0}
  }{
    \edef\nmp@xmag{\fpeval{max(abs(\xmin),abs(\xmax))}}
    \edef\nmp@xlabelexp{\fpeval{
      \nmp@xmag > 0 ? 3*floor(ln(\nmp@xmag)/ln(10)/3) : 0
    }}
    \ifnum\nmp@xlabelexp=0\relax
      \__numodelplot_xlabel_emit_unscaled:
    \else
      \edef\nmp@xlabelexpuse{\nmp@xlabelexp}
      \exp_args:NnV \__numodelplot_set_prefixmode:nn
        {\nmp@xlabelexp} \nmp@savedxlabelunit
      \__numodelplot_xlabel_emit_scaled:
      \edef\nmp@xlabelexpneg{\fpeval{-\nmp@xlabelexp}}
      \pznmpAppendXScaleTicks{\nmp@xlabelexpneg}
      \__numodelplot_suppress_xscale_label:
    \fi
  }
  \global\let\xlabelqty\undefined
  \global\let\xlabelunit\undefined
}

\NewDocumentCommand{\ylabelbuild}{}{
  \let\nmp@savedylabelqty\ylabelqty
  \let\nmp@savedylabelunit\ylabelunit
  \pzuIfUnitNonEngTF{\ylabelunit}{
    \__numodelplot_ylabel_emit_unscaled:
    \def\nmp@ylabelexp{0}
  }{
    \edef\nmp@ymag{\fpeval{max(abs(\ymin),abs(\ymax))}}
    \edef\nmp@ylabelexp{\fpeval{
      \nmp@ymag > 0 ? 3*floor(ln(\nmp@ymag)/ln(10)/3) : 0
    }}
    \ifnum\nmp@ylabelexp=0\relax
      \__numodelplot_ylabel_emit_unscaled:
    \else
      \edef\nmp@ylabelexpuse{\nmp@ylabelexp}
      \exp_args:NnV \__numodelplot_set_prefixmode:nn
        {\nmp@ylabelexp} \nmp@savedylabelunit
      \__numodelplot_ylabel_emit_scaled:
      \edef\nmp@ylabelexpneg{\fpeval{-\nmp@ylabelexp}}
      \pznmpAppendYScaleTicks{\nmp@ylabelexpneg}
      \__numodelplot_suppress_yscale_label:
    \fi
  }
  \global\let\ylabelqty\undefined
  \global\let\ylabelunit\undefined
}
\ExplSyntaxOff

\newcommand{\calcplotdims}{%
\edef\nmp@xlog{\fpeval{floor(ln(\xmax-\xmin)/ln(10))}}%
\edef\nmp@ylog{\fpeval{floor(ln(\ymax-\ymin)/ln(10))}}%
\edef\nmp@xS{\fpeval{(\xmax-\xmin) / (10^\nmp@xlog)}}%
\edef\nmp@yS{\fpeval{(\ymax-\ymin) / (10^\nmp@ylog)}}%
\edef\xcm{\fpeval{%
  \nmp@xS <= \xcmmax/10 ? 10 * \nmp@xS :%
  \nmp@xS <= \xcmmax/5  ?  5 * \nmp@xS :%
  \nmp@xS <= \xcmmax/4  ?  4 * \nmp@xS :%
  \nmp@xS <= \xcmmax/2  ?  2 * \nmp@xS :%
  \nmp@xS <= \xcmmax    ?      \nmp@xS :%
                        \nmp@xS / 2%
}}%
\edef\ycm{\fpeval{%
  \nmp@yS <= \ycmmax/10 ? 10 * \nmp@yS :%
  \nmp@yS <= \ycmmax/5  ?  5 * \nmp@yS :%
  \nmp@yS <= \ycmmax/4  ?  4 * \nmp@yS :%
  \nmp@yS <= \ycmmax/2  ?  2 * \nmp@yS :%
  \nmp@yS <= \ycmmax    ?      \nmp@yS :%
                        \nmp@yS / 2%
}}%
\edef\nmp@xtickdiv{\fpeval{\nmp@xS / \xcm}}%
\edef\nmp@ytickdiv{\fpeval{\nmp@yS / \ycm}}%
\edef\xtickdistance{\fpeval{\nmp@xtickdiv * 10^\nmp@xlog}}%
\edef\ytickdistance{\fpeval{\nmp@ytickdiv * 10^\nmp@ylog}}%
\edef\xmax{\fpeval{ceil(\xmax / \xtickdistance) * \xtickdistance}}%
\edef\ymax{\fpeval{ceil(\ymax / \ytickdistance) * \ytickdistance}}%
\edef\xmin{\fpeval{floor(\xmin / \xtickdistance) * \xtickdistance}}%
\edef\ymin{\fpeval{floor(\ymin / \ytickdistance) * \ytickdistance}}%
\edef\xcm{\fpeval{(\xmax - \xmin) / \xtickdistance}}%
\edef\ycm{\fpeval{(\ymax - \ymin) / \ytickdistance}}%
\edef\nmp@yStraddles{\fpeval{(\ymin<0 && \ymax>0) ? 1 : 0}}%
\edef\nmp@yAllNonPos{\fpeval{\ymax<=0 ? 1 : 0}}%
\edef\nmp@xStraddles{\fpeval{(\xmin<0 && \xmax>0) ? 1 : 0}}%
\edef\nmp@xAllNonPos{\fpeval{\xmax<=0 ? 1 : 0}}%
% The x-axis is drawn through the middle only when the data straddles
% y=0; there the x tick labels sit on top of the curves and need the
% white halo.  Edge placements (top / default bottom) leave the
% labels on the white margin, so no halo.
\ifnum\nmp@yStraddles=1
  \pgfplotsset{numodel/axis/.append style={axis x line*=middle, set layers,
    /pgfplots/numodel/xticklabel halo}}%
\else
  \ifnum\nmp@yAllNonPos=1
    \pgfplotsset{numodel/axis/.append style={axis x line*=top}}%
  \fi
\fi
% Likewise the y-axis only gets a halo when it is drawn through
% the middle (data straddles x=0).
\ifnum\nmp@xStraddles=1
  \pgfplotsset{numodel/axis/.append style={axis y line*=middle, set layers,
    /pgfplots/numodel/yticklabel halo}}%
\else
  \ifnum\nmp@xAllNonPos=1
    \pgfplotsset{numodel/axis/.append style={axis y line*=right}}%
  \fi
\fi
\edef\nmp@axisMoved{\fpeval{(\nmp@yStraddles || \nmp@yAllNonPos || \nmp@xStraddles || \nmp@xAllNonPos) ? 1 : 0}}%
\ifnum\nmp@axisMoved=1
\pgfplotsset{numodel/axis/.append style={%
  every x tick scale label/.style={at={(xticklabel cs:1.05)}, anchor=south west},%
  every y tick scale label/.style={at={(yticklabel cs:1.05)}, anchor=south east}%
}}%
\fi
\edef\nmp@fourQuad{\fpeval{(\nmp@xStraddles && \nmp@yStraddles) ? 1 : 0}}%
\ifnum\nmp@fourQuad=1
\pgfplotsset{numodel/axis/.append style={%
  xlabel style={at={(axis cs:\xmax,0)}, anchor=west, xshift=1em},%
  ylabel style={at={(axis cs:0,\ymax)}, anchor=south, yshift=1em, rotate=0}%
}}%
\fi
\edef\nmp@qIandII{\fpeval{(\nmp@xStraddles && \ymin>=0) ? 1 : 0}}%
\ifnum\nmp@qIandII=1
\pgfplotsset{numodel/axis/.append style={%
  ylabel style={at={(axis cs:0,\ymax)}, anchor=south, yshift=1em, rotate=0}%
}}%
\fi
\edef\nmp@qIIandIII{\fpeval{(\nmp@yStraddles && \nmp@xAllNonPos) ? 1 : 0}}%
\ifnum\nmp@qIIandIII=1
\pgfplotsset{numodel/axis/.append style={%
  xlabel style={at={(axis cs:\xmin,0)}, anchor=east, xshift=-1em}%
}}%
\fi
\edef\nmp@qIIIandIV{\fpeval{(\nmp@xStraddles && \nmp@yAllNonPos) ? 1 : 0}}%
\ifnum\nmp@qIIIandIV=1
\pgfplotsset{numodel/axis/.append style={%
  ylabel style={at={(axis cs:0,\ymin)}, anchor=north, yshift=-1em, rotate=0}%
}}%
\fi
\edef\nmp@qIVandI{\fpeval{(\nmp@yStraddles && \xmin>=0) ? 1 : 0}}%
\ifnum\nmp@qIVandI=1
\pgfplotsset{numodel/axis/.append style={%
  xlabel style={at={(axis cs:\xmax,0)}, anchor=west, xshift=1em}%
}}%
\fi
\ifdefined\xlabelqty\ifdefined\xlabelunit\xlabelbuild\fi\fi
\ifdefined\ylabelqty\ifdefined\ylabelunit\ylabelbuild\fi\fi
}

\newcommand{\drawplot}[1]{
  \calcplotdims
  \begin{tikzpicture}
    % Resolve halo-color=auto against the *current* surroundings; the
    % tikzpicture group keeps the \colorlet local to this plot.
    \nmp@resolvehalocolor
    \begin{axis}[
      numodel/grid,
      numodel/ticks,
      numodel/axis,
      xlabel=\xlabel,
      ylabel=\ylabel,
      xmin=\xmin, xmax=\xmax,
      ymin=\ymin, ymax=\ymax,
      xtick distance=\xtickdistance,
      ytick distance=\ytickdistance,
      width=\xcm cm,
      height=\ycm cm,
      scale only axis
    ]
      #1
    \end{axis}
  \end{tikzpicture}
}

\makeatother
%    \end{macrocode}
%
% \Finale
\endinput