File Coverage

lib/Astro/Montenbruck/MathUtils.pm
Criterion Covered Total %
statement 148 153 96.7
branch 30 38 78.9
condition 4 5 80.0
subroutine 28 30 93.3
pod 20 20 100.0
total 230 246 93.5


line stmt bran cond sub pod time code
1             package Astro::Montenbruck::MathUtils;
2              
3 13     13   416983 use 5.22.0;
  13         73  
4 13     13   64 use feature qw/signatures/;
  13         19  
  13         1718  
5 13     13   95 no warnings qw/experimental::signatures/;
  13         24  
  13         520  
6             # The line below disables wrong perlcritic warnings
7             ## no critic qw/Subroutines::ProhibitSubroutinePrototypes/
8              
9 13     13   74 use Exporter qw/import/;
  13         22  
  13         501  
10 13     13   67 use POSIX qw (floor ceil acos modf fmod);
  13         34  
  13         114  
11 13     13   12397 use List::Util qw/any reduce/;
  13         25  
  13         828  
12              
13 13     13   2796 use Math::Trig qw/:pi :radial deg2rad rad2deg/;
  13         79577  
  13         1993  
14 13     13   99 use constant { ARCS => 3600.0 * 180.0 / pi };
  13         34  
  13         18135  
15              
16             our %EXPORT_TAGS = (
17             all => [
18             qw/frac frac360 dms hms zdms ddd polynome sine
19             reduce_deg reduce_rad to_range opposite_deg opposite_rad
20             angle_s angle_c angle_c_rad diff_angle polar cart quad
21             ARCS/
22             ],
23             );
24             our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );
25             our $VERSION = 0.02;
26              
27 18340     18340 1 28125 sub frac($x) { ( modf($x) )[0] }
  18340         21062  
  18340         18299  
  18340         44184  
28              
29 0     0 1 0 sub frac360($x) { frac($x) * 360 }
  0         0  
  0         0  
  0         0  
30              
31 45     45 1 2013 sub dms ( $x, $places = 3 ) {
  45         49  
  45         62  
  45         47  
32 45 100       102 return $x if $places == 1;
33              
34 29         70 my ( $f, $i ) = modf($x);
35 29 50 66     96 $f = -$f if $i != 0 && $f < 0;
36              
37 29         58 ( $i, dms( $f * 60, $places - 1 ) );
38             }
39              
40 0     0 1 0 sub hms { dms @_ }
41              
42 12     12 1 10878 sub zdms($x) {
  12         17  
  12         14  
43 12         19 my ( $d, $m, $s ) = dms($x);
44 12         25 my $z = int( $d / 30 );
45 12         17 $d %= 30;
46              
47 12         28 $z, $d, $m, $s;
48             }
49              
50 110     110 1 6440 sub ddd(@args) {
  110         202  
  110         122  
51 110     241   660 my $b = any { $_ < 0 } @args;
  241         381  
52 110 100       417 my $sgn = $b ? -1 : 1;
53 110   100     239 my ( $d, $m, $s ) = map { abs( $args[$_] || 0 ) } ( 0 .. 2 );
  330         1008  
54 110         863 return $sgn * ( $d + ( $m + $s / 60.0 ) / 60.0 );
55             }
56              
57 83     83 1 3307 sub polynome ( $t, @terms ) {
  83         123  
  83         135  
  83         98  
58 83     174   548 reduce { $a * $t + $b } reverse @terms;
  174         621  
59             }
60              
61 74     74 1 4449 sub to_range ( $x, $limit ) {
  74         87  
  74         84  
  74         87  
62 74         162 $x = fmod( $x, $limit );
63 74 50       238 $x < 0 ? $x + $limit : $x;
64             }
65              
66             #sub reduce_deg($x) { to_range( $x, 360 ) }
67              
68 920     920 1 4559 sub reduce_deg($x) {
  920         1179  
  920         1073  
69 920         1621 my $res = Math::Trig::deg2deg($x);
70 920 100       5381 $res < 0 ? $res + 360 : $res;
71             }
72              
73             #sub reduce_rad($x) { to_range( $x, pi2 ) }
74              
75 1385     1385 1 7046 sub reduce_rad($x) {
  1385         2063  
  1385         1692  
76 1385         3431 my $res = Math::Trig::rad2rad($x);
77 1385 100       11075 $res < 0 ? $res + pi2 : $res;
78             }
79              
80 4600     4600 1 8199 sub sine($x) { sin( pi2 * frac($x) ) }
  4600         5406  
  4600         4995  
  4600         6005  
81              
82 4     4 1 4341 sub opposite_deg($x) { reduce_deg( $x + 180 ) }
  4         7  
  4         5  
  4         6  
