File Coverage

blib/lib/Math/PlanePath/ArchimedeanChords.pm
Criterion Covered Total %
statement 130 154 84.4
branch 14 24 58.3
condition n/a
subroutine 31 32 96.8
pod 4 4 100.0
total 179 214 83.6


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Kevin Ryde
2              
3             # This file is part of Math-PlanePath.
4             #
5             # Math-PlanePath is free software; you can redistribute it and/or modify it
6             # under the terms of the GNU General Public License as published by the Free
7             # Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-PlanePath is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-PlanePath. If not, see .
17              
18              
19             # Also possible would be circle involute spiral, unrolling string around
20             # centre of circumference 1, but is only very slightly different radius from
21             # an Archimedean spiral.
22              
23              
24             package Math::PlanePath::ArchimedeanChords;
25 1     1   7333 use 5.004;
  1         8  
26 1     1   5 use strict;
  1         1  
  1         30  
27 1     1   354 use Math::Libm 'hypot', 'asinh';
  1         4276  
  1         71  
28 1     1   445 use POSIX 'ceil';
  1         6489  
  1         5  
29             #use List::Util 'min', 'max';
30             *min = \&Math::PlanePath::_min;
31             *max = \&Math::PlanePath::_max;
32              
33 1     1   1145 use vars '$VERSION', '@ISA';
  1         2  
  1         88  
34             $VERSION = 128;
35 1     1   556 use Math::PlanePath;
  1         2  
  1         33  
36             @ISA = ('Math::PlanePath');
37              
38             use Math::PlanePath::Base::Generic
39 1         33 'is_infinite',
40 1     1   6 'round_nearest';
  1         1  
41 1     1   471 use Math::PlanePath::MultipleRings;
  1         2  
  1         28  
42              
43             # uncomment this to run the ### lines
44             # use Smart::Comments;
45              
46              
47 1     1   5 use constant figure => 'circle';
  1         3  
  1         43  
48 1     1   5 use constant n_start => 0;
  1         2  
  1         34  
49 1     1   5 use constant x_negative_at_n => 3;
  1         1  
  1         34  
50 1     1   5 use constant y_negative_at_n => 5;
  1         1  
  1         34  
51 1     1   5 use constant gcdxy_maximum => 1;
  1         1  
  1         36  
52 1     1   5 use constant dx_minimum => -1; # infimum when straight
  1         1  
  1         33  
53 1     1   4 use constant dx_maximum => 1; # at N=0
  1         2  
  1         47  
54 1     1   5 use constant dy_minimum => -1;
  1         1  
  1         40  
55 1     1   5 use constant dy_maximum => 1;
  1         2  
  1         46  
56 1     1   6 use constant dsumxy_minimum => -sqrt(2); # supremum when diagonal
  1         1  
  1         45  
57 1     1   5 use constant dsumxy_maximum => sqrt(2);
  1         2  
  1         36  
58 1     1   4 use constant ddiffxy_minimum => -sqrt(2); # supremum when diagonal
  1         1  
  1         46  
59 1     1   5 use constant ddiffxy_maximum => sqrt(2);
  1         1  
  1         38  
60 1     1   5 use constant turn_any_right => 0; # left always
  1         2  
  1         33  
61 1     1   4 use constant turn_any_straight => 0; # left always
  1         1  
  1         61  
62              
63              
64             #------------------------------------------------------------------------------
65              
66 1     1   5 use constant 1.02 _PI => 2*atan2(1,0);
  1         11  
  1         173  
