[lots of new elastica material
raph.levien@gmail.com**20070512202417] {
addfile ./conclusion.tex
hunk ./conclusion.tex 1
+\chapter{Conclusion and Results}
+
+\section{Summary of contributions}
+
+\begin{itemize}
+
+\item Accurately plot the domain for which the angle-constrained MEC
+  has solutions.
+
+\item Formulate a spline more generally with heterogeneous constraints
+  at the control points so that $G^2$ and $G^4$ points can be mixed.
+
+\item Introduce ``one-way'' $G^2$ constraints, simplifying smooth
+  straight-to-curved transitions.
+
+\item Formulate a variational solution of the MVC problem as a simple
+  differential equation relating curvature to arclength.
+
+\end{itemize}
+
+
+\begin{itemize}
+
+\item Efficient integration of third-order polynomial spirals (spiro
+  curves).
+
+\item Fast and robust fitting of Euler spirals to angle constraints.
+
+\item Selection of optimal Bezier curves to fit mathematically-defined
+  curves, including spiro splines.
+
+\end{itemize}
+
hunk ./curves.bib 489
+
+@book{Statics,
+  author = "Lawrence E. Goodman and William H. Warner",
+  title = "Statics",
+  publisher = "Wadsworth",
+  year = 1963
+}
+
hunk ./elastica.tex 156
-expressible in the form of equation \ref{var-constraints}. In terms of
-curvature $\kappa$, the endpoint constraint is a double integral,
-which does not. Cartesian coordinates would make the endpoint
+expressible in the form of equation \ref{var-constraints}.
+
+It's now easy to see why representing the curve in terms of tangent
+angle as a function of arclength yields most readily to the
+variational approach. In terms of curvature $\kappa$, the endpoint
+constraint is a double integral, which does not fit neatly into a
+formulation of constraints as a single integral of some
+function. Alternatively, Cartesian coordinates would make the endpoint
hunk ./elastica.tex 165
-function becomes more complex. Even so, Euler chose Cartesian
-coordinates for his 1744 solution to the elastica problem (see 
+function becomes more complex, requiring second derivatives as well as
+the first derivatives of the Euler-Lagrange equation. The first
+solution was indeed done using Cartesian coordinates, but this
+approach is most definitely simpler.
hunk ./elastica.tex 225
-\kappa)(d^2 kappa / ds^2)$ is the longitudinal force on the beam. ***
-Present equations of equilibrium, and show that variational solution
-is same as equilibrium.
+\kappa)(d^2 kappa / ds^2)$ is the longitudinal force on the beam.
+
+
+%Present equations of equilibrium, and show that variational solution
+%is same as equilibrium.
+
+\section{A family of solutions}
+\label{elastica-fam-sect}
+
+The equations derived in the previous section for the shape of the
+elastica don't have a single solution. Rather, they describe a family
+of solutions, characterized by a single scalar parameter. All such
+solutions are realizable by a phyiscal elastica; none are purely an
+artifact of the mathematical approach.
+
+For mathematical convenience, assign $\kappa = 1$ at the midpoint of
+the curve, and use $\lambda_1$ as the parameter, setting $\lambda_2 =
+0$. All other instances of the elastica can be transformed to meet
+these conditions by scaling and rotation. Since these assumptions
+remove the need to distinguish multiple $\lambda$ parameters, it is
+convenient to drop the ${}_1$ subscript and simply write $\lambda$.
+
+Figure \ref{figs/elastica-fam} shows a representative selection from
+this family, for various values of $\lambda$.
+
+\begin{figure*}
+\begin{center}
+\includegraphics[scale=.9]{figs/elastica-fam}
+\caption{\label{elastica-fam}The family of elastica solutions.}
+\end{center}
+\end{figure*}
+
+The remainder of this chapter is devoted to exploring these solutions
+and their physical meaning, but a few highlights are immediately
+apparent. All solutions have a mirror (even) symmetry at each
+curvature maximum. All solutions save $\lambda = 0.25$ are
+periodic. The solutions for values larger than $0.25$ have
+inflections, and have odd symmetry around each inflection point. For
+values smaller than $0.25$, there is no inflection; the curvature
+always remains the same sign.
+
+As $\lambda$ approaches $0$, the elastica tends to a circle. As
+$\lambda$ approaches $\infty$, the elastica tends to a sine wave. The
+solution at $\lambda = 0.5$ is special, and is known as the
+\emph{rectangular elastica}. The tangents at the inflection points of
+this curve perpendicular to the horizontal axis; the curve is tangent
+to the rectangle formed by inflections and curvature extrema.
+
+\begin{figure*}
+\begin{center}
+\includegraphics[scale=.9]{figs/elastica-fam-k}
+\caption{\label{elastica-fam-k}The family of pendulum solutions.}
+\end{center}
+\end{figure*}
+
+Figure \ref{elastica-fam-k} shows the curvature plots of the 
+instances of the elastica family shown in Figure \ref{elastica-fam},
+plotted with curvature on the vertical axis and arclength on the
+horizontal. As is usual, inflection points on the curve correspond to
+the curvature crossing (or touching) the horizontal axis.
+
+\section{Equilibrium of moments}
+
+\begin{figure*}
+\begin{center}
+\includegraphics[scale=.9]{figs/moment}
+\caption{\label{moment}Equilibrium with moments.}
+\end{center}
+\end{figure*}
+
+Analyzing the elastica using the theory of elasticity yields three
+physical quantities which must be in equilibrium along its length:
+longitudinal force (tension), shearing force, and bending moment. The
+tension and shearing forces are useful when working out the effect of
+constraints, but it is possible to derive the equation for the shape
+of the elastica from moments alone. Indeed, computing the moments is
+by far the easiest derivation of the elastica equation.
+
+This section uses of very elementary principles of statics to carry
+out the derivation of the elastica from equilibrium of moments. A
+classic textbook introducing these concepts is \emph{Statics,} by
+Goodman and Warner\cite{Statics}.
+
+Consider the elastica in Figure \ref{moment}. Two forces, $\mathbf{F}$ and
+$-\mathbf{F}$, push the ends of the elastica together. Because these forces are
+equal in magnitude, opposite in direction, and have parallel lines of
+action, they form a \emph{couple}\cite[p. 65]{Statics}.
+
+At any point along the elastica, then, the \emph{resultant moment} is
+simply the force of the couple times the distance from the line of
+action: $M = Fy$, where $F$ is the magnitude of the force
+vector $\mathbf{F}$.
+
+The simple theory of elasticity says that the relation between bending
+moment and curvature is linear. Thus, collecting all the constants
+together into a single $\lambda$ yields this remarkably simple formulation
+of the shape of a general elastica:
+
+\begin{equation}
+\kappa = \lambda y
+\end{equation}
+
+This equation is readily shown to be equivalent to the variational
+solution above. Take the derivative with respect to arclength $s$, and
+note that $y' = \sin \theta$. Then,
+
+\begin{equation}
+\kappa' = c \sin \theta,
+\end{equation}
+
+and this equation is obviously equivalent to Equation
+\ref{elastica-eq}, assuming without loss of generality that the line
+of action is parallel to the horizontal axis.
+
+All elastica solutions are equivalent to the system shown in Figure
+\ref{moment}, but in non-inflectional cases the line of action of the
+couple is virtual, in that it doesn't touch the curve.
hunk ./elastica.tex 356
-vertical.
+the pendulum in its vertical resting state.
hunk ./elastica.tex 373
-momentum. According to Love\cite{Love}, the systematic application of the
-kinetic analogue was worked out by W. Hess in 1885.
+momentum.
+
+% According to Love\cite{Love}, the systematic application of the
+%kinetic analogue was worked out by W. Hess in 1885.
hunk ./elastica.tex 386
-\section{A family of solutions}
-\label{elastica-fam-sect}
-
-The solutions to the elastica equation \ref{elastica-eq} form a family
-characterized by a single scalar parameter, once the trivial
-transforms of rotation and scaling are factored out. Without loss of
-generality, we will assign the bottom of the swing (i.e. maximum
-velocity) to $t = 0$ and assume a unit length. The remaining parameter
-is, then, the ratio of the total energy of the system (which remains
-unchanged over time) 
+The solutions of the motion of the pendulum, as do the solutions of
+the elastica equation \ref{elastica-eq}, form a family characterized
+by a single scalar parameter, once the trivial transforms of rotation
+and scaling are factored out. Without loss of generality, let $t = 0$
+a\t the bottom of the swing (i.e. maximum velocity) and let the
+pendulum have unit length. The remaining parameter is, then, the ratio
+of the total energy of the system (which remains unchanged over time)
hunk ./elastica.tex 422
-Figure \ref{elastica-fam-k} shows various solutions of the elastica
-equation. For each value of $\lambda$, velocity ($\theta'$) is plotted
+Figure \ref{elastica-fam-k}, the curvature plots for solutions of the
+elastica equation, thus also shows the velocity of the pendulum as a
+function of time for various solutions of the motion of the
+pendulum. For each value of $\lambda$, velocity ($\theta'$) is plotted
hunk ./elastica.tex 429
-\begin{figure*}
-\begin{center}
-\includegraphics[scale=.9]{figs/elastica-fam-k}
-\caption{\label{elastica-fam-k}The family of pendulum solutions.}
-\end{center}
-\end{figure*}
-
hunk ./elastica.tex 439
-corresponds to curvature. Thus, only values of $\lambda$ above
-$\frac{1}{4}$ have inflection points. The elastica curves corresponding
-to the same values of $\lambda$ plotted in Figure \ref{elastica-fam-k}
-are shown in Figure \ref{elastica-fam}.
+corresponds to curvature. Thus, values of $\lambda$ above
+$\frac{1}{4}$ have inflection points.
hunk ./elastica.tex 442
-\begin{figure*}
-\begin{center}
-\includegraphics[scale=.9]{figs/elastica-fam}
-\caption{\label{elastica-fam}The family of elastica solutions.}
-\end{center}
-\end{figure*}
hunk ./elastica.tex 443
-\section{Mechanical equilibrium}
+\section{Finite element approach}
hunk ./elastica.tex 449
+\begin{figure*}[tbh]
+\begin{center}
+\includegraphics{figs/chain}
+\caption{\label{chain}Finite element approximation of elastica as
+  chain of struts and pivots.}
+\end{center}
+\end{figure*}
+
hunk ./elastica.tex 459
-shear, normal, and tension forces, and flexural couples. The equations
-are relatively straightforward, but the relevant physics is not
-exactly trivial. This section, therefore, will present a simplified
-finite-element model with only bending and tension forces.
+shear, normal, and longitudinal forces, and flexural couples. The
+equations are relatively straightforward, but the relevant physics is
+not exactly trivial. This section, therefore, presents a
+simplified finite-element model with only bending and tension forces.
hunk ./elastica.tex 472
-struts and pivots with springs. Each strut generates a tension force
-at both ends, in the same direction as the strut's length. In
+struts and pivots with torsion springs. Such a chain, with spring
+torque shown schematically at each pivot, is shown in Figure
+\ref{chain}.
+
+Each strut generates a tension
+force at both ends, in the same direction as the strut's length. In
hunk ./elastica.tex 479
-$P$. By convention, compression force is denoted with a negative sign.
+$T$. By convention, compression force is denoted with a negative sign.
hunk ./elastica.tex 491
-the ``moment force'' by $M$. At the endpoints, the moment force is
-normal to the strut joining the center to the endpoint. According to
-Newton's third law, the sum of all forces arising from the element
-must be zero, so there is a balancing force at the center point equal
-to the vector sum of the forces at the endpoints, pointing in the
-opposite direction.
+the the force generated by the moment (or spring torque) by $M$. At
+the endpoints, the force generated by the moment is normal to the
+strut joining the center to the endpoint. According to Newton's third
+law, the sum of all forces arising from the element must be zero, so
+there is a balancing force at the center point equal to the vector sum
+of the forces at the endpoints, pointing in the opposite direction.
hunk ./elastica.tex 499
-is to say that the moment force is proportional to the turning
+is to say that the moment torque is proportional to the turning
hunk ./elastica.tex 518
-single pivot point: the tension force from each adjoining strut ($P_0$
-and $-P_1$), the moment forces from the previous and next pivots in
+single pivot point: the tension force from each adjoining strut ($T_0$
+and $-T_1$), the moment forces from the previous and next pivots in
hunk ./elastica.tex 529
--P_1 + \cos\theta_1\;P_0 - \sin\theta_1(M_0 - M_1) = 0
+-T_1 + \cos\theta_1\;T_0 - \sin\theta_1(M_0 - M_1) = 0
hunk ./elastica.tex 535
-\sin\theta_1\;P_0 + \cos\theta_1(M_0 - M_1) + M_2 - M_1 = 0
+\sin\theta_1\;T_0 + \cos\theta_1(M_0 - M_1) + M_2 - M_1 = 0
hunk ./elastica.tex 545
--P_1 + P_0 - c\theta_1(\theta_0 - \theta_1) = 0
+-T_1 + T_0 - c\theta_1(\theta_0 - \theta_1) = 0
hunk ./elastica.tex 548
-\theta_1\;P_0 + c(\theta_0 + \theta_2 - 2\theta_1) = 0
+\theta_1\;T_0 + c(\theta_0 + \theta_2 - 2\theta_1) = 0
hunk ./elastica.tex 557
--\frac{dP}{ds} + c\kappa\frac{d\kappa}{ds} = 0 \\
+-\frac{dT}{ds} + c\kappa\frac{d\kappa}{ds} = 0 \\
hunk ./elastica.tex 561
-\kappa P + c\frac{d^2\kappa}{ds^2} = 0
+\kappa T + c\frac{d^2\kappa}{ds^2} = 0
hunk ./elastica.tex 567
-P = \frac{c}{2}\kappa^2 + C
+T = \frac{c}{2}\kappa^2 + C
hunk ./elastica.tex 575
-regardless of the value of this coefficient.
+regardless of the value of this coefficient, as long as it is
+consistent across the length of the elastica.
hunk ./elastica.tex 588
-in terms of curvature and tension (longitudinal force).
+in terms of curvature and tension (longitudinal force). In particular,
+if $C$ is zero, then the tension vanishes at inflection points.
hunk ./elastica.tex 642
+\section{Elastica with general constraints}
+
+Working towards the application of an elastica to splines, this
+section considers the effect of general constraints on an
+elastica. A goal is to develop physical intuition, so the constraints
+are expressed in terms of forces.
+
+There are four basic types of constraints. Each constraint, can either
+rotate freely or fix the tangent angle of the elastica. Further, it
+can either allow the elastica to slide freely, or can fix the relative
+arclength of the elastica. A segment between two constraints that fix
+relative arclength has fixed arclength, even if there are sliding
+constraints in between. In all these types of constraints, the
+position of the constraint is fixed, and the elastica is constrained
+to go through that point.
+
+A position constraint, freely rotating, and with the elastica free to
+slide through the point, exerts only a normal force on the
+elastica. If the tangent angle is fixed, the constraint exerts a shear
+force on the elastica. Similarly, if the relative arclength is fixed,
+the constraint exerts a longitudinal force.
+
+\newtheorem{result}{Result}
+
+\begin{result}
+A normal force acting on an elastica preserves $G^2$
+  continuity but creates a discontinuity in the first derivative of curvature.
+\end{result}
+
+\begin{result}
+A shear force (torque) acting on an elastica preserves $G^1$
+continuity but creates a discontinuity in curvature.
+\end{result}
+
+\begin{result}
+An elastica segment with zero longitudinal force at the endpoints
+takes the form of a rectangular elastica. Nonzero longitudinal force
+induces values of lambda other than $.5$.
+\end{result}
+
+\begin{figure*}
+\begin{center}
+\includegraphics[scale=.8]{figs/rectel}
+\caption{\label{rectel-fig}Forces acting on rectangular elastica.}
+\end{center}
+\end{figure*}
+
+Considering a full cycle of the rectangular elastica, as shown in
+Figure \ref{rectel-fig}, it is clear that the longitudinal forces are
+zero at the endpoints, because the curve is perpendicular to the force
+$\mathbf{F}$ acting on it.
+
+Assume that the minimum energy curve for given angle constraints at
+the endpoints was found \emph{not} to be a rectangular elastica, in
+other words that the curve has lower energy than the rectangular
+elastica section matching the angle constraints. Then,
+find the chord across one cycle of the rectangular elastica, so that
+the angles of the elastica with respect to that chord match the angle
+constraints given. Now, scale and rotate the curve so that its
+endpoints match the endpoints of the chord across the rectangular
+elastica. Replace that segment of the rectangular elastica with the
+scaled and rotated curve. By the assumption, the resulting curve will
+have lower energy than the rectangular elastica. But then, the
+rectangular elastica is not the minimum energy curve for the system as
+shown in Figure \ref{rectel-fig}, which is a contradiction.
+
+\section{Angle constraints}
+
+The general constraint approach of the previous section yields a
+specific result relevant to using the elastica as a spline.
+
+\begin{result}
+An elastica segment with an angle constraint at each endpoint, but
+allowed to slide freely through those constraints, is a section cut from
+the rectangular elastica.
+\end{result}
+
+Equivalently, the minimum energy curve (of arbitrary length) with
+given endpoints and given angles at the endpoints is a section cut from
+the rectangular elastica.
+
hunk ./elastica.tex 726
-\caption{\label{mecrange-fig}mecrange}
+\caption{\label{mecrange-fig}Angle-constrained minimum energy curves.}
hunk ./elastica.tex 730
+Figure \ref{mecrange-fig} shows the full parameter space of
+angle-constrained minimum energy curves. For sufficiently large
+angles, no such solution exists; the envelope for which the MEC
+converges is shows as well.
+
hunk ./elastica.tex 738
-\caption{\label{mecrangek-fig}mecrangek}
+\caption{\label{mecrangek-fig}Curvature plots, angle-constrained
+  minimum energy curves.}
hunk ./elastica.tex 743
+Similarly, Figure \ref{mecrangek-fig} shows the curvature plots of all
+angle-constrained MEC's within their parameter space. It is worth
+noting that, for reasonably small angles, the curvature is very nearly
+linear with arclength.
hunk ./elastica.tex 748
-\section{The curve of least energy}
+Note also that the curvatures at the endpoints are fully determined by
+the angles at the endpoints.
hunk ./elastica.tex 751
-We now consider a specific set of constraints. Assume endpoint
-constraints through which the elastica can slide (friction-free), but
-which set a particular angle. Because of the sliding nature of the
-endpoint constraints, the length of the elastica becomes an additional
-parameter to be varied.
+\section{Elastica as interpolating spline}
hunk ./elastica.tex 753
-Forsythe and Lee presented a solution to this case\cite{ForsytheLee},
-as did Horn\cite{Horn} several years later. Simply, friction-free
-sliding at the endpoints selects the $\lambda = 0.5$ case from the
-family considered in Section \ref{elastica-fam-sect}. This case also
-corresponds to $C = 0$ in the curvature-only formulation of the
-elastica equation (Equation \ref{fl-elastica-eq}).
+The classical formulation of an interpolating spline is the mimimum
+energy curve passing through a sequence of control points.
+Equivalently, it is an elastica passing through a sequence of
+constraints, each of which rotates freely and allows the elastica to
+slide freely through it.
hunk ./elastica.tex 759
-Horn considered primarily the case where the endpoint angle constraints are
-both perpendicular to the chord between the endpoints, and facing in
-the same direction (i.e., even symmetry and 180 degrees of total arc).
+Because each of these constraints exerts only a normal force, the
+resulting curve is $G^2$ continuous. Since the elastica is unbent past
+either endpoint, the curvature at the endpoints is zero (this is also
+known as the ``natural'' endpoint condition).
hunk ./elastica.tex 764
-\begin{figure*}
-\begin{center}
-\includegraphics[scale=1]{figs/horn}
-\caption{\label{horn-fig}The curve of least energy.}
-\end{center}
-\end{figure*}
+Each segment of the spline, bounded by two control points, is an
+instance of an angle-constrained MEC as described in the previous
+section.
+
+In order to find an interpolating spline, it is necessary and
+sufficient to find a global solution assigning a tangent angle to each
+point so that the resulting curvatures are continuous across each
+point (and zero at the endpoints).
+
+\section{SIMEC}
hunk ./elastica.tex 775
-Horn solution. Arbitrary angle solution is cut from 180-degree
-solution. 180-deg soln has zero curvature at endpoints, but this is
-not true for arbitrary angles. Divergence when angles exceed 180
-degrees (or should this be in convergence, below).
+Scale Invariant Minimum Energy Curve (SIMEC) minimizes the functional:
+
+\begin{equation}
+SIMEC[\kappa(s)] = l\int_0^l \kappa(s)^2\ ds
+\end{equation}
hunk ./elastica.tex 781
-\section{Elastica with length constraints}
+Equivalent to an elastica with relative arclength constraints at
+endpoints, with total arclength chosen to be a minimum of SIMEC functional.
hunk ./elastica.tex 784
+Gluing two SIMEC's together is not a SIMEC, as longitudinal force at
+glue point may be non-zero.
hunk ./elastica.tex 787
-\section{Closed-form solutions of the elastica equation}
+SIMEC is round, and doesn't run away like MEC (length is penalized).
hunk ./elastica.tex 789
+\begin{result}
+\end{result}
hunk ./figs/Makefile 9
-	strutfig.pdf pivotfig.pdf chainfig.pdf watchspring.pdf
+	strutfig.pdf pivotfig.pdf chainfig.pdf watchspring.pdf \
+	moment.pdf rectel.pdf chain.pdf bernoulli91.pdf
hunk ./figs/Makefile 58
+moment.pdf:	elastica.py
+	python elastica.py $@ > /tmp/foo.ps
+	python ps2pdfbbox.py /tmp/foo.ps $@
+
+rectel.pdf:	elastica.py
+	python elastica.py $@ > /tmp/foo.ps
+	python ps2pdfbbox.py /tmp/foo.ps $@
+
+chain.pdf:	elastica.py
+	python elastica.py $@ > /tmp/foo.ps
+	python ps2pdfbbox.py /tmp/foo.ps $@
+
+bernoulli91.pdf:	elastica.py
+	python elastica.py $@ > /tmp/foo.ps
+	python ps2pdfbbox.py /tmp/foo.ps $@
+
hunk ./figs/elastica.py 93
-def plot(xyk, x0, y0, sx, sy):
+def plot(xys, x0, y0, sx, sy):
hunk ./figs/elastica.py 95
-    for i in range(len(xyk)):
-	x, y, k = xyk[i]
-	print x0 + x * sx, y0 + y * sy, cmd
+    for i in range(len(xys)):
+	xy = xys[i]
+	print x0 + xy[0] * sx, y0 + xy[1] * sy, cmd
hunk ./figs/elastica.py 113
+def moment_arrow(x0, y0, r, th0, th1):
+    wid = 2
+    len = min(3, r * (th1 - th0) * .005)
+    print x0 + r * cos(th0 * pi / 180), y0 + r * sin(th0 * pi / 180), 'moveto'
+    print x0, y0, r, th0, th1, 'arc'
+    print 'stroke'
+    print 'gsave'
+    print x0 + r * cos(th0 * pi / 180), y0 + r * sin(th0 * pi / 180), 'translate'
+    print th0 - 90, 'rotate'
+    print -len, -wid, 'moveto 0 0 lineto', -len, wid, 'lineto stroke'
+    print 'grestore'
+    print 'gsave'
+    print x0 + r * cos(th1 * pi / 180), y0 + r * sin(th1 * pi / 180), 'translate'
+    print th1 + 90, 'rotate'
+    print -len, -wid, 'moveto 0 0 lineto', -len, wid, 'lineto stroke'
+    print 'grestore'
+
+def draw_moment():
+    import elastfe
+    print '%!PS-Adobe-3.0 EPSF'
+    s = 3.12676
+    xyk, cost = run_elastica(-s, s, 1, 0, -.4)
+    print '% k = ', xyk[0][2]
+    x0, y0 = 300, 500
+    sx = -100
+    sy = -100
+    f = 50
+    plot(xyk, x0, y0, sx, sy)
+    print 'gsave 2 setlinewidth'
+    xl = x0 + xyk[-1][0] * sx
+    yl = y0 + xyk[-1][1] * sy
+    elastfe.arrow(xl - f, yl, f, 0)
+    xr = x0 + xyk[0][0] * sx
+    yr = y0 + xyk[0][1] * sy
+    elastfe.arrow(xr + f, yr, f, 180)
+    print 'grestore'
+    print '/Times-Bold 12 selectfont'
+    print xl - 30, yl + 4, 'moveto (-F) show'
+    print xr + 25, yl + 4, 'moveto (F) show'
+    print '.75 setlinewidth'
+    print x0, yl, 'moveto', 0, 300, 'rlineto stroke'
+    print xl, yl, 'moveto', xr, yr, 'lineto stroke'
+    print '/Times-Italic 12 selectfont'
+    if 0:
+	print x0 + 4, yl + 290, 'moveto (y) show'
+	print xr - 10, yl - 8, 'moveto (x) show'
+    xs = x0 + xyk[26][0] * sx
+    ys = y0 + xyk[26][1] * sy
+    print 'gsave [2] 0 setdash'
+    print xs, yl, 'moveto', xs, ys, 'lineto stroke'
+    print 'grestore'
+    print xs - 8, (yl + ys) * .5, 'moveto (y) show'
+    moment_arrow(xs, ys, 10, -70, 200)
+    print xs + 12, ys + 2, 'moveto'
+    #print '/Times-Bold 12 selectfont (M) show'
+    print '/Times-Roman 12 selectfont (M = F) show'
+    #print '/Times-Bold 12 selectfont (F) show'
+    print '/Times-Italic 12 selectfont (y) show'
+    print 'showpage'
+
+def draw_rectel():
+    import elastfe
+    print '%!PS-Adobe-3.0 EPSF'
+    s = 2.63
+    xyk, cost = run_elastica(-3 * s, s, 1, 0, -.5)
+    print '% k = ', xyk[0][2]
+    x0, y0 = 300, 500
+    sx = -50
+    sy = -50
+    f = 50
+    plot(xyk, x0, y0, sx, sy)
+    print 'gsave 2 setlinewidth'
+    xl = x0 + xyk[-1][0] * sx
+    yl = y0 + xyk[-1][1] * sy
+    elastfe.arrow(xl - f, yl, f, 0)
+    xr = x0 + xyk[0][0] * sx
+    yr = y0 + xyk[0][1] * sy
+    elastfe.arrow(xr + f, yr, f, 180)
+    print 'grestore'
+    print '/Times-Bold 12 selectfont'
+    print xl - 30, yl + 4, 'moveto (-F) show'
+    print xr + 25, yl + 4, 'moveto (F) show'
+    print '.75 setlinewidth'
+    print 'showpage'
+
+
+
+def run_elastica_half_th(sp, k0, lam1, lam2, n = 10):
+    def mec_derivs(x, ys):
+	th, k = ys[2], ys[3]
+	dx, dy = cos(th), sin(th)
+	return [dx, dy, k, lam1 * dx + lam2 * dy, k * k]
+    ys = [0, 0, 0, k0, 0]
+    xytk = [ys[:4]]
+    h = float(sp) / n
+    s = 0
+    for i in range(n):
+	dydx = mec_derivs(s, ys)
+	rk4(ys, dydx, s, h, mec_derivs)
+	xytk.append(ys[:4])
+    return xytk
+
+def draw_chain():
+    print '%!PS-Adobe-3.0 EPSF'
+    print '100 400 translate'
+    print '-70 rotate'
+    print '/ss 2 def'
+    print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
+    xytks = run_elastica_half_th(4, 1, 0, -.5, 10)
+    x0, y0 = 0, 0
+    sc = 100
+    th0 = -90
+    plot(xytks, x0, y0, sc, sc)
+    for i in range(len(xytks)):
+	xy = xytks[i]
+	print 'gsave', x0 + sc * xy[0], y0 + sc * xy[1], 'translate circle fill'
+	print 180 * xy[2] / pi, 'rotate'
+	if i > 0 and i < len(xytks) - 1:
+	    if (xy[3] > 0):
+		moment_arrow(0, 0, 12, th0 - 60 * xy[3], th0 + 60 * xy[3])
+	    else:
+		moment_arrow(0, 0, 12, -th0 + 60 * xy[3], -th0 - 60 * xy[3])
+	print 'grestore'
+    print 'showpage'
+
+def draw_bernoulli91(Eix = 45):
+    print '%!PS-Adobe-3.0 EPSF'
+    xytk = run_elastica_half_th(2.63, 1, 0, -.5, 100)
+    scale = 100
+    fontsize = 21
+    D = (200, 300)
+    print '2.5 setlinewidth'
+    cmd = 'moveto'
+    for i in range(len(xytk)):
+	xy = xytk[i]
+	print D[0] + xy[1] * scale, D[1] + xy[0] * scale, cmd
+	if i == Eix: E = (D[0] + xy[1] * scale, D[1] + xy[0] * scale)
+	cmd = 'lineto'
+    print 'stroke'
+    A = (D[0] + xytk[-1][1] * scale, D[1] + xytk[-1][0] * scale)
+    AC = A[0] - D[0]
+    # The evolute
+    cmd = 'moveto'
+    for i in range(len(xytk)):
+	xy = xytk[i]
+	r = 1 / xy[3]
+	x = xy[1] + r * cos(xy[2])
+	y = xy[0] - r * sin(xy[2])
+	print D[0] + x * scale, D[1] + y * scale, cmd
+	if i == 0: G = (D[0] + x * scale, D[1] + y * scale)
+	if i == Eix: L = (D[0] + x * scale, D[1] + y * scale)
+	cmd = 'lineto'
+	if y < -2:
+	    M = (D[0] + x * scale, D[1] + y * scale)
+	    break
+    print 'stroke'
+    print '1.5 setlinewidth'
+    print 'gsave [4] 0 setdash'
+    print L[0], L[1], 'moveto', E[0], E[1], 'lineto'
+    print A[0], E[1], 'lineto stroke'
+    print '[6] 0 setdash'
+    print D[0], D[1], 'moveto', D[0], A[1], 'lineto stroke'
+    print 'grestore'
+
+    # axes and other guide lines
+    N = (A[0], D[1] - 2 * scale)
+    B = (A[0], A[1] + AC)
+    print N[0], N[1], 'moveto', B[0], B[1], 'lineto stroke'
+    print D[0], D[1], 'moveto', A[0], D[1], 'lineto stroke' # DH
+    print L[0], L[1], 'moveto', A[0], L[1], 'lineto stroke' # LI
+    print A[0], A[1], 'moveto', D[0], A[1], 'lineto stroke' # AC
+    print B[0], B[1], 'moveto', A[0], A[1], AC, 90, 180, 'arc stroke'
+    Q = (.5 * (G[0] + A[0]), D[1] + .5 * (A[0] - G[0]))
+    print A[0], D[1], 'moveto'
+    print Q[0], D[1], Q[1] - D[1], 0, 180, 'arc stroke'
+    print G[0], G[1], 'moveto', Q[0], Q[1], 'lineto'
+    print A[0], D[1], 'lineto stroke'
+
+    # labels
+    print '/Times-Roman', fontsize, 'selectfont'
+    print A[0] + 3, A[1], 'moveto (A) show'
+    print B[0] - .9 * fontsize, B[1] - .8 * fontsize, 'moveto (B) show'
+    print D[0] - .9 * fontsize, A[1] - fontsize / 3, 'moveto (C) show'
+    print D[0] - .9 * fontsize, D[1] - fontsize / 3, 'moveto (D) show'
+    print E[0] - .6 * fontsize, E[1] + 1, 'moveto (E) show'
+    print A[0] + 3, E[1] - 2, 'moveto (F) show'
+    print G[0] - .9 * fontsize, G[1] + 2, 'moveto (G) show'
+    print A[0] + 3, G[1], 'moveto (H) show'
+    print A[0] + 3, L[1] - fontsize / 3, 'moveto (I) show'
+    print L[0] - .7 * fontsize, L[1] - fontsize / 3, 'moveto (L) show'
+    print M[0] - 1.2 * fontsize, M[1] + 12, 'moveto (M) show'
+    print N[0] + 3, N[1], 'moveto (N) show'
+    print Q[0] - .4 * fontsize, Q[1] + 5, 'moveto (Q) show'
+    print 'showpage'
+
hunk ./figs/elastica.py 347
-	eps_prologue(150, 290, 450, 501)
-	xyk = compute_elastica(-2.63, 2.63, 1, .5)
-	plot(xyk, 300, 500, -100, -100)
-	print '1 0 0 setrgbcolor'
+	if 0:
+	    eps_prologue(150, 290, 450, 501)
+	    xyk = compute_elastica(-2.63, 2.63, 1, .5)
+	    plot(xyk, 300, 500, -100, -100)
+	    print '1 0 0 setrgbcolor'
hunk ./figs/elastica.py 356
+    elif figure == 'moment.pdf':
+	draw_moment()
+
+    elif figure == 'rectel.pdf':
+	draw_rectel()
+
+    elif figure == 'chain.pdf':
+	draw_chain()
+
+    elif figure == 'bernoulli91.pdf':
+	if len(sys.argv) > 2:
+	    draw_bernoulli91(int(sys.argv[2]))
+	else:
+	    draw_bernoulli91()
+
+
hunk ./frame.tex 30
+
+\include{conclusion}
hunk ./hist_elast.tex 60
-\caption{\label{bernoulli-clxx}Drawing from James Bernoulli's 1692 Med. CLXX.}
+\hspace{.5in}
+\raisebox{.35in}{\includegraphics[width=1.25in]{figs/bernoulli91.pdf}}
+\caption{\label{bernoulli-clxx}Drawing from James Bernoulli's 1692 Med. CLXX, with modern reconstruction.}
hunk ./hist_elast.tex 68
-\ref{elastica}). In his Meditatione CLXX dated that year, titled
+\ref{elastica-chapter}). In his Meditatione CLXX dated that year, titled
hunk ./hist_elast.tex 76
-Radius circuli $AG = a$, $AED$ lamina elastica ab appenso pondere in
+Radius circuli $AB = a$, $AED$ lamina elastica ab appenso pondere in
}
