| 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 |  | 1903 | use strict; | 
|  | 16 |  |  |  |  | 38 |  | 
|  | 16 |  |  |  |  | 561 |  | 
| 50 | 16 |  |  | 16 |  | 86 | use warnings; | 
|  | 16 |  |  |  |  | 33 |  | 
|  | 16 |  |  |  |  | 805 |  | 
| 51 |  |  |  |  |  |  |  | 
| 52 |  |  |  |  |  |  | our $VERSION = '0.129_01'; | 
| 53 |  |  |  |  |  |  |  | 
| 54 | 16 |  |  | 16 |  | 108 | use base qw{Astro::Coord::ECI}; | 
|  | 16 |  |  |  |  | 33 |  | 
|  | 16 |  |  |  |  | 2020 |  | 
| 55 |  |  |  |  |  |  |  | 
| 56 | 16 |  |  | 16 |  | 111 | use Astro::Coord::ECI::Utils qw{ @CARP_NOT :mainstream }; | 
|  | 16 |  |  |  |  | 49 |  | 
|  | 16 |  |  |  |  | 7522 |  | 
| 57 | 16 |  |  | 16 |  | 130 | use Carp; | 
|  | 16 |  |  |  |  | 46 |  | 
|  | 16 |  |  |  |  | 1066 |  | 
| 58 | 16 |  |  | 16 |  | 101 | use POSIX qw{ ceil floor strftime }; | 
|  | 16 |  |  |  |  | 30 |  | 
|  | 16 |  |  |  |  | 137 |  | 
| 59 |  |  |  |  |  |  |  | 
| 60 | 16 |  |  | 16 |  | 1587 | use constant MEAN_MAGNITUDE => -26.8; | 
|  | 16 |  |  |  |  | 41 |  | 
|  | 16 |  |  |  |  | 6773 |  | 
| 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 | 1719 | my ($class, @args) = @_; | 
| 109 | 134 | 50 |  |  |  | 343 | ref $class and $class = ref $class; | 
| 110 | 134 | 100 | 66 |  |  | 694 | if ( $Singleton && $weaken && __classisa( $class, __PACKAGE__ ) ) { | 
|  |  |  | 66 |  |  |  |  | 
| 111 | 132 |  |  |  |  | 229 | my $self; | 
| 112 | 132 | 100 |  |  |  | 327 | if ( $self = $object{$class} ) { | 
| 113 | 97 | 100 |  |  |  | 224 | $self->set( @args ) if @args; | 
| 114 | 97 |  |  |  |  | 277 | return $self; | 
| 115 |  |  |  |  |  |  | } else { | 
| 116 | 35 |  |  |  |  | 247 | $self = $object{$class} = $class->SUPER::new (%static, @args); | 
| 117 | 35 |  |  |  |  | 174 | $weaken->( $object{$class} ); | 
| 118 | 35 |  |  |  |  | 150 | return $self; | 
| 119 |  |  |  |  |  |  | } | 
| 120 |  |  |  |  |  |  | } else { | 
| 121 | 2 |  |  |  |  | 10 | 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 |  |  |  |  | 7 | my $horizon = $station->__get_almanac_horizon(); | 
| 170 |  |  |  |  |  |  |  | 
| 171 | 2 |  |  |  |  | 17 | 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 |  | 51 | and return @{ $events[$inx++] }; | 
|  | 8 |  |  |  |  | 70 |  | 
| 185 | 2 |  |  |  |  | 8 | return; | 
| 186 | 2 |  |  |  |  | 11 | }; | 
| 187 |  |  |  |  |  |  | } | 
| 188 |  |  |  |  |  |  |  | 
| 189 | 16 |  |  | 16 |  | 6289 | use Astro::Coord::ECI::Mixin qw{ almanac }; | 
|  | 16 |  |  |  |  | 42 |  | 
|  | 16 |  |  |  |  | 1045 |  | 
| 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 |  | 108 | use Astro::Coord::ECI::Mixin qw{ almanac_hash }; | 
|  | 16 |  |  |  |  | 32 |  | 
|  | 16 |  |  |  |  | 9166 |  | 
| 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 | 464 | 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 | 2117 | my $self = shift; | 
| 261 | 1125 | 50 |  |  |  | 2342 | 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 |  |  |  |  | 2330 | 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 | 13910 | my ( $self, @args ) = @_; | 
| 278 | 6742 |  |  |  |  | 10288 | my @rslt; | 
| 279 | 6742 |  |  |  |  | 11329 | foreach my $name ( @args ) { | 
| 280 |  |  |  |  |  |  | push @rslt, 'sun' eq $name ? $self : | 
| 281 |  |  |  |  |  |  | $attrib{$name} ? | 
| 282 | 6742 | 0 |  |  |  | 23355 | ref $self ? $self->{$name} : $static{$name} : | 
|  |  | 50 |  |  |  |  |  | 
|  |  | 50 |  |  |  |  |  | 
| 283 |  |  |  |  |  |  | $self->SUPER::get( $name ); | 
| 284 |  |  |  |  |  |  | } | 
| 285 | 6742 | 100 |  |  |  | 19118 | 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 |  | 128 | use constant NEXT_QUARTER_INCREMENT => 86400 * 85;	# 85 days. | 
|  | 16 |  |  |  |  | 34 |  | 
|  | 16 |  |  |  |  | 5204 |  | 
| 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 | 965 | my ( $self, $quarter ) = @_; | 
| 420 |  |  |  |  |  |  | $self->{iterate_for_quarters} | 
| 421 | 7 | 50 |  |  |  | 22 | and goto &Astro::Coord::ECI::Mixin::next_quarter; | 
| 422 | 7 |  |  |  |  | 22 | my $time = $self->universal(); | 
| 423 | 7 |  |  |  |  | 33 | my $year = ( gmtime( $time ) )[5] + 1900; | 
| 424 | 7 |  |  |  |  | 16 | my $season; | 
| 425 | 7 | 50 |  |  |  | 38 | 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 |  |  |  |  | 21 | 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 |  |  |  |  | 42 | $quarter = ceil( ( rad2deg( $lon ) + 359 ) / 90 ) % 4; | 
| 442 | 7 |  |  |  |  | 29 | $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 |  |  |  | 24 | if ( $time - $season > SECSPERDAY * 180 ) { | 
| 446 | 1 |  |  |  |  | 2 | $year++; | 
| 447 | 1 |  |  |  |  | 4 | $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 |  |  |  | 18 | 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 |  |  |  |  | 20 | $season = ceil( $season );	# Make sure we're AFTER. | 
| 462 | 7 |  |  |  |  | 24 | $self->universal( $season ); | 
| 463 | 7 | 100 |  |  |  | 38 | 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 |  | 156 | use Astro::Coord::ECI::Mixin qw{ next_quarter_hash }; | 
|  | 16 |  |  |  |  | 61 |  | 
|  | 16 |  |  |  |  | 14169 |  | 
| 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 | 979 | sub period {return 31558149.7632}	# 365.256363 * 86400 | 
| 498 |  |  |  |  |  |  |  | 
| 499 |  |  |  |  |  |  | sub __horizon_name_tplt { | 
| 500 | 4 |  |  | 4 |  | 19 | my ( $self ) = @_; | 
| 501 | 4 | 50 |  |  |  | 23 | return $self->__object_is_self_named() ? | 
| 502 |  |  |  |  |  |  | [ '%sset', '%srise' ] : | 
| 503 |  |  |  |  |  |  | $self->SUPER::__horizon_name_tplt(); | 
| 504 |  |  |  |  |  |  | } | 
| 505 |  |  |  |  |  |  |  | 
| 506 |  |  |  |  |  |  | sub __quarter_name { | 
| 507 | 2 |  |  | 2 |  | 6 | my ( $self, $event, $tplt ) = @_; | 
| 508 | 2 | 50 | 33 |  |  | 12 | $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 |  |  |  |  | 4 | 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 |  |  |  |  | 9 | return $self->__event_name( $event, $tplt ); | 
| 521 |  |  |  |  |  |  | } | 
| 522 |  |  |  |  |  |  |  | 
| 523 |  |  |  |  |  |  | sub __transit_name_tplt { | 
| 524 | 4 |  |  | 4 |  | 9 | my ( $self ) = @_; | 
| 525 | 4 | 50 |  |  |  | 25 | 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 |  | 14 | my ( $self, $event, $tplt ) = @_; | 
| 532 | 4 | 50 | 33 |  |  | 22 | $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 | 26 | my ( $self, $year, $season ) = @_; | 
| 565 | 11 | 50 |  |  |  | 95 | 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 |  |  |  |  | 48 | my $JDE0 = ( ( ( $d->[0] * $Y + $d->[1] ) * $Y + $d->[2] ) * $Y + | 
| 583 |  |  |  |  |  |  | $d->[3] ) * $Y + $d->[4]; | 
| 584 | 11 |  |  |  |  | 21 | my $T = ( $JDE0 - 2451545.0 ) / 36525; | 
| 585 |  |  |  |  |  |  | $self->{debug} | 
| 586 | 11 | 50 |  |  |  | 31 | and print "Debug - T = $T\n"; | 
| 587 | 11 |  |  |  |  | 30 | my $W = mod2pi( deg2rad( 35999.373 * $T - 2.47 ) ); | 
| 588 | 11 |  |  |  |  | 35 | 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 |  |  |  |  | 17 | my $S = 0; | 
| 592 | 11 |  |  |  |  | 124 | 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 |  |  |  |  | 542 | $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 |  |  |  |  | 22 | my $JDE = 0.00001 * $S / $delta_lambda + $JDE0; | 
| 624 |  |  |  |  |  |  | $self->{debug} | 
| 625 | 11 | 50 |  |  |  | 25 | and print "Debug - JDE = $JDE\n"; | 
| 626 | 11 |  |  |  |  | 17 | 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 |  |  |  | 25 | and print "Debug - dynamical date is ", scalar gmtime( $time ), "\n"; | 
| 632 | 11 |  |  |  |  | 77 | 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 | 223 | my ( $self, @args ) = @_; | 
| 644 | 75 |  |  |  |  | 200 | while ( @args ) { | 
| 645 | 188 |  |  |  |  | 407 | my ( $name, $val ) = splice @args, 0, 2; | 
| 646 | 188 | 50 |  |  |  | 536 | if ( 'sun' eq $name ) { | 
|  |  | 100 |  |  |  |  |  | 
| 647 |  |  |  |  |  |  | } elsif ( $attrib{$name} ) { | 
| 648 | 37 | 50 |  |  |  | 151 | if ( ref $self ) { | 
| 649 | 37 |  |  |  |  | 136 | $self->{$name} = $val; | 
| 650 |  |  |  |  |  |  | } else { | 
| 651 | 0 |  |  |  |  | 0 | $static{$name} = $val; | 
| 652 |  |  |  |  |  |  | } | 
| 653 |  |  |  |  |  |  | } else { | 
| 654 | 151 |  |  |  |  | 393 | $self->SUPER::set( $name, $val ); | 
| 655 |  |  |  |  |  |  | } | 
| 656 |  |  |  |  |  |  | } | 
| 657 | 75 |  |  |  |  | 167 | 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 |  |  |  |  | 42 |  | 
|  | 16 |  |  |  |  | 70 |  | 
| 683 | 16 |  |  | 16 |  | 108 | use constant SUN_C1_1 => deg2rad (-0.004817); | 
|  | 16 |  |  |  |  | 36 |  | 
|  | 16 |  |  |  |  | 57 |  | 
| 684 | 16 |  |  | 16 |  | 105 | use constant SUN_C1_2 => deg2rad (-0.000014); | 
|  | 16 |  |  |  |  | 33 |  | 
|  | 16 |  |  |  |  | 56 |  | 
| 685 | 16 |  |  | 16 |  | 110 | use constant SUN_C2_0 => deg2rad (0.019993); | 
|  | 16 |  |  |  |  | 31 |  | 
|  | 16 |  |  |  |  | 77 |  | 
| 686 | 16 |  |  | 16 |  | 103 | use constant SUN_C2_1 => deg2rad (0.000101); | 
|  | 16 |  |  |  |  | 50 |  | 
|  | 16 |  |  |  |  | 81 |  | 
| 687 | 16 |  |  | 16 |  | 98 | use constant SUN_C3_0 => deg2rad (0.000289); | 
|  | 16 |  |  |  |  | 37 |  | 
|  | 16 |  |  |  |  | 58 |  | 
| 688 | 16 |  |  | 16 |  | 105 | use constant SUN_LON_2000 => deg2rad (- 0.01397); | 
|  | 16 |  |  |  |  | 33 |  | 
|  | 16 |  |  |  |  | 76 |  | 
| 689 |  |  |  |  |  |  |  | 
| 690 |  |  |  |  |  |  | sub time_set { | 
| 691 | 11257 |  |  | 11257 | 1 | 16741 | my $self = shift; | 
| 692 | 11257 |  |  |  |  | 22433 | my $time = $self->dynamical; | 
| 693 |  |  |  |  |  |  |  | 
| 694 |  |  |  |  |  |  | #	The following algorithm is from Meeus, chapter 25, page, 163 ff. | 
| 695 |  |  |  |  |  |  |  | 
| 696 | 11257 |  |  |  |  | 26218 | my $T = jcent2000 ($time);				# Meeus (25.1) | 
| 697 | 11257 |  |  |  |  | 27363 | my $L0 = mod2pi(deg2rad((.0003032 * $T + 36000.76983) * $T	# Meeus (25.2) | 
| 698 |  |  |  |  |  |  | + 280.46646)); | 
| 699 | 11257 |  |  |  |  | 25542 | my $M = mod2pi(deg2rad(((-.0001537) * $T + 35999.05029)	# Meeus (25.3) | 
| 700 |  |  |  |  |  |  | * $T + 357.52911)); | 
| 701 | 11257 |  |  |  |  | 20770 | my $e = (-0.0000001267 * $T - 0.000042037) * $T + 0.016708634;# Meeus (25.4) | 
| 702 | 11257 |  |  |  |  | 31251 | 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 |  |  |  |  | 19789 | my $O = $self->{_sun_geometric_longitude} = $L0 + $C; | 
| 706 | 11257 |  |  |  |  | 23516 | my $omega = mod2pi (deg2rad (125.04 - 1934.136 * $T)); | 
| 707 | 11257 |  |  |  |  | 27567 | my $lambda = mod2pi ($O - deg2rad (0.00569 + 0.00478 * sin ($omega))); | 
| 708 | 11257 |  |  |  |  | 18018 | my $nu = $M + $C; | 
| 709 | 11257 |  |  |  |  | 24130 | my $R = (1.000_001_018 * (1 - $e * $e)) / (1 + $e * cos ($nu)) | 
| 710 |  |  |  |  |  |  | * AU; | 
| 711 | 11257 | 50 |  |  |  | 23294 | $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 |  |  |  |  | 30392 | $self->ecliptic (0, $lambda, $R); | 
| 725 |  |  |  |  |  |  | ## $self->set (equinox_dynamical => $time); | 
| 726 | 11257 |  |  |  |  | 32868 | $self->equinox_dynamical ($time); | 
| 727 | 11257 |  |  |  |  | 20893 | 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 : |