'''
Length of an Archimedial spiral

Motivation:  How much toilet paper is on a roll?  One way to measure it
would be to roll it out.  This is perhaps the most accurate method.  But
it would be nice to be able to estimate it from the roll's dimensions.
The function ArchimedianSpiralArcLength below will help you do it.

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

The polar equation of this spiral is 

    r = a*theta

where theta is the angle and a is a constant.  For a spiral with 
multiple revolutions, the distance between the revolutions (pitch)
is 2*pi*a.

The arc length s is gotten from the integral from theta1 to theta2 
of 

    sqrt(r*r + (dr/dtheta)^2) dtheta

Substituting the equation for a spiral, we get

    A = sqrt(theta*theta + 1)
    s = a/2*[theta*A + ln(theta + A)]

This is the formula for the total arc length from theta = 0 to theta.

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

Copyright (C) 2006 Don Peterson

This program is free software; you can redistribute it
and/or modify it under the terms of the GNU General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be
useful, but WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the GNU General Public License for more
details.

You should have received a copy of the GNU General Public
License along with this program; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330,
Boston, MA 02111-1307  USA
'''

from __future__ import division
from math import sqrt, log, pi

def ArchimedianSpiralArcLength(a, theta, degrees=0):
    '''Calculate the arc length of an Archimedian spiral from angle
    0 to theta2.  theta is in radians unless degrees is true.  
    '''
    if a <= 0:
        raise "a must be >= 0"
    if degrees:
        c = pi/180
        theta *= c
    A = sqrt(theta*theta + 1)
    return a/2*(theta*A + log(theta + A))

if __name__ == "__main__":
    print '''Test of reasonableness for spiral arc length formula.  If we go
to a large radius of a spiral (many turns out), the arc length
over 2*pi radians should be about equal to the circumference of a
circle with that radius.  This approximation should get better
for larger theta.  Here, we'll print the ratio of the calculated
arc length to the ratio of the corresponding circle.  You'll see
that the ratio approaches 1 for larger theta.  The polar formula
for the spiral is r = a*theta.  The distance between successive
arcs of the spiral is constant at 2*pi*r.
'''
    a = 0.05
    print "a =", a, "\n"
    print " Theta,"
    print "radians    Ratio"
    print "--------  --------"
    thetas = [1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6]
    def p(xx):
        return "%9.3f" % xx
    for theta in thetas:
        r = a*theta  # Circular radius
        s1 = ArchimedianSpiralArcLength(a, theta)
        s2 = ArchimedianSpiralArcLength(a, theta + 2*pi)
        s = s2 - s1
        ratio = 2*pi*r/s
        #print "t =", p(theta), "r =", p(r), "s =", p(s), "ratio =", p(ratio)
        print "%7.0e   %f" % (theta, ratio)