1=head1 NAME 2 3Astro::Coord::ECI::Utils - Utility routines for astronomical calculations 4 5=head1 SYNOPSIS 6 7 use Astro::Coord::ECI::Utils qw{ :all }; 8 my $now = time (); 9 print "The current Julian day is ", julianday ($now); 10 11=head1 DESCRIPTION 12 13This module was written to provide a home for all the constants and 14utility subroutines used by B<Astro::Coord::ECI> and its descendants. 15What ended up here was anything that was essentially a subroutine, not 16a method. 17 18Because figuring out how to convert to and from Perl time bids fair to 19become complicated, this module is also responsible for figuring out how 20that is done, and exporting whatever is needful to export. See C<:time> 21below for the gory details. 22 23This package exports nothing by default. But all the constants, 24variables, and subroutines documented below are exportable, and the 25following export tags may be used: 26 27=over 28 29=item :all 30 31This imports everything exportable into your name space. 32 33=item :greg_time 34 35This imports all time routines except the discouraged routines 36C<time_gm()> and C<time_local()>. 37 38=item :mainstream 39 40This imports everything not deprecated into your name space. 41 42=item :params 43 44This imports the parameter validation routines C<__classisa> and 45C<__instance>. 46 47=item :ref 48 49This imports all the C<*_REF> constants. 50 51=item :time 52 53This imports the time routines into your name space. If 54L<Time::y2038|Time::y2038> is available, it will be loaded, and both 55this tag and C<:all> will import C<gmtime>, C<localtime>, C<time_gm>, 56C<time_local>, C<greg_time_gm>, and C<greg_time_local> into your name 57space. Otherwise, C<Time::Local|Time::Local> will be loaded, and both 58this tag and C<:all> will import C<time_gm>, C<time_local>, 59C<greg_time_gm>, and C<greg_time_local> into your name space. 60 61=item :vector 62 63This imports the vector arithmetic routines. It includes anything whose 64name begins with C<'vector_'>. 65 66=back 67 68Under Perl 5.6 you may find that, if you use any of the above tags, you 69need to specify them first in your import list. 70 71=head2 The following constants are exportable: 72 73 AU = number of kilometers in an astronomical unit 74 JD_OF_EPOCH = the Julian Day of Perl epoch 0 75 LIGHTYEAR = number of kilometers in a light year 76 PARSEC = number of kilometers in a parsec 77 PERL2000 = January 1 2000, 12 noon universal, in Perl time 78 PI = the circle ratio, computed as atan2 (0, -1) 79 PIOVER2 = half the circle ratio 80 SECSPERDAY = the number of seconds in a day 81 SECS_PER_SIDERIAL_DAY = seconds in a siderial day 82 SPEED_OF_LIGHT = speed of light in kilometers per second 83 TWOPI = twice the circle ratio 84 85 ARRAY_REF = 'ARRAY' 86 CODE_REF = 'CODE' 87 HASH_REF = 'HASH' 88 SCALAR_REF = 'SCALAR' 89 90=head2 The following global variables are exportable: 91 92=head3 $DATETIMEFORMAT 93 94This variable represents the POSIX::strftime format used to convert 95times to strings. The default value is '%a %b %d %Y %H:%M:%S' to be 96consistent with the behavior of gmtime (or, to be precise, the 97behavior of ctime as documented on my system). 98 99=head3 $JD_GREGORIAN 100 101This variable represents the Julian Day of the switch from Julian to 102Gregorian calendars. This is used by date2jd(), jd2date(), and the 103routines which depend on them, for deciding whether the date is to be 104interpreted as in the Julian or Gregorian calendar. Its initial setting 105is 2299160.5, which represents midnight October 15 1582 in the Gregorian 106calendar, which is the date that calendar was first adopted. This is 107slightly different than the value of 2299161 (noon of the same day) used 108by Jean Meeus. 109 110If you are interested in historical calculations, you may wish to reset 111this appropriately. If you use date2jd to calculate the new value, be 112aware of the effect the current setting of $JD_GREGORIAN has on the 113interpretation of the date you give. 114 115=head2 In addition, the following subroutines are exportable: 116 117=over 4 118 119=cut 120 121package Astro::Coord::ECI::Utils; 122 123use strict; 124use warnings; 125 126our $VERSION = '0.122'; 127our @ISA = qw{Exporter}; 128 129use Carp; 130## use Config; 131## use Data::Dumper; 132use POSIX qw{floor strftime}; 133use Scalar::Util qw{ blessed }; 134 135my @greg_time_routines; 136 137BEGIN { 138 139 # NOTE WELL 140 # 141 # The logic here should be consistent with the optional-module text 142 # emitted by inc/Astro/Coord/ECI/Recommend.pm. 143 # 144 145 eval { 146 require Time::y2038; 147 Time::y2038->import( qw{ gmtime localtime } ); 148 149 # sub time_gm 150 *time_gm = sub { 151 my @date = @_; 152 $date[5] = _year_adjust_y2038( $date[5] ); 153 return Time::y2038::timegm( @date ); 154 }; 155 # greg_time_local 156 *greg_time_gm = sub { 157 my @date = @_; 158 $date[5] -= 1900; 159 return Time::y2038::timegm( @date ); 160 }; 161 162 # sub time_local 163 *time_local = sub { 164 my @date = @_; 165 $date[5] = _year_adjust_y2038( $date[5] ); 166 return Time::y2038::timelocal( @date ); 167 }; 168 # sub greg_time_local 169 *greg_time_local = sub { 170 my @date = @_; 171 $date[5] -= 1900; 172 return Time::y2038::timelocal( @date ); 173 }; 174 175 @greg_time_routines = qw{ 176 gmtime localtime 177 greg_time_gm greg_time_local 178 __tle_year_to_Gregorian_year 179 }; 180 181 1; 182 } or do { 183 require Time::Local; 184 185 # sub time_gm 186 *time_gm = Time::Local->can( 'timegm' ); 187 # sub greg_time_gm 188 *greg_time_gm = Time::Local->can( 'timegm_modern' ) || sub { 189 my @date = @_; 190 $date[5] = _year_adjust_greg( $date[5] ); 191 return Time::Local::timegm( @date ); 192 }; 193 194 # sub time_local 195 *time_local = Time::Local->can( 'timelocal' ); 196 # sub greg_time_local 197 *greg_time_local = Time::Local->can( 'timelocal_modern' ) || sub { 198 my @date = @_; 199 $date[5] = _year_adjust_greg( $date[5] ); 200 return Time::Local::timelocal( @date ); 201 }; 202 203 @greg_time_routines = qw{ 204 greg_time_gm greg_time_local 205 __tle_year_to_Gregorian_year 206 }; 207 208 }; 209} 210 211# This subroutine is used to convert year numbers to Perl years in 212# accordance with the documentation in the 5.24.0 version of 213# Time::Local. It is intended to be called by the Time::y2038 code, 214# which expects Perl years. 215 216{ 217 # The following code is lifted verbatim from Time::Local 1.25. 218 # Because that code bases the window used for expanding two-digit 219 # years on the local year as of the time the module was loaded, I do 220 # too. 221 222 my $ThisYear = ( localtime() )[5]; 223 my $Breakpoint = ( $ThisYear + 50 ) % 100; 224 my $NextCentury = $ThisYear - $ThisYear % 100; 225 $NextCentury += 100 if $Breakpoint < 50; 226 my $Century = $NextCentury - 100; 227 228 # The above code is lifted verbatim from Time::Local 1.25. 229 230 use constant NOT_GREG => 231 '%d not interpreted as Gregorian year by Time::Local::timegm'; 232 233 # Adujst the year so that the Time::y2038 implementation of 234 # time_gm() and time_local() mirrors the Time::Local timegm() and 235 # timelocal() behavior. Kinda sorta. 236 sub _year_adjust_y2038 { 237 my ( $year ) = @_; 238 239 $year < 0 240 and return $year; 241 242 $year >= 1000 243 and return $year - 1900; 244 245 # The following line of code is lifted verbatim from Time::Local 246 # 1.25. 247 $year += ( $year > $Breakpoint ) ? $Century : $NextCentury; 248 249 return $year; 250 } 251} 252 253# Adjust a Gregorian year so that Time::Local timegm() and timelocal() 254# return epochs in that year. 255sub _year_adjust_greg { 256 my ( $year ) = @_; 257 return $year >= 1000 ? $year : $year - 1900; 258} 259 260our @CARP_NOT = qw{ 261 Astro::Coord::ECI 262 Astro::Coord::ECI::Mixin 263 Astro::Coord::ECI::Moon 264 Astro::Coord::ECI::Star 265 Astro::Coord::ECI::Sun 266 Astro::Coord::ECI::TLE 267 Astro::Coord::ECI::TLE::Set 268 Astro::Coord::ECI::Utils 269}; 270 271our @EXPORT; 272my @all_external = ( qw{ 273 AU $DATETIMEFORMAT $JD_GREGORIAN JD_OF_EPOCH LIGHTYEAR PARSEC 274 PERL2000 PI PIOVER2 SECSPERDAY SECS_PER_SIDERIAL_DAY 275 SPEED_OF_LIGHT TWOPI 276 ARRAY_REF CODE_REF HASH_REF SCALAR_REF 277 acos add_magnitudes asin 278 atmospheric_extinction date2epoch date2jd 279 decode_space_track_json_time deg2rad distsq dynamical_delta 280 embodies epoch2datetime find_first_true 281 fold_case format_space_track_json_time intensity_to_magnitude 282 jcent2000 jd2date jd2datetime jday2000 julianday 283 keplers_equation load_module looks_like_number max min mod2pi 284 omega 285 position_angle 286 rad2deg rad2dms rad2hms tan theta0 thetag vector_cross_product 287 vector_dot_product vector_magnitude vector_unitize __classisa 288 __default_station __instance __subroutine_deprecation 289 __sprintf 290 }, 291 qw{ time_gm time_local }, @greg_time_routines ); 292our @EXPORT_OK = ( 293 qw{ @CARP_NOT }, # Package-private, undocumented 294 @all_external, 295); 296 297my %deprecated_export = map { $_ => 1 } qw{ 298}; 299 300our %EXPORT_TAGS = ( 301 all => \@all_external, 302 greg_time => \@greg_time_routines, 303 mainstream => [ grep { ! $deprecated_export{$_} } @all_external ], 304 params => [ qw{ __classisa __instance } ], 305 ref => [ grep { m/ [[:upper:]]+ _REF \z /smx } @all_external ], 306 time => [ qw{ time_gm time_local }, @greg_time_routines ], 307 vector => [ grep { m/ \A vector_ /smx } @all_external ], 308); 309 310use constant AU => 149597870; # 1 astronomical unit, per 311 # Meeus, Appendix I pg 407. 312use constant LIGHTYEAR => 9.4607e12; # 1 light-year, per Meeus, 313 # Appendix I pg 407. 314use constant PARSEC => 30.8568e12; # 1 parsec, per Meeus, 315 # Appendix I pg 407. 316use constant PERL2000 => greg_time_gm( 0, 0, 12, 1, 0, 2000 ); 317use constant PI => atan2 (0, -1); 318use constant PIOVER2 => PI / 2; 319use constant SECSPERDAY => 86400; 320use constant SECS_PER_SIDERIAL_DAY => 86164.0905; # Appendix I, page 408. 321use constant SPEED_OF_LIGHT => 299792.458; # KM/sec, per NIST. 322### use constant SOLAR_RADIUS => 1392000 / 2; # Meeus, Appendix I, page 407. 323use constant TWOPI => PI * 2; 324 325use constant ARRAY_REF => ref []; 326use constant CODE_REF => ref sub {}; 327use constant HASH_REF => ref {}; 328use constant SCALAR_REF => ref \0; 329 330=item $angle = acos ($value) 331 332This subroutine calculates the arc in radians whose cosine is the given 333value. 334 335=cut 336 337sub acos { 338 abs ($_[0]) > 1 and confess <<eod; 339Programming error - Trying to take the arc cosine of a number greater 340 than 1. 341eod 342 return atan2 (sqrt (1 - $_[0] * $_[0]), $_[0]) 343} 344 345=item $mag = add_magnitudes( $mag1, $mag2, ... ); 346 347This subroutine computes the total magnitude of a list of individual 348magnitudes. The algorithm comes from Jean Meeus' "Astronomical 349Algorithms", Second Edition, Chapter 56, Page 393. 350 351=cut 352 353sub add_magnitudes { 354 my @arg = @_; 355 my $sum = 0; 356 foreach my $mag ( @arg ) { 357 $sum += 10 ** ( -0.4 * $mag ); 358 } 359 return -2.5 * log( $sum ) / log( 10 ); 360} 361 362=item $angle = asin ($value) 363 364This subroutine calculates the arc in radians whose sine is the given 365value. 366 367=cut 368 369sub asin {return atan2 ($_[0], sqrt (1 - $_[0] * $_[0]))} 370 371=item $magnitude = atmospheric_extinction ($elevation, $height); 372 373This subroutine calculates the typical atmospheric extinction in 374magnitudes at the given elevation above the horizon in radians and the 375given height above sea level in kilometers. 376 377The algorithm comes from Daniel W. E. Green's article "Magnitude 378Corrections for Atmospheric Extinction", which was published in 379the July 1992 issue of "International Comet Quarterly", and is 380available online at 381L<http://www.icq.eps.harvard.edu/ICQExtinct.html>. The text of 382this article makes it clear that the actual value of the 383atmospheric extinction can vary greatly from the typical 384values given even in the absence of cloud cover. 385 386=cut 387 388# Note that the "constant" 0.120 in Aaer (aerosol scattering) is 389# based on a compromise value A0 = 0.050 in Green's equation 3 390# (not exhibited here), which can vary from 0.035 in the winter to 391# 0.065 in the summer. This makes a difference of a couple tenths 392# at 20 degrees elevation, but a couple magnitudes at the 393# horizon. Green also remarks that the 1.5 denominator in the 394# same equation (a.k.a. the scale height) can be up to twice 395# that. 396 397sub atmospheric_extinction { 398 my ($elevation, $height) = @_; 399 my $cosZ = cos (PIOVER2 - $elevation); 400 my $X = 1/($cosZ + 0.025 * exp (-11 * $cosZ)); # Green 1 401 my $Aray = 0.1451 * exp (-$height / 7.996); # Green 2 402 my $Aaer = 0.120 * exp (-$height / 1.5); # Green 4 403 return ($Aray + $Aaer + 0.016) * $X; # Green 5, 6 404} 405 406=item $jd = date2jd ($sec, $min, $hr, $day, $mon, $yr) 407 408This subroutine converts the given date to the corresponding Julian day. 409The inputs are a Perl date and time; $mon is in the range 0 - 41011, and $yr is from 1900, with earlier years being negative. The year 1 411BC is represented as -1900. 412 413If less than 6 arguments are provided, zeroes will be prepended to the 414argument list as needed. 415 416The date is presumed to be in the Gregorian calendar. If the resultant 417Julian Day is before $JD_GREGORIAN, the date is reinterpreted as being 418from the Julian calendar. 419 420The only validation is that the month be between 0 and 11 inclusive, and 421that the year be not less than -6612 (4713 BC). Fractional days are 422accepted. 423 424The algorithm is from Jean Meeus' "Astronomical Algorithms", second 425edition, chapter 7 ("Julian Day"), pages 60ff, but the month is 426zero-based, not 1-based, and years are 1900-based. 427 428=cut 429 430our $DATETIMEFORMAT; 431our $JD_GREGORIAN; 432BEGIN { 433 $DATETIMEFORMAT = '%a %b %d %Y %H:%M:%S'; 434 $JD_GREGORIAN = 2299160.5; 435} 436 437sub date2jd { 438 my @args = @_; 439 unshift @args, 0 while @args < 6; 440 my ($sec, $min, $hr, $day, $mon, $yr) = @args; 441 $mon++; # Algorithm expects month 1-12. 442 $yr += 1900; # Algorithm expects zero-based year. 443 ($yr < -4712) and croak "Error - Invalid year $yr"; 444 ($mon < 1 || $mon > 12) and croak "Error - Invalid month $mon"; 445 if ($mon < 3) { 446 --$yr; 447 $mon += 12; 448 } 449 my $A = floor ($yr / 100); 450 my $B = 2 - $A + floor ($A / 4); 451 my $jd = floor (365.25 * ($yr + 4716)) + 452 floor (30.6001 * ($mon + 1)) + $day + $B - 1524.5 + 453 ((($sec || 0) / 60 + ($min || 0)) / 60 + ($hr || 0)) / 24; 454 $jd < $JD_GREGORIAN and 455 $jd = floor (365.25 * ($yr + 4716)) + 456 floor (30.6001 * ($mon + 1)) + $day - 1524.5 + 457 ((($sec || 0) / 60 + ($min || 0)) / 60 + ($hr || 0)) / 24; 458 return $jd; 459} 460 461use constant JD_OF_EPOCH => date2jd (gmtime (0)); 462 463=item $epoch = date2epoch ($sec, $min, $hr, $day, $mon, $yr) 464 465This is a convenience routine that converts the given date to seconds 466since the epoch, going through date2jd() to do so. The arguments are the 467same as those of date2jd(). 468 469If less than 6 arguments are provided, zeroes will be prepended to the 470argument list as needed. 471 472The functionality is the same as B<Time::Local::timegm>, but this 473function lacks timegm's limited date range under Perls before 5.12.0. If 474you have Perl 5.12.0 or better, the core L<Time::Local|Time::Local> 475C<timegm()> will probably do what you want. If you have an earlier 476Perl, L<Time::y2038|Time::y2038> C<timegm()> may do what you want. 477 478=cut 479 480sub date2epoch { 481 my @args = @_; 482 unshift @args, 0 while @args < 6; 483 my ($sec, $min, $hr, $day, $mon, $yr) = @args; 484 return (date2jd ($day, $mon, $yr) - JD_OF_EPOCH) * SECSPERDAY + 485 (($hr || 0) * 60 + ($min || 0)) * 60 + ($sec || 0); 486} 487 488=item $time = decode_space_track_json_time( $string ) 489 490This subroutine decodes a time in the format Space Track uses in their 491JSON code. This is ISO-8601-ish, but with a possible fractional part and 492no zone. 493 494=cut 495 496sub decode_space_track_json_time { 497 my ( $string ) = @_; 498 $string =~ m{ \A \s* 499 ( [0-9]+ ) [^0-9]+ ( [0-9]+ ) [^0-9]+ ( [0-9]+ ) [^0-9]+ 500 ( [0-9]+ ) [^0-9]+ ( [0-9]+ ) [^0-9]+ ( [0-9]+ ) (?: ( [.] [0-9]* ) )? \s* \z }smx 501 or return; 502 my @time = ( $1, $2, $3, $4, $5, $6 ); 503 my $frac = $7; 504 $time[0] = __tle_year_to_Gregorian_year( $time[0] ); 505 $time[1] -= 1; 506 my $rslt = greg_time_gm( reverse @time ); 507 defined $frac 508 and $frac ne '.' 509 and $rslt += $frac; 510 return $rslt; 511} 512 513# my ( $self, $station, @args ) = __default_station( @_ ) 514# 515# This exportable subroutine checks whether the second argument embodies 516# an Astro::Coord::ECI object. If so, the arguments are returned 517# verbatim. If not, the 'station' attribute of the invocant is inserted 518# after the first argument and the result returned. If the 'station' 519# attribute of the invocant has not been set, an exception is thrown. 520 521sub __default_station { 522 my ( $self, @args ) = @_; 523 if ( ! embodies( $args[0], 'Astro::Coord::ECI' ) ) { 524 my $sta = $self->get( 'station' ) 525 or croak 'Station attribute not set'; 526 unshift @args, $sta; 527 } 528 return ( $self, @args ); 529} 530 531=item $rad = deg2rad ($degr) 532 533This subroutine converts degrees to radians. If the argument is 534C<undef>, C<undef> will be returned. 535 536=cut 537 538sub deg2rad { return defined $_[0] ? $_[0] * PI / 180 : undef } 539 540=item $value = distsq (\@coord1, \@coord2) 541 542This subroutine calculates the square of the distance between the two 543sets of Cartesian coordinates. We do not take the square root here 544because of cases (e.g. the law of cosines) where we would just have 545to square the result again. 546 547B<Notice> that the subroutine does B<not> assume three-dimensional 548coordinates. If @coord1 and @coord2 have six entries, you will get a 549six-dimensional distance. 550 551=cut 552 553sub distsq { 554 my ( $x, $y ) = @_; 555 ARRAY_REF eq ref $x 556 and ARRAY_REF eq ref $y 557 and @{ $x } == @{ $y } 558 or confess <<'EOD'; 559Programming error - Both arguments to distsq must be references to 560 arrays of the same length. 561EOD 562 563 my $sum = 0; 564 my $size = @$x; 565 for (my $inx = 0; $inx < $size; $inx++) { 566 my $delta = $x->[$inx] - $y->[$inx]; 567 $sum += $delta * $delta; 568 } 569 return $sum 570} 571 572=item $seconds = dynamical_delta ($time); 573 574This method returns the difference between dynamical and universal time 575at the given universal time. That is, 576 577 $dynamical = $time + dynamical_delta ($time) 578 579if $time is universal time. 580 581The algorithm is from Jean Meeus' "Astronomical Algorithms", 2nd 582Edition, Chapter 10, page 78. Meeus notes that this is actually an 583observed quantity, and the algorithm is an approximation. 584 585=cut 586 587sub dynamical_delta { 588 my ($time) = @_; 589 my $year = (gmtime $time)[5] + 1900; 590 my $t = ($year - 2000) / 100; 591 my $correction = .37 * ($year - 2100); # Meeus' correction to (10.2) 592 return (25.3 * $t + 102) * $t + 102 # Meeus (10.2) 593 + $correction; # Meeus' correction. 594} 595 596=item $boolean = embodies ($thingy, $class) 597 598This subroutine represents a safe way to call the 'represents' method on 599$thingy. You get back true if and only if $thingy->can('represents') 600does not throw an exception and returns true, and 601$thingy->represents($class) returns true. Otherwise it returns false. 602Any exception is trapped and dismissed. 603 604This subroutine is called 'embodies' because it was too confusing to 605call it 'represents', both for the author and for the Perl interpreter. 606 607=cut 608 609sub embodies { 610 my ($thingy, $class) = @_; 611 my $code = eval {$thingy->can('represents')}; 612 return $code ? $code->($thingy, $class) : undef; 613} 614 615=item ($sec, $min, $hr, $day, $mon, $yr, $wday, $yday, 0) = epoch2datetime ($epoch) 616 617This convenience subroutine converts the given time in seconds from the 618system epoch to the corresponding date and time. It is implemented in 619terms of jd2date (), with the year and month returned from that 620subroutine. The day is a whole number, with the fractional part 621converted to hours, minutes, and seconds. The $wday is the day of the 622week, with Sunday being 0. The $yday is the day of the year, with 623January 1 being 0. The trailing 0 is the summer time (or daylight saving 624time) indicator which is always 0 to be consistent with gmtime. 625 626If called in scalar context, it returns the date formatted by 627POSIX::strftime, using the format string in $DATETIMEFORMAT. 628 629The functionality is the same as the core C<gmtime()>, but this function 630lacks gmtime's limited date range under Perls before 5.12.0. If you have 631Perl 5.12.0 or better, the core C<gmtime()> will probably do what you 632want. If you have an earlier Perl, L<Time::y2038|Time::y2038> 633C<gmtime()> may do what you want. 634 635The input must convert to a non-negative Julian date. The exact lower 636limit depends on the system, but is computed by -(JD_OF_EPOCH * 86400). 637For Unix systems with an epoch of January 1 1970, this is -210866760000. 638 639Additional algorithms for day of week and day of year come from Jean 640Meeus' "Astronomical Algorithms", 2nd Edition, Chapter 7 (Julian Day), 641page 65. 642 643=cut 644 645sub epoch2datetime { 646 my ($time) = @_; 647 my $day = floor ($time / SECSPERDAY); 648 my $sec = $time - $day * SECSPERDAY; 649 ($day, my $mon, my $yr, undef, my $leap) = jd2date ( 650 my $jd = $day + JD_OF_EPOCH); 651 $day = floor ($day + .5); 652 my $min = floor ($sec / 60); 653 $sec = $sec - $min * 60; 654 my $hr = floor ($min / 60); 655 $min = $min - $hr * 60; 656 my $wday = ($jd + 1.5) % 7; 657 my $yd = floor (275 * ($mon + 1) / 9) - (2 - $leap) * floor (($mon + 658 10) / 12) + $day - 31; 659 wantarray and return ($sec, $min, $hr, $day, $mon, $yr, $wday, $yd, 660 0); 661 return strftime ($DATETIMEFORMAT, $sec, $min, $hr, $day, $mon, $yr, 662 $wday, $yd, 0); 663} 664 665=item $time = find_first_true ($start, $end, \&test, $limit); 666 667This function finds the first time between $start and $end for which 668test ($time) is true. The resolution is $limit, which defaults to 1 if 669not specified. If the times are reversed (i.e. the start time is after 670the end time) the time returned is the last time test ($time) is true. 671 672The C<test()> function is called with the Perl time as its only 673argument. It is assumed to be false for the first part of the interval, 674and true for the rest. If this assumption is violated, the result of 675this subroutine should be considered meaningless. 676 677The calculation is done by, essentially, a binary search; the interval 678is repeatedly split, the function is evaluated at the midpoint, and a 679new interval selected based on whether the result is true or false. 680 681Actually, nothing in this function says the independent variable has to 682be time. 683 684=cut 685 686sub find_first_true { 687 my ($begin, $end, $test, $limit) = @_; 688 $limit ||= 1; 689 defined $begin 690 or confess 'Programming error - $begin undefined'; 691 defined $end 692 or confess 'Programming error - $end undefined'; 693 if ($limit >= 1) { 694 if ($begin <= $end) { 695 $begin = floor ($begin); 696 $end = floor ($end) == $end ? $end : floor ($end) + 1; 697 } else { 698 $end = floor ($end); 699 $begin = floor ($begin) == $begin ? $begin : floor ($begin) + 1; 700 } 701 } 702 my $iterator = ( 703 $end > $begin ? 704 sub {$end - $begin > $limit} : 705 sub {$begin - $end > $limit} 706 ); 707 while ($iterator->()) { 708 my $mid = $limit >= 1 ? 709 floor (($begin + $end) / 2) : ($begin + $end) / 2; 710 ($begin, $end) = ($test->($mid)) ? 711 ($begin, $mid) : ($mid, $end); 712 } 713 return $end; 714} 715 716=item $folded = fold_case( $text ); 717 718This function folds the case of its input, kinda sorta. It maps to 719C<CORE::fc> if that is available, otherwise it maps to C<CORE::lc>. 720 721=cut 722 723*fold_case = CORE->can( 'fc' ) || sub ($) { return lc $_[0] }; 724 725=item $fmtd = format_space_track_json_time( time() ) 726 727This function takes as input a Perl time, and returns that time 728in a format consistent with the Space Track JSON data. This is 729ISO-8601-ish, in Universal time, but without the zone indicated. 730 731=cut 732 733sub format_space_track_json_time { 734 my ( $time ) = @_; 735 defined $time 736 and $time =~ m/ \S /smx 737 or return; 738 my @parts = gmtime floor( $time + .5 ); 739 $parts[4] += 1; 740 $parts[5] += 1900; 741 return sprintf '%04d-%02d-%02d %02d:%02d:%02d', reverse 742 @parts[ 0 .. 5 ]; 743} 744 745=item $epoch = greg_time_gm( $sec, $min, $hr, $day, $mon, $yr ); 746 747This exportable subroutine is a wrapper for either 748C<Time::y2038::timegm()> (if that module is installed), 749C<Time::Local::timegm_modern()> (if that is available), or 750C<Time::Local::timegm()> (if not.) 751 752This subroutine interprets years as Gregorian years. 753 754The difference between this and C<time_gm()> is that C<time_gm()> 755interprets the year the way C<Time::Local::timegm()> does. For that 756reason, this subroutine is preferred over c<time_gm()>. 757 758This wrapper is needed because the routines have subtly different 759signatures. L<Time::y2038|Time::y2038> C<timegm()> interprets years 760strictly as Perl years. L<Time::Local|Time::Local> C<timegm_modern()> 761interprets them strictly as Gregorian years. L<Time::Local|Time::Local> 762C<timegm()> interprets them differently depending on the value of the 763year. Years greater than or equal to 1000 are Gregorian years, but all 764others are Perl years, except for the range 0 to 99 inclusive, which are 765within 50 years of the current year. 766 767If you are doing historical calculations, see 768L<Historical Calculations|Astro::Coord::ECI::Sun/Historical Calculations> 769in the L<Astro::Coord::ECI::Sun|Astro::Coord::ECI::Sun> documentation 770for a discussion of input and output time conversion. 771 772=item $epoch = greg_time_local( $sec, $min, $hr, $day, $mon, $yr ); 773 774This exportable subroutine is a wrapper for either 775C<Time::y2038::timelocal()> (if that module is installed), 776C<Time::Local::timelocal_modern()> (if that is available), or 777C<Time::Local::timelocal()> (if not.) 778 779This subroutine interprets years as Gregorian years. 780 781The difference between this and c<time_local()> is that C<time_local()> 782interprets the year the way C<Time::Local::timelocal()> does. For that 783reason, this subroutine is preferred over c<time_local()>. 784 785This wrapper is needed for the same reason C<greg_time_gm()> is 786needed. 787 788If you are doing historical calculations, see 789L<Historical Calculations|Astro::Coord::ECI::Sun/Historical Calculations> 790in the L<Astro::Coord::ECI::Sun|Astro::Coord::ECI::Sun> documentation 791for a discussion of input and output time conversion. 792 793=item $difference = intensity_to_magnitude ($ratio) 794 795This function converts a ratio of light intensities to a difference in 796stellar magnitudes. The algorithm comes from Jean Meeus' "Astronomical 797Algorithms", Second Edition, Chapter 56, Page 395. 798 799Note that, because of the way magnitudes work (a more negative number 800represents a brighter star) you get back a positive result for an 801intensity ratio less than 1, and a negative result for an intensity 802ratio greater than 1. 803 804=cut 805 806{ # Begin local symbol block 807 my $intensity_to_mag_factor; # Calculate only if needed. 808 sub intensity_to_magnitude { 809 return - ($intensity_to_mag_factor ||= 2.5 / log (10)) * log ($_[0]); 810 } 811} 812 813=item ($day, $mon, $yr, $greg, $leap) = jd2date ($jd) 814 815This subroutine converts the given Julian day to the corresponding date. 816The returns are year - 1900, month (0 to 11), day (which may have a 817fractional part), a Gregorian calendar indicator which is true if the 818date is in the Gregorian calendar and false if it is in the Julian 819calendar, and a leap (or bissextile) year indicator which is true if the 820year is a leap year and false otherwise. The year 1 BC is returned as 821-1900 (i.e. as year 0), and so on. The date will probably have a 822fractional part (e.g. 2006 1 1.5 for noon January first 2006). 823 824If the $jd is before $JD_GREGORIAN, the date will be in the Julian 825calendar; otherwise it will be in the Gregorian calendar. 826 827The input may not be less than 0. 828 829The algorithm is from Jean Meeus' "Astronomical Algorithms", second 830edition, chapter 7 ("Julian Day"), pages 63ff, but the month is 831zero-based, not 1-based, and the year is 1900-based. 832 833=cut 834 835sub jd2date { 836 my ($time) = @_; 837 my $mod_jd = $time + 0.5; 838 my $Z = floor ($mod_jd); 839 my $F = $mod_jd - $Z; 840 my $A = (my $julian = $Z < $JD_GREGORIAN) ? $Z : do { 841 my $alpha = floor (($Z - 1867216.25)/36524.25); 842 $Z + 1 + $alpha - floor ($alpha / 4); 843 }; 844 my $B = $A + 1524; 845 my $C = floor (($B - 122.1) / 365.25); 846 my $D = floor (365.25 * $C); 847 my $E = floor (($B - $D) / 30.6001); 848 my $day = $B - $D - floor (30.6001 * $E) + $F; 849 my $mon = $E < 14 ? $E - 1 : $E - 13; 850 my $yr = $mon > 2 ? $C - 4716 : $C - 4715; 851 return ($day, $mon - 1, $yr - 1900, !$julian, ($julian ? !($yr % 4) : ( 852 $yr % 400 ? $yr % 100 ? !($yr % 4) : 0 : 1))); 853## % 400 ? 1 : $yr % 100 ? 0 : !($yr % 4)))); 854} 855 856=item ($sec, $min, $hr, $day, $mon, $yr, $wday, $yday, 0) = jd2datetime ($jd) 857 858This convenience subroutine converts the given Julian day to the 859corresponding date and time. All this really does is converts its 860argument to seconds since the system epoch, and pass off to 861epoch2datetime(). 862 863The input may not be less than 0. 864 865=cut 866 867sub jd2datetime { 868 my ($time) = @_; 869 return epoch2datetime(($time - JD_OF_EPOCH) * SECSPERDAY); 870} 871 872=item $century = jcent2000 ($time); 873 874Several of the algorithms in Jean Meeus' "Astronomical Algorithms" 875are expressed in terms of the number of Julian centuries from epoch 876J2000.0 (e.g equations 12.1, 22.1). This subroutine encapsulates 877that calculation. 878 879=cut 880 881sub jcent2000 {return jday2000 ($_[0]) / 36525} 882 883=item $jd = jday2000 ($time); 884 885This subroutine converts a Perl date to the number of Julian days 886(and fractions thereof) since Julian 2000.0. This quantity is used 887in a number of the algorithms in Jean Meeus' "Astronomical 888Algorithms". 889 890The computation makes use of information from Jean Meeus' "Astronomical 891Algorithms", 2nd Edition, Chapter 7, page 62. 892 893=cut 894 895sub jday2000 {return ($_[0] - PERL2000) / SECSPERDAY} 896 897=item $jd = julianday ($time); 898 899This subroutine converts a Perl date to a Julian day number. 900 901The computation makes use of information from Jean Meeus' "Astronomical 902Algorithms", 2nd Edition, Chapter 7, page 62. 903 904=cut 905 906sub julianday {return jday2000($_[0]) + 2_451_545.0} 907 908=item $ea = keplers_equation( $M, $e, $prec ); 909 910This subroutine solves Kepler's equation for the given mean anomaly 911C<$M> in radians, eccentricity C<$e> and precision C<$prec> in radians. 912It returns the eccentric anomaly in radians, to the given precision. 913 914The C<$prec> argument is optional, and defaults to the radian equivalent 915of 0.001 degrees. 916 917The algorithm is Meeus' equation 30.7, with John M. Steele's amendment 918for large values for the correction, given on page 205 of Meeus' book, 919 920This subroutine is B<not> used in the computation of satellite orbits, 921since the models have their own implementation. 922 923=cut 924 925sub keplers_equation { 926 my ( $mean_anomaly, $eccentricity, $precision ) = @_; 927 defined $precision 928 or $precision = 1.74533e-5; # Radians, = 0.001 degrees 929 my $curr = $mean_anomaly; 930 my $prev; 931 # Meeus' equation 30.7, page 199. 932 { 933 $prev = $curr; 934 my $delta = ( $mean_anomaly + $eccentricity * sin( $curr 935 ) - $curr ) / ( 1 - $eccentricity * cos $curr ); 936 # Steele's correction, page 205 937 $curr = $curr + max( -.5, min( .5, $delta ) ); 938 $precision < abs( $curr - $prev ) 939 and redo; 940 } 941 return $curr; 942} 943 944=item $rslt = load_module ($module_name) 945 946This convenience method loads the named module (using 'require'), 947throwing an exception if the load fails. If the load succeeds, it 948returns the result of the 'require' built-in, which is required to be 949true for a successful load. Results are cached, and subsequent attempts 950to load the same module simply give the cached results. 951 952=cut 953 954{ # Local symbol block. Oh, for 5.10 and state variables. 955 my %error; 956 my %rslt; 957 sub load_module { 958 my ($module) = @_; 959 exists $error{$module} and croak $error{$module}; 960 exists $rslt{$module} and return $rslt{$module}; 961 # I considered Module::Load here, but it appears not to support 962 # .pmc files. No, it's not an issue at the moment, but it may be 963 # if Perl 6 becomes a reality. 964 $rslt{$module} = eval "require $module" 965 or croak( $error{$module} = $@ ); 966 return $rslt{$module}; 967 } 968} # End local symbol block. 969 970=item $boolean = looks_like_number ($string); 971 972This subroutine returns true if the input looks like a number. It uses 973Scalar::Util::looks_like_number if that is available, otherwise it uses 974its own code, which is lifted verbatim from Scalar::Util 1.19, which in 975turn leans heavily on perlfaq4. 976 977=cut 978 979unless (eval {require Scalar::Util; Scalar::Util->import 980 ('looks_like_number'); 1}) { 981 no warnings qw{once}; 982 *looks_like_number = sub { 983 local $_ = shift; 984 985 # checks from perlfaq4 986 return 0 if !defined($_) || ref($_); 987 return 1 if (/^[+-]?[0-9]+$/); # is a +/- integer 988 return 1 if (/^([+-]?)(?=[0-9]|\.[0-9])[0-9]*(\.[0-9]*)?([Ee]([+-]?[0-9]+))?$/); # a C float 989 return 1 if ($] >= 5.008 and /^(Inf(inity)?|NaN)$/i) 990 or ($] >= 5.006001 and /^Inf$/i); 991 992 return 0; 993 }; 994} 995 996=item $maximum = max (...); 997 998This subroutine returns the maximum of its arguments. If List::Util can 999be loaded and 'max' imported, that's what you get. Otherwise you get a 1000pure Perl implementation. 1001 1002=cut 1003 1004unless (eval {require List::Util; List::Util->import ('max'); 1}) { 1005 no warnings qw{once}; 1006 *max = sub { 1007 my $rslt; 1008 foreach (@_) { 1009 defined $_ or next; 1010 if (!defined $rslt || $_ > $rslt) { 1011 $rslt = $_; 1012 } 1013 } 1014 $rslt; 1015 }; 1016} 1017 1018=item $minimum = min (...); 1019 1020This subroutine returns the minimum of its arguments. If List::Util can 1021be loaded and 'min' imported, that's what you get. Otherwise you get a 1022pure Perl implementation. 1023 1024=cut 1025 1026unless (eval {require List::Util; List::Util->import ('min'); 1}) { 1027 no warnings qw{once}; 1028 *min = sub { 1029 my $rslt; 1030 foreach (@_) { 1031 defined $_ or next; 1032 if (!defined $rslt || $_ < $rslt) { 1033 $rslt = $_; 1034 } 1035 } 1036 $rslt; 1037 }; 1038} 1039 1040=item $theta = mod2pi ($theta) 1041 1042This subroutine reduces the given angle in radians to the range 0 <= 1043$theta < TWOPI. 1044 1045=cut 1046 1047sub mod2pi {return $_[0] - floor ($_[0] / TWOPI) * TWOPI} 1048 1049=item $radians = omega ($time); 1050 1051This subroutine calculates the ecliptic longitude of the ascending node 1052of the Moon's mean orbit at the given B<dynamical> time. 1053 1054The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd 1055Edition, Chapter 22, pages 143ff. 1056 1057=cut 1058 1059sub omega { 1060 my $T = jcent2000 (shift); # Meeus (22.1) 1061 return mod2pi (deg2rad ((($T / 450000 + .0020708) * $T - 1062 1934.136261) * $T + 125.04452)); 1063} 1064 1065=item $pa = position_angle( $alpha1, $delta1, $alpha2, $delta2 ); 1066 1067This low-level subroutine calculates the position angle in right 1068ascension of the second body with respect to the first, given the first 1069body's right ascension and declination and the second body's right 1070ascension and declination in that order, B<in radians>. 1071 1072The return is the position angle B<in radians>, in the range 1073C<< -PI <= $pa < PI >>. 1074 1075The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd 1076Edition, page 116, but his algorithm is for the position angle of the 1077first body with respect to the second (i.e. the roles of the two bodies 1078are reversed). The order of arguments for this subroutine is consistent 1079with The IDL Astronomy User's Library at 1080L<https://idlastro.gsfc.nasa.gov/>, function C<posang()>. 1081 1082This is exposed because in principal you could calculate the position 1083angle in any spherical coordinate system, you would just need to get the 1084order of arguments right (e.g. azimuth, elevation or longitude, 1085latitude). 1086 1087=cut 1088 1089sub position_angle { 1090 my ( $alpha1, $delta1, $alpha2, $delta2 ) = @_; 1091 my $delta_alpha = $alpha2 - $alpha1; 1092 return atan2( sin( $delta_alpha ), 1093 cos( $delta1 ) * tan( $delta2 ) - 1094 sin( $delta1 ) * cos( $delta_alpha ) ); 1095} 1096 1097=item $degrees = rad2deg ($radians) 1098 1099This subroutine converts the given angle in radians to its equivalent 1100in degrees. If the argument is C<undef>, C<undef> will be returned. 1101 1102=cut 1103 1104sub rad2deg { return defined $_[0] ? $_[0] / PI * 180 : undef } 1105 1106=item $deg_min_sec = rad2dms( $radians, $decimals ) 1107 1108This subroutine converts the given angle in radians to its equivalent in 1109degrees, minutes and seconds, represented as a string. Degrees will be 1110suppressed if zero, and minutes will be suppressed if both degrees and 1111minutes are zero. If degrees are present the delimiter will be a degree 1112sign (C<"\N{DEGREE SIGN}>, a.k.a. C<"\N{U+00B0}">). The delimiters for 1113minutes and seconds are C<'> and C<"> respectively, with the C<"> 1114appearing before the decimal point. 1115 1116The optional C<$decimals> argument specifies the number of decimal 1117places in the seconds value, and defaults to C<3>. 1118 1119=cut 1120 1121sub rad2dms { 1122 my ( $rad, $dp ) = @_; 1123 defined $rad 1124 or return $rad; 1125 defined $dp 1126 or $dp = 3; 1127 my $sec = rad2deg( $rad ) * 3600; 1128 ( $sec, my $sgn ) = $sec < 0 ? ( - $sec, '-' ) : ( $sec, '' ); 1129 my $frc = sprintf "%.${dp}f", $sec; 1130 $frc =~ s/ [^.]* //smx; 1131 $sec = floor( $sec ); 1132 my $min = floor( $sec / 60 ); 1133 $sec %= 60; 1134 my $deg = floor( $min / 60 ); 1135 $min %= 60; 1136 $deg or $min 1137 or return sprintf q<%s%d"%s>, $sgn, $sec, $frc; 1138 $deg 1139 or return sprintf q<%s%d'%02d"%s>, $sgn, $sec, $frc; 1140 return sprintf qq<%s%d°%02d'%02d"%s>, 1141 $sgn, $deg, $min, $sec, $frc; 1142} 1143 1144=item $hr_min_sec = rad2hms( $radians, $decimals ) 1145 1146This subroutine converts the given angle in radians to its equivalent in 1147hours, minutes and seconds (presumably of right ascension), represented 1148as a string. Hours will be suppressed if zero, and minutes will be 1149suppressed if both hours and minutes are zero. The delimiters for hours, 1150minutes, and seconds are C<'h'>, C<'m'>, and C<'s'> respectively, with 1151the C<'s'> appearing before the decimal point. 1152 1153The optional C<$decimals> argument specifies the number of decimal 1154places in the seconds value, and defaults to C<3>. 1155 1156=cut 1157 1158sub rad2hms { 1159 my ( $rad, $dp ) = @_; 1160 defined $rad 1161 or return $rad; 1162 defined $dp 1163 or $dp = 3; 1164 my $sec = $rad * 12 / PI * 3600; 1165 ( $sec, my $sgn ) = $sec < 0 ? ( - $sec, '-' ) : ( $sec, '' ); 1166 my $frc = sprintf "%.${dp}f", $sec; 1167 $frc =~ s/ [^.]* //smx; 1168 $sec = floor( $sec ); 1169 my $min = floor( $sec / 60 ); 1170 $sec %= 60; 1171 my $hr = floor( $min / 60 ); 1172 $min %= 60; 1173 $hr or $min 1174 or return sprintf q<%s%ds%s>, $sgn, $sec, $frc; 1175 $hr 1176 or return sprintf q<%s%dm%02ds%s>, $sgn, $sec, $frc; 1177 return sprintf qq<%s%dh%02dm%02ds%s>, 1178 $sgn, $hr, $min, $sec, $frc; 1179} 1180 1181=item $value = tan ($angle) 1182 1183This subroutine computes the tangent of the given angle in radians. 1184 1185=cut 1186 1187sub tan {return sin ($_[0]) / cos ($_[0])} 1188 1189=item $value = theta0 ($time); 1190 1191This subroutine returns the Greenwich hour angle of the mean equinox at 11920 hours universal on the day whose time is given (i.e. the argument is 1193a standard Perl time). 1194 1195=cut 1196 1197sub theta0 { 1198 my ($time) = @_; 1199 my @t = gmtime $time; 1200 $t[5] += 1900; 1201 return thetag( greg_time_gm( 0, 0, 0, @t[3 .. 5] ) ); 1202} 1203 1204=item $value = thetag ($time); 1205 1206This subroutine returns the Greenwich hour angle of the mean equinox at 1207the given time. 1208 1209The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd 1210Edition, equation 12.4, page 88. 1211 1212=cut 1213 1214# Meeus, pg 88, equation 12.4, converted to radians and Perl dates. 1215 1216sub thetag { 1217 my ($time) = @_; 1218 my $T = jcent2000 ($time); 1219 return mod2pi (4.89496121273579 + 6.30038809898496 * 1220 jday2000 ($time)) 1221 + (6.77070812713916e-06 - 4.5087296615715e-10 * $T) * $T * $T; 1222} 1223 1224# time_gm and time_local are actually created at the top of the module. 1225 1226=item $epoch = time_gm( $sec, $min, $hr, $day, $mon, $yr ); 1227 1228This exportable subroutine is a wrapper for either 1229C<Time::y2038::timegm()> (if that module is installed), or 1230C<Time::Local::timegm()> (if not.) 1231 1232This subroutine interprets years the same way C<Time::Local::timegm()> 1233does. 1234 1235This wrapper is needed because the routines have subtly different 1236signatures. L<Time::y2038|Time::y2038> C<timegm()> interprets years 1237strictly as Perl years. L<Time::Local|Time::Local> C<timegm()> 1238interprets years differently depending on the value of the year; greater 1239than 999 as Gregorian years, but other years are Perl years, except for 1240the years 0 to 99 inclusive, which are interpreted as within 50 years of 1241the current year. 1242 1243This subroutine is discouraged in favor of C<greg_time_gm()>, which 1244always interprets years as Gregorian years. 1245 1246=item $epoch = time_local( $sec, $min, $hr, $day, $mon, $yr ); 1247 1248This exportable subroutine is a wrapper for either 1249C<Time::y2038::timelocal()> (if that module is installed), or 1250C<Time::Local::timelocal()> (if not.) 1251 1252This subroutine interprets years the same way 1253C<Time::Local::timelocal()> does. 1254 1255This wrapper is needed for the same reason C<time_gm()> is 1256needed. 1257 1258This subroutine is discouraged in favor of C<greg_time_local()>, which 1259always interprets years as Gregorian years. 1260 1261=item $a = vector_cross_product( $b, $c ); 1262 1263This subroutine computes and returns the vector cross product of $b and 1264$c. Vectors are represented by array references. The cross product is 1265only defined if both arrays have 3 elements. 1266 1267=cut 1268 1269sub vector_cross_product { 1270 my ( $b, $c ) = @_; 1271 @{ $b } == 3 and @{ $c } == 3 1272 or confess 'Programming error - vector_cross_product arguments', 1273 ' must be references to arrays of length 3'; 1274 return [ 1275 $b->[1] * $c->[2] - $b->[2] * $c->[1], 1276 $b->[2] * $c->[0] - $b->[0] * $c->[2], 1277 $b->[0] * $c->[1] - $b->[1] * $c->[0], 1278 ]; 1279} 1280 1281=item $a = vector_dot_product( $b, $c ); 1282 1283This subroutine computes and returns the vector dot product of $b and 1284$c. Vectors are represented by array references. The dot product is only 1285defined if both arrays have the same number of elements. 1286 1287=cut 1288 1289sub vector_dot_product { 1290 my ( $b, $c ) = @_; 1291 @{ $b } == @{ $c } 1292 or confess 'Programming error - vector_dot_product arguments ', 1293 'must be references to arrays of the same length'; 1294 my $prod = 0; 1295 my $size = @{ $b } - 1; 1296 foreach my $inx ( 0 .. $size ) { 1297 $prod += $b->[$inx] * $c->[$inx]; 1298 } 1299 return $prod; 1300} 1301 1302=item $m = vector_magnitude( $x ); 1303 1304This subroutine computes and returns the magnitude of vector $x. The 1305vector is represented by an array reference. 1306 1307=cut 1308 1309sub vector_magnitude { 1310 my ( $x ) = @_; 1311 ARRAY_REF eq ref $x 1312 or confess 'Programming error - vector_magnitude argument ', 1313 'must be a reference to an array'; 1314 my $mag = 0; 1315 my $size = @{ $x } - 1; 1316 foreach my $inx ( 0 .. $size ) { 1317 $mag += $x->[$inx] * $x->[$inx]; 1318 } 1319 return sqrt $mag; 1320} 1321 1322=item $u = vector_unitize( $x ); 1323 1324This subroutine computes and returns a unit vector pointing in the same 1325direction as $x. The vectors are represented by array references. 1326 1327=cut 1328 1329sub vector_unitize { 1330 my ( $x ) = @_; 1331 ARRAY_REF eq ref $x 1332 or confess 'Programming error - vector_unitize argument ', 1333 'must be a reference to an array'; 1334 my $mag = vector_magnitude( $x ); 1335 return [ map { $_ / $mag } @{ $x } ]; 1336} 1337 1338# __classisa( 'Foo', 'Bar' ); 1339# 1340# Returns true if Foo->isa( 'Bar' ) is true, and false otherwise. 1341# In particular, returns false if the first argument is a 1342# reference. 1343 1344sub __classisa { 1345 my ( $invocant, $class ) = @_; 1346 ref $invocant and return; 1347 defined $invocant or return; 1348 return $invocant->isa( $class ); 1349} 1350 1351# __instance( $foo, 'Bar' ); 1352# 1353# Returns true if $foo is an instance of 'Bar', and false 1354# otherwise. In particular, returns false if $foo is not a 1355# reference, or if it is not blessed. 1356 1357sub __instance { 1358 my ( $object, $class ) = @_; 1359 ref $object or return; 1360 blessed( $object ) or return; 1361 return $object->isa( $class ); 1362} 1363 1364# $epoch_time = __parse_time_iso_8601 1365# 1366# Parse ISO 8601 date/time. I do not intend to expose this, since 1367# it will probably go away when the satpass script is dropped. It 1368# would actually be in that script except for the fact that it can 1369# be more easily tested here, and because of the possibility that 1370# it will be used in App::Satpass2. 1371{ 1372 1373 my %special_day_offset = ( 1374 yesterday => -SECSPERDAY(), 1375 today => 0, 1376 tomorrow => SECSPERDAY(), 1377 ); 1378 1379 sub __parse_time_iso_8601 { 1380 my ( $string ) = @_; 1381 1382 my @zone; 1383 $string =~ s/ \s* (?: ( Z ) | 1384 ( [+-] ) ( [0-9]{2} ) :? ( [0-9]{2} )? ) \z //smx 1385 and @zone = ( $1, $2, $3, $4 ); 1386 my @date; 1387 1388 # ISO 8601 date 1389 if ( $string =~ m{ \A 1390 ( [0-9]{4} [^0-9]? | [0-9]{2} [^0-9] ) # year: $1 1391 (?: ( [0-9]{1,2} ) [^0-9]? # month: $2 1392 (?: ( [0-9]{1,2} ) (?: \s* | [^0-9]? ) # day: $3 1393 (?: ( [0-9]{1,2} ) [^0-9]? # hour: $4 1394 (?: ( [0-9]{1,2} ) [^0-9]? # minute: $5 1395 (?: ( [0-9]{1,2} ) [^0-9]? # second: $6 1396 ( [0-9]* ) # fract: $7 1397 )? 1398 )? 1399 )? 1400 )? 1401 )? 1402 \z 1403 }smx ) { 1404 @date = ( $1, $2, $3, $4, $5, $6, $7, undef ); 1405 1406 # special-case 'yesterday', 'today', and 'tomorrow'. 1407 } elsif ( $string =~ m< \A 1408 ( yesterday | today | tomorrow ) # day: $1 1409 (?: [^0-9]* ( [0-9]{1,2} ) [^0-9]? # hour: $2 1410 (?: ( [0-9]{1,2} ) [^0-9]? # minute: $3 1411 (?: ( [0-9]{1,2} ) [^0-9]? # second: $4 1412 ( [0-9]* ) # fract: $5 1413 )? 1414 )? 1415 )? 1416 \z >smx ) { 1417 my @today = @zone ? gmtime : localtime; 1418 @date = ( $today[5] + 1900, $today[4] + 1, $today[3], $2, $3, 1419 $4, $5, $special_day_offset{$1} ); 1420 1421 } else { 1422 1423 return; 1424 1425 } 1426 1427 my $offset = pop @date || 0; 1428 if ( @zone && !$zone[0] ) { 1429 my ( undef, $sign, $hr, $min ) = @zone; 1430 $offset -= $sign . ( ( $hr * 60 + ( $min || 0 ) ) * 60 ) 1431 } 1432 1433 foreach ( @date ) { 1434 defined $_ and s/ [^0-9]+ //smxg; 1435 } 1436 1437 $date[0] = __tle_year_to_Gregorian_year( $date[0] ); 1438 1439 defined $date[1] and --$date[1]; 1440 defined $date[2] or $date[2] = 1; 1441 my $frc = pop @date; 1442 1443 foreach ( @date ) { 1444 defined $_ or $_ = 0; 1445 } 1446 1447 my $time; 1448 if ( @zone ) { 1449 $time = greg_time_gm( reverse @date ); 1450 } else { 1451 $time = greg_time_local( reverse @date ); 1452 } 1453 1454 if ( defined $frc && $frc ne '') { 1455 my $denom = 1 . ( 0 x length $frc ); 1456 $time += $frc / $denom; 1457 } 1458 1459 return $time + $offset; 1460 } 1461 1462} 1463 1464sub __sprintf { 1465 my ( $tplt, @args ) = @_; 1466 defined $tplt 1467 or return undef; ## no critic (ProhibitExplicitReturnUndef) 1468 no if $] gt '5.021002', qw{ warnings redundant }; 1469 return sprintf $tplt, @args; 1470} 1471 1472{ 1473 my %deprecate = ( 1474 ); 1475 1476 sub __subroutine_deprecation { 1477 ( my $sub = ( caller 1 )[3] ) =~ s/ .* :: //smx; 1478 my $info = $deprecate{$sub} or return; 1479 $info->{level} 1480 or return; 1481 my $msg = "Subroutine $sub() deprecated in favor of @{[ 1482 $info->{method} || $sub ]}() method"; 1483 $info->{level} >= 3 1484 and croak $msg; 1485 carp $msg; 1486 $info->{level} == 1 1487 and $info->{level} = 0; 1488 return; 1489 } 1490} 1491 1492=item $year = __tle_year_to_Gregorian_year( $year ) 1493 1494The TLE data contain the year in two-digit form. NORAD decided to deal 1495with Y2K by decreeing that year numbers lower than 57 (the launch of 1496Sputnik 1) are converted to Gregorian by adding 2000. Years numbers of 149757 or greater are still converted to Gregorian by adding 1900. This 1498subroutine encapsulates this logic. Years of 100 or greater are returned 1499unmodified. 1500 1501This subroutine is B<private> to this package, and can be changed or 1502revoked without notice. 1503 1504=cut 1505 1506sub __tle_year_to_Gregorian_year { 1507 my ( $year ) = @_; 1508 return $year < 57 ? $year + 2000 : 1509 $year < 100 ? $year + 1900 : $year; 1510} 1511 15121; 1513 1514__END__ 1515 1516=back 1517 1518=head1 ACKNOWLEDGMENTS 1519 1520The author wishes to acknowledge Jean Meeus, whose book "Astronomical 1521Algorithms" (second edition) published by Willmann-Bell Inc provided 1522several of the algorithms implemented herein. Willmann-Bell ceased to be 1523a separate entity in 2021, but their publications, including Dr. Meeus' 1524book, are still available through Sky and Telescope's Willmann-Bell 1525imprint at L<https://shopatsky.com/collections/willmann-bell>. 1526 1527=head1 BUGS 1528 1529Support is by the author. Please file bug reports at 1530L<https://rt.cpan.org/Public/Dist/Display.html?Name=Astro-satpass>, 1531L<https://github.com/trwyant/perl-Astro-Coord-ECI/issues>, or in 1532electronic mail to the author. 1533 1534=head1 AUTHOR 1535 1536Thomas R. Wyant, III (F<wyant at cpan dot org>) 1537 1538=head1 COPYRIGHT AND LICENSE 1539 1540Copyright (C) 2005-2021 by Thomas R. Wyant, III 1541 1542This program is free software; you can redistribute it and/or modify it 1543under the same terms as Perl 5.10.0. For more details, see the full text 1544of the licenses in the directory LICENSES. 1545 1546This program is distributed in the hope that it will be useful, but 1547without any warranty; without even the implied warranty of 1548merchantability or fitness for a particular purpose. 1549 1550=cut 1551 1552# ex: set textwidth=72 : 1553