File Coverage

blib/lib/Math/PlanePath/MultipleRings.pm
Criterion Covered Total %
statement 252 403 62.5
branch 138 280 49.2
condition 27 41 65.8
subroutine 30 52 57.6
pod 27 27 100.0
total 474 803 59.0


line stmt bran cond sub pod time code
1             # Copyright 2010, 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
6             # it under the terms of the GNU General Public License as published by the
7             # Free 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              
20             # math-image --path=MultipleRings --lines
21             #
22             # math-image --wx --path=MultipleRings,ring_shape=polygon,step=5 --scale=50 --figure=ring --all
23              
24             #
25             # FIXME: $y equal across bottom side centre ?
26              
27              
28             package Math::PlanePath::MultipleRings;
29 15     15   2999 use 5.004;
  15         53  
30 15     15   82 use strict;
  15         30  
  15         426  
31 15     15   78 use Carp 'croak';
  15         36  
  15         1159  
32             #use List::Util 'min','max';
33             *min = \&Math::PlanePath::_min;
34             *max = \&Math::PlanePath::_max;
35              
36             # Math::Trig has asin_real() too, but it just runs the blob of code in
37             # Math::Complex -- prefer libm
38 15     15   868 use Math::Libm 'asin', 'hypot';
  15         6183  
  15         864  
39              
40 15     15   93 use vars '$VERSION', '@ISA';
  15         29  
  15         949  
41             @ISA = ('Math::PlanePath');
42 15     15   1636 use Math::PlanePath;
  15         30  
  15         842  
43             *_sqrtint = \&Math::PlanePath::_sqrtint;
44             $VERSION = 128;
45              
46             use Math::PlanePath::Base::Generic
47 15     15   97 'is_infinite';
  15         28  
  15         625  
48 15     15   1629 use Math::PlanePath::SacksSpiral;
  15         35  
  15         553  
49              
50             # uncomment this to run the ### lines
51             # use Smart::Comments;
52              
53              
54 15     15   87 use constant 1.02; # for leading underscore
  15         237  
  15         597  
55 15     15   90 use constant _PI => 2*atan2(1,0);
  15         35  
  15         888  
56              
57 15     15   103 use constant figure => 'circle';
  15         37  
  15         723  
58 15     15   85 use constant n_frac_discontinuity => 0;
  15         27  
  15         723  
59 15     15   97 use constant gcdxy_minimum => 0;
  15         31  
  15         1480  
60              
61 15         20058 use constant parameter_info_array =>
62             [{ name => 'step',
63             display => 'Step',
64             share_key => 'step_6_min3',
65             type => 'integer',
66             minimum => 0,
67             default => 6,
68             width => 3,
69             description => 'How much longer each ring is than the preceding.',
70             },
71              
72             { name => 'ring_shape',
73             display => 'Ring Shape',
74             type => 'enum',
75             default => 'circle',
76             choices => ['circle','polygon'],
77             choices_display => ['Circle','Polygon'],
78             description => 'The shape of each ring, either a circle or a polygon of "step" many sides.',
79             },
80 15     15   101 ];
  15         33  
