~tex-sx/tex-sx/development

% \iffalse meta-comment

%<*internal>

\def\nameofplainTeX{plain}

\ifx\fmtname\nameofplainTeX\else

  \expandafter\begingroup

\fi

%</internal>

%<*install>

\input docstrip.tex

\keepsilent

\askforoverwritefalse

\preamble

Copyright (C) 2012 by Andrew Stacey

-------------------------------------------

This file 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 2005/12/01 or later.

\endpreamble

\generate{\file{tikzlibraryhobby.code.tex} {\from{hobby.dtx}{tikzlibrary}}}

\generate{\file{pgflibraryhobby.code.tex} {\from{hobby.dtx}{pgflibrary}}}

\generate{\file{hobby.code.tex}

{\from{hobby.dtx}{hobby}}}

\generate{\file{pml3array.sty}

{\from{hobby.dtx}{array}}}

%</install>

%<install>\endbatchfile

%<*internal>

\generate{

  \file{\jobname.ins}{\from{\jobname.dtx}{install}}

}

\ifx\fmtname\nameofplainTeX

  \expandafter\endbatchfile

\else

  \expandafter\endgroup

\fi

%</internal>

%<*driver>

\documentclass{l3doc}

\usepackage[T1]{fontenc}

\usepackage{csquotes}

\usepackage{lmodern}

\usepackage{tikz}

\usepackage{amsmath}

\usetikzlibrary{hobby,decorations.pathreplacing}

\usepackage[margin=3cm]{geometry}

\EnableCrossrefs

\CodelineIndex

\RecordChanges

\begin{document}

  \DocInput{\jobname.dtx}

\end{document}

%</driver>

% \fi

%

% \title{Hobby's Algorithm in TikZ/PGF}

% \author{Andrew Stacey}

% \date{2012-05-15}

% \maketitle

%

% \tikzset{

%   show curve controls/.style={

%    decoration={

%      show path construction,

%      curveto code={

%      \draw [blue, dashed]

%          (\tikzinputsegmentfirst)    -- (\tikzinputsegmentsupporta)

%          node [at end, draw, solid, red, inner sep=2pt]{};

%        \draw [blue, dashed]

%          (\tikzinputsegmentsupportb) -- (\tikzinputsegmentlast)

%          node [at start, draw, solid, red, inner sep=2pt]{};

%      }

%    },decorate

%  },

% }

%

% \begin{figure}

% \centering

% \begin{tikzpicture}[scale=.5]

% \draw[scale=.1,postaction=show curve controls,line width=1mm,red] (0,0)

% .. controls (26.76463,-1.84543) and (51.4094,14.58441) .. (60,40)

% .. controls (67.09875,61.00188) and (59.76253,84.57518) .. (40,90)

% .. controls (25.35715,94.01947) and (10.48064,84.5022) .. (10,70)

% .. controls (9.62895,58.80421) and (18.80421,49.62895) .. (30,50);

% 

% \fill[green] (0,0) circle[radius=2pt]

%  (6,4) circle[radius=2pt]

%  (4,9) circle[radius=2pt]

%  (1,7) circle[radius=2pt]

%  (3,5) circle[radius=2pt];

% \draw[postaction=show curve controls,thick] (0,0) to[curve through={(6,4) .. (4,9) .. (1,7)}] (3,5);

% \begin{scope}[xshift=10cm]

% \draw[scale=.1,postaction=show curve controls,line width=1mm,red] (0,0)

% .. controls (5.18756,-26.8353) and (60.36073,-18.40036) .. (60,40)

% .. controls (59.87714,59.889) and (57.33896,81.64203) .. (40,90)

% .. controls (22.39987,98.48387) and (4.72404,84.46368) .. (10,70)

% .. controls (13.38637,60.7165) and (26.35591,59.1351) .. (30,50)

% .. controls (39.19409,26.95198) and (-4.10555,21.23804) .. (0,0); % 

% \fill[green] (0,0) circle[radius=2pt]

%  (6,4) circle[radius=2pt]

%  (4,9) circle[radius=2pt]

%  (1,7) circle[radius=2pt]

%  (3,5) circle[radius=2pt];

% \draw[postaction=show curve controls,thick] (0,0) to[closed,curve through={(6,4) .. (4,9) .. (1,7)}] (3,5);

% \end{scope}

% \end{tikzpicture}

% \caption{Hobby's algorithm in TikZ overlaying the output of MetaPost}

% \end{figure}

%

% \section{Usage}

%

% The package is provided in form of a TikZ library.

% It can be loaded with

% \begin{verbatim}

% \usetikzlibrary{hobby}

% \end{verbatim}

%

% The TikZ library installs a \Verb+to path+ which draws a smooth curve through the given points:

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw (0,0) to[curve through={(6,4) .. (4,9) .. (1,7)}] (3,5);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}[scale=.5]

% \draw (0,0) to[curve through={(6,4) .. (4,9) .. (1,7)}] (3,5);

% \end{tikzpicture}

% \end{center}

%

% The path can be open, as above, or closed:

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw (0,0) to[closed,curve through={(6,4) .. (4,9) .. (1,7)}] (3,5);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}[scale=.5]

% \draw (0,0) to[closed,curve through={(6,4) .. (4,9) .. (1,7)}] (3,5);

% \end{tikzpicture}

% \end{center}

% There is also the facility to subvert TikZ's path processor and define curves simply using the \Verb+..+ separator between points.

% Note that this relies on something a little special in TikZ: the syntax \Verb+(0,0) .. (2,3)+ is currently detected and processed but there is no action assigned to that syntax.

% If a later version of TikZ assigns some action to that syntax, this package should make its override optional.

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw (-3,0) -- (0,0) .. (6,4) .. (4,9) .. (1,7) .. (3,5) -- ++(2,0);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}[scale=.5]

% \draw (-3,0) -- (0,0) .. (6,4) .. (4,9) .. (1,7) .. (3,5) -- ++(2,0);

% \end{tikzpicture}

% \end{center}

%

% The algorithm can deal with open or closed paths, it is possible to vary the ``tensions'' between the specified points of the paths, and for an open path it is possible to specify the incoming and outgoing angles either directly or via certain ``curl'' parameters.

% See the Examples section for more examples.

%

% The algorithm is actually implemented in \LaTeX3 with (almost\footnote{At the moment, \LaTeX3 lacks a \Verb+atan2+ function so \Verb+PGFMath+ is used to remedy that.}) no reference to TikZ or PGF.

% The TikZ library is simply a wrapper that takes the user's input, converts it into the right format for the \LaTeX3 code, and then calls that code to generate the path.

% 

