% \iffalse
% -------------------------------------------------------------------
%
% Copyright 2002--2004, Daniel H. Luecking
%
% Splines.mp may be distributed and/or modified under the conditions of the
% LaTeX Project Public License, either version 1.3 of this license or (at
% your option) any later version. The latest version of this license is in
%    <http://www.latex-project.org/lppl.txt>
% and version 1.3 or later is part of all distributions of LaTeX version
% 2003/12/01 or later.
%
% Splines has maintenance status "author-maintained". The Current Maintainer
% is Daniel H. Luecking. The Base Interpreter is MetaPost (or Metafont).
%
%<*driver>
\ProvidesFile{splines.dtx}
  [2005/02/05 v0.2 Metafont/post macros to compute splines.]%
\documentclass[draft]{ltxdoc}
\usepackage{docmfp}

\addtolength{\textwidth}{.5878pt}

\def\mytt{\upshape\mdseries\ttfamily}
\renewcommand\marg[1]{{\mytt \{#1\}}}
\renewcommand\oarg[1]{{\mytt  [#1]}}
\renewcommand\parg[1]{{\mytt  (#1)}}
\renewcommand{\meta}[1]{{$\langle$\rmfamily\itshape#1$\rangle$}}
\DeclareRobustCommand\cs[1]{{\mytt\char`\\#1}}
\def\prog#1{{\mdseries\scshape #1}}




\def\MF{\prog{meta\-font}}
\def\MP{\prog{meta\-post}}

\def\CMF{\prog{Meta\-font}}
\def\CMP{\prog{Meta\-post}}
\def\opt#1{{\sffamily\upshape#1}}
\def\mfc#1{{\mytt#1}}
\let\env\mfc
\let\file\mfc
\let\gbc\mfc
\renewcommand\{{{\mytt\char`\{}}
\renewcommand\}{{\mytt\char`\}}}
\renewcommand\|{${}\mathrel{|}{}$}

\makeatletter
\newcommand\bsl{{\mytt\@backslashchar}}
% Stupid lists!
\def\@listi{\leftmargin\leftmargini
  \parsep \z@ \@plus\p@ \@minus\z@
  \topsep 4\p@ \@plus\p@ \@minus2\p@
  \itemsep\parsep}
\let\@listI\@listi \@listi
\renewcommand\labelitemi{\normalfont\bfseries \textendash}
\renewcommand\labelitemii{\textasteriskcentered}
\renewcommand\labelitemiii{\textperiodcentered}
\leftmargini\parindent
% Stupid index!
\def\usage#1{\textrm{#1}}
\def\index@prologue{\section*{Index}\markboth{Index}{Index}Numbers
refer to the page where the corresponding entry is described.}
\def\IndexParms{%
  \parindent \z@ \columnsep 15pt
  \parskip 0pt plus 1pt
  \rightskip 5pt plus2em \mathsurround \z@
  \parfillskip=-5pt \small
  % less hanging:
  \def\@idxitem{\par\hangindent 20pt}%
  \def\subitem{\@idxitem\hspace*{15pt}}%
  \def\subsubitem{\@idxitem\hspace*{25pt}}%
  \def\indexspace{\par\vspace{10pt plus 2pt minus 3pt}}}
\renewcommand\routinestring{}
\renewcommand\variablestring{\space(var.)}
% Why does every command have to be indexed twice?
\renewcommand\SpecialMfpIndex[3]{\@bsphack
  \index{%
    \string#1\actualchar
    \string\verb\quotechar*\verbatimchar\string#1\verbatimchar
    #2 \encapchar usage}%
  \@esphack}
\def\close@crossref{\SpecialEscapechar{:}}
\makeatother
\def\VariableIndex#1{\SpecialMfpIndex{#1}{\variablestring}{}}
\def\RoutineIndex #1{\SpecialMfpIndex{#1}{}{}}

\title{Macros to Compute Splines\thanks{This file has version number
        \fileversion, last revised \filedate.}}
\author{Dan Luecking}
\date{\filedate}
\SpecialEscapechar{:}
\def\bslash{:}
\DisableCrossrefs
\CodelineIndex
\AlsoImplementation

\begin{document}
  \DeleteShortVerb{\|}
  \DocInput{splines.dtx}
\end{document}
%</driver>
%\fi
%
% \CheckSum{25}
% \CharacterTable
%  {Upper-case    \A\B\C\D\E\F\G\H\I\J\K\L\M\N\O\P\Q\R\S\T\U\V\W\X\Y\Z
%   Lower-case    \a\b\c\d\e\f\g\h\i\j\k\l\m\n\o\p\q\r\s\t\u\v\w\x\y\z
%   Digits        \0\1\2\3\4\5\6\7\8\9
%   Exclamation   \!     Double quote  \"     Hash (number) \#
%   Dollar        \$     Percent       \%     Ampersand     \&
%   Acute accent  \'     Left paren    \(     Right paren   \)
%   Asterisk      \*     Plus          \+     Comma         \,
%   Minus         \-     Point         \.     Solidus       \/
%   Colon         \:     Semicolon     \;     Less than     \<
%   Equals        \=     Greater than  \>     Question mark \?
%   Commercial at \@     Left bracket  \[     Backslash     \\
%   Right bracket \]     Circumflex    \^     Underscore    \_
%   Grave accent  \`     Left brace    \{     Vertical bar  \|
%   Right brace   \}     Tilde         \~}
%
% \catcode`\_=12
% \GetFileInfo{splines.dtx}
% \maketitle
%
% \StopEventually{\PrintIndex}
% A cubic spline through a set of points is a curve obtained by joining
% each point to the next with a cubic parametrized curve, where adjoining
% cubics must have matching first and second derivative at their common
% point.
%
% It is possible for \MP{} to compute the necessary controls.
% Unfortunately, the controls are not uniquely determined unless the curve
% is required to be closed. For open curves, there is need for two
% additional equations at the end points. A `relaxed spline' is produced
% if we require that the second derivative is $(0,0)$ at those points.
% For a closed curve, the equality of the first and second derivatives at
% the common beginning/ending point gives the needed additional equations.
%
% Note that this equates \emph{time} derivatives, so this works best when
% points are relatively evenly spaced and so the speed is relatively
% uniform. If points are differently spaced then the relatively slower
% speed between closely spaced points allows sharper turns without large
% second derivatives. Curves produced tend to have a more natural look,
% and relaxed splines are most suitable for smoothing data that is
% obtained by taking observations at evenly spaced times. Still, the
% technique is somewhat unstable when points are closely spaced, for
% example when a small change in the position of one point can produce a
% large change in its direction when viewed from another point.
%
% Start with version control information.
%    \begin{macrocode}
%<*package>
if known splines_fileversion: endinput fi;
string splines_fileversion;
splines_fileversion := "2006/09/25, v0.2a";
message "Loading splines.mp " & splines_fileversion;

%    \end{macrocode}
%
% \DescribeRoutine{list_to_array}
% Now for a command that takes a variable name \gbc{(suffix arr)} and
% copies a list of pairs to \mfc{arr1}, \mfc{arr2}, etc..
% The suffix must be declared by the calling program so that \mfc{arr}
% is numeric but \mfc{arr[\,]} are pairs, for example by ``\mfc{save arr; pair
% arr[\,];}''.
%    \begin{macrocode}
def list_to_array (suffix arr) (text list) =
  arr := 0;
  for _itm = list :
    arr[incr arr] := _itm;
  endfor
enddef;

%    \end{macrocode}
%
% \DescribeRoutine{compute_spline}
% In this command we generate the equations common to all cubic
% splines: the equality of derivatives at all interior points.
% This command accepts three suffixes: \gbc{points}, \gbc{pr}, and
% \gbc{po} which should represent previously declared arrays of pairs.
% \gbc{points} is the array of points to be connected, and must be known.
% Arrays \gbc{pr} and \gbc{po} \emph{must} be unknown and will hold the
% computed control points. See the code of \mfc{mkrelaxedspline} for
% an example of how this was arranged for the variables \mfc{rs_pr} and
% \mfc{rs_po}.
%
% Contrary to the previous (unreleased) version, \mfc{compute_spline}
% takes a boolean argument and appends the same equations at the first (=
% last) point.
%
% \DescribeRoutine{mkrelaxedspline}
% The first of these macros appends the necessary additional equations to
% get zero second derivatives at the ends.
% \DescribeRoutine{mkclosedspline}
% The second simply calls \mfc{compute_spline} with the boolean set to
% true. Both return the computed path. In theory the knowledgeable user
% can call \gbc{compute_spline (false)}, append a choice of equations for
% the ends, and then call \gbc{mksplinepath (false)}.
%
% \DescribeRoutine{dospline}
% This version accepts a list of pairs and produces a spline through
% them. It simply stores the list in an array and calls the appropriate
% version that operates on an array.
%    \begin{macrocode}
def compute_spline (expr closed) (suffix points, pr, po) =
  % interior equations:
  for j= 2 upto points - 1 :
    % equate first derivatives:
    po[j] + pr[j] = 2 points[j];
    %   and second derivatives:
    pr[j+1] + 2 pr[j] = 2 po[j] + po[j-1];
  endfor
  % for a closed curve, the first and last are also interior:
  if closed:
   po 1 + pr 1 = 2 points 1;
   po[points] + pr[points] = 2 points[points];
   pr 2 + 2 pr 1 = 2 po 1 + po[points];
   pr 1 + 2 pr[points] = 2 po[points] + po[points-1];
  fi
enddef;

vardef mksplinepath (expr closed) (suffix points, pr, po) =
  points1..controls po1 and
  for j = 2 upto points if not closed: -1 fi:
    pr[j]..points[j]..controls po[j] and
  endfor
    if closed: pr 1..cycle else: pr[points]..points[points] fi
enddef;

vardef mkrelaxedspline (suffix pnts) =
  save rs_pr, rs_po;
  pair rs_po[], rs_pr[];
  % Equate second derivative to zero at both end points
  rs_pr 2 + pnts 1 = 2 rs_po 1 ;
  pnts[pnts] + rs_po[pnts-1] = 2 rs_pr[pnts];
  compute_spline (false) (pnts, rs_pr, rs_po);
  mksplinepath   (false) (pnts, rs_pr, rs_po)
enddef;

vardef mkclosedspline (suffix pnts) =
  save cs_pr, cs_po;
  pair cs_pr[], cs_po[];
  compute_spline (true) (pnts, cs_pr, cs_po);
  mksplinepath   (true) (pnts, cs_pr, cs_po)
enddef;

vardef dospline (expr closed) (text the_list) =
  save _sp; pair _sp[];
  list_to_array (_sp) (the_list);
  if closed :
    mkclosedspline (_sp)
  else:
    mkrelaxedspline (_sp)
  fi
enddef;

%    \end{macrocode}
% The above computations produce a $2$-dimensional spline. A $1$-dimensional
% cubic spline would be a function $f(t)$ with numeric values rather
% than pair values. Such are often used to interpolate functions. That is,
% given pairs $(x\sb j,y\sb{j})$, and assuming they lie on the graph of
% some function (generally unknown), fill in the graph with $y = f(x)$
% where $f$ is a cubic function of $x$ in each interval $x\sb j < x < x\sb
% {j+1}$, making sure that the resulting graph is as smooth as possible at
% the points $x\sb j$.
%
% The requirements on our $2$-dimensional path are the following:
% \begin{enumerate}
%   \item The $j$th link should connect $(x\sb{j},y\sb{j})$ to $(x\sb{j+1},
%       y\sb{j+1})$.
%   \item The $x$-part of that link should increase linearly from $x\sb{j}$ to
%     $x\sb{j+1}$ as $t$ goes from $0$ to $1$.
%   \item The $y$-part should be a cubic $y = f(x)$.
%   \item The $x$-derivatives $df/dx$ and $d^2f/dx^2$ should match at the
% connecting points.
% \end{enumerate}
%
% Two necessary equations for converting between $x$ and $t$ coordinates
% are:
% \begin{equation}\label{first}
%         x = x\sb{j} + t \Delta x\sb{j}
% \end{equation}
% (where $\Delta x\sb{j} = x\sb{j+1} - x\sb{j}$) and
% \begin{equation}\label{second}
%       \frac{df}{dt} = \frac{dx}{dt}\frac{df}{dx} = (x\sb{j+1} -
%       x\sb{j}) \frac{df}{dx}.
% \end{equation}
% Thus we want to choose controls so that (\ref{first}) is maintained and
% so that $x$-derivatives match. It turns out that this requires controls
% at
% \begin{equation}
% \begin{array}{c}
%   (x\sb{j}, y\sb{j}) - (\Delta x\sb{j-1}, s\sb{j} \Delta x\sb{j-1})/3\\
%   (x\sb{j}, y\sb{j}) + (\Delta x\sb{j}  , s\sb{j} \Delta x\sb{j}  )/3
% \end{array}
% \end{equation}
% where $s\sb{j}$ is the slope (derivative) at $x\sb{j}$. These control
% points will produce matching slopes regardless of the values
% chosen for the $s\sb j$. To get matching second derivatives we need the
% same conditions as in parametric splines. But those equations simplify
% to the form:
% \begin{displaymath}
%  s\sb{j+1} dx\sb{j} - 2s\sb{j} (dx\sb{j} + dx\sb{j-1}) +
%       s\sb{j-1}dx\sb{j-1}
%           =  3y\sb{j+1} - 3y\sb{j-1}.
% \end{displaymath}
% As with 2-D splines there can be almost any equations at the end points.
% For a relaxed spline we equate the second derivatives to 0. To get a
% periodic function, we equate the slope and second derivative at
% beginning to those at the end. This makes it possible to put a shifted
% copy of the graph with starting point at the end of the original and
% have the same smoothness at that connection as at the other points.
%
% \DescribeRoutine{compute_fcnspline}
% This issues the equation for the slopes (array \gbc{sl} of
% \emph{unknown} numerics). The array \gbc{points} contains the $(x,y)$
% values and \gbc{dx} is a temporary numeric array which will be
% overwritten if known.
%
% \DescribeRoutine{mkfcnsplinepath}
% This simply assembles the path from the information computed by the
% above equations (and the extra equations given in the calling command).
%
% \DescribeRoutine{mkrelaxedfcnspline}
% This sets up arrays for the  \gbc{dx} and \gbc{sl} parameters of
% \mfc{compute_fcnspline}, emit the necessary endpoint equations (zero
% second derivatives) and calls the previous two routines.
%
% \DescribeRoutine{mkperiodicfcnspline}
% This does the same as the previous command, but the endpoint equations
% make the first and second derivatives at the start equal to those at
% the end.
%
% \DescribeRoutine{fcnspline}
% Finally, this command copies a list of pairs into an array and calls the
% appropriate command to process them.
%    \begin{macrocode}
def compute_fcnspline (suffix points, dx, sl) =
  % Get delta_x:
  for j = 1 upto points - 1: dx[j] := xpart (points[j+1]-points[j]);
  endfor
  for j=2 upto points - 1:
    sl[j + 1] * dx[j] + 2sl[j]*(dx[j] + dx[j-1]) + sl[j-1]*dx[j-1]
        = 3*ypart(points[j+1] - points[j-1]);
  endfor
enddef;

vardef mkfcnsplinepath (suffix points, dx, sl) =
  points1..controls (points1 + (1, sl1)*dx1/3) and
    for j = 2 upto points - 1:
      (points[j] - (1, sl[j])*dx[j-1]/3) ..points[j]..
        controls (points[j] + (1,sl[j])*dx[j]/3) and
    endfor
    (points[points] - (1,sl[points])*dx[points-1]/3)..points[points]
enddef;

vardef mkperiodicfcnspline (suffix pnts) =
  save _sl, _dx; numeric _dx[], _sl[];
  compute_fcnspline (pnts, _dx, _sl);
  % periodicity equations:
  _sl 1 = _sl[pnts];
  _sl 2 * _dx 1 + 2 _sl 1 * _dx 1 + 2 _sl[pnts] * _dx[pnts-1]
      + _sl[pnts-1] * _dx[pnts-1]
          = 3 * ypart(pnts[2] - pnts[pnts-1]);
  mkfcnsplinepath  (pnts, _dx, _sl)
enddef;

vardef mkrelaxedfcnspline (suffix pnts) =
  save _sl, _dx; numeric _dx[], _sl[];
  compute_fcnspline (pnts, _dx, _sl);
  % relaxation equations.
  _sl 2 * _dx 1 + 2 _sl1 * _dx 1 = 3 * ypart(pnts2 - pnts1);
  _sl[pnts-1] * _dx[pnts-1] + 2 _sl[pnts] * _dx[pnts-1]
        = 3 * ypart(pnts[pnts] - pnts[pnts-1]);
  mkfcnsplinepath  (pnts, _dx, _sl)
enddef;

vardef fcnspline (expr periodic) (text the_list) =
  save _fs; pair _fs[];
  list_to_array (_fs) (the_list);
  if periodic:
    mkperiodicfcnspline (_fs)
  else:
    mkrelaxedfcnspline (_fs)
  fi
enddef;

%</package>
%    \end{macrocode}
% \clearpage
%\Finale