81              
82             sub turn_any_left {
83 0     0 1 0 my ($self) = @_;
84             # step == 0 is always straight ahead
85 0         0 return ($self->{'step'} != 0);
86             }
87             sub turn_any_right {
88 0     0 1 0 my ($self) = @_;
89             # step=0 is always straight ahead
90             # step=1 is never right
91 0         0 return ($self->{'step'} >= 2);
92             }
93             {
94             my @_UNDOCUMENTED__turn_any_right_at_n
95             = (undef, # 0
96             undef, # 1
97             131, # 2
98             44, # 3
99             23, # 4
100             29, # 5
101             17, # 6
102             20, # 7
103             23); # 8
104             sub _UNDOCUMENTED__turn_any_right_at_n {
105 0     0   0 my ($self) = @_;
106 0 0       0 $self->turn_any_right or return undef;
107 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
108             # step=8 24, 9, 10, 11
109             return $self->n_start - 1 + ($self->{'step'} < 9 ? 3*$self->{'step'}
110 0 0       0 : $self->{'step'});
111             }
112             return $self->n_start
113             + ($self->{'step'} <= $#_UNDOCUMENTED__turn_any_right_at_n
114             ? $_UNDOCUMENTED__turn_any_right_at_n[$self->{'step'}]
115 0 0       0 : $self->{'step'} - 1);
116             }
117             }
118              
119             sub turn_any_straight {
120 0     0 1 0 my ($self) = @_;
121             # step=0 straight line
122             # step=1 straight at N=2
123             # step=2 straight at N=2
124             return ($self->{'step'} <= 2 ? 1
125 0 0       0 : $self->{'ring_shape'} eq 'circle' ? 0 # never straight
    0          
126             : 1); # ring_shape=polygon sides straight
127             }
128              
129              
130             #------------------------------------------------------------------------------
131             # Electricity transmission cable in sixes, with one at centre ?
132             # 7 poppy
133             # 19 hyacinth
134             # 37 marigold
135             # 61 cowslip
136             # 127 bluebonnet
137              
138             # An n-gon of points many vertices has each angle
139             # alpha = 2*pi/points
140             # The radius r to a vertex, using a line perpendicular to the line segment
141             # sin(alpha/2) = (1/2)/r
142             # r = 0.5 / sin(pi/points)
143             # And with points = d*step, starting from d=1
144             # r = 0.5 / sin(pi/(d*step))
145              
146             # step==0 is a straight line y==0 x=0,1,2,..., anything else whole plane
147             sub x_negative {
148 4     4 1 7 my ($self) = @_;
149 4         13 return ($self->{'step'} > 0);
150             }
151             *y_negative = \&x_negative;
152              
153             sub y_maximum {
154 0     0 1 0 my ($self) = @_;
155 0 0       0 return ($self->{'step'} == 0 ? 0 # step=0 always Y=0
156             : undef);
157             }
158              
159             sub x_negative_at_n {
160 0     0 1 0 my ($self) = @_;
161             return ($self->{'step'} == 0 ? undef # no negatives
162             : $self->{'step'} == 1 ? 3
163 0 0       0 : $self->n_start + int($self->{'step'}/4) + 1);
    0          
164             }
165             sub y_negative_at_n {
166 0     0 1 0 my ($self) = @_;
167             return ($self->{'step'} == 0 ? undef # no negatives
168             : $self->{'step'} <= 2 ? 6
169 0 0       0 : $self->n_start + int($self->{'step'}/2) + 1);
    0          
170             }
171              
172             sub sumxy_minimum {
173 0     0 1 0 my ($self) = @_;
174 0 0       0 return ($self->{'step'} == 0 ? 0 : undef);
175             }
176             sub sumabsxy_minimum {
177 0     0 1 0 my ($self) = @_;
178             # first point N=1 innermost ring
179 0         0 my ($x,$y) = $self->n_to_xy($self->n_start);
180 0         0 return $x;
181             }
182             *diffxy_minimum = \&sumxy_minimum;
183              
184             # step=0 X=0,Y=0 AbsDiff=0
185             # step=3 N=88 X=Y=5.3579957587697 ring of 24 is a multiple of 8
186              
187             sub rsquared_minimum {
188 0     0 1 0 my ($self) = @_;
189 0         0 my $step = $self->{'step'};
190 0 0       0 if ($step <= 1) {
191             # step=0 along X axis starting X=0,Y=0
192             # step=1 start at origin
193 0         0 return 0;
194             }
195              
196             # step=3 *--___
197             # circle | --__ o 0.5/r = sin60 = sqrt(3)/2
198             # | o __* / | \ r = 1/sqrt(3)
199             # | ___-- / | \ r^2 = 1/3
200             # *-- *---------*
201             # 1/2
202             # polygon
203             # o 0.5/r = sin60 = sqrt(3)/2
204             # / | \ r = 1/sqrt(3)
205             # / | \ r^2 = 1/3
206             # *---------*
207             # 1/2
208             #
209 0 0       0 if ($step == 3) {
210 0 0       0 return ($self->{'ring_shape'} eq 'polygon' ? 3/4 : 1/3);
211             }
212 0 0       0 if ($step == 4) {
213             # radius = sqrt(2)/2, rsquared=1/2
214 0         0 return 0.5;
215             }
216              
217             # _numsides_to_r() returns 1, no need for a special case here
218             # if ($step == 6) {
219             # # hexagon
220             # return 1;
221             # }
222              
223 0         0 my $r;
224 0 0 0     0 if ($step >= 6 || $self->{'ring_shape'} eq 'polygon') {
225 0         0 $r = _numsides_to_r($step,_PI);
226             } else {
227 0         0 $r = $self->{'base_r'} + 1;
228             }
229 0         0 return $r*$r;
230             }
231              
232              
233             #------------------------------------------------------------------------------
234             # dx_minimum() etc
235              
236             # step <= 6
237             # R=base_r+d
238             # theta = 2*$n * $pi / ($d * $step)
239             # = 2pi/(d*step)
240             # dX -> R*sin(theta)
241             # -> R*theta
242             # = (base_r+d)*2pi/(d*step)
243             # -> 2pi/step
244             #
245             # step=5 across first ring
246             # N=6 at X=base_r+2, Y=0
247             # N=5 at R=base_r+1 theta = 2pi/5
248             # X=(base_r+1)*cos(theta)
249             # dX = base_r+2 - (base_r+1)*cos(theta)
250             #
251             # step=6 across first ring
252             # base_r = 0.5/sin(_PI/6) - 1
253             # = 0.5/0.5 - 1
254             # = 0
255             # N=7 at X=base_r+2, Y=0
256             # N=6 at R=base_r+1 theta = 2pi/6
257             # X=(base_r+1)*cos(theta)
258             # dX = base_r+2 - (base_r+1)*cos(theta)
259             # = base_r+2 - (base_r+1)*0.5
260             # = 1.5*base_r + 1.5
261             # = 1.5
262             #
263             # step > 6
264             # R = 0.5 / sin($pi / ($d*$step))
265             # diff = 0.5 / sin($pi / ($d*$step)) - 0.5 / sin($pi / (($d-1)*$step))
266             # -> 0.5 / ($pi / ($d*$step)) - 0.5 / ($pi / (($d-1)*$step))
267             # = 0.5 * ($d*$step) / $pi - 0.5 * (($d-1)*$step) / $pi
268             # = step*0.5/pi * ($d - ($d-1))
269             # = step*0.5/pi
270             # and extra from N=step to N=step+1
271             # * (1-cos(2pi/step))
272             #
273             sub dx_minimum {
274 0     0 1 0 my ($self) = @_;
275 0 0       0 if ($self->{'step'} == 0) {
276 0         0 return 1; # horizontal only
277             }
278              
279 0 0       0 if ($self->{'step'} > 6) {
280 0         0 return -1; # supremum, unless polygon and step even
281             }
282 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
283             # step=3,4,5
284 0         0 return (-2*_PI()) / $self->{'step'};
285             } else {
286 0         0 return (-2*_PI()) / $self->{'step'};
287             }
288             }
289              
290             sub dx_maximum {
291 0     0 1 0 my ($self) = @_;
292             return ($self->{'step'} == 0
293             ? 1 # horizontal only
294              
295             : $self->{'step'} == 5
296             ? $self->{'base_r'}+2 - ($self->{'base_r'}+1)*cos(2*_PI()/5)
297              
298             : $self->{'step'} == 6
299             ? 1.5
300              
301             : $self->{'step'} <= 6
302             ? (2*_PI()) / $self->{'step'}
303              
304             # step > 6, between rings
305             : (0.5/_PI()) * $self->{'step'}
306 0 0       0 * (2-cos(2*_PI()/$self->{'step'})));
    0          
    0          
    0          
307             }
308              
309             sub dy_minimum {
310 0     0 1 0 my ($self) = @_;
311             return ($self->{'step'} == 0 ? 0 # horizontal only
312 0 0       0 : $self->{'step'} <= 6 ? (-2*_PI) / $self->{'step'}
    0          
313             : -1); # supremum
314             }
315             sub dy_maximum {
316 0     0 1 0 my ($self) = @_;
317             return ($self->{'step'} == 0 ? 0 # horizontal only
318 0 0       0 : $self->{'step'} <= 6 ? (2*_PI) / $self->{'step'}
    0          
319             : 1); # supremum
320             }
321             sub _UNDOCUMENTED__dxdy_list {
322 0     0   0 my ($self) = @_;
323 0 0       0 return ($self->{'step'} == 0 ? (1,0) # E only
324             : ()); # unlimited
325             }
326              
327             sub absdx_minimum {
328 0     0 1 0 my ($self) = @_;
329 0         0 my $step = $self->{'step'};
330 0 0       0 if ($step == 0) {
331 0         0 return 1; # horizontal dX=1 always
332             }
333 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
334 0 0       0 if ($step % 2) {
335 0         0 return 0; # polygons with odd num sides have left vertical dX=0
336             } else {
337 0         0 return sin(_PI/2 /$step);
338             }
339              
340             # if ($self->{'step'} % 2 == 1) {
341             #
342             # return 0;
343             # } else {
344             # return abs($self->dx_minimum);
345             # }
346             }
347 0         0 return 0;
348             }
349             sub absdy_minimum {
350 0     0 1 0 my ($self) = @_;
351 0         0 my $step = $self->{'step'};
352 0 0       0 if ($step == 0) {
353 0         0 return 0; # horizontal dX=1 always
354             }
355 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
356 0 0       0 if ($step == 3) {
357 0         0 return 0.5; # sin(30 degrees) innermost polygon
358             }
359 0         0 my $frac = ($step+2) % 4;
360 0 0       0 if ($frac == 3) { $frac = 1; }
  0         0  
361 0         0 return sin(_PI/2 * $frac/$step);
362             }
363 0         0 return 0;
364             }
365              
366             sub dsumxy_minimum {
367 0     0 1 0 my ($self) = @_;
368 0 0       0 return ($self->{'step'} == 0
369             ? 1 # horizontal only
370             : -1); # infimum
371             }
372 15     15   131 use constant dsumxy_maximum => 1;
  15         39  
  15         44972  
