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