1Math::PlanePath::ArchimUesdeeranCCohnotrrdisb(u3t)ed PerMlatDho:c:uPmleannteaPtaitohn::ArchimedeanChords(3)
2
3
4

NAME

6       Math::PlanePath::ArchimedeanChords -- radial spiral chords
7

SYNOPSIS

9        use Math::PlanePath::ArchimedeanChords;
10        my $path = Math::PlanePath::ArchimedeanChords->new;
11        my ($x, $y) = $path->n_to_xy (123);
12

DESCRIPTION

14       This path puts points at unit chord steps along an Archimedean spiral.
15       The spiral goes outwards by 1 unit each revolution and the points are
16       spaced 1 apart.
17
18           R = theta/(2*pi)
19
20       The result is roughly
21
22                                31
23                          32          30         ...                3
24                    33                   29
25                             14
26              34       15          13          28    50             2
27                                         12
28                    16        3
29           35                       2             27    49          1
30                           4                11
31                 17
32           36           5        0     1          26    48     <- Y=0
33                                            10
34                 18
35           37              6                      25    47         -1
36                                          9
37                    19        7     8          24    46
38              38                                                   -2
39                       20                23
40                 39          21    22             45
41                                                                   -3
42                       40                   44
43                          41    42    43
44
45
46                                 ^
47              -3    -2    -1    X=0    1     2     3     4
48
49       X,Y positions returned are fractional.  Each revolution is about 2*pi
50       longer than the previous, so the effect is a kind of 6.28 increment
51       looping.
52
53       Because the spacing is by unit chords, adjacent unit circles centred on
54       each N position touch but don't overlap.  The spiral spacing of 1 unit
55       per revolution means they don't overlap radially either.
56
57       The unit chords here are a little like the "TheodorusSpiral".  But the
58       "TheodorusSpiral" goes by unit steps at a fixed right-angle and
59       approximates an Archimedean spiral (of 3.14 radial spacing).  Whereas
60       this "ArchimedeanChords" is an actual Archimedean spiral (of radial
61       spacing 1), with unit steps angling along that.
62

FUNCTIONS

64       See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
65       classes.
66
67       "$path = Math::PlanePath::ArchimedeanChords->new ()"
68           Create and return a new path object.
69
70       "($x,$y) = $path->n_to_xy ($n)"
71           Return the X,Y coordinates of point number $n on the path.
72
73           $n can be any value "$n >= 0" and fractions give positions on the
74           chord between the integer points (ie. straight line between the
75           points).  "$n==0" is the origin 0,0.
76
77           For "$n < 0" the return is an empty list, it being considered there
78           are no negative points in the spiral.
79
80       "$n = $path->xy_to_n ($x,$y)"
81           Return an integer point number for coordinates "$x,$y".  Each
82           integer N is considered the centre of a circle of diameter 1 and an
83           "$x,$y" within that circle returns N.
84
85           The unit spacing of the spiral means those circles don't overlap,
86           but they also don't cover the plane and if "$x,$y" is not within
87           one then the return is "undef".
88
89           The current implementation is a bit slow.
90
91       "$n = $path->n_start ()"
92           Return 0, the first $n on the path.
93
94       "$str = $path->figure ()"
95           Return "circle".
96

FORMULAS