67              
68             # Starting at polar angle position t in radians,
69             #
70             # r = t / 2pi
71             #
72             # x = r * cos(t) = t * cos(t) / 2pi
73             # y = r * sin(t) = t * sin(t) / 2pi
74             #
75             # Want a polar angle amount u to move by a chord distance of 1. Hypot
76             # square distance from t to t+u is
77             #
78             # dist(u) = ( (t+u)/2pi*cos(t+u) - t/2pi*cos(t) )^2 # X
79             # + ( (t+u)/2pi*sin(t+u) - t/2pi*sin(t) )^2 # Y
80             # = [ (t+u)^2*cos^2(t+u) - 2*(t+u)*t*cos(t+u)*cos(t) + t^2*cos^2(t)
81             # + (t+u)^2*sin^2(t+u) - 2*(t+u)*t*sin(t+u)*sin(t) + t^2*sin^2(t)
82             # ] / (4*pi^2)
83             #
84             # and from sin^2 + cos^2 = 1
85             # and addition cosA*cosB + sinA*sinB = cos(A-B)
86             #
87             # = [ (t+u)^2 - 2*(t+u)*t*cos((t+u)-t) + t^2 ] /4pi^2
88             # = [ (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) ] / (4*pi^2)
89             #
90             # then double angle cos(u) = 1 - 2*sin^2(u/2) to go to the sine since if u
91             # is small then cos(u) near 1.0 might lose accuracy
92             #
93             # dist(u) = [(t+u)^2 + t^2 - 2*t*(t+u)*(1 - 2*sin^2(u/2))] / (4*pi^2)
94             # = [(t+u)^2 + t^2 - 2*t*(t+u) + 2*t*(t+u)*2*sin^2(u/2)] / (4*pi^2)
95             # = [((t+u) - t)^2 + 4*t*(t+u)*sin^2(u/2)] / (4*pi^2)
96             # = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2)
97             #
98             # Seek d(u) = 1 by letting f(u)=4*pi^2*(d(u)-1) and seeking f(u)=0
99             #
100             # f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) - 4*pi^2
101             #
102             # Derivative f'(u) for the slope, starting from the cosine form,
103             #
104             # f(u) = (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) - 4*pi^2
105             #
106             # f'(u) = 2*(t+u) - 2*t*[ cos(u) - (t+u)*sin(u) ]
107             # = 2*(t+u) - 2*t*[ 1 - 2*sin^2(u/2) - (t+u)*sin(u) ]
108             # = 2*t + 2*u - 2*t + 2*t*2*sin^2(u/2) + 2*t*(t+u)*sin(u)
109             # = 2*[ u + 2*t*sin^2(u/2) + t*(t+u)*sin(u) ]
110             # = 2*[ u + t * [2*sin^2(u/2) + (t+u)*sin(u) ] ]
111             #
112             # Newton's method
113             # */ <- f(x) high
114             # */|
115             # * / |
116             # * / |
117             # ---------*------------------
118             # +---+ <- subtract
119             #
120             # f(x) / sub = f'(x)
121             # sub = f(x) / f'(x)
122             #
123             #
124             # _chord_angle_inc() takes $t is a polar angle around the Archimedean spiral.
125             # Returns an increment polar angle $u which may be added to $t to move around
126             # the spiral by a chord length 1 unit.
127             #
128             # The loop is Newton's method, $f=f(u), $slope=f'(u) so $u-$f/$slope is a
129             # better $u, ie. f($u) closer to 0. Stop when the subtract becomes small,
130             # usually only about 3 iterations.
131             #
132             sub _chord_angle_inc {
133 599     599   765 my ($t) = @_;
134             # ### _chord_angle_inc(): $t
135              
136 599         691 my $u = 2*_PI/$t; # estimate
137              
138 599         732 foreach (0 .. 10) {
139 2097         2497 my $shu = sin($u/2); $shu *= $shu; # sin^2(u/2)
  2097         2223  
140 2097         2585 my $tu = ($t+$u);
141              
142 2097         2694 my $f = $u*$u + 4*$t*$tu*$shu - 4*_PI*_PI;
143 2097         2826 my $slope = 2 * ( $u + $t*(2*$shu + $tu*sin($u)));
144              
145             # unsimplified versions ...
146             # $f = ($t+$u)**2 + $t**2 - 2*$t*($t+$u)*cos($u) - 4*_PI*_PI;
147             # $slope = 2*($t+$u) - 2*$t*( cos($u) - ($t+$u)*sin($u) );
148              
149 2097         2286 my $sub = $f/$slope;
150 2097         2210 $u -= $sub;
151              
152             # printf ("f=%.6f slope=%.6f sub=%.20f u=%.6f\n", $f, $slope, $sub, $u);
153 2097 100       3318 last if (abs($sub) < 1e-15);
154             }
155             # printf ("return u=%.6f\n", $u);
156 599         793 return $u;
157             }
158              
159 1     1   6 use constant 1.02; # for leading underscore
  1         25  
  1         22  