83              
84 4     4 1 1578 sub opposite_rad($x) { reduce_rad( $x + pi ) }
  4         7  
  4         5  
  4         9  
85              
86 27     27 1 6378 sub angle_c ( $a, $b ) {
  27         33  
  27         33  
  27         30  
87 27         46 my $x = abs( $a - $b );
88 27 100       101 $x > 180 ? 360 - $x : $x;
89             }
90              
91 4     4 1 1669 sub angle_c_rad ( $a, $b ) {
  4         7  
  4         5  
  4         5  
92 4         7 my $x = abs( $a - $b );
93 4 50       12 $x > pi ? pi2 - $x : $x;
94             }
95              
96             sub angle_s {
97 1     1 1 2039 my ( $x1, $y1, $x2, $y2 ) = map { deg2rad $_ } @_;
  4         24  
98 1         23 rad2deg(
99             acos( sin($y1) * sin($y2) + cos($y1) * cos($y2) * cos( $x1 - $x2 ) ) );
100             }
101              
102 18     18 1 2596 sub diff_angle($a, $b, $mode = 'degrees') {
  18         22  
  18         21  
  18         29  
  18         19  
103 18         31 my $m = lc $mode;
104 18 0       39 my $whole = $m eq 'degrees' ? 360
    50          
105             : $m eq 'radians' ? pi2
106             : undef;
107 18 50       35 die "Expected 'degrees' or 'radians' mode" unless $whole;
108 18 50       37 my $half = $m eq 'degrees' ? 180 : pi;
109 18 100       38 my $x = $b < $a ? $b + $whole : $b;
110 18         23 $x -= $a;
111 18 100       46 return $x - $whole if $x > $half;
112 9         19 return $x;
113             }
114              
115              
116 220     220 1 7055 sub cart( $r, $theta, $phi ) {
  220         417  
  220         414  
  220         497  
  220         400  
117 220         594 my $rcst = $r * cos($theta);
118 220         908 $rcst * cos($phi), $rcst * sin($phi), $r * sin($theta);
119             }
120              
121             # in previous versions was named 'polar'
122 721     721 1 7330 sub polar ( $x, $y, $z ) {
  721         929  
  721         925  
  721         961  
  721         877  
123 721         1105 my $rho = $x * $x + $y * $y;
124 721         1135 my $r = sqrt( $rho + $z * $z );
125 721         2360 my $phi = atan2( $y, $x );
126 721 100       1811 $phi += pi2 if $phi < 0;
127 721         1187 $rho = sqrt($rho);
128 721         1221 my $theta = atan2( $z, $rho );
129 721         1856 $r, $theta, $phi;
130             }
131              
132             sub quad {
133 312     312 1 5897 my ( $y_minus, $y_0, $y_plus ) = @_;
134 312         506 my $nz = 0;
135 312         555 my $a = 0.5 * ( $y_minus + $y_plus ) - $y_0;
136 312         537 my $b = 0.5 * ( $y_plus - $y_minus );
137 312         456 my $c = $y_0;
138              
139 312         632 my $xe = -$b / ( 2 * $a );
140 312         628 my $ye = ( $a * $xe + $b ) * $xe + $c;
141 312         670 my $dis = $b * $b - 4 * $a * $c; # discriminant of y = axx+bx+c
142 312         476 my @zeroes;
143 312 100       854 if ( $dis >= 0 ) {
144              
145             # parabola intersects x-axis
146 305         520 my $dx = 0.5 * sqrt($dis) / abs($a);
147 305         825 @zeroes[ 0, 1 ] = ( $xe - $dx, $xe + $dx );
148 305 100       791 $nz++ if abs( $zeroes[0] ) <= 1;
149 305 100       699 $nz++ if abs( $zeroes[1] ) <= 1;
150 305 100       709 $zeroes[0] = $zeroes[1] if $zeroes[0] < -1;
151             }
152 312         1162 $nz, $xe, $ye, @zeroes;
153             }
154              
155              
156              
157             1; # End of Astro::Montenbruck::Core::MathUtils
158              
159             __END__
160              
161             =pod
162              
163             =encoding UTF-8
164              
165             =head1 NAME
166              
167             Astro::Montenbruck::Core::MathUtils - Core mathematical routines used by Astro::Montenbruck modules.
168              
169             =head1 VERSION
170              
171             Version 0.01
172              
173              
174             =head1 SYNOPSIS
175              
176             use Astro::Montenbruck::Core::MathUtils qw/dms/;
177              
178             my ($d, $m, $s) = dms(55.75); # (55, 45, 0)
179             ...
180              
181             =head1 EXPORT
182              
183             =over
184              
185             =item * L</frac($x)>
186              
187             =item * L</frac360($x)>
188              
189             =item * L</dms($x)>
190              
191             =item * L</hms($x)>
192              
193             =item * L</zdms($x)>
194              
195             =item * L</ddd($deg[, $min[, $sec]])>
196              
197             =item * L</polynome($t, @terms)>
198              
199             =item * L</to_range($x, $range)>
200              
201             =item * L</reduce_deg($x)>
202              
203             =item * L</reduce_rad($x)>
204              
205             =item * L</opposite_deg($x)>
206              
207             =item * L</opposite_rad($x)>
208              
209             =item * L</angle_c($x, $y)>
210              
211             =item * L</angle_c_rad($x, $y)>
212              
213             =item * L</angle_c_rad($x, $y)>
214              
215             =item * L</angle_s($x1, $y1, $x2, $y2)>
216              
217             =item * L</diff_angle($a, $b, $mode='degrees')>
218              
219             =item * L</diff_angle($a, $b, $mode='degrees')>
220              
221             =item * L</cart($r, $theta, $phi)>
222              
223             =item * L</polar($x, $y, $z)>
224              
225             =back
226              
227              
228             =head1 SUBROUTINES
229              
230              
231             =head2 frac($x)
232              
233             Fractional part of a decimal number.
234              
235              
236             =head2 frac360($x)
237              
238             Range function, similar to L<to_range($x, $range)>, used with polinomial function for better accuracy.
239              
240              
241             =head2 dms($x)
242              
243             Given decimal hours (or degrees), return nearest hours (or degrees), int,
244             minutes, int, and seconds, float.
245              
246             =head3 Positional arguments:
247              
248             =over
249              
250             =item * decimal value, 0..360 for angular mode, 0..24 for time
251              
252             =back
253              
254             =head3 Named arguments:
255              
256             =over
257              
258             =item * B<places> (optional) amount of required sexagesimal values to be returned (1-3);
259             default = 3 (degrees/hours, minutes, seconds)
260              
261             =back
262              
263             =head3 Returns:
264              
265             =over
266              
267             =item * array of degrees (int), minutes (int), seconds (float)
268              
269             =back
270              
271              
272             =head2 hms($x)
273              
274             Alias for L</dms>
275              
276             =head2 zdms($x)
277              
278             Converts decimal degrees to zodiac sign number (zero based), zodiac degrees, minutes and seconds.
279              
280             =head3 Positional arguments:
281              
282             =over
283              
284             =item * decimal value, 0..360 for angular mode, 0..24 for time
285              
286             =back
287              
288             =head3 Returns:
289              
290             =over
291              
292             =item * array of zodiac sign (0-11), degrees (int), minutes (int), seconds (float)
293              
294             =back
295              
296              
297             =head2 ddd($deg[, $min[, $sec]])
298              
299             Converts sexagesimal values to decimal.
300              
301             =head3 Arguments
302              
303             =over
304              
305             1 to 3 sexagesimal values, such as: degrees, minutes and
306             seconds, or degrees and minutes, or just degrees:
307              
308             =over
309              
310             =item * C<ddd(11)>
311              
312             =item * C<ddd(11, 46)>
313              
314             =item * C<ddd(11, 46, 20)>
315              
316             =back
317              
318             If any non-zero argument is negative, the result is negative.
319              
320             =over
321              
322             =item * C<ddd(-11, 46, 0) = -11.766666666666667>
323              
324             =item * C<ddd(11, -46, 0) = 11.766666666666667>
325              
326             =back
327              
328             Negative sign in wrong position is ignored.
329              
330             =back
331              
332             =head3 Returns:
333              
334             =over
335              
336             =item * decimal (degrees or hours)
337              
338             =back
339              
340              
341             =head2 polynome($t, @terms)
342              
343             Calculates polynome: $a1 + $a2*$t + $a3*$t*$t + $a4*$t*$t*$t...
344              
345             =head3 Arguments
346              
347             =over
348              
349             =item * $t coefficient, in astronomical routines usually time in centuries
350              
351             =item * any number of decimal values
352              
353             =back
354              
355             =head3 Returns:
356              
357             =over
358              
359             =item * decimal number
360              
361             =back
362              
363              
364              
365             =head2 to_range($x, $range)
366              
367             Reduces $x to 0 >= $x < $range
368              
369             =head3 Arguments
370              
371             =over
372              
373             =item * number to reduce
374              
375             =item * limit (non-inclusive), e.g: 360 for degrees, 24 for hours
376              
377             =back
378              
379             =head3 Returns
380              
381             =over
382              
383             =item * number
384              
385             =back
386              
387              
388              
389             =head2 reduce_deg($x)
390              
391             Reduces $x to 0 >= $x < 360
392              
393              
394             =head2 reduce_rad($x)
395              
396             Reduces $x to 0 >= $x < pi2
397              
398             =head2 opposite_deg($x)
399              
400             Returns opposite degree.
401              
402              
403             =head2 opposite_rad($x)
404              
405             Returns opposite radian.
406              
407             =head2 angle_c($x, $y)
408              
409             Calculate shortest arc in dergees between $x and $y.
410              
411             =head2 angle_c_rad($x, $y)
412              
413             Calculates shortest arc in radians between $x and $y.
414              
415             =head2 angle_s($x1, $y1, $x2, $y2)
416              
417             Calculates arc between 2 points on a sphere.
418             Expected arguments: 2 pairs of coordinates (X, Y) of the 2 points.
419              
420             The coordinates may be ecliptic, equatorial or horizontal.
421              
422             =head2 diff_angle($a, $b, $mode='degrees')
423              
424             Return angle C<$b - $a>, accounting for circular values.
425              
426             Parameters $a and $b should be in the range 0..pi*2 or 0..360, depending on
427             optional B<$mode> argument. The result will be in the range I<-pi..pi> or I<-180..180>.
428             This allows us to directly compare angles which cross through 0:
429             I<359 degress... 0 degrees... 1 degree...> etc.
430              
431             =head3 Positional Arguments
432              
433             =over
434              
435             =item * B<$a> first angle, in radians or degrees
436              
437             =item * B<$b> second angle, in radians or degrees
438              
439             =back
440              
441             =head3 Named Arguments
442              
443             =over
444              
445             =item * B<$mode> C<"degrees"> (default) or C<"radians">, case insensitive.
446              
447             =back
448              
449              
450             =head2 sine($x)
451              
452             Calculate sin(phi); phi in units of 1 revolution = 360 degrees
453              
454             =head2 cart($r, $theta, $phi)
455              
456             Conversion of polar coordinates (r,theta,phi) into cartesian (x,y,z).
457              
458             =head3 Arguments
459              
460             =over
461              
462             =item * B<$r>, distance from the origin;
463              
464             =item * B<$theta> (in radians) corresponding to [-90 deg, +90 deg];
465              
466             =item * B<$phi> (in radians) corresponding to [-360 deg, +360 deg])
467              
468             =back
469              
470             =head3 Returns
471              
472             Rectangular coordinates:
473              
474             =over
475              
476             =item * B<$x>, X
477              
478             =item * B<$y>, Y
479              
480             =item * B<$z>, Z
481              
482              
483             =back
484              
485             =head2 polar($x, $y, $z)
486              
487             Conversion of cartesian coordinates (x,y,z) into polar (r,theta,phi).
488              
489             =head3 Arguments
490              
491             =over
492              
493             =item * B<$x>, X
494              
495             =item * B<$y>, Y
496              
497             =item * B<$z>, Z
498              
499             =back
500              
501             =head3 Returns
502              
503             =over
504              
505             =item * B<$r>, distance from the origin;
506              
507             =item * B<$theta> (in radians) corresponding to [-90 deg, +90 deg];
508              
509             =item * B<$phi> (in radians) corresponding to [-360 deg, +360 deg])
510              
511             =back
512              
513             =head2 quad($y_minus, $y_0, $y_plus)
514              
515             Quadratic interpolation
516              
517             Finds a parabola through 3 points C<(-1 , y_minus), (0, Y_0), (1, y_plus)>,
518             that do not lie on a straight line.
519              
520             =head3 Arguments
521              
522             Three y-values:
523              
524             =over
525              
526             =item * B<$y_minus> value of function at x = -1
527              
528             =item * B<$y_0> value of function at x = 0
529              
530             =item * B<$y_plus> value of function at x = 1
531              
532             =back
533              
534             =head3 Returns
535              
536             =over
537              
538             =item * B<$xe>, abscissa of extremum (may be outside C<[-1, 1]>)
539              
540             =item * B<$ye>, Value of function at xe
541              
542             =item * B<$root1>, first root found
543              
544             =item * B<$root2>, second root found
545              
546             =item * B<$n_root>, number of roots within the interval C<[-1, +1]>
547              
548             =back
549              
550              
551             =head1 AUTHOR
552              
553             Sergey Krushinsky, C<< <krushi at cpan.org> >>
554              
555             =head1 COPYRIGHT & LICENSE
556             Copyright (C) 2009-2022 by Sergey Krushinsky
557              
558             This library is free software; you can redistribute it and/or modify
559             it under the same terms as Perl itself.