File Coverage

blib/lib/Astro/Coord/ECI/Sun.pm
Criterion Covered Total %
statement 158 202 78.2
branch 41 78 52.5
condition 10 18 55.5
subroutine 34 36 94.4
pod 11 11 100.0
total 254 345 73.6


line stmt bran cond sub pod time code
1             =head1 NAME
2              
3             Astro::Coord::ECI::Sun - Compute the position of the Sun.
4              
5             =head1 SYNOPSIS
6              
7             use Astro::Coord::ECI;
8             use Astro::Coord::ECI::Sun;
9             use Astro::Coord::ECI::Utils qw{deg2rad};
10            
11             # 1600 Pennsylvania Ave, Washington DC USA
12             # latitude 38.899 N, longitude 77.038 W,
13             # altitude 16.68 meters above sea level
14             my $lat = deg2rad (38.899); # Radians
15             my $long = deg2rad (-77.038); # Radians
16             my $alt = 16.68 / 1000; # Kilometers
17             my $sun = Astro::Coord::ECI::Sun->new ();
18             my $sta = Astro::Coord::ECI->
19             universal (time ())->
20             geodetic ($lat, $long, $alt);
21             my ($time, $rise) = $sta->next_elevation ($sun);
22             print "Sun @{[$rise ? 'rise' : 'set']} is ",
23             scalar gmtime $time, " UT\n";
24              
25             Although this example computes the Sun rise or set in Washington D.C.
26             USA, the time is displayed in Universal Time. This is because I did not
27             want to complicate the example by adding machinery to convert the time
28             to the correct zone for Washington D.C. (which is UT - 5 except when
29             Summer Time is in effect, when it is UT - 4).
30              
31             =head1 DESCRIPTION
32              
33             This module implements the position of the Sun as a function of time, as
34             described in Jean Meeus' "Astronomical Algorithms," second edition. It
35             is a subclass of B, with the id, name, and diameter
36             attributes initialized appropriately, and the time_set() method
37             overridden to compute the position of the Sun at the given time.
38              
39             =head2 Methods
40              
41             The following methods should be considered public:
42              
43             =over
44              
45             =cut
46              
47             package Astro::Coord::ECI::Sun;
48              
49 16     16   1925 use strict;
  16         35  
  16         527  
50 16     16   90 use warnings;
  16         34  
  16         825  
51              
52             our $VERSION = '0.129';
53              
54 16     16   105 use base qw{Astro::Coord::ECI};
  16         35  
  16         1964  
55              
56 16     16   108 use Astro::Coord::ECI::Utils qw{ @CARP_NOT :mainstream };
  16         45  
  16         7542  
57 16     16   127 use Carp;
  16         34  
  16         1102  
58 16     16   103 use POSIX qw{ ceil floor strftime };
  16         33  
  16         136  
59              
60 16     16   1489 use constant MEAN_MAGNITUDE => -26.8;
  16         35  
  16         6640  