% \section{Implementing Hobby's Algorithm}

% 

% We start with a list of \(n+1\) points, \(z_0, \dotsc, z_n\).

% The base code assumes that these are already stored in two arrays\footnote{Arrays are thinly disguised property lists}: the \(x\)--coordinates in \Verb+\l_hobby_points_x_array+ and the \(y\)--coordinates in \Verb+\l_hobby_points_y_array+.

% As our arrays are \(0\)--indexed, the actual number of points is one more than this.

% For a closed curve, we have \(z_n = z_0\)\footnote{Note that there is a difference between a closed curve and an open curve whose endpoints happen to overlap}.

% For closed curves it will be convenient to add an additional point at \(z_1\): thus \(z_{n+1} = z_1\).

% This makes \(z_n\) an internal point and makes the algorithms for closed paths and open paths agree longer than they would otherwise.

% The number of apparent points is stored as \Verb+\l_hobby_npoints_int+.

% Thus for an open path, \Verb+\l_hobby_npoints_int+ is \(n\), whilst for a closed path, it is \(n+1\)\footnote{In fact, we allow for the case where the user specifies a closed path but with \(z_n \ne z_0\).

% In that case, we assume that the user meant to repeat \(z_0\).

% This adds another point to the list.}.

% Following Hobby, let us write \(n'\) for \(n\) if the path is open and \(n+1\) if closed.

%

% From this we compute the distances and angles between successive points, storing these again as arrays.

% These are \Verb+\l_hobby_distances_array+ and \Verb+\l_hobby_angles_array+.

% The term indexed by \(k\) is the distance (or angle) of the line between the \(k\)th point and the \(k+1\)th point.

% For the internal nodes\footnote{Hobby calls the specified points \emph{knots}}, we store the difference in the angles in \Verb+\l_hobby_psi_array+.

% The \(k\)th value on this is the angle subtended at the \(k\)th node.

% This is thus indexed from \(1\) to \(n'-1\).

% 

% The bulk of the work consists in setting up a linear system to compute the angles of the control points.

% At a node, say \(z_i\), we have various pieces of information:

%

% \begin{enumerate}

% \item The angle of the incoming curve, \(\phi_i\), relative to the straight line from \(z_{i-1}\) to \(z_i\)

% \item The angle of the outgoing curve, \(\theta_i\), relative to the straight line from \(z_i\) to \(z_{i+1}\)

% \item The tension of the incoming curve, \(\overline{\tau}_i\)

% \item The tension of the outgoing curve, \(\tau_i\)

% \item The speed of the incoming curve, \(\sigma_i\)

% \item The speed of the outgoing curve, \(\rho_i\)

% \end{enumerate}

%

% The tensions are known at the start.

% The speeds are computed from the angles.

% Thus the key thing to compute is the angles.

% This is done by imposing a ``mock curvature'' condition.%

% The formula for the mock curvature is:

%^^A

% \[

%   \hat{k}(\theta,\phi,\tau,\overline{\tau}) = \tau^2 \left( \frac{2(\theta + \phi)}{\overline{\tau}} - 6\theta\right)

% \]

%^^A

% and the condition that the mock curvatures have to satisfy is that at each \emph{internal} node, the curvatures must match:

% %

% \[

%   \hat{k}(\phi_i,\theta_{i-1},\overline{\tau}_i,\tau_{i-1})/d_{i-1} = \hat{k}(\theta_i,\phi_{i+1},\tau_i,\overline{\tau}_{i+1})/d_i.

% \]

%^^A

% Substituting in yields:

%^^A

% \[

%   \frac{\overline{\tau}_i^2}{d_{i-1}} \left( \frac{2(\phi_i + \theta_{i-1})}{\tau_{i-1}} - 6\phi_i\right) = \frac{\tau_i^2}{d_i} \left( \frac{2(\theta_i + \phi_{i+1})}{\overline{\tau}_{i+1}} - 6\theta_i \right).

% \]

%^^A

% Let us rearrange that to the following:

%^^A

% \begin{align*}

%   d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 &\theta_{i-1} \\

%^^A

% +

% d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (1 - 3 \tau_{i-1}) &\phi_i \\

%^^A

% -

%  d_{i-1} \tau_{i-1} \tau_i^2 (1 - 3 \overline{\tau}_{i+1}) &\theta_i \\

%^^A

% -

%  d_{i-1} \tau_{i-1} \tau_i^2 &\phi_{i+1} \\

%^^A

% =

% 0

% \end{align*}

%^^A

% For both open and closed paths this holds for \(i=1\) to \(i=n' - 1\).

% 

% We also have the condition that \(\theta_i + \phi_i = -\psi_i\) where \(\psi_i\) is the angle subtended at a node by the lines to the adjacent nodes.

% This holds for the internal nodes\footnote{Recall that by dint of repetition, all nodes are effectively internal for a closed path}.

% Therefore for \(i=1\) to \(n'-1\) the above simplifies to the following:

% %

% \begin{align*}

%   d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 &\theta_{i-1} \\

% +

% (d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (3 \tau_{i-1} - 1)

% +

%  d_{i-1} \tau_{i-1} \tau_i^2 (3 \overline{\tau}_{i+1} - 1)) &\theta_i \\

% +

%  d_{i-1} \tau_{i-1} \tau_i^2 & \theta_{i+1} \\

% =

% - d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (3 \tau_{i-1} - 1) &\psi_i \\

% - d_{i-1} \tau_{i-1} \tau_i^2& \psi_{i+1}

% \end{align*}

% 

% For an open path we have two more equations.

% One involves \(\theta_0\).

% The other is the above for \(i = n'-1 = n-1\) with additional information regarding \(\psi_n\).

% It may be that one or either of \(\theta_0\) or \(\psi_n\) is specified in advance.

% In that case, the first equation is simply setting \(\theta_0\) to that value and the last equation involves substituting the value for \(\psi_n\) into the above.

% If not, they are given by formulae involving ``curl'' parameters \(\chi_0\) and \(\chi_n\) and result in the equations:

% %

% \begin{align*}

% \theta_0 &= \frac{\tau_0^3 + \chi_0 \overline{\tau}_1^3(3 \tau_0 - 1)}{\tau_0^3(3 \overline{\tau}_1 - 1) + \chi_0 \overline{\tau}_1^3} \phi_1 \\

% \phi_n &= \frac{\overline{\tau}_n^3 + \chi_n \tau_{n-1}^3(3 \overline{\tau}_n - 1)}{\overline{\tau}_n^3(3 \tau_{n-1} - 1) + \chi_n \tau_{n-1}^3} \theta_{n-1}

% \end{align*}

%^^A

% Using \(\phi_1 = - \psi_1 - \theta_1\), the first rearranges to:

%^^A

% \[

% (\tau_0^3(3 \overline{\tau}_1 - 1) + \chi_0 \overline{\tau}_1^3) \theta_0 + (\tau_0^3 + \chi_0 \overline{\tau}_1^3(3 \tau_0 - 1)) \theta_1 = - (\tau_0^3 + \chi_0 \overline{\tau}_1^3(3 \tau_0 - 1)) \psi_1.

% \]

%^^A

% The second should be substituted in to the general equation with \(i = n-1\).

% This yields:

%^^A  

% \begin{align*}

%   d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 &\theta_{n-2} \\

% +

% (d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 (3 \tau_{n-2} - 1)

% +

%  d_{n-2} \tau_{n-2} \tau_{n-1}^2 (3 \overline{\tau}_{n} - 1) \\

% - d_{n-2} \tau_{n-2} \tau_{n-1}^2  \frac{\overline{\tau}_n^3 + \chi_n \tau_{n-1}^3(3 \overline{\tau}_n - 1)}{\overline{\tau}_n^3(3 \tau_{n-1} - 1) + \chi_n \tau_{n-1}^3}) & \theta_{n-1} \\

% =

% - d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 (3 \tau_{n-2} - 1) &\psi_{n-1}

% \end{align*}

%^^A

% This gives \(n'\) equations in \(n'\) unknowns (\(\theta_0\) to \(\theta_{n-1}\)).

% The coefficient matrix is tridiagonal.

% It is more natural to index the entries from \(0\).

% Let us write \(A_i\) for the subdiagonal, \(B_i\) for the main diagonal, and \(C_i\) for the superdiagonal.

% Let us write \(D_i\) for the target vector.

% Then for an open path we have the following formulae:

%^^A

% \begin{align*}

% A_i &= d_i \overline{\tau}_{i+1} \overline{\tau}^2_i \\

% B_0 &= \begin{cases}

% 1 & \text{if}\; \theta_0\; \text{given} \\

% \tau_0^3(3 \overline{\tau}_1 - 1) + \chi_0 \overline{\tau}^3_1 & \text{otherwise}

% \end{cases} \\

% B_i &= d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (3 \tau_{i-1} -1) + d_{i-1} \tau_{i-1} \tau_i^2(3 \overline{\tau}_{i+1} - 1) \\

% B_{n-1} &= \begin{cases} d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 (3 \tau_{n-2} - 1) + d_{n-2} \tau_{n-2} \tau_{n-1}^2(3 \overline{\tau}_{n} - 1) & \text{if}\; \phi_n\; \text{given} \\

%  d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 (3 \tau_{n-2} - 1) + d_{n-2} \tau_{n-2} \tau_{n-1}^2(3 \overline{\tau}_{n} - 1)

% \\

% - d_{n-2} \tau_{n-2} \tau_{n-1}^2  \frac{\overline{\tau}_n^3 + \chi_n \tau_{n-1}^3(3 \overline{\tau}_n - 1)}{\overline{\tau}_n^3(3 \tau_{n-1} - 1) + \chi_n \tau_{n-1}^3}) & \text{otherwise}

% \end{cases} \\

% C_0 &= \begin{cases}

% 0 & \text{if}\; \theta_0\; \text{given} \\

% \tau_0^3 + \chi_0 \overline{\tau}_1^3(3\tau_0 - 1) & \text{otherwise}

% \end{cases} \\

% C_i &= d_{i-1} \tau_{i-1} \tau_i^2 \\

% D_0 &= \begin{cases}

% \overline{\theta}_0 & \text{if}\; \theta_0\; \text{given} \\

%  - (\tau_0^3 + \chi_0 \overline{\tau}_1^3(3 \tau_0 - 1)) \psi_1 & \text{otherwise}

% \end{cases} \\

% D_i &= - d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (3 \tau_{i-1} - 1) \psi_i

% - d_{i-1} \tau_{i-1} \tau_i^2 \psi_{i+1} \\

% D_{n-1} &= \begin{cases}

% - d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 (3 \tau_{n-2} - 1) \psi_{n-1} - d_{n-2} \tau_{n-2} \tau_{n-1}^2 \overline{\phi}_n & \text{if}\; \phi_n\; \text{given} \\

% - d_{n-1} \overline{\tau}_{n} \overline{\tau}_{n-1}^2 (3 \tau_{n-2} - 1) \psi_{n-1} & \text{otherwise}

% \end{cases}

% \end{align*}

%

% For a closed path, we have \(n\) equations in \(n+2\) unknowns (\(\theta_0\) to \(\theta_{n+1}\)).

% However, we have not included all the information.

% Since we have repeated points, we need to identify \(\theta_0\) with \(\theta_n\) and \(\theta_1\) with \(\theta_{n+1}\).

% To get a system with \(n'\) equations in \(n'\) unknowns, we add the equation \(\theta_0 - \theta_n = 0\) and substitute in \(\theta_{n+1} = \theta_1\).

% The resulting matrix is not quite tridiagonal but has extra entries on the off-corners.

% However, it can be written in the form \(M + u v^\top\) with \(M\) tridiagonal.

% There is some freedom in choosing \(u\) and \(v\).

% For simplest computation, we take \(u = e_0 + e_{n'-1}\).

% This means that \(v = d_{n'-2} \tau_{n'-2} \tau_{n'-1}^2 e_1 - e_{n'-1}\).

% With the same notation as above, the matrix \(M\) is given by the following formulae:

%^^A

% \begin{align*}

% A_i &= d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 \\

%^^A

% B_0 &= 1 \\

%^^A

% B_i &= d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (3 \tau_{i-1} -1) + d_{i-1} \tau_{i-1} \tau_i^2(3 \overline{\tau}_{i+1} - 1) \\

%^^A

% B_{n'-1} &= d_{n'-1} \overline{\tau}_{n'} \overline{\tau}_{n'-1}^2 (3 \tau_{n'-2} -1) + d_{n'-2} \tau_{n'-2} \tau_{n'-1}^2(3 \overline{\tau}_{n'} - 1) + 1\\

%^^A

% C_0 &= - d_{n'-2} \tau_{n'-2} \tau_{n'-1}^2 \\

%^^A

% C_i &= d_{i-1} \tau_{i-1} \tau_i^2 \\

%^^A

% D_0 &= 0 \\

%^^A

% D_i &= - d_i \overline{\tau}_{i+1} \overline{\tau}_i^2 (3 \tau_{i-1} - 1) \psi_i

% - d_{i-1} \tau_{i-1} \tau_i^2 \psi_{i+1} \\

%^^A

% D_{n'-1} &= - d_{n'-1} \overline{\tau}_{n'} \overline{\tau}_{n'-1}^2 (3 \tau_{n'-2} - 1) \psi_{n'-1}

% - d_{n'-2} \tau_{n'-2} \tau_{n'-1}^2 \psi_1

% \end{align*}

%

% The next step in the implementation is to compute these coefficients and store them in appropriate arrays.

% Having done that, we need to solve the resulting tridiagonal system.

% This is done by looping through the arrays doing the following substitutions (starting at \(i = 1\)):

% %

% \begin{align*}

% B_i' &= B_{i-1}' B_i - A_i C_{i-1}' \\

% C_i' &= B_{i-1}' C_i \\

% D_i' &= B_{i-1}' D_i - A_i D_{i-1}'

% \end{align*}

%^^A

% followed by back-substitution:

%^^A

% \begin{align*}

% \theta_{n-1} &= D_{n-1}'/B_{n-1}' \\

% \theta_i &= (D_i' - C_i' \theta_{i+1})/B_i'

% \end{align*}

%^^A

% For a closed path, we run this both with the vector \(D\) and the vector \(u = e_0 + e_{n'-1}\).

% Then to get the real answer, we use the Sherman--{}Morrison formula:

%^^A

% \[

% (M + u v^\top)^{-1} D = M^{-1} D - \frac{M^{-1} u v^\top M^{-1} D}{1 + v^\top M^{-1} u}.

% \]

%^^A

% 

% This leaves us with the values for \(\theta_i\).

% We now substitute these into Hobby's formulae for the lengths:

% %

% \begin{align*}

% \rho_i &= \frac{2 + \alpha_i}{1 + (1 - c) \cos \theta_i + c \cos \phi_{i+1}} \\

% \sigma_{i+1} &= \frac{2 - \alpha_i}{1 + (1 - c) \cos \phi_{i+1} + c \cos \theta_i} \\

% \alpha_i &= a (\sin \theta_i - b \sin \phi_{i+1})(\sin \phi_{i+1} - b \sin \theta_i)(\cos \theta_i - \cos \phi_{i+1})

% \end{align*}

% %

% where \(a = \sqrt{2}\), \(b = 1/16\), and \(c = (3 - \sqrt{5})/2\).

% 

% Now \(\theta_i\) is the angle relative to the line from \(z_i\) to \(z_{i+1}\), so to get the true angle we need to add back that angle.

% Fortunately, we stored those angles at the start.

% So the control points are:

% %

% \begin{gather*}

%   (\rho_i \cos (\theta_i + \omega_i), \rho_i \sin (\theta_i + \omega_i)) + z_i \\

% (-\sigma_i \cos(\theta_i + \omega_i), -\sigma_i \sin(\theta_i + \omega_i)) + z_{i+1}

% \end{gather*}

% 

% \section{Examples}

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[curve through={(1,.5) .. (2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[curve through={(1,.5) .. (2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{center}

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[out angle=0,in angle=180,curve through={(1,.5) .. (2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[out angle=0,in angle=180,curve through={(1,.5) .. (2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{center}

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[curve through={(1,.5) .. ([tension in=2]2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[curve through={(1,.5) .. ([tension in=2]2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{center}

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[curve through={(1,.5) .. ([tension=2]2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) to[curve through={(1,.5) .. ([tension=2]2,0) .. (3,.5)}] (4,0);

% \end{tikzpicture}

% \end{center}

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) .. (1,.5) .. (2,0) .. (3,.5) .. (4,0);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw[postaction=show curve controls]

% (0,0) .. (1,.5) .. (2,0) .. (3,.5) .. (4,0);

% \end{tikzpicture}

% \end{center}

%

% \begin{verbatim}

% \begin{tikzpicture}[scale=.5]

% \draw (0,0) .. (6,4) .. (4,9) .. (1,7) .. (3,5) .. cycle;

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}[scale=.5]

% \draw (0,0) .. (6,4) .. (4,9) .. (1,7) .. (3,5) .. cycle;

% \end{tikzpicture}

% \end{center}

%

% \section{Edge Cases}

%

% Angles are constrained to lie in the interval \((-\pi,pi]\).

% This can introduce edge cases as there is a point where we have to compare an angle with \(-\pi\) and if it is equal, add \(2 \pi\).

% This will occur if the path ``doubles back'' on itself as in the next example.

% By nudging the repeated point slightly, the behaviour changes drastically.

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw (0,0) .. (1,0) .. (0,0) .. (0,-1);

% \draw[xshift=2cm] (0,0) .. (1,0) .. (0,0.1) .. (0,-1);

% \draw[xshift=4cm] (0,0) .. (1,0) .. (0,-0.1) .. (0,-1);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw (0,0) .. (1,0) .. (0,0) .. (0,-1);

% \draw[xshift=2cm] (0,0) .. (1,0) .. (0,0.1) .. (0,-1);

% \draw[xshift=4cm] (0,0) .. (1,0) .. (0,-0.1) .. (0,-1);

% \end{tikzpicture}

% \end{center}

%

% Due to the precision of the computations, it is not possible to always get this test correct.

% The simplest solution is to nudge the repeated point in one direction or the other.

% Experimenting shows that the ``nudge factor'' can be extremely small (note that it will be proportional to the distance between the specified points).

% It is best to nudge it in the direction most normal to the line between the specified points as the goal is to nudge the difference of the angles.

% An alternative solution is to add an additional point for the curve to go through.

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw (0,0) .. (1,0) .. (0,0) .. (0,-1);

% \draw[xshift=2cm] (0,0) .. (1,0) .. (0,0.002) .. (0,-1);

% \draw[xshift=4cm] (0,0) .. (1,0) .. (0,-0.002) .. (0,-1);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw (0,0) .. (1,0) .. (0,0) .. (0,-1);

% \draw[xshift=2cm] (0,0) .. (1,0) .. (0,0.002) .. (0,-1);

% \draw[xshift=4cm] (0,0) .. (1,0) .. (0,-0.002) .. (0,-1);

% \end{tikzpicture}

% \end{center}

%

% Lastly, it is possible to add an \Verb+excess angle+ key to a coordinate.

% This will add the corresponding multiple of \(2\pi\) to the angle difference.

%

% \begin{verbatim}

% \begin{tikzpicture}

% \draw (0,0) .. (1,0) .. (0,0) .. (0,-1);

% \draw[xshift=2cm] (0,0) .. ([excess angle=1]1,0) .. (0,0) .. (0,-1);

% \draw[xshift=4cm] (0,0) .. ([excess angle=-1]1,0) .. (0,0) .. (0,-1);

% \end{tikzpicture}

% \end{verbatim}

%

% \begin{center}

% \begin{tikzpicture}

% \draw (0,0) .. (1,0) .. (0,0) .. (0,-1);

% \draw[xshift=2cm] (0,0) .. ([excess angle=1]1,0) .. (0,0) .. (0,-1);

% \draw[xshift=4cm] (0,0) .. ([excess angle=-1]1,0) .. (0,0) .. (0,-1);

% \end{tikzpicture}

% \end{center}

%

% Although this is intended to be an integer, no check is done and so some quite odd curves can result from changing this parameter.

%

% \StopEventually{\PrintChanges}

% \section{Implementation}

%

% \subsection{Main Code}

%

% \iffalse

%<*hobby>

% \fi

%

% We use \LaTeX3 syntax so need to load the requisite packages

%    \begin{macrocode}

\RequirePackage{expl3}

\RequirePackage{xparse}

\RequirePackage{pml3array}

\ExplSyntaxOn

%    \end{macrocode}

% \subsubsection{Initialisation}

%

% We declare all our variables.

%

% Now we define our objects for use in generating the path.

%

% \begin{macro}{\l_hobby_closed_bool}

% \Verb+\l_hobby_closed_bool+ is \Verb+true+ if the path is closed.

%    \begin{macrocode}

\bool_new:N \l_hobby_closed_bool

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_disjoint_bool}

% \Verb+\l_hobby_disjoint_bool+ is \Verb+true+ if the path should start with a \Verb+moveto+ command.

%    \begin{macrocode}

\bool_new:N \l_hobby_disjoint_bool

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_points_array}

% \Verb+\l_hobby_points_array+ is an array holding the specified points on the path.

% In the \LaTeX3 code, a ``point'' is a token list of the form \Verb+x = <number>, y = <number>+.

% This gives us the greatest flexibility in passing points back and forth between the \LaTeX3 code and any calling code.

% The array is indexed by integers beginning with \(0\).

% In the documentation, we will use the notation \(z_k\) to refer to the \(k\)th point.

%    \begin{macrocode}

\array_new:N \l_hobby_points_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_points_x_array}

% \Verb+\l_hobby_points_x_array+ is an array holding the \(x\)--{}coordinates of the specified points.

%    \begin{macrocode}

\array_new:N \l_hobby_points_x_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_points_y_array}

% \Verb+\l_hobby_points_y_array+ is an array holding the \(y\)--{}coordinates of the specified points.

%    \begin{macrocode}

\array_new:N \l_hobby_points_y_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_angles_array}

% \Verb+\l_hobby_angles_array+ is an array holding the angles of the lines between the points.

% Specifically, the angle indexed by \(k\) is the angle in radians of the line from \(z_k\) to \(z_{k+1}\).

%    \begin{macrocode}

\array_new:N \l_hobby_angles_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_distances_array}

% \Verb+\l_hobby_distances_array+ is an array holding the distances between the points.

% Specifically, the distance indexed by \(k\), which we will write as \(d_k\), is the length of the line from \(z_k\) to \(z_{k+1}\).

%    \begin{macrocode}

\array_new:N \l_hobby_distances_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tension_out_array}

% \Verb+\l_hobby_tension_out_array+ is an array holding the tension for the path as it leaves each point.

% This is a parameter that controls how much the curve ``flexes'' as it leaves the point.

% In the following, this will be written \(\tau_k\).

%    \begin{macrocode}

\array_new:N \l_hobby_tension_out_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tension_in_array}

% \Verb+\l_hobby_tension_in_array+ is an array holding the tension for the path as it arrives at each point.

% This is a parameter that controls how much the curve ``flexes'' as it gets to the point.

% In the following, this will be written \(\overline{\tau}_k\).

%    \begin{macrocode}

\array_new:N \l_hobby_tension_in_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tria_array}

% \Verb+\l_hobby_tria_array+ is an array holding the subdiagonal of the linear system that has to be solved to find the angles of the control points.

% In the following, this will be denoted by \(A_i\).

% The first index is \(1\).

%    \begin{macrocode}

\array_new:N \l_hobby_tria_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_trib_array}

% \Verb+\l_hobby_trib_array+ is an array holding the diagonal of the linear system that has to be solved to find the angles of the control points.

% In the following, this will be denoted by \(B_i\).

% The first index is \(0\).

%    \begin{macrocode}

\array_new:N \l_hobby_trib_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tric_array}

% \Verb+\l_hobby_tric_array+ is an array holding the superdiagonal of the linear system that has to be solved to find the angles of the control points.

% In the following, this will be denoted by \(C_i\).

% The first index is \(0\).

%    \begin{macrocode}

\array_new:N \l_hobby_tric_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_trid_array}

% \Verb+\l_hobby_trid_array+ is an array holding the target vector of the linear system that has to be solved to find the angles of the control points.

% In the following, this will be denoted by \(D_i\).

% The first index is \(1\).

%    \begin{macrocode}

\array_new:N \l_hobby_trid_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_triu_array}

% \Verb+\l_hobby_triu_array+ is an array holding the perturbation of the linear system for closed paths.

% The coefficient matrix for a \emph{open} path is tridiagonal and that means that Gaussian elimination runs faster than expected (\(O(n)\) instead of \(O(n^3)\)).

% The matrix for a closed path is not tridiagonal but is not far off.

% It can be solved by perturbing it to a tridiagonal matrix and then modifying the result.

% This array represents a utility vector in that perturbation. 

% In the following, the vector will be denoted by \(u\).

% The first index is \(1\).

%    \begin{macrocode}

\array_new:N \l_hobby_triu_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_excessangle_array}

% \Verb+\l_hobby_excessangle_array+ is an array that allows the user to say that the algorithm should add a multiple of \(2 \pi\) to the angle differences.

% This is because these angles are wrapped to the interval \((-\pi,\pi]\) but the wrapping might go wrong near the end points due to computation accuracy.

% The first index is \(1\).

%    \begin{macrocode}

\array_new:N \l_hobby_excessangle_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_psi_array}

% \Verb+\l_hobby_psi_array+ is an array holding the difference of the angles of the lines entering and exiting a point.

% That is, \(\psi_k\) is the angle between the lines joining \(z_k\) to \(z_{k-1}\) and \(z_{k+1}\).

% The first index is \(1\).

%    \begin{macrocode}

\array_new:N \l_hobby_psi_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_theta_array}

% \Verb+\l_hobby_theta_array+ is an array holding the angles of the outgoing control points for the generated path.

% These are measured relative to the line joining the point to the next point on the path.

% The first index is \(0\).

%    \begin{macrocode}

\array_new:N \l_hobby_theta_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_phi_array}

% \Verb+\l_hobby_phi_array+ is an array holding the angles of the incoming control points for the generated path.

% These are measured relative to the line joining the point to the previous point on the path.

% The first index is \(1\).

%    \begin{macrocode}

\array_new:N \l_hobby_phi_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_sigma_array}

% \Verb+\l_hobby_sigma_array+ is an array holding the lengths of the outgoing control points for the generated path.

% The units are such that the length of the line to the next specified point is one unit.

%    \begin{macrocode}

\array_new:N \l_hobby_sigma_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_rho_array}

% \Verb+\l_hobby_rho_array+ is an array holding the lengths of the incoming control points for the generated path.

% The units are such that the length of the line to the previous specified point is one unit.

%    \begin{macrocode}

\array_new:N \l_hobby_rho_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_controla_array}

% \Verb+\l_hobby_controla_array+ is an array holding the coordinates of the first control points on the curves.

% The format is the same as for \Verb+\l_hobby_points_array+.

%    \begin{macrocode}

\array_new:N \l_hobby_controla_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_controlb_array}

% \Verb+\l_hobby_controlb_array+ is an array holding the coordinates of the second control points on the curves.

% The format is the same as for \Verb+\l_hobby_points_array+.

%    \begin{macrocode}

\array_new:N \l_hobby_controlb_array

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_triv_fp}

% \Verb+\l_hobby_triv_fp+ is a number which is used when doing the perturbation of the solution of the linear system for a closed curve.

% There is actually a vector, \(v\), that this corresponds to but that vector only has one component that needs computation.

%    \begin{macrocode}

\fp_new:N \l_hobby_triv_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tempa_fp}

% \Verb+\l_hobby_tempa_fp+ is a temporary variable of type \Verb+fp+.

%    \begin{macrocode}

\fp_new:N \l_hobby_tempa_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tempb_fp}

% \Verb+\l_hobby_tempb_fp+ is a temporary variable of type \Verb+fp+.

%    \begin{macrocode}

\fp_new:N \l_hobby_tempb_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tempc_fp}

% \Verb+\l_hobby_tempc_fp+ is a temporary variable of type \Verb+fp+.

%    \begin{macrocode}

\fp_new:N \l_hobby_tempc_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_tempd_fp}

% \Verb+\l_hobby_tempd_fp+ is a temporary variable of type \Verb+fp+.

%    \begin{macrocode}

\fp_new:N \l_hobby_tempd_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_temps_fp}

% \Verb+\l_hobby_temps_fp+ is a temporary variable of type \Verb+fp+.

%    \begin{macrocode}

\fp_new:N \l_hobby_temps_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_incurl_fp}

% \Verb+\l_hobby_incurl_fp+ is the ``curl'' at the end of an open path.

% This is used if the angle at the end is not specified.

%    \begin{macrocode}

\fp_new:N \l_hobby_incurl_fp

\fp_set:Nn \l_hobby_incurl_fp {1}

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_outcurl_fp}

% \Verb+\l_hobby_outcurl_fp+ is the ``curl'' at the start of an open path.

% This is used if the angle at the start is not specified.

%    \begin{macrocode}

\fp_new:N \l_hobby_outcurl_fp

\fp_set:Nn \l_hobby_outcurl_fp {1}

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_inang_fp}

% \Verb+\l_hobby_inang_fp+ is the angle at the end of an open path.

% If this is not specified, it will be computed automatically.

% It is set to \Verb+\c_undefined_fp+ to allow easy detection of when it has been specified.

%    \begin{macrocode}

\fp_new:N \l_hobby_inang_fp

\fp_set_eq:NN \l_hobby_inang_fp \c_undefined_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_outang_fp}

% \Verb+\l_hobby_outang_fp+ is the angle at the start of an open path.

% If this is not specified, it will be computed automatically.

% It is set to \Verb+\c_undefined_fp+ to allow easy detection of when it has been specified.

%    \begin{macrocode}

\fp_new:N \l_hobby_outang_fp

\fp_set_eq:NN \l_hobby_outang_fp \c_undefined_fp

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\l_hobby_npoints_int}

% \Verb+\l_hobby_npoints_int+ is one less than the number of points on the curve.

% As our list of points starts at \(0\), this is the index of the last point.

% In the algorithm for a closed curve, some points are repeated whereupon this is incremented so that it is always the index of the last point. 

%    \begin{macrocode}

\int_new:N \l_hobby_npoints_int

%    \end{macrocode}

% \end{macro}

%

% A ``point'' is a key-value list setting the x-value, the y-value, and the tensions at that point.

% Using keys makes it easier to pass points from the algorithm code to the calling code and vice versa without either knowing too much about the other.

%    \begin{macrocode}

\keys_define:nn {hobby / read in all} {

  x .fp_set:N = \l_hobby_tempa_fp,

  y .fp_set:N = \l_hobby_tempb_fp,

  tension~out .fp_set:N = \l_hobby_tempc_fp,

  tension~in .fp_set:N = \l_hobby_tempd_fp,

  excess~angle .fp_set:N = \l_hobby_temps_fp,

  tension .meta:n = { tension~out=#1, tension~in=#1 },

  tension~out .default:n = 1,

  tension~in .default:n = 1,

  excess~angle .default:n = 0,

}

%    \end{macrocode}

% There are certain other parameters than can be set for a give curve.

%    \begin{macrocode}

\keys_define:nn { hobby / read in params} {

  in~angle .fp_set:N = \l_hobby_inang_fp,

  out~angle .fp_set:N = \l_hobby_outang_fp,

  in~curl .fp_set:N = \l_hobby_incurl_fp,

  out~curl .fp_set:N = \l_hobby_outcurl_fp,

  closed .bool_set:N = \l_hobby_closed_bool,

  closed .default:n = true,

  disjoint .bool_set:N = \l_hobby_disjoint_bool,

  disjoint .default:n = true,

}

%    \end{macrocode}

% \begin{macro}{\hobby_distangle:n}

% Computes the distance and angle between successive points.

% The argument given is the index of the current point.

% Assumptions: the points are in \Verb+\l_hobby_points_x_array+ and \Verb+\l_hobby_points_y_array+ and the index of the last point is \Verb+\l_hobby_npoints_int+.

%    \begin{macrocode}

\cs_set:Nn \hobby_distangle:n {

  \fp_set:Nn \l_hobby_tempa_fp

    {

      \array_get:Nn \l_hobby_points_x_array { #1 + 1 }

      - \array_get:Nn \l_hobby_points_x_array { #1 }

    }

  \fp_set:Nn \l_hobby_tempb_fp

    {

      \array_get:Nn \l_hobby_points_y_array { #1 + 1 }

      - \array_get:Nn \l_hobby_points_y_array { #1 }

    }

  \fp_atantwo:NNN \l_hobby_tempc_fp \l_hobby_tempb_fp   \l_hobby_tempa_fp

  \array_put:Nnx \l_hobby_angles_array {#1} {\fp_to_tl:N \l_hobby_tempc_fp}

  \array_put:Nnx \l_hobby_distances_array {#1}

    {

      \fp_eval:n

        { ( \l_hobby_tempa_fp ** 2 + \l_hobby_tempb_fp ** 2) ** .5 }

    }

  }

%    \end{macrocode}

% \end{macro}

%

%

% \begin{macro}{\fp_atantwo:NNN}

% Computes the angle of the point specified by the latter two arguments, storing the answer in the first.

% The inverse tangent function is not yet implemented in \LaTeX3 so for now we use the \Verb+pgfmath+ function.

% When this is implemented in \LaTeX3 this should be replaced.

%    \begin{macrocode}

\cs_new:Nn \fp_atantwo:NNN {

  \pgfmathparse{rad(atan2(\fp_use:N #3,\fp_use:N #2))}

  \exp_args:NNo \fp_set:Nn #1 {\pgfmathresult}

}

%    \end{macrocode}

% \end{macro}

%

%

% \begin{macro}{\hobby_ctrllen:NN}

%   Computes the length of the control point vector from the two angles.

%    \begin{macrocode}

\cs_new:Nn \hobby_ctrllen:NN {

  (2 - (cos(#1) - cos(#2))

       * (sin(#2) - sin(#1) / 16)

       * (sin(#1) - sin(#2) / 16)

       * 2 ** .5)

  / ((cos(#1) - cos(#2)) * (3 - 5 ** .5) / 2 + cos(#2) + 1)

}

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\hobby_append_point_copy:n}

%   This function adds a copy of the point (numbered by its argument) to

%   the end of the list of points, copying all the relevant data

%   (coordinates, tension, etc.).

%    \begin{macrocode}

\cs_new_protected:Npn \hobby_append_point_copy:n #1

  {

    \int_incr:N \l_hobby_npoints_int

    \hobby_append_point_copy_aux:Nn \l_hobby_points_array {#1}

    \hobby_append_point_copy_aux:Nn \l_hobby_points_x_array {#1}

    \hobby_append_point_copy_aux:Nn \l_hobby_points_y_array {#1}

    \hobby_append_point_copy_aux:Nn \l_hobby_tension_in_array {#1}

    \hobby_append_point_copy_aux:Nn \l_hobby_tension_out_array {#1}

    \hobby_append_point_copy_aux:Nn \l_hobby_excessangle_array {#1}

  }

\cs_new_protected:Npn \hobby_append_point_copy_aux:Nn #1#2

  { \array_put:Nnx #1 { \l_hobby_npoints_int } { \array_get:Nn #1 {#2} } }

%    \end{macrocode}

% \end{macro}

%

% \begin{macro}{\hobby_genpath:}

% This is the curve generation function.

% We assume at the start that we have an array containing all the points that the curve must go through, and the various curve parameters have been initialised.

% So these must be set up by a wrapper function which then calls this one.

% The list of required information is:

% \begin{enumerate}

% \item \Verb+\l_hobby_points_x_array+

% \item \Verb+\l_hobby_points_y_array+

% \item \Verb+\l_hobby_npoints_int+

% \item \Verb+\l_hobby_tension_out_array+

% \item \Verb+\l_hobby_tension_in_array+

% \item \Verb+\l_hobby_excessangle_array+

% \item \Verb+\l_hobby_incurl_fp+

% \item \Verb+\l_hobby_outcurl_fp+

% \item \Verb+\l_hobby_inang_fp+

% \item \Verb+\l_hobby_outang_fp+

% \item \Verb+\l_hobby_closed_bool+

% \end{enumerate}

%

%    \begin{macrocode}

\cs_new:Nn \hobby_genpath:

  {

%    \end{macrocode}

% For much of the time, we can pretend that a closed path is the same as

% an open path.  To do this, we need to make the end node an internal

% node by repeating the \(z_1\) node as the \(z_{n+1}\)th node.  We also

% check that the last (\(z_n\)) and first (\(z_0\)) nodes are the same,

% otherwise we repeat the \(z_0\) node as well.  If the first and last

% points are different, we need to duplicate the first point, with all

% its data.  Once we are sure that the first and last points are

% identical, we need to duplicate the first-but-one point (and all of

% its data).

%    \begin{macrocode}

    \bool_if:NT \l_hobby_closed_bool

      {

        \bool_if:nF

          {

            \fp_compare_p:n

              {

                \array_get:Nn \l_hobby_points_x_array {0}

                = \array_get:Nn \l_hobby_points_x_array { \l_hobby_npoints_int }

              }

            &&

            \fp_compare_p:n

              {

                \array_get:Nn \l_hobby_points_y_array {0}

                = \array_get:Nn \l_hobby_points_y_array { \l_hobby_npoints_int }

              }

          }

          { \hobby_append_point_copy:n {0} }

        \hobby_append_point_copy:n {1}

      }

%    \end{macrocode}

%

% Our first step is to go through the list of points and compute the distances and angles between successive points.

% Thus \(d_i\) is the distance from \(z_i\) to \(z_{i+1}\) and the angle is the angle of the line from \(z_i\) to \(z_{i+1}\).

%

%    \begin{macrocode}

    \prg_stepwise_function:nnnN

      {0} {1} {\l_hobby_npoints_int - 1}

      \hobby_distangle:n

%    \end{macrocode}

%

% For the majority of the code, we're only really interested in the differences of the angles.

% So for each internal point we compute the differences in the angles.

%    \begin{macrocode}

  \prg_stepwise_inline:nnnn {1} {1} {\l_hobby_npoints_int - 1}

    {

      \fp_set:Nn \l_hobby_tempa_fp

        {

          \array_get:Nn \l_hobby_angles_array {##1}

          - \array_get:Nn \l_hobby_angles_array { ##1 - 1 }

        }

%    \end{macrocode}

% We want to ensure that these angles lie in the range \((-\pi,\pi]\).

% So if the angle is bigger than \(\pi\), we subtract \(2 \pi\).

% (It shouldn't be that we can get bigger than \(3 \pi\) - check this.)

% Similarly, we check to see if the angle is less than \(-\pi\).

%    \begin{macrocode}

      \fp_compare:nTF { \l_hobby_tempa_fp > pi }

        { \fp_sub:Nn \l_hobby_tempa_fp { 2 pi } }

        {

          \fp_compare:nT { \l_hobby_tempa_fp <= - pi }

            { \fp_add:Nn \l_hobby_tempa_fp { 2 pi } }

        }

%    \end{macrocode}

% At the moment, the mixing of \Verb+pgfmath+ and \LaTeX3 that goes into computing the differences means that we can never get \Verb+\c_pi_fp+ so testing for equality is meaningless.

% However, when \Verb+atan2+ is implemented in \LaTeX3 then this might make sense again.

%

% The wrapping routine might not get it right at the edges so we add in the override.

%    \begin{macrocode}

      \array_get:NnNTF \l_hobby_excessangle_array {##1} \l_tmpa_tl

        { \fp_add:Nn \l_hobby_tempa_fp { 2 pi * \l_tmpa_tl } }

        { }

%    \end{macrocode}

%    \begin{macrocode}

      \array_put:Nnx \l_hobby_psi_array {##1}

        { \fp_to_tl:N \l_hobby_tempa_fp }

    }

%    \end{macrocode}

%

% Next, we generate the matrix.

% We start with the subdiagonal.

% This is indexed from \(1\) to \(n-1\).

%    \begin{macrocode}

  \prg_stepwise_inline:nnnn {1} {1} {\l_hobby_npoints_int - 1} {

    \array_put:Nnx \l_hobby_tria_array {##1}

      {

        \fp_eval:n

          {

            ( \array_get:Nn \l_hobby_tension_in_array {##1} ) ** 2

            * \array_get:Nn \l_hobby_distances_array {##1}

            * \array_get:Nn \l_hobby_tension_in_array { ##1 + 1 }

          }

      }

  }

%    \end{macrocode}

%

% Next, we attack main diagonal.

% We might need to adjust the first and last terms, but we'll do that in a minute.

%    \begin{macrocode}

  \prg_stepwise_inline:nnnn {1} {1} {\l_hobby_npoints_int - 1}

    {

      \array_put:Nnx \l_hobby_trib_array {##1}

        {

          \fp_eval:n

            {

              (3 * \array_get:Nn \l_hobby_tension_in_array {##1 + 1} - 1)

              * \array_get:Nn \l_hobby_tension_out_array {##1} ** 2

              * \array_get:Nn \l_hobby_tension_out_array {##1 - 1}

              * \array_get:Nn \l_hobby_distances_array {##1 - 1}

              +

              (3 * \array_get:Nn \l_hobby_tension_out_array {##1 - 1} - 1)

              * \array_get:Nn \l_hobby_tension_in_array {##1} ** 2

              * \array_get:Nn \l_hobby_tension_in_array {##1 + 1}

              * \array_get:Nn \l_hobby_distances_array {##1}

            }

        }

    }

%    \end{macrocode}

%

% Next, the superdiagonal.

%    \begin{macrocode}

  \prg_stepwise_inline:nnnn {1} {1} {\l_hobby_npoints_int - 2} {

    \array_put:Nnx \l_hobby_tric_array {##1}

      {

        \fp_eval:n

          {

            \array_get:Nn \l_hobby_tension_in_array {##1} ** 2

            * \array_get:Nn \l_hobby_tension_in_array {##1 - 1}

            * \array_get:Nn \l_hobby_distances_array {##1 - 1}

          }

      }

  }

%    \end{macrocode}

%

% Lastly (before the adjustments), the target vector.

%    \begin{macrocode}

  \prg_stepwise_inline:nnnn {1} {1} {\l_hobby_npoints_int - 2} {

  \array_get:NnN \l_hobby_psi_array {##1 + 1} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_tension_out_array {##1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_tension_out_array {##1 - 1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_distances_array {##1 - 1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_neg:N \l_hobby_tempa_fp

  \array_get:NnN \l_hobby_tension_out_array {##1 - 1} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {3}

  \fp_sub:Nn \l_hobby_tempb_fp {1}

  \array_get:NnN \l_hobby_psi_array {##1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_tension_in_array {##1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_tension_in_array {##1 + 1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_distances_array {##1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_sub:Nn \l_hobby_tempa_fp {\l_hobby_tempb_fp}

  \array_put:Nnx \l_hobby_trid_array {##1} {\fp_to_tl:N \l_hobby_tempa_fp}

}

%    \end{macrocode}

%

% Next, there are some adjustments at the ends.

% These differ depending on whether the path is open or closed.

%    \begin{macrocode}

\bool_if:NTF \l_hobby_closed_bool {

%    \end{macrocode}

% Closed path

%    \begin{macrocode}

\array_get:NnN \l_hobby_distances_array {\l_hobby_npoints_int - 2} \l_tmpa_tl

\fp_set:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int - 2} \l_tmpa_tl

\fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int - 1} \l_tmpa_tl

\fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\fp_neg:N \l_hobby_tempa_fp

\array_put:Nnx \l_hobby_tric_array {0} {\fp_to_tl:N \l_hobby_tempa_fp}

\array_put:Nnn \l_hobby_trib_array {0} {1}

\array_put:Nnn \l_hobby_trid_array {0} {0}

\array_get:NnN \l_hobby_trib_array {\l_hobby_npoints_int - 1} \l_tmpa_tl

\fp_set:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\fp_add:Nn \l_hobby_tempa_fp {1}

\array_put:Nnx \l_hobby_trib_array {\l_hobby_npoints_int - 1} {\fp_to_tl:N \l_hobby_tempa_fp}

  \array_get:NnN \l_hobby_psi_array {1} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

 \array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int -1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int -2} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

\array_get:NnN \l_hobby_distances_array {\l_hobby_npoints_int - 2} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_neg:N \l_hobby_tempa_fp

 \array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int - 2} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {3}

  \fp_sub:Nn \l_hobby_tempb_fp {1}

\array_get:NnN \l_hobby_psi_array {\l_hobby_npoints_int - 1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

\array_get:NnN \l_hobby_tension_in_array {\l_hobby_npoints_int - 1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

 \array_get:NnN \l_hobby_tension_in_array {\l_hobby_npoints_int} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

\array_get:NnN \l_hobby_distances_array {\l_hobby_npoints_int -1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_sub:Nn \l_hobby_tempa_fp {\l_hobby_tempb_fp}

 \array_put:Nnx \l_hobby_trid_array {\l_hobby_npoints_int - 1} {\fp_to_tl:N \l_hobby_tempa_fp}

%    \end{macrocode}

% We also need to populate the \(u\)-vector

%    \begin{macrocode}

  \array_put:Nnn \l_hobby_triu_array {0} {1}

\array_put:Nnn \l_hobby_triu_array {\l_hobby_npoints_int - 1} {1}

  \prg_stepwise_inline:nnnn {1} {1} {\l_hobby_npoints_int - 2} {

  \array_put:Nnn \l_hobby_triu_array {##1} {0}

  }

%    \end{macrocode}

% And define the significant entry in the \(v\)-vector.

%    \begin{macrocode}

\array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int -1} \l_tmpa_tl

\fp_set:Nn \l_hobby_triv_fp {\l_tmpa_tl}

\fp_mul:Nn \l_hobby_triv_fp {\l_tmpa_tl}

 \array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int -2} \l_tmpa_tl

\fp_mul:Nn \l_hobby_triv_fp {\l_tmpa_tl}

 \array_get:NnN \l_hobby_distances_array {\l_hobby_npoints_int -2} \l_tmpa_tl

\fp_mul:Nn \l_hobby_triv_fp {\l_tmpa_tl}

}

{

%    \end{macrocode}

% Open path.

% First, we test to see if \(\theta_0\) has been specified.

%    \begin{macrocode}

\fp_if_undefined:NTF \l_hobby_outang_fp

{

  \array_get:NnN \l_hobby_tension_in_array {1} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_hobby_incurl_fp}

  \array_get:NnN \l_hobby_tension_in_array {1} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {3}

  \fp_sub:Nn \l_hobby_tempb_fp {1}

  \array_get:NnN \l_hobby_tension_out_array {0} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_add:Nn \l_hobby_tempa_fp {\l_hobby_tempb_fp}

  \array_put:Nnx \l_hobby_trib_array {0}  {\fp_to_tl:N \l_hobby_tempa_fp}

  \array_get:NnN \l_hobby_tension_out_array {0} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \array_get:NnN \l_hobby_tension_out_array {0} \l_tmpa_tl

  \fp_set:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {3}

  \fp_sub:Nn \l_hobby_tempb_fp {1}

  \array_get:NnN \l_hobby_tension_in_array {1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_tmpa_tl}

  \fp_mul:Nn \l_hobby_tempb_fp {\l_hobby_incurl_fp}

  \fp_add:Nn \l_hobby_tempa_fp {\l_hobby_tempb_fp}

  \array_put:Nnx \l_hobby_tric_array {0} {\fp_to_tl:N \l_hobby_tempa_fp}

  \fp_neg:N \l_hobby_tempa_fp

  \array_get:NnN \l_hobby_psi_array {1} \l_tmpa_tl

  \fp_mul:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \array_put:Nnx \l_hobby_trid_array {0} {\fp_to_tl:N \l_hobby_tempa_fp}

}

{

  \array_put:Nnn \l_hobby_trib_array {0} {1}

  \array_put:Nnn \l_hobby_tric_array {0} {0}

  \fp_set_eq:NN \l_hobby_tempa_fp \l_hobby_outang_fp

  \array_get:NnN \l_hobby_angles_array {0} \l_tmpa_tl

  \fp_add:Nn \l_hobby_tempa_fp {\l_tmpa_tl}

  \array_put:Nnx \l_hobby_trid_array {0} {\fp_to_tl:N \l_hobby_tempa_fp}

}

%    \end{macrocode}

%

% Next, if \(\phi_n\) has been given.

%    \begin{macrocode}

\fp_if_undefined:NTF \l_hobby_inang_fp

{

 \array_get:NnN \l_hobby_tension_out_array {\l_hobby_npoints_int - 1} \l_tmpa_tl