373              
374             # FIXME: for step=1 is there a supremum at 9 or thereabouts?
375             # and for other step<6 too?
376             # 2*dXmax * sqrt(2) ?
377             sub ddiffxy_minimum {
378 0     0 1 0 my ($self) = @_;
379             return ($self->{'step'} == 0 ? 1 # horizontal only
380 0 0       0 : $self->{'step'} <= 6 ? $self->dx_minimum * sqrt(2)
    0          
381             : -1); # infimum
382             }
383             sub ddiffxy_maximum {
384 0     0 1 0 my ($self) = @_;
385             return ($self->{'step'} == 0 ? 1 # horizontal only
386 0 0       0 : $self->{'step'} <= 6 ? $self->dx_maximum * sqrt(2)
    0          
387             : 1); # supremum
388             }
389              
390             #------------------------------------------------------------------------------
391             # dir_maximum_dxdy()
392              
393             # polygon step many sides
394             # start at vertical angle 1/4 plus 0.5/step, then k*1/step each side
395             # a = 1/4 + (k+1/2)/step
396             # = (1 + 4(k+1/2)/step) / 4
397             # = ((4*k+2)/step + 1) / 4
398             #
399             # maximum want 1 > a >= 1-1/step
400             # 1/4 + (k+1/2)/step >= 1-1/step
401             # (k+1/2)/step >= 3/4-1/step
402             # k+1/2 >= 3*step/4-1
403             # k >= 3*step/4-3/2
404             # k >= (3*step-6)/4
405             # k = ceil((3*step-6)/4)
406             # = floor((3*step-6)/4 + 3/4)
407             # = floor((3*step-3)/4)
408             # high side
409             # 1/4 + (k+1/2)/step < 1
410             # (k+1/2)/step < 3/4
411             # k+1/2 < 3*step/4
412             # k < (3*step-2)/4
413             # k = floor((3*step-2)/4 - 1/4)
414             # = floor((3*step-3)/4)
415             #
416             # so
417             # a = 1/4 + (floor((3*step-3)/4) + 1/2)/step
418             # = (1 + 4*(floor((3*step-3)/4) + 1/2)/step) / 4
419             # = ((floor((3*step-3)/4)*4 + 2)/step + 1) / 4
420             # step=4 a = 7/8
421             # step=5 a = 19/20
422             # step=6 a = 5/6
423             # step=7 a = 25/28
424             # step=8 a = 15/16
425             # step=10 a = 9/10
426             # return (int((3*$step-3)/4) * 4 + 2)/$step + 1;
427             # is full circle less 4,3,2,1 as step-2 mod 4
428             #
429             # sub dir4_maximum {
430             # my ($self) = @_;
431             # if ($self->{'step'} == 0) {
432             # return 0; # horizontal only
433             # }
434             # my $step = $self->{'step'};
435             # if ($self->{'ring_shape'} eq 'polygon') {
436             # return (($step-2)%4 - 4)/$step + 4;
437             # }
438             # return 4; # supremum, full circle
439             # }
440              
441             # want a >= 1
442             # 1/4 + (k+1/2)/step >= 1
443             # (k+1/2)/step >= 3/4
444             # k+1/2 >= 3*step/4
445             # k >= 3*step/4 - 1/2
446             # k >= (3*step-2)/4
447             # k = ceil((3*step-2)/4)
448             # = floor((3*step-2)/4 + 3/4)
449             # = floor((3*step+1)/4)
450             # min_a = 1/4 + (floor((3*step+1)/4) + 1/2)/step - 1
451             # = (1 + 4*(floor((3*step+1)/4) + 1/2)/step ) / 4
452             # = ((4*floor((3*step+1)/4) + 2)/step + 1) / 4 - 1
453             # = ((floor((3*step+1)/4)*4 + 2)/step - 3) / 4
454             # return (int((3*$step+1)/4) * 4 + 2)/$step - 3;
455             # is 0,1,2,3 as step-2 mod 4
456             # return (($step-2) % 4) / $step;
457             #
458             # but last of ring across to first of next may be shallower
459             #
460             # sub dir4_minimum {
461             # my ($self) = @_;
462             # my $step = $self->{'step'};
463             # if ($self->{'ring_shape'} eq 'polygon') {
464             # if ($step % 4 != 2) { # polygon step=2mod4 includes horizontal ...
465             # my ($dx,$dy) = $self->n_to_dxdy($self->{'step'});
466             # return min (atan2($dy,$dx) * (2/_PI),
467             # (($step-2) % 4) / $step);
468             # }
469             #
470             # }
471             # return 0; # horizontal
472             # }
473              
474             sub dir_minimum_dxdy {
475 0     0 1 0 my ($self) = @_;
476 0         0 my $step = $self->{'step'};
477 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
478 0 0       0 return $self->n_to_dxdy($step == 9
479             ? 9
480             : int((3*$step+5)/4));
481             }
482 0         0 return (1,0); # horizontal
483             }
484             sub dir_maximum_dxdy {
485 0     0 1 0 my ($self) = @_;
486 0 0       0 if ($self->{'step'} == 0) {
487 0         0 return (1,0); # step=0 horizontal always
488             }
489              
490 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
491 0         0 my $step = $self->{'step'};
492 0         0 return $self->n_to_dxdy(int((3*$step+1)/4)); # 1 before the minimum
493              
494             # # just before 3/4 way around, then half back ....
495             # # sides side
496             # # ----- ----
497             # # 3 1
498             # # 4 2
499             # # 5 3
500             # # 6 3
501             # # 7 4
502             # # 8 5
503             # # 9 6
504             # # 10 6
505             # return _circlefrac_to_xy (1, int((3*$step-3)/4), $step, _PI);
506             }
507              
508 0         0 return (0,0); # supremum, full circle
509             }
510              
511             #------------------------------------------------------------------------------
512              
513             sub new {
514             ### MultipleRings new() ...
515 146     146 1 15154 my $self = shift->SUPER::new(@_);
516              
517 146         259 my $step = $self->{'step'};
518 146 50       339 $step = $self->{'step'} = (! defined $step ? 6 # default
    100          
519             : $step < 0 ? 0 # minimum
520             : $step);
521             ### $step
522              
523 146   100     437 my $ring_shape = ($self->{'ring_shape'} ||= 'circle');
524 146 50 66     271 if (! ($ring_shape eq 'circle' || $ring_shape eq 'polygon')) {
525 0         0 croak "Unrecognised ring_shape option: ", $ring_shape;
526             }
527 146 100       243 if ($step < 3) {
528             # polygon shape only for step >= 3
529 79         116 $ring_shape = $self->{'ring_shape'} = 'circle';
530             }
531              
532 146 100       352 if ($ring_shape eq 'polygon') {
    100          
    100          
    100          
    100          
533             ### polygon ...
534 4 50       9 if ($step == 6) {
    0          
535             ### 0.5/sin(PI/6)=1 exactly ...
536 4         8 $self->{'base_r'} = 1;
537             } elsif ($step == 3) {
538             ### 0.5/sin(PI/3)=sqrt(3)/3 ...
539 0         0 $self->{'base_r'} = sqrt(3)/3;
540             } else {
541 0         0 $self->{'base_r'} = 0.5/sin(_PI/$step);
542             }
543              
544             } elsif ($step == 6) {
545             ### 0.5/sin(PI/6) = 1 exactly ...
546 18         32 $self->{'base_r'} = 0;
547              
548             } elsif ($step == 4) {
549             ### 0.5/sin(PI/4) = sqrt(2)/2 ...
550 13         21 $self->{'base_r'} = sqrt(2)/2 - 1;
551              
552             } elsif ($step == 3) {
553             ### 0.5/sin(PI/3) = sqrt(3)/3 ...
554 12         20 $self->{'base_r'} = sqrt(3)/3 - 1;
555              
556             } elsif ($step < 6) {
557             ### sin: $step>1 && sin(_PI/$step)
558 83   66     190 $self->{'base_r'} = ($step > 1 && 0.5/sin(_PI/$step)) - 1;
559             }
560             ### base r: $self->{'base_r'}
561              
562 146         230 return $self;
563             }
564              
565             # with N decremented
566             # d = [ 1, 2, 3, 4, 5 ]
567             # N = [ 0, 1, 3, 6, 10 ]
568             #
569             # N = (1/2 d^2 - 1/2 d)
570             # = (1/2*$d**2 - 1/2*$d)
571             # = ((0.5*$d - 0.5)*$d)
572             # = 0.5*$d*($d-1)
573             #
574             # d = 1/2 + sqrt(2 * $n + 1/4)
575             # = 0.5 + sqrt(2*$n + 0.25)
576             # = [ 1 + 2*sqrt(2n + 1/4) ] / 2
577             # = [ 1 + sqrt(8n + 1) ] / 2
578             #
579             # (d+1)d/2 - d(d-1)/2
580             # = [ (d^2 + d) - (d^2-d) ] / 2
581             # = [ d^2 + d - d^2 + d ] / 2
582             # = 2d/2 = d
583             #
584             # radius
585             # step > 6 1 / (2 * sin(pi / ($d*$step))
586             # step <= 6 Rbase + d
587             #
588             # usual polygon formula R = a / 2*sin(pi/n)
589             # cf inner radius r = a / 2*tan(pi/n)
590             # along chord
591             #
592             # polygon horizontal when a=1
593             # 1/4 + (k+1/2)/step = 1
594             # (k+1/2)/step = 3/4
595             # k+1/2 = 3*step/4
596             # k = 3*step/4 - 1/2
597             # k = ()/4
598             # 4*k = 3*step-2
599             # and when a=1/2
600             # 1/4 + (k+1/2)/step = 1/2
601             # (k+1/2)/step = 1/4
602             # k+1/2 = step/4
603             # 4*k+2 = step
604              
605             # 1/2 / R = sin(2pi/sides)
606             # 1/2 / (R^2 - 1/4) = tan(2pi/sides)
607             # f(x) = 1/2 / R - sin(2pi/sides) = $f
608             # f'(x) = -1/2 / R^2 - cos(2pi/sides) = $slope
609             # $r-$f/$slope better approx
610             # (1/2 / R - sin(2pi/sides)) / (-1/2 / R^2 - cos(2pi/sides))
611             # = (R/2 - R^2 sin(2pi/sides)) / (-1/2 - R^2 * cos(2pi/sides))
612              
613             sub n_to_xy {
614 179     179 1 947 my ($self, $n) = @_;
615             ### MultipleRings n_to_xy(): "n=$n step=$self->{'step'} shape=$self->{'ring_shape'}"
616              
617             # "$n<1" separate test from decrement so as to warn on undef
618             # don't have anything sensible for infinity, and _PI / infinity would
619             # throw a div by zero
620 179 50       295 if ($n < 1) { return; }
  0         0  
621 179 50       875 if (is_infinite($n)) { return ($n,$n); }
  0         0  
622 179         1571 $n -= 1;
623              
624             ### decremented n: $n
625 179         940 my $step = $self->{'step'};
626 179 100       278 if (! $step) {
627             ### step==0 goes along X axis ...
628 13         33 return ($n, 0);
629             }
630              
631 166         382 my $d = int((_sqrtint(8*$n/$step + 1) + 1) / 2);
632              
633             ### d frac: (sqrt(int(8*$n) + 1) + 1) / 2
634             ### d int: "$d"
635             ### base: ($d*($d-1)/2).''
636             ### next base: (($d+1)*$d/2).''
637             ### assert: $n >= ($d*($d-1)/2)
638             ### assert: $n < ($step * ($d+1) * $d / 2)
639              
640 166         1433 $n -= $d*($d-1)/2 * $step;
641             ### n remainder: "$n"
642             ### assert: $n >= 0
643             ### assert: $n < $d*$step
644              
645 166         1808 my $zero = $n * 0;
646 166 100       780 if (ref $n) {
647 2 100       10 if ($n->isa('Math::BigInt')) {
    50          
648 1         5 $n = Math::PlanePath::SacksSpiral::_bigfloat()->new($n);
649             } elsif ($n->isa('Math::BigRat')) {
650 0         0 $n = $n->as_float;
651             }
652 2 50       176 if ($n->isa('Math::BigFloat')) {
653             ### bigfloat ...
654 2         27 $d = Math::BigFloat->new($d);
655             }
656             }
657 166         431 my $pi = _pi($n);
658             ### $pi
659              
660             # my $base_r = $self->{'base_r'};
661             # $base_r = Math::BigFloat->new($base_r);
662              
663             {
664 166         1176 my $numsides;
  166         225  
665             my $r;
666 166 100       255 if ($self->{'ring_shape'} eq 'circle') {
667             ### circle ...
668 162         217 $numsides = $d * $step;
669 162 100       873 if ($step > 6) {
670 20         32 $r = 0.5 / sin($pi / $numsides);
671             } else {
672 142         151 my $base_r;
673 142 100       263 if ($step == 6) {
    100          
    100          
    100          
674 17         23 $base_r = 0; # exactly
675             } elsif ($step == 4) {
676             ### 0.5/sin(PI/4)=sqrt(2)/2 ...
677 21         34 $base_r = sqrt(0.5 + $zero) - 1; # sqrt() instead of sin()
678             } elsif ($step == 3) {
679             ### 0.5/sin(PI/3)=sqrt(3)/3 ...
680 19         31 $base_r = sqrt(3 + $zero)/3 - 1; # sqrt() instead of sin()
681             } elsif ($step == 1) {
682 51         60 $base_r = -1; # so initial d=1 at $r=0
683             } else {
684 34         55 $base_r = 0.5/sin($pi/$step) - 1;
685             }
686 142         173 $r = $base_r + $d;
687             }
688             } else {
689             ### polygon ...
690 4         7 $numsides = $step;
691 4         10 my $base_r = _numsides_to_r($step,$pi);
692 4 50       7 if ($step > 6) {
693 0         0 $r = $base_r*$d;
694             } else {
695 4         9 $r = $base_r + ($d-1)/cos($pi/$step);
696             }
697 4         7 $n /= $d;
698             }
699             ### n with frac: $n
700              
701             # numsides even N > numsides/2
702             # numsides odd N >= (numsides+1)/2 = ceil(numsides/2)
703 166         687 my $y_neg;
704 166 100       277 if (2*$n >= $numsides) {
705 51         597 $n = $numsides - $n;
706 51         321 $y_neg = 1;
707             }
708              
709 166         585 my $x_neg;
710             my $xy_transpose;
711 166 100       293 if ($numsides % 2 == 0) {
712 120 100       1578 if (4*$n >= $numsides) {
713 48         524 $n = $numsides/2 - $n;
714 48         1058 $x_neg = 1;
715             }
716 120 100 100     753 if ($numsides % 4 == 0 && 8*$n >= $numsides) {
717 19         36 $n = $numsides/4 - $n;
718 19         23 $xy_transpose = 1;
719             }
720             }
721              
722 166         1578 my $side = int($n);
723 166         346 $n -= $side;
724             ### $side
725              
726 166         469 my ($x, $y) = _circlefrac_to_xy($r, $side, $numsides, $pi);
727              
728 166 100       485548 if ($n) {
729             # fractional n offset into side ...
730              
731 25         29 my ($to_x, $to_y);
732 25         30 $side += 1;
733 25 100 66     73 if (2*$side == $numsides+1) {
    100          
734             # vertical at left, so X unchanged Y negate
735 3         4 $to_x = $x;
736 3         4 $to_y = - $y;
737              
738             } elsif (4*$side == $numsides+2 || 4*$side == 3*$numsides-2) {
739             # horizontal at top or bottom, so Y unchanged X negate
740 10         13 $to_x = - $x;
741 10         14 $to_y = $y;
742              
743             } else {
744 12         15 ($to_x, $to_y) = _circlefrac_to_xy($r, $side, $numsides, $pi);
745             }
746              
747             ### $side
748             ### $r
749             ### from: "$x, $y"
750             ### to: "$to_x, $to_y"
751              
752             # If vertical or horizontal then don't apply the proportions since the
753             # two parts $x*$n and $to_x*(1-$n) can round off giving the sum != to
754             # the original $x.
755 25 100       55 if ($to_x != $x) {
756 22         32 $x = $x*(1-$n) + $to_x*$n;
757             }
758 25 100       40 if ($to_y != $y) {
759 14         68 $y = $y*(1-$n) + $to_y*$n;
760             }
761             }
762              
763 166 100       306 if ($xy_transpose) {
764 19         36 ($x,$y) = ($y,$x);
765             }
766 166 100       228 if ($x_neg) {
767 48         65 $x = -$x;
768             }
769 166 100       271 if ($y_neg) {
770 51         60 $y = -$y;
771             }
772              
773             ### final: "x=$x y=$y"
774 166         502 return ($x, $y);
775             }
776              
777             # {
778             # # && $d != 0 # watch out for overflow making d==0 ??
779             # #
780             # my $d_step = $d*$step;
781             # my $r = ($step > 6
782             # ? 0.5 / sin($pi / $d_step)
783             # : $base_r + $d);
784             # ### r: "$r"
785             #
786             # my $n2 = 2*$n;
787             #
788             # if ($n2 == int($n2)) {
789             # if (($n2 % $d_step) == 0) {
790             # ### theta=0 or theta=pi, exactly on X axis ...
791             # return ($n ? -$r : $r, # n remainder 0 means +ve X axis, non-zero -ve
792             # 0);
793             # }
794             # if (($d_step % 2) == 0) {
795             # my $n2sub = $n2 - $d_step/2;
796             # if (($n2sub % $d_step) == 0) {
797             # ### theta=pi/2 or theta=3pi/2, exactly on Y axis ...
798             # return (0,
799             # $n2sub ? -$r : $r);
800             # }
801             # }
802             # }
803             #
804             # my $theta = $n2 * $pi / $d_step;
805             #
806             # ### theta frac: (($n - $d*($d-1)/2)/$d).''
807             # ### theta: "$theta"
808             #
809             # return ($r * cos($theta),
810             # $r * sin($theta));
811             # }
812             }
813              
814             # $side is 0 to $numsides-1
815             sub _circlefrac_to_xy {
816 178     178   270 my ($r, $side, $numsides, $pi) = @_;
817             ### _circlefrac_to_xy(): "r=$r side=$side numsides=$numsides pi=$pi"
818              
819 178 50       263 if (2*$side == $numsides) {
820             ### 180-degrees, so X=R, Y=0 ...
821 0         0 return (-$r, 0);
822              
823             }
824 178 100       1092 if (4*$side == $numsides) {
825             ### 90-degrees, so X=0, Y=R ...
826 4         9 return (0, $r);
827             }
828 174 100       1032 if (6*$side == $numsides) {
829             ### 60-degrees, so X=R/2, Y=sqrt(3)/2*R ...
830 7         19 return ($r / 2,
831             $r * sqrt(3 + $r*0) / 2);
832             }
833 167 100       1003 if (8*$side == $numsides) {
834             ### 45-degrees, so X=Y=R/sqrt(2) ...
835 1         4 my $x = $r / sqrt(2 + $r*0);
836 1         2 return ($x, $x);
837             }
838              
839             # my $two_pi = (ref $r && $r->isa('Math::BigFloat')
840             # ? 2*Math::BigFloat->bpi;
841             # : 2*_PI);
842             #
843             # if (2*$side == $numsides+1) {
844             # ### first below X axis ...
845             # my $theta = 2*$pi * ($side-1)/$numsides;
846             # return ($r * cos($theta),
847             # - $r * sin($theta));
848             # }
849             # if (4*$side == $numsides+1) {
850             # ### first past Y axis ...
851             # my $theta = 2*$pi * ($side-1)/$numsides;
852             # return (- $r * cos($theta),
853             # $r * sin($theta));
854             # }
855              
856             ### general case ...
857 166         1068 my $theta = 2 * $pi * $side/$numsides;
858 166         2847 return (cos($theta) * $r,
859             sin($theta) * $r);
860             }
861              
862             # my $numsides = $step;
863             # if ($self->{'ring_shape'} eq 'polygon') {
864             # $n /= $d;
865             # my $base_r = _numsides_to_r($step,$pi);
866             # if ($step > 6) {
867             # $r = $base_r*$d;
868             # } else {
869             # $r = $base_r + ($d-1)/cos($pi/$step);
870             # }
871             # } else {
872             # $numsides *= $d;
873             # if ($step > 6) {
874             # $r = _numsides_to_r($numsides,$pi);
875             # } else {
876             # $r = _numsides_to_r($step,$pi) + $d;
877             # }
878             # }
879             # my $side = int($n);
880             # $n -= $side;
881              
882             sub _numsides_to_r {
883 4     4   8 my ($numsides, $pi) = @_;
884 4 50       7 if ($numsides == 3) { return sqrt(0.75 + $pi*0); }
  0         0  
885 4 50       8 if ($numsides == 4) { return sqrt(0.5 + $pi*0); }
  0         0  
886 4 50       7 if ($numsides == 6) { return 1 + $pi*0; }
  4         9  
887 0         0 return 0.5 / sin($pi/$numsides);
888             }
889              
890              
891             # for step=4
892             # R = sqrt(2)/2 + d
893             # R^2 = (sqrt(2)/2 + d)^2
894             # = 2/4 + 2*sqrt(2)/2*d + d^2
895             # = 1/2 + d*sqrt(2) + d^2
896             # not an integer
897             #
898             sub n_to_rsquared {
899 107     107 1 7577 my ($self, $n) = @_;
900             ### MultipleRings n_to_rsquared(): "n=$n"
901 107 50       211 if ($n < 1) { return undef; }
  0         0  
902 107 50       214 if (is_infinite($n)) { return $n; }
  0         0  
903              
904 107 100       187 if (defined (my $r = _n_to_radius_exact($self,$n))) {
905 55         112 return $r*$r;
906             }
907 52 100       85 if ($self->{'step'} == 1) {
908             # $n < 4 covered by _n_to_radius_exact()
909              
910 26 100 66     72 if ($n >= 4 && $n < 7) {
911             # triangle numsides=3
912             # N=4 at X=2, Y=0
913             # N=5 at X=-1, Y=sqrt(3)
914             # N=4+f at X=2-3*f Y=f*sqrt(3)
915             # R^2 = (2-3f)^2 + 3*f^2
916             # = 4-12f+9*f^2 + 3*f^2
917             # = 4-12f+12*f^2
918             # = 4*(1 - 3f + 3*f^2)
919             # = 4 - 6*(2*f) + 3*(2*f)^2
920             # f=1/2 is R^2 = 1
921             # N=5+f at X=-1 Y = sqrt(3)*(1-2*f)
922             # R^2 = 1 + 3*(1-2*f)^2
923             # = 1 + 3 - 3*4*f + 3*4*f^2
924             # = 4 - 12*f + 12*f^2
925             # = 4 - 12*(f - f^2)
926             # = 4 - 12*f*(1 - f)
927              
928 12         19 $n -= int($n);
929 12         29 return 4 - 12*$n*(1-$n);
930             }
931              
932 14 100 66     37 if ($n >= 7 && $n < 11) {
933             ### square numsides=4 ...
934             # X=3-3*f Y=3*f
935             # R^2 = (3-3*f)^2 + (3*f)^2
936             # = 9*[ (1-f)^2 + f^2) ]
937             # = 9*[ 1 - 2f + f^2 + f^2) ]
938             # = 9*[ 1 - 2f + 2f^2 ]
939             # = 9*[ 1 - 2(f - f^2) ]
940             # = 9 - 18*f*(1 - f)
941             # eg f=1/2 R^2 = (sqrt(2)/2*3)^2 = 2/4*9 = 9/2
942              
943 8         14 $n -= int($n);
944 8         20 return 9 - 18*$n*(1-$n);
945             }
946              
947 6 50 33     18 if ($n >= 16 && $n < 22) {
948             ### hexagon numsides=6 ...
949             # X=5 Y=0 to X=5*1/2 Y=5*sqrt(3)/2
950             # R^2 = (5 - 5/2*f)^2 + (5*sqrt(3)/2*f)^2
951             # = 25 - 25*f + 25*f^2
952             # = 25 - 25*f*(1-f)
953             # eg f=1/2 R^2 = 18.75
954             # or f=1/5 R^2 = 21 exactly, though 1/5 not exact in binary floats
955              
956 6         12 $n -= int($n);
957 6         15 return 25 - 25*$n*(1-$n);
958             }
959              
960             # other numsides don't have sin(pi/numsides) an integer or sqrt so
961             # aren't an exact R^2
962             }
963              
964             # ENHANCE-ME: step=1 various exact values for ring of 4 and ring of 6
965              
966 26         60 return $self->SUPER::n_to_rsquared($n);
967             }
968             sub n_to_radius {
969 43     43 1 2729 my ($self, $n) = @_;
970             ### n_to_radius(): $n
971              
972 43 50       92 if ($n < 1) { return undef; }
  0         0  
973 43 50       90 if (is_infinite($n)) { return $n; }
  0         0  
974              
975 43 100       78 if (defined (my $r = _n_to_radius_exact($self,$n))) {
976 30         53 return $r;
977             }
978 13         38 return sqrt($self->n_to_rsquared($n));
979             # return $self->SUPER::n_to_radius($n);
980             }
981              
982             # step=6 shape=polygon exact integer for some of second ring too
983             # sub n_to_trsquared {
984             # my ($self, $n) = @_;
985             # ### MultipleRings n_to_rsquared(): "n=$n"
986             # }
987              
988             sub _n_to_radius_exact {
989 150     150   208 my ($self, $n) = @_;
990             ### _n_to_radius_exact(): "n=$n step=$self->{'step'}"
991              
992 150 50       229 if ($n < 1) { return undef; }
  0         0  
993 150 50       241 if (is_infinite($n)) { return $n; }
  0         0  
994              
995 150         245 my $step = $self->{'step'};
996 150 100       227 if ($step == 0) {
997 13         37 return $n - 1; # step=0 goes along X axis starting X=0,Y=0
998             }
999              
1000 137 100       215 if ($step == 1) {
    100          
1001 89 100       134 if ($n < 4) {
1002 26 100       46 if ($n < 2) {
1003 4         9 return 0; # 0,0 only, no jump across to next ring
1004             }
1005 22         28 $n -= int($n);
1006 22         68 return abs(1-2*$n);
1007             }
1008 63 100       114 if ($n == int($n)) {
1009             ### step=1 radius=integer steps for integer N ...
1010 22         41 return _n0_to_d($self,$n-1) - 1;
1011             }
1012 41         46 my $two_n = 2*$n;
1013 41 50 66     160 if ($two_n == 9 || $two_n == 11 || $two_n == 13) {
      66        
1014             # N=4.5 at X=1/2 Y=sqrt(3)/2 R^2 = 1/4 + 3/4 = 1 exactly
1015             # N=5.5 at X=-1, Y=0 so R^2 = 1 exactly
1016             # N=6.5 same as N=4.5
1017 2         7 return 1;
1018             }
1019              
1020             } elsif ($step == 6) {
1021 22 50       43 if ($n == int($n)) {
1022             # step=6 circle all integer N has exact integer radius
1023             # step=6 polygon only innermost ring N<=6 exact integer radius
1024 22 50 66     49 if ($self->{'ring_shape'} eq 'circle'
1025             || $n <= 6) { # ring_shape=polygon
1026 22         34 return _n0_to_d($self,$n-1);
1027             }
1028             }
1029             }
1030              
1031             ### no exact radius ...
1032 65         127 return undef;
1033             }
1034             sub _n0_to_d {
1035 44     44   61 my ($self, $n) = @_;
1036 44         166 return int((_sqrtint(8*$n/$self->{'step'} + 1) + 1) / 2);
1037             }
1038             sub _d_to_n0base {
1039 51     51   65 my ($self, $d) = @_;
1040 51         100 return $d*($d-1)/2 * $self->{'step'};
1041             }
1042              
1043             # From above
1044             # r = 0.5 / sin(pi/(d*step))
1045             #
1046             # sin(pi/(d*step)) = 0.5/r
1047             # pi/(d*step) = asin(1/(2*r))
1048             # 1/d * pi/step = asin(1/(2*r))
1049             # d = pi/(step*asin(1/(2*r)))
1050             #
1051             # r1 = 0.5 / sin(pi/(d*step))
1052             # r2 = 0.5 / sin(pi/((d+1)*step))
1053             # r2 - r1 = 0.5 / sin(pi/(d*step)) - 0.5 / sin(pi/((d+1)*step))
1054             # r2-r1 >= 1 when step>=7 ?
1055              
1056             sub _xy_to_d {
1057 51     51   75 my ($self, $x, $y) = @_;
1058             ### _xy_to_d(): "x=$x y=$y"
1059              
1060 51         107 my $r = hypot ($x, $y);
1061 51 50       102 if ($r < 0.5) {
1062             ### r smaller than 0.5 ring, treat as d=1
1063             # 1/(2*r) could be div-by-zero
1064             # or 1/(2*r) > 1 would be asin()==-nan
1065 51         108 return 1;
1066             }
1067 0         0 my $two_r = 2*$r;
1068 0 0       0 if (is_infinite($two_r)) {
1069             ### 1/inf is a divide by zero, avoid that ...
1070 0         0 return $two_r;
1071             }
1072             ### $r
1073              
1074 0         0 my $step = $self->{'step'};
1075 0 0       0 if ($self->{'ring_shape'} eq 'polygon') {
1076 0         0 my $theta_frac = _xy_to_angle_frac($x,$y);
1077 0         0 $theta_frac -= int($theta_frac*$step) / $step; # modulo 1/step
1078              
1079 0         0 my $r = hypot ($x, $y);
1080 0         0 my $alpha = 2*_PI/$step;
1081 0         0 my $theta = 2*_PI * $theta_frac;
1082             ### $r
1083             ### x=r*cos(theta): $r*cos($theta)
1084             ### y=r*sin(theta): $r*sin($theta)
1085              
1086 0         0 my $p = $r*cos($theta) + $r*sin($theta) * sin($alpha/2)/cos($alpha/2);
1087             ### $p
1088             ### base_r: $self->{'base_r'}
1089             ### p - base_r: $p - $self->{'base_r'}
1090              
1091 0 0       0 if ($step >= 6) {
1092 0         0 return $p / $self->{'base_r'};
1093             } else {
1094 0         0 return ($p - $self->{'base_r'}) * cos(_PI/$step) + 1;
1095             }
1096             }
1097              
1098 0 0       0 if ($step > 6) {
1099             ### d frac by asin: _PI / ($step * asin(1/$two_r))
1100 0         0 return _PI / ($step * asin(1/$two_r));
1101             } else {
1102             # $step <= 6
1103             ### d frac by base: $r - $self->{'base_r'}
1104 0         0 return $r - $self->{'base_r'};
1105             }
1106             }
1107              
1108             sub xy_to_n {
1109 56     56 1 138 my ($self, $x, $y) = @_;
1110             ### MultipleRings xy_to_n(): "$x, $y step=$self->{'step'} shape=$self->{'ring_shape'}"
1111              
1112 56         60 my $n;
1113 56         83 my $step = $self->{'step'};
1114 56 100       87 if ($step == 0) {
1115             # step==0
1116 5         11 $n = int ($x + 1.5);
1117              
1118             } else {
1119 51         83 my $theta_frac = _xy_to_angle_frac($x,$y);
1120             ### $theta_frac
1121             ### assert: (0 <= $theta_frac && $theta_frac < 1) || $theta_frac!=$theta_frac
1122              
1123 51         56 my $d;
1124 51 50       81 if ($self->{'ring_shape'} eq 'polygon') {
1125 0         0 $n = int($theta_frac*$step);
1126 0         0 $theta_frac -= $n/$step;
1127             ### theta modulo 1/step: $theta_frac
1128             ### $n
1129              
1130 0         0 my $r = hypot ($x, $y);
1131 0         0 my $alpha = 2*_PI/$step;
1132 0         0 my $theta = 2*_PI * $theta_frac;
1133             ### $r
1134             ### so x=r*cos(theta): $r*cos($theta)
1135             ### so y=r*sin(theta): $r*sin($theta)
1136              
1137 0         0 my $pi = _PI;
1138 0         0 my $p = $r*cos($theta) + $r*sin($theta) * sin($alpha/2)/cos($alpha/2);
1139 0         0 my $base_r = Math::PlanePath::MultipleRings::_numsides_to_r($step,$pi);
1140             ### $p
1141             ### $base_r
1142              
1143 0 0       0 if ($step > 6) {
1144 0         0 $d = $p / $base_r;
1145             } else {
1146 0         0 $d = ($p - $base_r) * cos($pi/$step) + 1;
1147             }
1148             ### d frac: $d
1149 0         0 $d = int($d+0.5);
1150             ### $d
1151             ### cf _xy_to_d(): _xy_to_d($self,$x,$y)
1152              
1153 0 0       0 my $f = ($p == 0 ? 0 : $r*sin($theta) / ($p*sin($alpha)));
1154 0         0 $n = int(($n+$f)*$d + 0.5);
1155              
1156             ### e: $r*sin($theta) * sin($alpha/2)/cos($alpha/2)
1157             ### $f
1158             ### $n
1159              
1160             } else {
1161 51         83 $d = int(_xy_to_d($self,$x,$y) + 0.5);
1162             ### $d
1163 51         80 $n = int (0.5 + $theta_frac * $d*$step);
1164 51 50       84 if ($n >= $d*$step) { $n = 0; }
  0         0  
1165             }
1166              
1167             ### n within ring: $n
1168             ### n ring start: _d_to_n0base($self,$d) + 1
1169              
1170 51         74 $n += _d_to_n0base($self,$d) + 1;
1171             ### $d
1172             ### d base: 0.5*$d*($d-1)
1173             ### d base M: $step * 0.5*$d*($d-1)
1174             ### $theta_frac
1175             ### theta offset: $theta_frac*$d
1176             ### $n
1177             }
1178              
1179             ### trial n: $n
1180 56 50       94 if (my ($nx, $ny) = $self->n_to_xy($n)) {
1181             ### nxy: "nx=$nx ny=$ny hypot=".hypot($x-$nx,$y-$ny)
1182             ### cf orig xy: "x=$x y=$y"
1183 56 100       158 if (hypot($x-$nx, $y-$ny) <= 0.5) {
1184 17         71 return $n;
1185             }
1186             }
1187 39         131 return undef;
1188             }
1189              
1190             # ENHANCE-ME: step>=3 small rectangles around 0,0 don't cover any pixels
1191             #
1192             # not exact
1193             sub rect_to_n_range {
1194 22     22 1 69 my ($self, $x1,$y1, $x2,$y2) = @_;
1195             ### MultipleRings rect_to_n_range(): "$x1,$y1, $x2,$y2 step=$self->{'step'}"
1196              
1197 22   66     50 my $zero = ($x1<0) != ($x2<0) || ($y1<0) != ($y2<0);
1198 22         34 my $step = $self->{'step'};
1199              
1200 22         60 my ($r_lo, $r_hi) = Math::PlanePath::SacksSpiral::_rect_to_radius_range
1201             ($x1,$y1, $x2,$y2);
1202             ### $r_lo
1203             ### $r_hi
1204 22 50       46 if (is_infinite($r_hi)) {
1205 0         0 return (1,$r_hi);
1206             }
1207 22 100       42 if ($r_hi < 1) { $r_hi = 1; }
  11         12  
1208 22 50       43 if ($self->{'ring_shape'} eq 'polygon') {
1209 0         0 $r_hi /= cos(_PI/$self->{'step'});
1210             ### poly increase r_hi: $r_hi
1211             }
1212              
1213 22         27 my ($d_lo, $d_hi);
1214 22 50       29 if ($self->{'ring_shape'} eq 'polygon') {
1215 0 0       0 if ($step >= 6) {
1216 0         0 $d_lo = $r_lo / $self->{'base_r'};
1217 0         0 $d_hi = $r_hi / $self->{'base_r'};
1218             } else {
1219 0         0 $d_lo = ($r_lo - $self->{'base_r'}) * cos(_PI/$step) + 1;
1220 0         0 $d_hi = ($r_hi - $self->{'base_r'}) * cos(_PI/$step) + 1;
1221             }
1222             } else {
1223 22 100       33 if ($step > 6) {
1224 8 50       14 $d_lo = ($r_lo > 0
1225             ? _PI / ($step * asin(0.5/$r_lo))
1226             : 0);
1227 8         26 $d_hi = _PI / ($step * asin(0.5/$r_hi));
1228             } else {
1229 14         20 $d_lo = $r_lo - $self->{'base_r'};
1230 14         18 $d_hi = $r_hi - $self->{'base_r'};
1231             }
1232             }
1233             ### $d_lo
1234             ### $d_hi
1235              
1236 22         31 $d_lo = int($d_lo - 1);
1237 22         30 $d_hi = int($d_hi + 2);
1238 22 50       34 if ($d_lo < 1) { $d_lo = 1; }
  22         29  
1239              
1240 22 100       41 if ($step) {
1241             # start of ring is N= 0.5*$d*($d-1) * $step + 1
1242             ### n_lo: 0.5*$d_lo*($d_lo-1) * $step + 1
1243             ### n_hi: 0.5*$d_hi*($d_hi+1) * $step
1244 20         58 return ($d_lo*($d_lo-1)/2 * $step + 1,
1245             $d_hi*($d_hi+1)/2 * $step);
1246             } else {
1247             # $step == 0
1248 2         6 return ($d_lo, $d_hi);
1249             }
1250              
1251              
1252              
1253              
1254              
1255             # # if x1,x2 pos and neg then 0 is covered and it's the minimum
1256             # # ENHANCE-ME: might be able to be a little tighter on $d_lo
1257             # my $d_lo = ($zero
1258             # ? 1
1259             # : max (1, -2 + int (_xy_to_d ($self,
1260             # min($x1,$x2),
1261             # min($y1,$y2)))));
1262             # my $d_hi = 1 + int (_xy_to_d ($self,
1263             # max($x1,$x2),
1264             # max($y1,$y2)));
1265             # ### $d_lo
1266             # ### $d_hi
1267             # if ((my $step = $self->{'step'})) {
1268             # # start of ring is N= 0.5*$d*($d-1) * $step + 1
1269             # ### n_lo: 0.5*$d_lo*($d_lo-1) * $step + 1
1270             # ### n_hi: 0.5*$d_hi*($d_hi+1) * $step
1271             # return ($d_lo*($d_lo-1)/2 * $step + 1,
1272             # $d_hi*($d_hi+1)/2 * $step);
1273             # } else {
1274             # # $step == 0
1275             # return ($d_lo, $d_hi);
1276             # }
1277             }
1278              
1279             #------------------------------------------------------------------------------
1280             # generic
1281              
1282             # _xy_to_angle_frac() returns the angle of X,Y as a fraction 0 <= angle < 1
1283             # measured anti-clockwise around from the X axis.
1284             #
1285             sub _xy_to_angle_frac {
1286 120     120   557 my ($x, $y) = @_;
1287              
1288             # perlfunc.pod warns atan2(0,0) is implementation dependent. The C99 spec
1289             # is atan2(+/-0, -0) returns +/-pi, both of which would come out 0.5 here.
1290             # Prefer 0 for any +/-0,+/-0.
1291 120 100 100     350 if ($x == 0 && $y == 0) {
1292 53         97 return 0;
1293             }
1294              
1295 67         140 my $frac = atan2($y,$x) * (0.5 / _PI);
1296             ### $frac
1297 67 100       126 if ($frac < 0) { $frac += 1; }
  16 50       23  
1298 0         0 elsif ($frac >= 1) { $frac -= 1; }
1299 67         111 return $frac;
1300             }
1301              
1302             # return pi=3.14159 etc, inheriting precision etc from $n if it's a BigFloat
1303             # or other overload
1304             sub _pi {
1305 168     168   2163 my ($n) = @_;
1306 168 100       265 if (ref $n) {
1307 3 50       11 if ($n->isa('Math::BigFloat')) {
1308 3         25 my $digits;
1309 3 100       13 if (defined($digits = $n->accuracy)) {
    50          
    50          
    50          
1310             ### n accuracy ...
1311             } elsif (defined($digits = $n->precision)) {
1312             ### n precision ...
1313 0         0 $digits = -$digits + 1;
1314             } elsif (defined($digits = Math::BigFloat->accuracy)) {
1315             ### global accuracy ...
1316             } elsif (defined($digits = Math::BigFloat->precision)) {
1317             ### global precision ...
1318 0         0 $digits = -$digits + 1;
1319             } else {
1320             ### div_scale ...
1321 1         72 $digits = Math::BigFloat->div_scale+1;
1322             }
1323             ### $digits
1324 3         48 $digits = max (1, $digits);
1325 3         11 return Math::BigFloat->bpi($digits);
1326             }
1327             ### other overload n class: ref $n
1328 0         0 my $zero = $n * 0;
1329 0         0 return 2*atan2($zero,1+$zero);
1330             }
1331 165         267 return _PI;
1332             }
1333              
1334             1;
1335             __END__