61              
62             my %attrib = (
63             iterate_for_quarters => 1,
64             );
65              
66             my %static = (
67             id => 'Sun',
68             name => 'Sun',
69             diameter => 1392000,
70             iterate_for_quarters => undef,
71             );
72              
73             my $weaken = eval {
74             require Scalar::Util;
75             Scalar::Util->can('weaken');
76             };
77              
78             our $Singleton = $weaken;
79              
80             my %object; # By class
81              
82             =item $sun = Astro::Coord::ECI::Sun->new();
83              
84             This method instantiates an object to represent the coordinates of the
85             Sun. This is a subclass of L, with
86             the id and name attributes set to 'Sun', and the diameter attribute set
87             to 1392000 km per Jean Meeus' "Astronomical Algorithms", 2nd Edition,
88             Appendix I, page 407.
89              
90             Any arguments are passed to the set() method once the object has been
91             instantiated. Yes, you can override the "hard-wired" id, name, and so
92             forth in this way.
93              
94             If C<$Astro::Coord::ECI::Sun::Singleton> is true, you get a singleton
95             object; that is, only one object is instantiated and subsequent calls to
96             C just return that object. If higher-accuracy subclasses are ever
97             implemented, there will be one singleton for each class.
98              
99             The singleton logic only works if L exports
100             C. If it does not, the setting of
101             C<$Astro::Coord::ECI::Sun::Singleton> is silently ignored. The default
102             is true if L can be loaded and exports
103             C, and false otherwise.
104              
105             =cut
106              
107             sub new {
108 134     134 1 1949 my ($class, @args) = @_;
109 134 50       311 ref $class and $class = ref $class;
110 134 100 66     688 if ( $Singleton && $weaken && __classisa( $class, __PACKAGE__ ) ) {
      66        
111 132         214 my $self;
112 132 100       362 if ( $self = $object{$class} ) {
113 97 100       225 $self->set( @args ) if @args;
114 97         293 return $self;
115             } else {
116 35         380 $self = $object{$class} = $class->SUPER::new (%static, @args);
117 35         240 $weaken->( $object{$class} );
118 35         156 return $self;
119             }
120             } else {
121 2         9 return $class->SUPER::new (%static, @args);
122             }
123             }
124              
125             =item @almanac = $sun->almanac( $station, $start, $end );
126              
127             This method produces almanac data for the Sun for the given observing
128             station, between the given start and end times. The station is assumed
129             to be Earth-Fixed - that is, you can't do this for something in orbit.
130              
131             The C<$station> argument may be omitted if the C attribute has
132             been set. That is, this method can also be called as
133              
134             @almanac = $sun->almanac( $start, $end )
135              
136             The start time defaults to the current time setting of the $sun
137             object, and the end time defaults to a day after the start time.
138              
139             The almanac data consists of a list of list references. Each list
140             reference points to a list containing the following elements:
141              
142             [0] => time
143             [1] => event (string)
144             [2] => detail (integer)
145             [3] => description (string)
146              
147             The @almanac list is returned sorted by time.
148              
149             The following events, details, and descriptions are at least
150             potentially returned:
151              
152             horizon: 0 = Sunset, 1 = Sunrise;
153             transit: 0 = local midnight, 1 = local noon;
154             twilight: 0 = end twilight, 1 = begin twilight;
155             quarter: 0 = spring equinox, 1 = summer solstice,
156             2 = fall equinox, 3 = winter solstice.
157              
158             Twilight is calculated based on the current value of the twilight
159             attribute of the $sun object. This attribute is inherited from
160             L, and documented there.
161              
162             =cut
163              
164             sub __almanac_event_type_iterator {
165 2     2   6 my ( $self, $station ) = @_;
166              
167 2         4 my $inx = 0;
168              
169 2         8 my $horizon = $station->__get_almanac_horizon();
170              
171 2         13 my @events = (
172             [ $station, next_elevation => [ $self, $horizon, 1 ],
173             horizon => '__horizon_name' ],
174             [ $station, next_meridian => [ $self ],
175             transit => '__transit_name' ],
176             [ $station, next_elevation =>
177             [ $self, $self->get( 'twilight' ) + $horizon, 0 ],
178             twilight => '__twilight_name' ],
179             [ $self, next_quarter => [], quarter => '__quarter_name', ],
180             );
181              
182             return sub {
183             $inx < @events
184 10 100   10   38 and return @{ $events[$inx++] };
  8         53  
185 2         9 return;
186 2         13 };
187             }
188              
189 16     16   5961 use Astro::Coord::ECI::Mixin qw{ almanac };
  16         47  
  16         973  
190              
191             =item @almanac = $sun->almanac_hash( $station, $start, $end );
192              
193             This convenience method wraps $sun->almanac(), but returns a list of
194             hash references, sort of like Astro::Coord::ECI::TLE->pass()
195             does. The hashes contain the following keys:
196              
197             {almanac} => {
198             {event} => the event type;
199             {detail} => the event detail (typically 0 or 1);
200             {description} => the event description;
201             }
202             {body} => the original object ($sun);
203             {station} => the observing station;
204             {time} => the time the quarter occurred.
205              
206             The {time}, {event}, {detail}, and {description} keys correspond to
207             elements 0 through 3 of the list returned by almanac().
208              
209             =cut
210              
211 16     16   119 use Astro::Coord::ECI::Mixin qw{ almanac_hash };
  16         34  
  16         9255  
