Math::PlanePath::TheodorusSpiral -- right-angle unit step spiral
use Math::PlanePath::TheodorusSpiral; my $path = Math::PlanePath::TheodorusSpiral->new; my ($x, $y) = $path->n_to_xy (123);
This path puts points on the spiral of Theodorus, also called the square root spiral.
61 6 60 27 26 25 24 5 28 23 59 29 22 58 4 30 21 57 3 31 20 4 56 2 32 5 3 19 6 2 55 1 33 18 7 0 1 54 <- Y=0 34 17 8 53 -1 35 16 9 52 -2 36 15 10 14 51 -3 37 11 12 13 50 -4 38 49 39 48 -5 40 47 41 46 -6 42 43 44 45 ^ -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7
Each step is a unit distance at right angles to the previous radial spoke. So for example,
3 -- Y=1+1/sqrt(2) \ \ ..2 Y=1 .. | . | 0-----1 <- Y=0 ^ X=0 X=1
1 to 2 is a unit step at right angles to the 0 to 1 radial. Then 2 to 3 steps at a right angle to radial 0 to 2 (which is 45 degrees), etc. The distance 0 to 2 is sqrt(2), the distance 0 to 3 is sqrt(3), and in general r(N) = sqrt(N) since each step is a right triangle with radius(N+1)^2 = radius(N)^2 + 1. The resulting shape is very close to an Archimedean spiral with successive loops increasing in radius by pi = 3.14159 or thereabouts each time.
X,Y positions returned are fractional and each integer N position is exactly 1 away from the previous. Fractional N values give positions on the straight line between the integer points. (An analytic continuation for a rounded curve between points is possible, but not currently implemented.)
Each loop is just under 2*pi^2 = 19.7392 many N points longer than the previous. This means quadratic values 9.8696*k^2 for integer k are an almost straight line. Quadratics close to 9.87 (or a square multiple of that) nearly line up. For example the 22-polygonal numbers have 10*k^2 and at low values are nearly straight, but then spiral away.
See "FUNCTIONS" in Math::PlanePath for the behaviour common to all path classes.
The code is currently implemented by adding unit steps in X,Y coordinates, so it's not particularly fast. The last X,Y is saved in the object anticipating successively higher N (not necessarily consecutive), and previous positions 1000 apart are saved for re-use or to go back to.
$path = Math::PlanePath::TheodorusSpiral->new ()
Create and return a new theodorus spiral object.
($x,$y) = $path->n_to_xy ($n)
Return the X,Y coordinates of point number $n on the path.
$n
$n can be any value $n >= 0 and fractions give positions on the spiral in between the integer points.
$n >= 0
For $n < 0 the return is an empty list, it being currently considered there are no negative points in the spiral. (The analytic continuation by Davis would be a possibility, though the resulting "inner spiral" makes positive and negative points overlap a bit. A spiral starting at X=-1 would fit in between the positive points.)
$n < 0
$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.
$x,$y
The unit steps of the spiral means those unit circles don't overlap, but the loops are 3.14 apart so there's gaps in between. If $x,$y is not within one of the unit circles then the return is undef.
undef
$n = $path->figure ()
Return "circle".
For a given X,Y the radius r=hypot(X,Y) determines the N position as N=r^2. An N point up to 0.5 away radially might above cover X,Y, so the range of N to consider is
Nlo = (r-.5)^2 Nhi = (r+.5)^2
A simple search through those N's checking for which, if any, covers X,Y is then done. The number of N's there is Nhi-Nlo = 2*r+1 which is about 1/3 of a loop around the spiral (2*r/2*pi*r ~= 1/3). Actually 0.51 is used to guard against floating point round-off, which is then about 4*.51 = 2.04*r.
The angle of the X,Y position determines which part of the spiral is intersected, but using that doesn't seem particularly easy. The angle for a given N is an arctan sum and don't yet have a good closed-form for that to invert, or some Newton's method, or whatever.
For rect_to_n_range() the corner furthest from the origin determines the high N. For that corner
rect_to_n_range()
rhi=hypot(xhi,yhi) Nhi = (rhi+.5)^2
The extra .5 is since a unit circle figure centred as much as .5 further out might intersect the xhi,yhi. The worst case for this estimate is when Nhi doesn't intersect the xhi,yhi corner but is just before it, counter-clockwise. It's then a full revolution bigger than it need be (depending where the other corners fall).
Similarly for the corner nearest the origin,
rlo = hypot(xlo,ylo) Nlo = (rlo-.5)^2, or 0 if origin covered by rectangle
The worst case is when this Nlo doesn't intersect the xlo,ylo corner but is just after it counter-clockwise, so Nlo is a full revolution smaller than it need be (depending where the other corners fall).
Math::PlanePath, Math::PlanePath::ArchimedeanChords, Math::PlanePath::SacksSpiral, Math::PlanePath::MultipleRings
http://user42.tuxfamily.org/math-planepath/index.html
Copyright 2010, 2011 Kevin Ryde
This file is part of Math-PlanePath.
Math-PlanePath 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 3, or (at your option) any later version.
Math-PlanePath 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 Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
To install Math::PlanePath, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Math::PlanePath
CPAN shell
perl -MCPAN -e shell install Math::PlanePath
For more information on module installation, please visit the detailed CPAN module installation guide.