160 1     1   4 use constant _SAVE => 500;
  1         2  
  1         749  
161              
162             my @save_n = (1);
163             my @save_t = (2*_PI);
164             my $next_save = $save_n[0] + _SAVE;
165              
166             sub new {
167             ### ArchimedeanChords new() ...
168 3     3 1 993 return shift->SUPER::new (i => $save_n[0],
169             t => $save_t[0],
170             @_);
171             }
172              
173             sub n_to_xy {
174 13     13 1 26 my ($self, $n) = @_;
175              
176 13 50       26 if ($n < 0) { return; }
  0         0  
177 13 50       110 if (is_infinite($n)) { return ($n,$n); }
  0         0  
178              
179 13 100       28 if ($n <= 1) {
180 12         24 return ($n, 0); # exactly Y=0
181             }
182              
183             {
184             # ENHANCE-ME: look at the N+1 position for the frac directly, without
185             # the full call for N+1
186 1         2 my $int = int($n);
  1         3  
187 1 50       3 if ($n != $int) {
188 0         0 my $frac = $n - $int; # inherit possible BigFloat/BigRat
189 0         0 my ($x1,$y1) = $self->n_to_xy($int);
190 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
191 0         0 my $dx = $x2-$x1;
192 0         0 my $dy = $y2-$y1;
193 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
194             }
195             }
196              
197 1         5 my $i = $self->{'i'};
198 1         2 my $t = $self->{'t'};
199              
200 1 50       5 if ($i > $n) {
201 0         0 my $pos = min ($#save_n, int (($n - $save_n[0]) / _SAVE));
202 0         0 $i = $save_n[$pos];
203 0         0 $t = $save_t[$pos];
204             ### resume: "$i t=$t"
205             }
206              
207 1         3 while ($i < $n) {
208 599         779 $t += _chord_angle_inc($t);
209 599 100       1051 if (++$i == $next_save) {
210 1         6 push @save_n, $i;
211 1         3 push @save_t, $t;
212 1         3 $next_save += _SAVE;
213             }
214             }
215              
216 1         4 $self->{'i'} = $i;
217 1         2 $self->{'t'} = $t;
218              
219 1         3 my $r = $t * (1 / (2*_PI));
220 1         7 return ($r*cos($t),
221             $r*sin($t));
222             }
223              
224             sub _xy_to_nearest_r {
225 56     56   2116 my ($x, $y) = @_;
226 56         118 my $frac = Math::PlanePath::MultipleRings::_xy_to_angle_frac($x,$y);
227             ### assert: 0 <= $frac && $frac < 1
228              
229             # if $frac > 0.5 then 0.5-$frac is negative and int() rounds towards zero
230             # giving $r==$frac
231 56         174 return int(hypot($x,$y) + 0.5 - $frac) + $frac;
232             }
233              
234             sub xy_to_n {
235 12     12 1 212 my ($self, $x, $y) = @_;
236             ### ArchimedeanChords xy_to_n(): "$x, $y"
237              
238 12         23 my $r = _xy_to_nearest_r($x,$y);
239 12         18 my $r_limit = 1.001 * $r;
240              
241             ### hypot: hypot($x,$y)
242             ### $r
243             ### $r_limit
244             ### save_t: "end index=$#save_t save_t[0]=".($save_t[0]//'undef')
245              
246 12 50       29 if (is_infinite($r_limit)) {
247             ### infinite range, r inf or too big ...
248 0         0 return undef;
249             }
250              
251 12         18 my $theta = 0.999 * 2*_PI*$r;
252 12         15 my $n_lo = 0;
253 12         23 foreach my $i (1 .. $#save_t) {
254 12 50       94 if ($save_t[$i] > $theta) {
255 12         21 $n_lo = $save_n[$i-1];
256 12 50       47 if ($n_lo == 1) { $n_lo = 0; } # for finding X=0,Y=0
  12         16  
257 12         18 last;
258             }
259             }
260             ### $n_lo
261              
262             # loop with for(;;) since $n_lo..$n_hi limited to IV range
263 12         16 for (my $n = $n_lo; ; $n += 1) {
264 12         24 my ($nx,$ny) = $self->n_to_xy($n);
265             # #### $n
266             # #### $nx
267             # #### $ny
268             # #### hypot: hypot ($x-$nx,$y-$ny)
269 12 50       37 if (hypot($x-$nx,$y-$ny) <= 0.5) {
270             ### hypot in range ...
271 12         55 return $n;
272             }
273 0 0         if (hypot($nx,$ny) >= $r_limit) {
274 0           last;
275             }
276             }
277              
278             ### n not found ...
279 0           return undef;
280             }
281              
282             # int (max (0, int(_PI*$r2) - 4*$r));
283             #
284             # my $r2 = $r * $r;
285             # my $n_lo = int (max (0, int(_PI*$r2) - 4*$r));
286             # my $n_hi = $n_lo + 7*$r + 2;
287             # ### $r2
288             # $n_lo == $n_lo-1 ||
289              
290              
291              
292              
293             # x,y has radius hypot(x,y), then want the next higher spiral arc which is r
294             # >= hypot(x,y)+0.5, with the 0.5 being the width of the circle figure on
295             # the spiral.
296             #
297             # The polar angle of x,y is a=atan2(y,x) and frac=a/2pi is the extra away
298             # from an integer radius for the spiral. So seek integer k with k+a/2pi >=
299             # h with h=hypot(x,y)+0.5.
300             #
301             # k + a/2pi >= h
302             # k >= h-a/2pi
303             # k = ceil(h-a/2pi)
304             # = ceil(hypot(x,y) + 0.5 - atan2(y,x)/2pi)
305             #
306             #
307             # circle radius i has circumference 2*pi*i and at most that many N on it
308             # rectangle corner at radius Rcorner = hypot(x,y)
309             #
310             # sum i=1 to i=Rlimit of 2*pi*i = 2*pi/2 * Rlimit*(Rlimit+1)
311             # = pi * Rlimit*(Rlimit+1)
312             # is an upper bound, though a fairly slack one
313             #
314             #
315             # cf. arc length along the spiral r=a*theta with a=1/2pi
316             # arclength = (1/2) * a * (theta*sqrt(1+theta^2) + asinh(theta))
317             # = (1/4*pi) * (theta*sqrt(1+theta^2) + asinh(theta))
318             # and theta = 2*pi*r
319             # = (1/4*pi) * (4*pi^2*r^2 * sqrt(1+1/theta^2) + asinh(theta))
320             # = pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
321             #
322             # and to compare to the circles formula
323             #
324             # = pi * (r*(r+1) * r/(r+1) * sqrt(1+1/r^2)
325             # + asinh(theta)/(4*pi^2))
326             #
327             # so it's smaller hence better upper bound. Only a little smaller than the
328             # squaring once get to 100 loops or so.
329             #
330             #
331             # not exact
332             sub rect_to_n_range {
333 0     0 1   my ($self, $x1,$y1, $x2,$y2) = @_;
334             ### rect_to_n_range() ...
335              
336 0           my $rhi = 0;
337 0           foreach my $x ($x1, $x2) {
338 0           foreach my $y ($y1, $y2) {
339 0           my $frac = atan2($y,$x) / (2*_PI); # -0.5 <= $frac <= 0.5
340 0           $frac += ($frac < 0); # 0 <= $frac < 1
341 0           $rhi = max ($rhi, ceil(hypot($x,$y)+0.5 - $frac) + $frac);
342             }
343             }
344             ### $rhi
345              
346             # arc length pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
347             # = pi*r^2*sqrt(1+1/r^2) + asinh(theta)/4pi
348 0           my $rhi2 = $rhi*$rhi;
349 0           return (0,
350             ceil (_PI * $rhi2 * sqrt(1+1/$rhi2)
351             + asinh(2*_PI*$rhi) / (4*_PI)));
352              
353              
354             # # each loop begins at N = pi*k^2 - 2 or thereabouts
355             # return (0,
356             # int(_PI*$rhi*($rhi+1) + 1));
357             }
358              
359             1;
360             __END__