212              
213             =item $coord2 = $coord->clone ();
214              
215             If singleton objects are enabled, this override of the superclass'
216             method simply returns the invocant. Otherwise it does a deep clone of an
217             object, producing a different but identical object.
218              
219             Prior to version 0.099_01 it always returned a clone. Yes,
220             this is a change in long-standing functionality, but a long-standing bug
221             is still a bug.
222              
223             =cut
224              
225             sub clone {
226 2     2 1 612 my ( $self ) = @_;
227 2 100 66     11 $Singleton
228             and $weaken
229             and return $self;
230 1         10 return $self->SUPER::clone();
231             }
232              
233             =item $elevation = $tle->correct_for_refraction( $elevation )
234              
235             This override of the superclass' method simply returns the elevation
236             passed to it. I have no algorithm for refraction at the surface of the
237             photosphere or anywhere else in the environs of the Sun, and explaining
238             why I make no correction at all seemed easier than explaining why I make
239             an incorrect correction.
240              
241             See the L C and
242             C documentation for whether this class'
243             C method is actually called by those methods.
244              
245             =cut
246              
247             sub correct_for_refraction {
248 0     0 1 0 my ( undef, $elevation ) = @_; # Invocant unused
249 0         0 return $elevation;
250             }
251              
252             =item $long = $sun->geometric_longitude ()
253              
254             This method returns the geometric longitude of the Sun in radians at
255             the last time set.
256              
257             =cut
258              
259             sub geometric_longitude {
260 1125     1125 1 1748 my $self = shift;
261 1125 50       2387 croak <{_sun_geometric_longitude};
262             Error - You must set the time of the Sun object before the geometric
263             longitude can be returned.
264             eod
265              
266 1125         2416 return $self->{_sun_geometric_longitude};
267             }
268              
269             =item $sun->get( ... )
270              
271             This method has been overridden to return the invocant as the C<'sun'>
272             attribute.
273              
274             =cut
275              
276             sub get {
277 6742     6742 1 13362 my ( $self, @args ) = @_;
278 6742         9649 my @rslt;
279 6742         11651 foreach my $name ( @args ) {
280             push @rslt, 'sun' eq $name ? $self :
281             $attrib{$name} ?
282 6742 0       25575 ref $self ? $self->{$name} : $static{$name} :
    50          
    50          
283             $self->SUPER::get( $name );
284             }
285 6742 100       21342 return wantarray ? @rslt : $rslt[0];
286             }
287              
288             =item ($point, $intens, $central) = $sun->magnitude ($theta, $omega);
289              
290             This method returns the magnitude of the Sun at a point $theta radians
291             from the center of its disk, given that the disk's angular radius
292             (B diameter) is $omega radians. The returned $point is the
293             magnitude at the given point (undef if $theta > $omega), $intens is the
294             ratio of the intensity at the given point to the central intensity (0
295             if $theta > $omega), and $central is the central magnitude.
296              
297             If this method is called in scalar context, it returns $point, the point
298             magnitude.
299              
300             If the $omega argument is omitted or undefined, it is calculated based
301             on the geocentric range to the Sun at the current time setting of the
302             object.
303              
304             If the $theta argument is omitted or undefined, the method returns
305             the average magnitude of the Sun, which is taken to be -26.8.
306              
307             The limb-darkening algorithm and the associated constants come from
308             L.
309              
310             For consistency's sake, an observing station can optionally be passed as
311             the first argument (i.e. before C<$theta>). This is currently ignored.
312              
313             =cut
314              
315             { # Begin local symbol block
316              
317             my $central_mag;
318             my @limb_darkening = (.3, .93, -.23);
319              
320             sub magnitude {
321 0     0 1 0 my ( $self, @args ) = @_;
322             # We want to accept a station as the second argument for
323             # consistency's sake, though we do not (at this point) use it.
324 0 0       0 embodies( $args[0], 'Astro::Coord::ECI' )
325             and shift @args;
326 0         0 my ( $theta, $omega ) = @args;
327 0 0       0 return MEAN_MAGNITUDE unless defined $theta;
328 0 0       0 unless (defined $omega) {
329 0         0 my @eci = $self->eci ();
330 0         0 $omega = $self->get ('diameter') / 2 /
331             sqrt (distsq (\@eci[0 .. 2], [0, 0, 0]));
332             }
333 0 0       0 unless (defined $central_mag) {
334 0         0 my $sum = 0;
335 0         0 my $quotient = 2;
336 0         0 foreach my $a (@limb_darkening) {
337 0         0 $sum += $a / $quotient++;
338             }
339 0         0 $central_mag = MEAN_MAGNITUDE - intensity_to_magnitude (2 * $sum);
340             }
341 0         0 my $intens = 0;
342 0         0 my $point;
343 0 0       0 if ($theta < $omega) {
344 0         0 my $costheta = cos ($theta);
345 0         0 my $cosomega = cos ($omega);
346 0         0 my $sinomega = sin ($omega);
347 0         0 my $cospsi = sqrt ($costheta * $costheta -
348             $cosomega * $cosomega) / $sinomega;
349 0         0 my $psiterm = 1;
350 0         0 foreach my $a (@limb_darkening) {
351 0         0 $intens += $a * $psiterm;
352 0         0 $psiterm *= $cospsi;
353             }
354 0         0 $point = $central_mag + intensity_to_magnitude ($intens);
355             }
356 0 0       0 return wantarray ? ($point, $intens, $central_mag) : $point;
357             }
358             } # End local symbol block.
359              
360             =item ($time, $quarter, $desc) = $sun->next_quarter($want);
361              
362             This method calculates the time of the next equinox or solstice after
363             the current time setting of the $sun object. The return is the time,
364             which equinox or solstice it is as a number from 0 (March equinox) to 3
365             (December solstice), and a string describing the equinox or solstice. If
366             called in scalar context, you just get the time.
367              
368             If the C attribute is not set or set to a location on or north
369             of the Equator, the descriptor strings are
370              
371             0 - Spring equinox
372             1 - Summer solstice
373             2 - Fall equinox
374             3 - Winter solstice
375              
376             If the C attribute is set to a location south of the Equator,
377             the descriptor strings are
378              
379             0 - Fall equinox
380             1 - Winter solstice
381             2 - Spring equinox
382             3 - Summer solstice
383              
384             The optional $want argument says which equinox or solstice you want, as
385             a number from 0 through 3.
386              
387             As a side effect, the time of the $sun object ends up set to the
388             returned time.
389              
390             As of version 0.088_01, the algorithm given in Jean Meeus' "Astronomical
391             Algorithms", 2nd Edition, Chapter 27 ("Equinoxes and Solstices"), pages
392             278ff is used. This should be good for the range -1000 to 3000
393             Gregorian, and good to within a minute or so within the range 1951 to
394             2050 Gregorian, but the longitude of the Sun at the calculated time may
395             be as much as 0.01 degree off the exact time for the event.
396              
397             If you take the United States Naval Observatory's times (given to the
398             nearest minute) as the standard, the maximum deviation from that
399             standard in the range 1700 to 2100 is 226 seconds. I have no information
400             on this algorithm's accuracy outside that range. I.
401              
402             In version 0.088 and before, this calculation was done by successive
403             approximation based on the position of the Sun, and was good to about 15
404             minutes.
405              
406             If you want the old iterative version back, set attribute
407             C to a true value.
408              
409             =cut
410              
411 16     16   127 use constant NEXT_QUARTER_INCREMENT => 86400 * 85; # 85 days.
  16         34  
  16         4959  
