++ed by:
Kevin Ryde
and 1 contributors

# SYNOPSIS

`````` use Math::PlanePath::ArchimedeanChords;
my \$path = Math::PlanePath::ArchimedeanChords->new;
my (\$x, \$y) = \$path->n_to_xy (123);``````

# DESCRIPTION

This path puts points at unit chord steps along an Archimedean spiral. The spiral goes outwards by 1 unit each revolution and the points are spaced 1 apart.

``    R = theta/(2*pi)``

The result is roughly

``````                         31
32          30         ...                3
33                   29
14
34       15          13          28    50             2
12
16        3
35                       2             27    49          1
4                11
17
36           5        0     1          26    48     <- Y=0
10
18
37              6                      25    47         -1
9
19        7     8          24    46
38                                                   -2
20                23
39          21    22             45
-3
40                   44
41    42    43

^
-3    -2    -1    X=0    1     2     3     4``````

X,Y positions returned are fractional. Each revolution is about 2*pi longer than the previous, so the effect is a kind of 6.28 increment looping.

Because the spacing is by unit chords, adjacent unit circles centred on each N position touch but don't overlap. The spiral spacing of 1 unit per revolution means they don't overlap radially either.

The unit chords here are a little like the `TheodorusSpiral`. But the `TheodorusSpiral` goes by unit steps at a fixed right-angle and approximates an Archimedean spiral (of 3.14 radial spacing). Whereas this `ArchimedeanChords` is an actual Archimedean spiral (of radial spacing 1), with unit steps angling along that.

# FUNCTIONS

See "FUNCTIONS" in Math::PlanePath for behaviour common to all path classes.

`\$path = Math::PlanePath::ArchimedeanChords->new ()`

Create and return a new path object.

`(\$x,\$y) = \$path->n_to_xy (\$n)`

Return the X,Y coordinates of point number `\$n` on the path.

`\$n` can be any value `\$n >= 0` and fractions give positions on the chord between the integer points (ie. straight line between the points). `\$n==0` is the origin 0,0.

For `\$n < 0` the return is an empty list, it being considered there are no negative points in the spiral.

`\$n = \$path->xy_to_n (\$x,\$y)`

Return an integer point number for coordinates `\$x,\$y`. Each integer N is considered the centre of a circle of diameter 1 and an `\$x,\$y` within that circle returns N.

The unit spacing of the spiral means those circles don't overlap, but they also don't cover the plane and if `\$x,\$y` is not within one then the return is `undef`.

The current implementation is a bit slow.

`\$n = \$path->n_start ()`

Return 0, the first `\$n` on the path.

`\$str = \$path->figure ()`

Return "circle".

# FORMULAS

## N to X,Y

The current code keeps a position as a polar angle t and calculates an increment u needed to move along by a unit chord. If dist(u) is the straight-line distance between t and t+u, then squared is the hypotenuse

``````    dist^2(u) =   ((t+u)/2pi*cos(t+u) - t/2pi*cos(t))^2     # X
+ ((t+u)/2pi*sin(t+u) - t/2pi*sin(t))^2     # Y``````

which simplifies to

``    dist^2(u) = [ (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) ] / (4*pi^2)``

Switch from cos to sin using the half angle cos(u) = 1 - 2*sin^2(u/2) in case if u is small then the cos(u) near 1.0 might lose floating point accuracy, and also as a slight simplification,

``    dist^2(u) = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2)``

Then want the u which has dist(u)=1 for a unit chord. The u*sin(u) part probably doesn't have a good closed form inverse, so the current code is a Newton/Raphson iteration on f(u) = dist^2(u)-1, seeking f(u)=0

``    f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) - 4*pi^2``

Derivative f'(u) for the slope from the cos form is

``    f'(u) = 2*(t+u) - 2*t*[ cos(u) - (t+u)*sin(u) ]``

And again switching from cos to sin in case u is small,

``    f'(u) = 2*[ u + t*[2*sin^2(u/2) + (t+u)*sin(u)] ]``

## X,Y to N

A given x,y point is at a fraction of a revolution

``````    frac = atan2(y,x) / 2pi     # -.5 <= frac <= .5
frac += (frac < 0)          # 0 <= frac < 1``````

And the nearest spiral arm measured radially from x,y is then

``    r = int(hypot(x,y) + .5 - frac) + frac``

Perl's `atan2` is the same as the C library and gives -pi <= angle <= pi, hence allowing for frac<0. It may also be "unspecified" for x=0,y=0, and give +/-pi for x=negzero, which has to be a special case so 0,0 gives r=0. The `int` rounds towards zero, so frac>.5 ends up as r=0.

So the N point just before or after that spiral position may cover the x,y, but how many N chords it takes to get around to there is 's not so easily calculated.

The current code looks in saved `n_to_xy()` positions for an N below the target, and searches up from there until past the target and thus not covering x,y. With `n_to_xy()` points saved 500 apart this means searching somewhere between 1 and 500 points.

One possibility for calculating a lower bound for N, instead of the saved positions, and both for `xy_to_n()` and `rect_to_n_range()`, would be to add up chords in circles. A circle of radius k fits pi/arcsin(1/2k) many unit chords, so

``````             k=floor(r)     pi
total = sum         ------------
k=0        arcsin(1/2k)``````

and this is less than the chords along the spiral. Is there a good polynomial over-estimate of arcsin, to become an under-estimate total, without giving away so much?

## Rectangle to N Range

For the `rect_to_n_range()` upper bound, the current code takes the arc length along with spiral with the usual formula

``    arc = 1/4pi * (theta*sqrt(1+theta^2) + asinh(theta))``

Written in terms of the r radius (theta = 2pi*r) as calculated from the biggest of the rectangle x,y corners,

``    arc = pi*r^2*sqrt(1+1/r^2) + asinh(2pi*r)/4pi``

The arc length is longer than chords, so N=ceil(arc) is an upper bound for the N range.

An upper bound can also be calculated simply from the circumferences of circles 1 to r, since a spiral loop from radius k to k+1 is shorter than a circle of radius k.

``````             k=ceil(r)
total = sum         2pi*k
k=1

= pi*r*(r+1)``````

This is bigger than the arc length, thus a poorer upper bound, but an easier calculation. (Incidentally, for smallish r have arc length <= pi*(r^2+1) which is a tighter bound and an easy calculation, but it only holds up to somewhere around r=10^7.)

http://user42.tuxfamily.org/math-planepath/index.html