1Math::PlanePath::ArchimUesdeeranCCohnotrrdisb(u3t)ed PerMlatDho:c:uPmleannteaPtaitohn::ArchimedeanChords(3)
2
6 Math::PlanePath::ArchimedeanChords -- radial spiral chords
7
9 use Math::PlanePath::ArchimedeanChords;
10 my $path = Math::PlanePath::ArchimedeanChords->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
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 R = theta/(2*pi)
19 The result is roughly
21 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 ^
47 -3 -2 -1 X=0 1 2 3 4
48 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 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 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
64 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
65 classes.
66 "$path = Math::PlanePath::ArchimedeanChords->new ()"
68 Create and return a new path object.
69 "($x,$y) = $path->n_to_xy ($n)"
71 Return the X,Y coordinates of point number $n on the path.
72 $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 For "$n < 0" the return is an empty list, it being considered there
78 are no negative points in the spiral.
79 "$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 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 The current implementation is a bit slow.
90 "$n = $path->n_start ()"
92 Return 0, the first $n on the path.
93 "$str = $path->figure ()"
95 Return "circle".
96
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 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 which simplifies to
108 dist^2(u) = [ (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) ] / (4*pi^2)
110 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 dist^2(u) = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2)
116 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 f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) - 4*pi^2
123 Derivative f'(u) for the slope from the cos form is
125 f'(u) = 2*(t+u) - 2*t*[ cos(u) - (t+u)*sin(u) ]
127 And again switching from cos to sin in case u is small,
129 f'(u) = 2*[ u + t*[2*sin^2(u/2) + (t+u)*sin(u)] ]
131 X,Y to N
133 A given x,y point is at a fraction of a revolution
134 frac = atan2(y,x) / 2pi # -.5 <= frac <= .5
136 frac += (frac < 0) # 0 <= frac < 1
137 And the nearest spiral arm measured radially from x,y is then
139 r = int(hypot(x,y) + .5 - frac) + frac
141 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 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 The current code looks in saved n_to_xy() positions for an N below the
153 target, and searches up from there until past the target and thus not
154 covering x,y. With n_to_xy() points saved 500 apart this means
155 searching somewhere between 1 and 500 points.
156 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(), would be
159 to add up chords in circles. A circle of radius k fits pi/arcsin(1/2k)
160 many unit chords, so
161 k=floor(r) pi
163 total = sum ------------
164 k=0 arcsin(1/2k)
165 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 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 arc = 1/4pi * (theta*sqrt(1+theta^2) + asinh(theta))
175 Written in terms of the r radius (theta = 2pi*r) as calculated from the
177 biggest of the rectangle x,y corners,
178 arc = pi*r^2*sqrt(1+1/r^2) + asinh(2pi*r)/4pi
180 The arc length is longer than chords, so N=ceil(arc) is an upper bound
182 for the N range.
183 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 k=ceil(r)
189 total = sum 2pi*k
190 k=1
191 = pi*r*(r+1)
193 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
200 Math::PlanePath, Math::PlanePath::TheodorusSpiral,
201 Math::PlanePath::SacksSpiral
202
204 <http://user42.tuxfamily.org/math-planepath/index.html>
205
207 Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019,
208 2020 Kevin Ryde
209 This file is part of Math-PlanePath.
211 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 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 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/>.
224perl v5.36.0 2023-01-M2a0th::PlanePath::ArchimedeanChords(3)