412              
413             *__next_quarter_coordinate = __PACKAGE__->can(
414             'ecliptic_longitude' );
415              
416             # use Astro::Coord::ECI::Mixin qw{ next_quarter };
417              
418             sub next_quarter {
419 7     7 1 1308 my ( $self, $quarter ) = @_;
420             $self->{iterate_for_quarters}
421 7 50       21 and goto &Astro::Coord::ECI::Mixin::next_quarter;
422 7         21 my $time = $self->universal();
423 7         31 my $year = ( gmtime( $time ) )[5] + 1900;
424 7         15 my $season;
425 7 50       18 if ( defined $quarter ) {
426             # I can't think of an edge case that makes the first calculation
427             # give a quarter that is too late. Wish that made me feel better
428             # than it in fact does.
429 0         0 $season = $self->season( $year, $quarter );
430 0 0       0 $season < $time
431             and $season = $self->season( $year + 1, $quarter );
432 0         0 $self->universal( $season );
433             } else {
434 7         19 my ( undef, $lon ) = $self->ecliptic();
435             # The fudged-in 359 is because I am worried about the limited
436             # accuracy of the longitude calculation causing me to pick the
437             # wrong quarter. So I essentially back the Sun up by about a
438             # day. That may indeed put me a quarter too early, but I have to
439             # check for that anyway because of the accuracy (or lack
440             # thereof) of the calculated longitude.
441 7         20 $quarter = ceil( ( rad2deg( $lon ) + 359 ) / 90 ) % 4;
442 7         26 $season = $self->season( $year, $quarter );
443             # The above calculation gives the wrong $season between the
444             # December equinox and the end of the year. We fix it up here:
445 7 100       20 if ( $time - $season > SECSPERDAY * 180 ) {
446 1         4 $year++;
447 1         5 $season = $self->season( $year, $quarter );
448             }
449             # If we're a quarter too early, add one and repeat the
450             # calculation. We shouldn't have to do this more than once,
451             # since our maximum error even with the fudge factor is a day.
452 7 100       22 if ( $season < $time ) {
453 3         5 $quarter++;
454 3 50       8 if ( $quarter > 3 ) {
455 0         0 $quarter -= 4;
456 0         0 $year++;
457             }
458 3         8 $season = $self->season( $year, $quarter );
459             }
460             }
461 7         19 $season = ceil( $season ); # Make sure we're AFTER.
462 7         24 $self->universal( $season );
463 7 100       35 return wantarray ? ( $season, $quarter, $self->__quarter_name(
464             $quarter ) ) : $season;
465             }
466              
467             =item $hash_reference = $sun->next_quarter_hash($want);
468              
469             This convenience method wraps $sun->next_quarter(), but returns the
470             data in a hash reference, sort of like Astro::Coord::ECI::TLE->pass()
471             does. The hash contains the following keys:
472              
473             {body} => the original object ($sun);
474             {almanac} => {
475             {event} => 'quarter',
476             {detail} => the quarter number (0 through 3);
477             {description} => the quarter description;
478             }
479             {time} => the time the quarter occurred.
480              
481             The {time}, {detail}, and {description} keys correspond to elements 0
482             through 2 of the list returned by next_quarter().
483              
484             =cut
485              
486 16     16   133 use Astro::Coord::ECI::Mixin qw{ next_quarter_hash };
  16         35  
  16         13487  