98   N to X,Y
99       The current code keeps a position as a polar angle t and calculates an
100       increment u needed to move along by a unit chord.  If dist(u) is the
101       straight-line distance between t and t+u, then squared is the
102       hypotenuse
103
104           dist^2(u) =   ((t+u)/2pi*cos(t+u) - t/2pi*cos(t))^2     # X
105                       + ((t+u)/2pi*sin(t+u) - t/2pi*sin(t))^2     # Y
106
107       which simplifies to
108
109           dist^2(u) = [ (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) ] / (4*pi^2)
110
111       Switch from cos to sin using the half angle cos(u) = 1 - 2*sin^2(u/2)
112       in case if u is small then the cos(u) near 1.0 might lose floating
113       point accuracy, and also as a slight simplification,
114
115           dist^2(u) = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2)
116
117       Then want the u which has dist(u)=1 for a unit chord.  The u*sin(u)
118       part probably doesn't have a good closed form inverse, so the current
119       code is a Newton/Raphson iteration on f(u) = dist^2(u)-1, seeking
120       f(u)=0
121
122           f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) - 4*pi^2
123
124       Derivative f'(u) for the slope from the cos form is
125
126           f'(u) = 2*(t+u) - 2*t*[ cos(u) - (t+u)*sin(u) ]
127
128       And again switching from cos to sin in case u is small,
129
130           f'(u) = 2*[ u + t*[2*sin^2(u/2) + (t+u)*sin(u)] ]
131
132   X,Y to N
133       A given x,y point is at a fraction of a revolution
134
135           frac = atan2(y,x) / 2pi     # -.5 <= frac <= .5
136           frac += (frac < 0)          # 0 <= frac < 1
137
138       And the nearest spiral arm measured radially from x,y is then
139
140           r = int(hypot(x,y) + .5 - frac) + frac
141
142       Perl's "atan2" is the same as the C library and gives -pi <= angle <=
143       pi, hence allowing for frac<0.  It may also be "unspecified" for
144       x=0,y=0, and give +/-pi for x=negzero, which has to be a special case
145       so 0,0 gives r=0.  The "int" rounds towards zero, so frac>.5 ends up as
146       r=0.
147
148       So the N point just before or after that spiral position may cover the
149       x,y, but how many N chords it takes to get around to there is 's not so
150       easily calculated.
151
152       The current code looks in saved "n_to_xy()" positions for an N below
153       the target, and searches up from there until past the target and thus
154       not covering x,y.  With "n_to_xy()" points saved 500 apart this means
155       searching somewhere between 1 and 500 points.
156
157       One possibility for calculating a lower bound for N, instead of the
158       saved positions, and both for "xy_to_n()" and "rect_to_n_range()",
159       would be to add up chords in circles.  A circle of radius k fits
160       pi/arcsin(1/2k) many unit chords, so
161
162                    k=floor(r)     pi
163           total = sum         ------------
164                    k=0        arcsin(1/2k)
165
166       and this is less than the chords along the spiral.  Is there a good
167       polynomial over-estimate of arcsin, to become an under-estimate total,
168       without giving away so much?
169
170   Rectangle to N Range
171       For the "rect_to_n_range()" upper bound, the current code takes the arc
172       length along with spiral with the usual formula
173
174           arc = 1/4pi * (theta*sqrt(1+theta^2) + asinh(theta))
175
176       Written in terms of the r radius (theta = 2pi*r) as calculated from the
177       biggest of the rectangle x,y corners,
178
179           arc = pi*r^2*sqrt(1+1/r^2) + asinh(2pi*r)/4pi
180
181       The arc length is longer than chords, so N=ceil(arc) is an upper bound
182       for the N range.
183
184       An upper bound can also be calculated simply from the circumferences of
185       circles 1 to r, since a spiral loop from radius k to k+1 is shorter
186       than a circle of radius k.
187
188                    k=ceil(r)
189           total = sum         2pi*k
190                    k=1
191
192                 = pi*r*(r+1)
193
194       This is bigger than the arc length, thus a poorer upper bound, but an
195       easier calculation.  (Incidentally, for smallish r have arc length <=
196       pi*(r^2+1) which is a tighter bound and an easy calculation, but it
197       only holds up to somewhere around r=10^7.)
198

SEE ALSO

200       Math::PlanePath, Math::PlanePath::TheodorusSpiral,
201       Math::PlanePath::SacksSpiral
202

HOME PAGE

204       <http://user42.tuxfamily.org/math-planepath/index.html>
205

LICENSE

207       Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019,
208       2020 Kevin Ryde
209
210       This file is part of Math-PlanePath.
211
212       Math-PlanePath is free software; you can redistribute it and/or modify
213       it under the terms of the GNU General Public License as published by
214       the Free Software Foundation; either version 3, or (at your option) any
215       later version.
216
217       Math-PlanePath is distributed in the hope that it will be useful, but
218       WITHOUT ANY WARRANTY; without even the implied warranty of
219       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
220       General Public License for more details.
221
222       You should have received a copy of the GNU General Public License along
223       with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.
224
225
226
227perl v5.34.0                      2022-01-M2a1th::PlanePath::ArchimedeanChords(3)
Impressum