487              
488             =item $period = $sun->period ()
489              
490             Although this method is attached to an object that represents the
491             Sun, what it actually returns is the sidereal period of the Earth,
492             per Appendix I (pg 408) of Jean Meeus' "Astronomical Algorithms,"
493             2nd edition.
494              
495             =cut
496              
497 246     246 1 963 sub period {return 31558149.7632} # 365.256363 * 86400
498              
499             sub __horizon_name_tplt {
500 4     4   11 my ( $self ) = @_;
501 4 50       25 return $self->__object_is_self_named() ?
502             [ '%sset', '%srise' ] :
503             $self->SUPER::__horizon_name_tplt();
504             }
505              
506             sub __quarter_name {
507 2     2   7 my ( $self, $event, $tplt ) = @_;
508 2 50 33     15 $tplt ||= $self->__object_is_self_named() ?
509             [
510             'Spring equinox', 'Summer solstice',
511             'Fall equinox', 'Winter solstice',
512             ] : [
513             '%s Spring equinox', '%s Summer solstice',
514             '%s Fall equinox', '%s Winter solstice',
515             ];
516 2         5 my $station;
517             $station = $self->get( 'station' )
518             and ( $station->geodetic() )[0] < 0
519 2 50 66     6 and $event = ( $event + @{ $tplt } / 2 ) % @{ $tplt };
  0         0  
  0         0  
520 2         10 return $self->__event_name( $event, $tplt );
521             }
522              
523             sub __transit_name_tplt {
524 4     4   10 my ( $self ) = @_;
525 4 50       13 return $self->__object_is_self_named() ?
526             [ 'local midnight', 'local noon' ] :
527             [ '%s transits nadir', '%s transits meridian' ];
528             }
529              
530             sub __twilight_name {
531 4     4   12 my ( $self, $event, $tplt ) = @_;
532 4 50 33     27 $tplt ||= $self->__object_is_self_named() ?
533             [ 'end twilight', 'begin twilight' ] :
534             [ 'end %s twilight', 'begin %s twilight' ];
535 4         16 return $self->__event_name( $event, $tplt );
536             }
537              
538             =begin comment
539              
540             =item $time = $self->season( $year, $season );
541              
542             This method calculates the time of the given season of the given
543             Gregorian year. The $season is an integer from 0 to 3, with 0 being the
544             astronomical Spring equinox (first point of Aries), and so on.
545              
546             The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
547             Edition, Chapter 27 ("Equinoxes and Solstices), pages 278ff.
548              
549             THIS METHOD IS UNSUPPORTED. I understand the temptation to call it if
550             all you want are the seasons, but if possible I would like to be able to
551             remove it if its use in next_quarter() turns out to e a bad idea. I am
552             not unwilling to support it; if you want me to, please contact me.
553              
554             Because it is unsupported, its name may change without warning if
555             L becomes smart enough to
556             realize that the =begin/end comment markers mean that this method is not
557             documented after all.
558              
559             =end comment
560              
561             =cut
562              
563             sub season {
564 11     11 1 50 my ( $self, $year, $season ) = @_;
565 11 50       108 my ( $Y, $d ) = $year < 1000 ? (
566             $year / 1000,
567             [ # Meeus table 27 A
568             [ -0.00071, 0.00111, 0.06134, 365242.13740, 1721139.29189 ],
569             [ 0.00025, 0.00907, -0.05323, 365241.72562, 1721233.25401 ],
570             [ 0.00074, -0.00297, -0.11677, 365242.49558, 1721325.70455 ],
571             [ -0.00006, -0.00933, -0.00769, 365242.88257, 1721414.39987 ],
572             ]->[ $season ],
573             ) : (
574             ( $year - 2000 ) / 1000,
575             [ # Meeus table 27 B
576             [ -0.00057, -0.00411, 0.05169, 365242.37404, 2451623.80984 ],
577             [ -0.00030, 0.00888, 0.00325, 365241.62603, 2451716.56767 ],
578             [ 0.00078, 0.00337, -0.11575, 365242.01767, 2451810.21715 ],
579             [ 0.00032, -0.00823, -0.06223, 365242.74049, 2451900.05952 ],
580             ]->[ $season ],
581             );
582 11         47 my $JDE0 = ( ( ( $d->[0] * $Y + $d->[1] ) * $Y + $d->[2] ) * $Y +
583             $d->[3] ) * $Y + $d->[4];
584 11         19 my $T = ( $JDE0 - 2451545.0 ) / 36525;
585             $self->{debug}
586 11 50       26 and print "Debug - T = $T\n";
587 11         33 my $W = mod2pi( deg2rad( 35999.373 * $T - 2.47 ) );
588 11         34 my $delta_lambda = 1 + 0.0334 * cos( $W ) + 0.0007 * cos( 2 * $W );
589             $self->{debug}
590 11 50       26 and print "Debug - delta lambda = $delta_lambda\n";
591 11         16 my $S = 0;
592 11         121 foreach my $term ( # Meeus table 27 C
593             [ 485, 324.96, 1934.136 ],
594             [ 203, 337.23, 32964.467 ],
595             [ 199, 342.08, 20.186 ],
596             [ 182, 27.85, 445267.112 ],
597             [ 156, 73.14, 45036.886 ],
598             [ 136, 171.52, 22518.443 ],
599             [ 77, 222.54, 65928.934 ],
600             [ 74, 296.72, 3034.906 ],
601             [ 70, 243.58, 9037.513 ],
602             [ 58, 119.81, 33718.147 ],
603             [ 52, 297.17, 150.678 ],
604             [ 50, 21.02, 2281.226 ],
605             [ 45, 247.54, 29929.562 ],
606             [ 44, 325.15, 31555.956 ],
607             [ 29, 60.93, 4443.417 ],
608             [ 18, 155.12, 67555.328 ],
609             [ 17, 288.79, 4562.452 ],
610             [ 16, 198.04, 62894.029 ],
611             [ 14, 199.76, 31436.921 ],
612             [ 12, 95.39, 14577.848 ],
613             [ 12, 287.11, 31931.756 ],
614             [ 12, 320.81, 34777.259 ],
615             [ 9, 227.73, 1222.114 ],
616             [ 8, 15.45, 16859.074 ],
617             ) {
618 264         556 $S += $term->[0] * cos( mod2pi( deg2rad( $term->[1] + $term->[2]
619             * $T ) ) );
620             }
621             $self->{debug}
622 11 50       56 and print "Debug - S = $S\n";
623 11         18 my $JDE = 0.00001 * $S / $delta_lambda + $JDE0;
624             $self->{debug}
625 11 50       23 and print "Debug - JDE = $JDE\n";
626 11         19 my $time = ( $JDE - JD_OF_EPOCH ) * SECSPERDAY;
627             # Note that gmtime() in the following needs the parens because it
628             # might have come from Time::y2038, which appears to take more than
629             # one argument -- even though, as I read it, its prototype is (;$).
630             $self->{debug}
631 11 50       21 and print "Debug - dynamical date is ", scalar gmtime( $time ), "\n";
632 11         27 return $time - dynamical_delta( $time ); # Not quite right.
633             }
634              
635             =item $sun->set( ... )
636              
637             This method has been overridden to silently ignore any attempt to set
638             the C<'sun'> attribute.
639              
640             =cut
641              
642             sub set {
643 75     75 1 215 my ( $self, @args ) = @_;
644 75         200 while ( @args ) {
645 188         472 my ( $name, $val ) = splice @args, 0, 2;
646 188 50       511 if ( 'sun' eq $name ) {
    100          
647             } elsif ( $attrib{$name} ) {
648 37 50       104 if ( ref $self ) {
649 37         132 $self->{$name} = $val;
650             } else {
651 0         0 $static{$name} = $val;
652             }
653             } else {
654 151         402 $self->SUPER::set( $name, $val );
655             }
656             }
657 75         159 return $self;
658             }
659              
660             =item $sun->time_set ()
661              
662             This method sets coordinates of the object to the coordinates of the
663             Sun at the object's currently-set universal time. The velocity
664             components are arbitrarily set to 0. The 'equinox_dynamical' attribute
665             is set to the object's currently-set dynamical time.
666              
667             Although there's no reason this method can't be called directly, it
668             exists to take advantage of the hook in the B
669             object, to allow the position of the Sun to be computed when the
670             object's time is set.
671              
672             The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
673             Edition, Chapter 25, pages 163ff.
674              
675             =cut
676              
677             # The following constants are used in the time_set() method,
678             # because Meeus' equations are in degrees, I was too lazy
679             # to hand-convert them to radians, but I didn't want to
680             # penalize the user for the conversion every time.
681              
682 16     16   133 use constant SUN_C1_0 => deg2rad (1.914602);
  16         37  
  16         76  
683 16     16   106 use constant SUN_C1_1 => deg2rad (-0.004817);
  16         36  
  16         67  
684 16     16   107 use constant SUN_C1_2 => deg2rad (-0.000014);
  16         41  
  16         71  
685 16     16   104 use constant SUN_C2_0 => deg2rad (0.019993);
  16         43  
  16         66  
686 16     16   100 use constant SUN_C2_1 => deg2rad (0.000101);
  16         41  
  16         62  
687 16     16   98 use constant SUN_C3_0 => deg2rad (0.000289);
  16         43  
  16         56  
688 16     16   103 use constant SUN_LON_2000 => deg2rad (- 0.01397);
  16         30  
  16         90  
689              
690             sub time_set {
691 11257     11257 1 16724 my $self = shift;
692 11257         23012 my $time = $self->dynamical;
693              
694             # The following algorithm is from Meeus, chapter 25, page, 163 ff.
695              
696 11257         26697 my $T = jcent2000 ($time); # Meeus (25.1)
697 11257         28045 my $L0 = mod2pi(deg2rad((.0003032 * $T + 36000.76983) * $T # Meeus (25.2)
698             + 280.46646));
699 11257         26261 my $M = mod2pi(deg2rad(((-.0001537) * $T + 35999.05029) # Meeus (25.3)
700             * $T + 357.52911));
701 11257         20817 my $e = (-0.0000001267 * $T - 0.000042037) * $T + 0.016708634;# Meeus (25.4)
702 11257         32573 my $C = ((SUN_C1_2 * $T + SUN_C1_1) * $T + SUN_C1_0) * sin ($M)
703             + (SUN_C2_1 * $T + SUN_C2_0) * sin (2 * $M)
704             + SUN_C3_0 * sin (3 * $M);
705 11257         20819 my $O = $self->{_sun_geometric_longitude} = $L0 + $C;
706 11257         24552 my $omega = mod2pi (deg2rad (125.04 - 1934.136 * $T));
707 11257         26548 my $lambda = mod2pi ($O - deg2rad (0.00569 + 0.00478 * sin ($omega)));
708 11257         19154 my $nu = $M + $C;
709 11257         26109 my $R = (1.000_001_018 * (1 - $e * $e)) / (1 + $e * cos ($nu))
710             * AU;
711 11257 50       25574 $self->{debug} and print <
712 0         0 Debug sun - @{[strftime '%d-%b-%Y %H:%M:%S', gmtime( $time )]}
713             T = $T
714 0         0 L0 = @{[rad2deg ($L0)]} degrees
715 0         0 M = @{[rad2deg ($M)]} degrees
716             e = $e
717 0         0 C = @{[rad2deg ($C)]} degrees
718 0         0 O = @{[rad2deg ($O)]} degrees
719 0         0 R = @{[$R / AU]} AU
720 0         0 omega = @{[rad2deg ($omega)]} degrees
721 0         0 lambda = @{[rad2deg ($lambda)]} degrees
722             eod
723              
724 11257         31606 $self->ecliptic (0, $lambda, $R);
725             ## $self->set (equinox_dynamical => $time);
726 11257         32926 $self->equinox_dynamical ($time);
727 11257         21470 return $self;
728             }
729              
730             # The Sun is normally positioned in inertial coordinates.
731              
732 37     37   157 sub __initial_inertial { return 1 }
733              
734             1;
735              
736             =back
737              
738             =head2 Attributes
739              
740             This class has the following public attributes. The description gives
741             the data type.
742              
743             =over
744              
745             =item iterate_for_quarters (Boolean)
746              
747             If this attribute is true, the C method uses the old
748             (pre-0.088_01) algorithm.
749              
750             If this attribute is false, the new algorithm is used.
751              
752             The default is C, i.e. false, because I believe the new algorithm
753             to be more accurate for reasonably-current times.
754              
755             This attribute is new with version 0.088_01.
756              
757             =back
758              
759             =head2 Historical Calculations
760              
761             This class was written for the purpose of calculating whether the Sun
762             was shining on a given point on the Earth (or in space) at a given time
763             in or reasonably close to the present. I can not say how accurate it is
764             at times far from the present. Those interested in such calculations may
765             want to consider using
766             L
767             instead.
768              
769             Historical calculations will need to use a 64-bit Perl, or at least one
770             with 64-bit integers, to represent times more than about 38 years from
771             the system epoch.
772              
773             In addition, you will need to be careful how you do input and output
774             conversions.
775              
776             L is a core module and an obvious choice, but
777             it only does Gregorian dates. Historical calculations prior to 1587
778             typically use the Julian calendar. For this you will need to go to
779             something like L,
780             or L which
781             does either Julian or Gregorian as needed.
782              
783             Should you decide to use L, you should be aware
784             that its C and C interpret the year argument
785             strangely: years in the range C<0 - 999> inclusive are not interpreted
786             as Gregorian years, though years outside that range are so interpreted.
787             Beginning with version 1.27 (released July 9 2018), additional
788             subroutines C and C were added.
789             These always interpret the year argument as a Gregorian year.
790              
791             =head1 ACKNOWLEDGMENTS
792              
793             The author wishes to acknowledge Jean Meeus, whose book "Astronomical
794             Algorithms" (second edition) formed the basis for this module.
795              
796             =head1 SEE ALSO
797              
798             The L
799             documentation for a discussion of how the pieces/parts of this
800             distribution go together and how to use them.
801              
802             L by Brett Hamilton, which contains a
803             function-based module to compute the current phase, distance and angular
804             diameter of the Moon, as well as the angular diameter and distance of
805             the Sun.
806              
807             L by Ron Hill and Jean Forget, which
808             contains a function-based module to compute sunrise and sunset for the
809             given day and location.
810              
811             L by Rob Fugina, which provides
812             functionality similar to B.
813              
814             =head1 SUPPORT
815              
816             Support is by the author. Please file bug reports at
817             L,
818             L, or in
819             electronic mail to the author.
820              
821             =head1 AUTHOR
822              
823             Thomas R. Wyant, III (F)
824              
825             =head1 COPYRIGHT AND LICENSE
826              
827             Copyright (C) 2005-2023 by Thomas R. Wyant, III
828              
829             This program is free software; you can redistribute it and/or modify it
830             under the same terms as Perl 5.10.0. For more details, see the full text
831             of the licenses in the directory LICENSES.
832              
833             This program is distributed in the hope that it will be useful, but
834             without any warranty; without even the implied warranty of
835             merchantability or fitness for a particular purpose.
836              
837             =cut
838              
839             # ex: set textwidth=72 :