1# -*- encoding: utf-8; indent-tabs-mode: nil -*- 2# 3# Perl DateTime extension for computing the sunrise/sunset on a given day 4# Copyright © 1999-2004, 2013-2014, 2020 Ron Hill and Jean Forget 5# 6# See the license in the embedded documentation below. 7# 8package DateTime::Event::Sunrise; 9 10use strict; 11use warnings; 12require Exporter; 13use POSIX qw(floor); 14use Math::Trig; 15use Carp; 16use DateTime; 17use DateTime::Set; 18use Params::Validate qw(:all); 19use Set::Infinite qw(inf $inf); 20use vars qw( $VERSION $RADEG $DEGRAD @ISA ); 21@ISA = qw( Exporter ); 22$VERSION = '0.0506'; 23$RADEG = ( 180 / pi ); 24$DEGRAD = ( pi / 180 ); 25my $INV360 = ( 1.0 / 360.0 ); 26 27# Julian day number for the 0th January 2000 (that is, 31st December 1999) 28my $jd_2000_Jan_0 = DateTime->new(year => 1999, month => 12, day => 31, time_zone => 'UTC')->jd; 29 30 31sub new { 32 my $class = shift; 33 34 if (@_ % 2 != 0) { 35 croak "Odd number of parameters"; 36 } 37 my %args = @_; 38 if (exists $args{iteration} && exists $args{precise}) { 39 croak "Parameter 'iteration' is deprecated, use only 'precise'"; 40 } 41 42 %args = validate( 43 @_, { 44 longitude => { type => SCALAR, optional => 1, default => 0 }, 45 latitude => { type => SCALAR, optional => 1, default => 0 }, 46 altitude => { 47 type => SCALAR, 48 default => '-0.833', 49 regex => qr/^(-?\d+(?:\.\d+)?)$/ 50 }, 51 iteration => { type => SCALAR, default => '0' }, 52 precise => { type => SCALAR, default => '0' }, 53 upper_limb => { type => SCALAR, default => '0' }, 54 silent => { type => SCALAR, default => '0' }, 55 trace => { type => GLOB | GLOBREF | SCALAR, default => '0' }, 56 } 57 ); 58 59 # Making old and new parameters synonymous 60 unless (exists $args{precise}) { 61 $args{precise} = $args{iteration}; 62 } 63 # TODO : get rid of the old parameters after this point 64 $args{iteration} = $args{precise}; 65 66 return bless \%args, $class; 67} 68 69 # 70 # 71 # FUNCTIONAL SEQUENCE for sunrise 72 # 73 # _GIVEN 74 # A sunrise object that was created by the new method 75 # 76 # _THEN 77 # 78 # setup subs for following/previous sunrise times 79 # 80 # 81 # _RETURN 82 # 83 # A new DateTime::Set recurrence object 84 # 85sub sunrise { 86 87 my $class = shift; 88 my $self = $class->new(@_); 89 return DateTime::Set->from_recurrence( 90 next => sub { 91 return $_[0] if $_[0]->is_infinite; 92 $self->_following_sunrise( $_[0] ); 93 }, 94 previous => sub { 95 return $_[0] if $_[0]->is_infinite; 96 $self->_previous_sunrise( $_[0] ); 97 } ); 98} 99 100 # 101 # 102 # FUNCTIONAL SEQUENCE for sunset 103 # 104 # _GIVEN 105 # 106 # A sunrise object that was created by the new method 107 # _THEN 108 # 109 # Setup subs for following/previous sunset times 110 # 111 # 112 # _RETURN 113 # 114 # A new DateTime::Set recurrence object 115 # 116sub sunset { 117 118 my $class = shift; 119 my $self = $class->new(@_); 120 return DateTime::Set->from_recurrence( 121 next => sub { 122 return $_[0] if $_[0]->is_infinite; 123 $self->_following_sunset( $_[0] ); 124 }, 125 previous => sub { 126 return $_[0] if $_[0]->is_infinite; 127 $self->_previous_sunset( $_[0] ); 128 } ); 129} 130 131 # 132 # 133 # FUNCTIONAL SEQUENCE for sunset_datetime 134 # 135 # _GIVEN 136 # 137 # A sunrise object 138 # A DateTime object 139 # 140 # _THEN 141 # 142 # Validate the DateTime object is valid 143 # Compute sunrise and sunset 144 # 145 # 146 # _RETURN 147 # 148 # DateTime object that contains the sunset time 149 # 150sub sunset_datetime { 151 152 my $self = shift; 153 my $dt = shift; 154 my $class = ref($dt); 155 156 if ( ! $dt->isa('DateTime') ) { 157 croak("Dates need to be DateTime objects"); 158 } 159 my ( undef, $tmp_set ) = _sunrise( $self, $dt, 0, 1 ); 160 return $tmp_set; 161} 162 163 # 164 # 165 # FUNCTIONAL SEQUENCE for sunrise_datetime 166 # 167 # _GIVEN 168 # 169 # A sunrise object 170 # A DateTime object 171 # 172 # _THEN 173 # 174 # Validate the DateTime object is valid 175 # Compute sunrise and sunset 176 # 177 # 178 # _RETURN 179 # 180 # DateTime object that contains the sunrise times 181 # 182sub sunrise_datetime { 183 184 my $self = shift; 185 my $dt = shift; 186 my $class = ref($dt); 187 188 if ( ! $dt->isa('DateTime') ) { 189 croak("Dates need to be DateTime objects"); 190 } 191 my ( $tmp_rise, undef ) = _sunrise( $self, $dt, 1, 0 ); 192 return $tmp_rise; 193} 194 195 # 196 # 197 # FUNCTIONAL SEQUENCE for sunrise_sunset_span 198 # 199 # _GIVEN 200 # 201 # A sunrise object 202 # A DateTime object 203 # 204 # _THEN 205 # 206 # Validate the DateTime object is valid 207 # Compute sunrise and sunset 208 # 209 # 210 # _RETURN 211 # 212 # DateTime Span object that contains the sunrise/sunset times 213 # 214sub sunrise_sunset_span { 215 216 my $self = shift; 217 my $dt = shift; 218 my $class = ref($dt); 219 220 if ( ! $dt->isa('DateTime') ) { 221 croak("Dates need to be DateTime objects"); 222 } 223 my ( $tmp_rise, $tmp_set ) = _sunrise( $self, $dt, 1, 1 ); 224 225 return DateTime::Span->from_datetimes( 226 start => $tmp_rise, 227 end => $tmp_set 228 ); 229} 230 231# 232# FUNCTIONAL SEQUENCE for is_polar_night 233# 234# _GIVEN 235# 236# A sunrise object 237# A DateTime object 238# 239# _THEN 240# 241# Validate the DateTime object is valid 242# Compute sunrise and sunset 243# 244# _RETURN 245# 246# A boolean flag telling whether the sun will stay under the horizon or not 247# 248sub is_polar_night { 249 250 my $self = shift; 251 my $dt = shift; 252 my $class = ref($dt); 253 254 if ( ! $dt->isa('DateTime') ) { 255 croak("Dates need to be DateTime objects"); 256 } 257 my ( undef, undef, $rise_season, $set_season ) = _sunrise( $self, $dt, 1, 1, 1 ); 258 return ($rise_season < 0 || $set_season < 0); 259} 260 261# 262# FUNCTIONAL SEQUENCE for is_polar_day 263# 264# _GIVEN 265# 266# A sunrise object 267# A DateTime object 268# 269# _THEN 270# 271# Validate the DateTime object is valid 272# Compute sunrise and sunset 273# 274# _RETURN 275# 276# A boolean flag telling whether the sun will stay above the horizon or not 277# 278sub is_polar_day { 279 280 my $self = shift; 281 my $dt = shift; 282 my $class = ref($dt); 283 284 if ( ! $dt->isa('DateTime') ) { 285 croak("Dates need to be DateTime objects"); 286 } 287 my ( undef, undef, $rise_season, $set_season ) = _sunrise( $self, $dt, 1, 1, 1 ); 288 return ($rise_season > 0 || $set_season > 0); 289} 290 291# 292# FUNCTIONAL SEQUENCE for is_day_and_night 293# 294# _GIVEN 295# 296# A sunrise object 297# A DateTime object 298# 299# _THEN 300# 301# Validate the DateTime object is valid 302# Compute sunrise and sunset 303# 304# _RETURN 305# 306# A boolean flag telling whether the sun will rise and set or not 307# 308sub is_day_and_night { 309 310 my $self = shift; 311 my $dt = shift; 312 my $class = ref($dt); 313 314 if ( ! $dt->isa('DateTime') ) { 315 croak("Dates need to be DateTime objects"); 316 } 317 my ( undef, undef, $rise_season, $set_season ) = _sunrise( $self, $dt, 1, 1, 1 ); 318 return ($rise_season == 0 && $set_season == 0); 319} 320 321 # 322 # 323 # FUNCTIONAL SEQUENCE for _following_sunrise 324 # 325 # _GIVEN 326 # 327 # A sunrise object 328 # A DateTime object 329 # 330 # _THEN 331 # 332 # Validate the DateTime object is valid 333 # Compute sunrise and return if it is greater 334 # than the original if not add one day and recompute 335 # 336 # _RETURN 337 # 338 # A new DateTime object that contains the sunrise time 339 # 340sub _following_sunrise { 341 342 my $self = shift; 343 my $dt = shift; 344 croak( "Dates need to be DateTime objects (" . ref($dt) . ")" ) 345 unless ( $dt->isa('DateTime') ); 346 my ( $tmp_rise, undef ) = _sunrise( $self, $dt, 1, 0 ); 347 return $tmp_rise if $tmp_rise > $dt; 348 my $d = DateTime::Duration->new( 349 days => 1, 350 ); 351 my $new_dt = $dt + $d; 352 ( $tmp_rise, undef ) = _sunrise( $self, $new_dt, 1, 0 ); 353 return $tmp_rise if $tmp_rise > $dt; 354 $new_dt = $new_dt + $d; 355 ( $tmp_rise, undef ) = _sunrise( $self, $new_dt, 1, 0 ); 356 return $tmp_rise; 357} 358 359 # 360 # 361 # FUNCTIONAL SEQUENCE for _previous_sunrise 362 # 363 # _GIVEN 364 # A sunrise object 365 # A DateTime object 366 # 367 # _THEN 368 # 369 # Validate the DateTime Object 370 # Compute sunrise and return if it is less than 371 # the original object if not subtract one day and recompute 372 # 373 # _RETURN 374 # 375 # A new DateTime Object that contains the sunrise time 376 # 377sub _previous_sunrise { 378 379 my $self = shift; 380 my $dt = shift; 381 croak( "Dates need to be DateTime objects (" . ref($dt) . ")" ) 382 unless ( $dt->isa('DateTime') ); 383 my ( $tmp_rise, undef ) = _sunrise( $self, $dt, 1, 0 ); 384 return $tmp_rise if $tmp_rise < $dt; 385 my $d = DateTime::Duration->new( 386 days => 1, 387 ); 388 my $new_dt = $dt - $d; 389 ( $tmp_rise, undef ) = _sunrise( $self, $new_dt, 1, 0 ); 390 return $tmp_rise if $tmp_rise < $dt; 391 $new_dt = $new_dt - $d; 392 ( $tmp_rise, undef ) = _sunrise( $self, $new_dt, 1, 0 ); 393 return $tmp_rise; 394} 395 396 # 397 # 398 # FUNCTIONAL SEQUENCE for _following_sunset 399 # 400 # _GIVEN 401 # A sunrise object 402 # A DateTime object 403 # 404 # _THEN 405 # 406 # Validate the DateTime object is valid 407 # Compute sunset and return if it is greater 408 # than the original if not add one day and recompute 409 # 410 # _RETURN 411 # 412 # A DateTime object with sunset time 413 # 414sub _following_sunset { 415 416 my $self = shift; 417 my $dt = shift; 418 croak( "Dates need to be DateTime objects (" . ref($dt) . ")" ) 419 unless ( ref($dt) eq 'DateTime' ); 420 my ( undef, $tmp_set ) = _sunrise( $self, $dt, 0, 1 ); 421 return $tmp_set if $tmp_set > $dt; 422 my $d = DateTime::Duration->new( 423 days => 1, 424 ); 425 my $new_dt = $dt + $d; 426 ( undef, $tmp_set ) = _sunrise( $self, $new_dt, 0, 1 ); 427 return $tmp_set if $tmp_set > $dt; 428 $new_dt = $new_dt + $d; 429 ( undef, $tmp_set ) = _sunrise( $self, $new_dt, 0, 1 ); 430 return $tmp_set; 431} 432 433 # 434 # 435 # FUNCTIONAL SEQUENCE for _previous_sunset 436 # 437 # _GIVEN 438 # A sunrise object 439 # A DateTime object 440 # 441 # _THEN 442 # 443 # Validate the DateTime Object 444 # Compute sunset and return if it is less than 445 # the original object if not subtract one day and recompute 446 # 447 # _RETURN 448 # 449 # A DateTime object with sunset time 450 # 451sub _previous_sunset { 452 453 my $self = shift; 454 my $dt = shift; 455 croak( "Dates need to be DateTime objects (" . ref($dt) . ")" ) 456 unless ( $dt->isa('DateTime') ); 457 my ( undef, $tmp_set ) = _sunrise( $self, $dt, 0, 1 ); 458 return $tmp_set if $tmp_set < $dt; 459 my $d = DateTime::Duration->new( 460 days => 1, 461 ); 462 my $new_dt = $dt - $d; 463 ( undef, $tmp_set ) = _sunrise( $self, $new_dt, 0, 1 ); 464 return $tmp_set if $tmp_set < $dt; 465 $new_dt = $new_dt - $d; 466 ( undef, $tmp_set ) = _sunrise( $self, $new_dt, 0, 1 ); 467 return $tmp_set; 468} 469 470# 471# 472# FUNCTIONAL SEQUENCE for _sunrise 473# 474# _GIVEN 475# A sunrise object and a DateTime object 476# three booleans, to control the iterative computations and the warning messages 477# 478# _THEN 479# 480# Check if precise is set to one if so 481# initially compute sunrise/sunset (using division 482# by 15.04107 instead of 15.0) then recompute rise/set time 483# using exact moment last computed. IF precise is set 484# to zero divide by 15.0 (only once) 485# UPDATE: actually, with the precise algorithm as currently implemented 486# the 15.0 value gives better results than the 15.04107 value. Results 487# cross-checked with the NOAA's solar calculator, with Astro::Coords + Astro::PAL 488# and with Stellarium 489# 490# If using the precise algorithm, the $want_sunrise and $want_sunset booleans control the computation 491# of the corresponding events, to eliminate computations that will be discarded upon return 492# from the sub (that is, "stored" into undef). 493# These booleans are not used for the basic algorithm. 494# 495# The $silent boolean, if provided, override the silent attribute of the sunrise object 496# to control the emission of warnings. 497# 498# _RETURN 499# 500# two DateTime objects with the date and time for sunrise and sunset 501# two season flags for sunrise and sunset respectively 502# 503sub _sunrise { 504 505 my ($self, $dt, $want_sunrise, $want_sunset, $silent) = @_; 506 my $cloned_dt = $dt->clone; 507 my $altit = $self->{altitude}; 508 my $precise = defined( $self->{precise} ) ? $self->{precise} : 0; 509 my $trace = defined( $self->{trace} ) ? $self->{trace} : 0; 510 unless (defined $silent) { 511 $silent = defined( $self->{silent} ) ? $self->{silent} : 0; 512 } 513 $cloned_dt->set_time_zone('floating'); 514 515 if (!$precise) { 516 if ($trace) { 517 printf $trace "\nBasic computation for %s, lon %.3f, lat %.3f, altitude %.3f, upper limb %d\n" 518 , $dt->ymd 519 , $self->{longitude} 520 , $self->{latitude} 521 , $self->{altitude} 522 , $self->{upper_limb}; 523 } 524 my $d = days_since_2000_Jan_0($cloned_dt) + 0.5 - $self->{longitude} / 360.0; 525 my $revsub = \&rev180; # normalizing angles around 0 degrees 526 my ( $h1, $h2, $season ) = _sunrise_sunset( $d 527 , $self->{longitude} 528 , $self->{latitude} 529 , $altit 530 , 15.0 531 , $self->{upper_limb} 532 , $silent 533 , $trace 534 , $revsub); 535 my ( $seconds_rise, $seconds_set ) = convert_hour( $h1, $h2 ); 536 my $rise_dur = DateTime::Duration->new( seconds => $seconds_rise ); 537 my $set_dur = DateTime::Duration->new( seconds => $seconds_set ); 538 my $tmp_dt1 = DateTime->new( 539 year => $dt->year, 540 month => $dt->month, 541 day => $dt->day, 542 hour => 0, 543 minute => 0, 544 time_zone => 'UTC' 545 ); 546 547 my $rise_time = $tmp_dt1 + $rise_dur; 548 my $set_time = $tmp_dt1 + $set_dur; 549 my $tz = $dt->time_zone; 550 $rise_time->set_time_zone($tz) unless $tz->is_floating; 551 $set_time ->set_time_zone($tz) unless $tz->is_floating; 552 return ( $rise_time, $set_time, $season, $season ); 553 } 554 555 else { 556 my $ang_speed = 15.0; 557 my $d = days_since_2000_Jan_0($cloned_dt) - $self->{longitude} / 360.0; # UTC decimal days at midnight LMT 558 my $tmp_dt1 = DateTime->new( 559 year => $dt->year, 560 month => $dt->month, 561 day => $dt->day, 562 hour => 0, 563 minute => 0, 564 time_zone => 'UTC' 565 ); 566 my $tz = $dt->time_zone; 567 my $rise_time; 568 my $set_time; 569 my $rise_season; 570 my $set_season; 571 my $revsub = sub { _rev_lon($_[0], $self->{longitude}) }; # normalizing angles around the local longitude 572 573 if ($want_sunrise) { 574 if ($trace) { 575 printf $trace "\nPrecise sunrise computation for %s, lon %.3f, lat %.3f, altitude %.3f, upper limb %d angular speed %.5f\n" 576 , $dt->ymd 577 , $self->{longitude} 578 , $self->{latitude} 579 , $self->{altitude} 580 , $self->{upper_limb} 581 , $ang_speed; 582 } 583 # This is the initial start 584 585 my $h1_lmt = 12; # LMT decimal hours, noon then the successive values of sunrise 586 my $h1_utc; # UTC decimal hours, noon LMT then the successive values of sunrise 587 for my $counter (1..9) { 588 # 9 is a arbitrary value to stop runaway loops. Normally, we should leave at the second or third iteration 589 my $h2_utc; 590 ($h2_utc, undef, $rise_season) = _sunrise_sunset( $d + $h1_lmt / 24 591 , $self->{longitude} 592 , $self->{latitude} 593 , $altit 594 , $ang_speed 595 , $self->{upper_limb} 596 , $silent 597 , $trace 598 , $revsub); 599 if ($rise_season != 0) { 600 $h1_utc = $h2_utc; 601 last; 602 } 603 $h1_utc = $h1_lmt - $self->{longitude} / 15; 604 if (equal($h1_utc, $h2_utc, 5)) { 605 # equal within 1e-5 hour, less than 0.04 second 606 $h1_utc = $h2_utc; 607 last; 608 } 609 $h1_utc = $h2_utc; 610 $h1_lmt = $h1_utc + $self->{longitude} / 15; 611 } 612 my $second_rise = _convert_1_hour($h1_utc); 613 # This is to fix the datetime object to use a duration 614 # instead of blindly setting the hour/min 615 my $rise_dur = DateTime::Duration->new( seconds => $second_rise ); 616 $rise_time = $tmp_dt1 + $rise_dur; 617 $rise_time->set_time_zone($tz) unless $tz->is_floating; 618 } 619 620 if ($want_sunset) { 621 if ($trace) { 622 printf $trace "\nPrecise sunset computation for %s, lon %.3f, lat %.3f, altitude %.3f, upper limb %d angular speed %.5f\n" 623 , $dt->ymd 624 , $self->{longitude} 625 , $self->{latitude} 626 , $self->{altitude} 627 , $self->{upper_limb} 628 , $ang_speed; 629 } 630 my $h3_lmt = 12; # LMT decimal hours, noon then the successive values of sunset 631 my $h3_utc; # UTC decimal hours, noon LMT then the successive values of sunset 632 for my $counter (1..9) { 633 # 9 is a arbitrary value to stop runaway loops. Normally, we should leave at the second or third iteration 634 my $h4_utc; 635 (undef, $h4_utc, $set_season) = _sunrise_sunset( $d + $h3_lmt / 24 636 , $self->{longitude} 637 , $self->{latitude} 638 , $altit 639 , $ang_speed 640 , $self->{upper_limb} 641 , $silent 642 , $trace 643 , $revsub); 644 if ($set_season != 0) { 645 $h3_utc = $h4_utc; 646 last; 647 } 648 $h3_utc = $h3_lmt - $self->{longitude} / 15; 649 if (equal($h3_utc, $h4_utc, 5)) { 650 # equal within 1e-5 hour, less than 0.04 second 651 $h3_utc = $h4_utc; 652 last; 653 } 654 $h3_utc = $h4_utc; 655 $h3_lmt = $h3_utc + $self->{longitude} / 15; 656 } 657 658 my $second_set = _convert_1_hour( $h3_utc ); 659 # This is to fix the datetime object to use a duration 660 # instead of blindly setting the hour/min 661 my $set_dur = DateTime::Duration->new( seconds => $second_set ); 662 $set_time = $tmp_dt1 + $set_dur; 663 $set_time ->set_time_zone($tz) unless $tz->is_floating; 664 } 665 666 return ( $rise_time, $set_time, $rise_season, $set_season ); 667 } 668 669} 670 671# 672# 673# FUNCTIONAL SEQUENCE for _sunrise_sunset 674# 675# _GIVEN 676# 677# days since Jan 0 2000, 678# the fractional part is the UTC time of day. E.g. 7458.66667 represents 2020-06-01T16:00:00 UTC 679# longitude, 680# latitude, 681# reference sun height, 682# all three in decimal degrees 683# angular speed 684# either Earth's spin 15.04107 degrees per hour or the combination or Earth's spin with Earth-Sun mean orbital speed 15 degrees per hour 685# "upper limb" flag 686# and "silent" flag 687# 688# _THEN 689# 690# Compute the sunrise/sunset times for that day 691# 692# _RETURN 693# 694# sunrise and sunset times as hours (GMT Time) 695# season flag: -1 for polar night, +1 for midnight sun, 0 for day and night 696# 697sub _sunrise_sunset { 698 699 my ( $d, $lon, $lat, $altit, $ang_spd, $upper_limb, $silent, $trace, $revsub ) = @_; 700 701 if ($trace) { 702 printf $trace "\n"; 703 } 704 705 # Compute local sidereal time of this moment 706 my $gmst0 = GMST0($d); 707 #my $sidtime = revolution($gmst0 + 180.0 + $lon); 708 my $sidtime = revolution($gmst0 + 180.0); 709 710 # Compute Sun's RA + Decl + distance at this moment 711 my ($sRA, $sdec, $sr) = sun_RA_dec($d, $lon, $trace); 712 713 # Compute time when Sun is at south - in hours (LMT then UTC) 714 my $tsouth_lmt = 12.0 - $revsub->( $sidtime - $sRA ) / 15; 715 my $tsouth = $tsouth_lmt - $lon / 15; 716 717 # Compute the Sun's apparent radius, degrees 718 my $sradius = 0.2666 / $sr; 719 720 if ($trace) { 721 printf $trace "For day $d (%s), GMST0 $gmst0 %s %s\n", _fmt_hr(24 * ($d - int($d)), $lon, 0), _fmt_angle($gmst0 ), _fmt_dur($gmst0 / 15); 722 printf $trace "For day $d (%s), sidereal time $sidtime, %s %s\n", _fmt_hr(24 * ($d - int($d)), $lon, 0), _fmt_angle($sidtime), _fmt_dur($sidtime / 15); 723 printf $trace "For day $d (%s), right asc $sRA %s %s\n", _fmt_hr(24 * ($d - int($d)), $lon, 0), _fmt_angle($sRA ), _fmt_dur($sRA / 15); 724 printf $trace "For day $d (%s), declination $sdec %s %s\n", _fmt_hr(24 * ($d - int($d)), $lon, 0), _fmt_angle($sdec ), _fmt_dur($sdec / 15); 725 printf $trace "For day $d (%s), solar noon at $tsouth (%s)\n", _fmt_hr(24 * ($d - int($d)), $lon, 0), _fmt_hr($tsouth, $lon, 0); 726 } 727 # Do correction to upper limb, if necessary 728 if ($upper_limb) { 729 $altit -= $sradius; 730 } 731 732 # Compute the diurnal arc that the Sun traverses to reach 733 # the specified height altit: 734 735 my $cost = (sind($altit) - sind($lat) * sind($sdec)) 736 / (cosd($lat) * cosd($sdec)); 737 738 if ($trace) { 739 print $trace "altit = $altit, sind(altit) = ", sind($altit), ", lat = $lat, sind(lat) = ", sind($lat), "\n"; 740 print $trace "sdec = $sdec, sind(sdec) = ", sind($sdec), ", lat = $lat, cosd(lat) = ", cosd($lat), "\n"; 741 print $trace "sdec = $sdec, cosd(sdec) = ", cosd($sdec), ", cost = $cost\n"; 742 } 743 744 my $t; 745 my $season = 0; 746 if ( $cost >= 1.0 ) { 747 unless ($silent) { 748 carp "Sun never rises!!\n"; 749 } 750 $t = 0.0; # Sun always below altit 751 $season = -1; 752 } 753 elsif ( $cost <= -1.0 ) { 754 unless ($silent) { 755 carp "Sun never sets!!\n"; 756 } 757 $t = 12.0; # Sun always above altit 758 $season = +1; 759 } 760 else { 761 my $arc = acosd($cost); # The diurnal arc 762 $t = $arc / $ang_spd; # Time to traverse the diurnal arc, hours 763 if ($trace) { 764 printf $trace "Diurnal arc $arc -> $t hours (%s)\n", _fmt_dur($t); 765 } 766 } 767 768 # Store rise and set times - in hours UT 769 770 my $hour_rise_ut = $tsouth - $t; 771 my $hour_set_ut = $tsouth + $t; 772 if ($trace) { 773 printf $trace "For day $d (%s), sunrise at $hour_rise_ut (%s)\n", _fmt_hr(24 * ($d - int($d)), $lon), _fmt_hr($hour_rise_ut, $lon); 774 printf $trace "For day $d (%s), sunset at $hour_set_ut (%s)\n", _fmt_hr(24 * ($d - int($d)), $lon), _fmt_hr($hour_set_ut , $lon); 775 } 776 return ( $hour_rise_ut, $hour_set_ut, $season ); 777 778} 779 780 # 781 # 782 # FUNCTIONAL SEQUENCE for GMST0 783 # 784 # _GIVEN 785 # Day number 786 # 787 # _THEN 788 # 789 # computes GMST0, the Greenwich Mean Sidereal Time 790 # at 0h UT (i.e. the sidereal time at the Greenwhich meridian at 791 # 0h UT). GMST is then the sidereal time at Greenwich at any 792 # time of the day. 793 # 794 # 795 # _RETURN 796 # 797 # Sidtime 798 # 799sub GMST0 { 800 my ($d) = @_; 801 my $sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) + ( 0.9856002585 + 4.70935E-5 ) * $d); 802 return $sidtim0; 803} 804 805 # 806 # 807 # FUNCTIONAL SEQUENCE for sunpos 808 # 809 # _GIVEN 810 # day number 811 # 812 # _THEN 813 # 814 # Computes the Sun's ecliptic longitude and distance 815 # at an instant given in d, number of days since 816 # 2000 Jan 0.0. 817 # 818 # 819 # _RETURN 820 # 821 # ecliptic longitude and distance 822 # ie. $True_solar_longitude, $Solar_distance 823 # 824sub sunpos { 825 826 my ($d) = @_; 827 828 # Mean anomaly of the Sun 829 # Mean longitude of perihelion 830 # Note: Sun's mean longitude = M + w 831 # Eccentricity of Earth's orbit 832 # Eccentric anomaly 833 # x, y coordinates in orbit 834 # True anomaly 835 836 # Compute mean elements 837 my $Mean_anomaly_of_sun = revolution( 356.0470 + 0.9856002585 * $d ); 838 my $Mean_longitude_of_perihelion = 282.9404 + 4.70935E-5 * $d; 839 my $Eccentricity_of_Earth_orbit = 0.016709 - 1.151E-9 * $d; 840 841 # Compute true longitude and radius vector 842 my $Eccentric_anomaly = $Mean_anomaly_of_sun 843 + $Eccentricity_of_Earth_orbit * $RADEG 844 * sind($Mean_anomaly_of_sun) 845 * ( 1.0 + $Eccentricity_of_Earth_orbit * cosd($Mean_anomaly_of_sun) ); 846 847 my $x = cosd($Eccentric_anomaly) - $Eccentricity_of_Earth_orbit; 848 849 my $y = sqrt( 1.0 - $Eccentricity_of_Earth_orbit * $Eccentricity_of_Earth_orbit ) 850 * sind($Eccentric_anomaly); 851 852 my $Solar_distance = sqrt( $x * $x + $y * $y ); # Solar distance 853 my $True_anomaly = atan2d( $y, $x ); # True anomaly 854 855 my $True_solar_longitude = 856 $True_anomaly + $Mean_longitude_of_perihelion; # True solar longitude 857 858 if ( $True_solar_longitude >= 360.0 ) { 859 $True_solar_longitude -= 360.0; # Make it 0..360 degrees 860 } 861 862 return ( $Solar_distance, $True_solar_longitude ); 863} 864 865 # 866 # 867 # FUNCTIONAL SEQUENCE for sun_RA_dec 868 # 869 # _GIVEN 870 # day number, $r and $lon (from sunpos) 871 # 872 # _THEN 873 # 874 # compute RA and dec 875 # 876 # 877 # _RETURN 878 # 879 # Sun's Right Ascension (RA), Declination (dec) and distance (r) 880 # 881 # 882sub sun_RA_dec { 883 884 my ($d, $lon_noon, $trace) = @_; 885 886 # Compute Sun's ecliptical coordinates 887 my ( $r, $lon ) = sunpos($d); 888 if ($trace) { 889 #my $datetime = DateTime->new(year => 1999, month => 12, day => 31)->add(days => $d)->ymd; 890 printf $trace "For day $d (%s), solar noon at ecliptic longitude $lon %s\n", _fmt_hr(24 * ($d - int($d)), $lon_noon), _fmt_angle($lon); 891 } 892 893 # Compute ecliptic rectangular coordinates (z=0) 894 my $x = $r * cosd($lon); 895 my $y = $r * sind($lon); 896 897 # Compute obliquity of ecliptic (inclination of Earth's axis) 898 my $obl_ecl = 23.4393 - 3.563E-7 * $d; 899 900 # Convert to equatorial rectangular coordinates - x is unchanged 901 my $z = $y * sind($obl_ecl); 902 $y = $y * cosd($obl_ecl); 903 904 # Convert to spherical coordinates 905 my $RA = atan2d( $y, $x ); 906 my $dec = atan2d( $z, sqrt( $x * $x + $y * $y ) ); 907 908 return ( $RA, $dec, $r ); 909 910} # sun_RA_dec 911 912 # 913 # 914 # FUNCTIONAL SEQUENCE for days_since_2000_Jan_0 915 # 916 # _GIVEN 917 # A Datetime object 918 # 919 # _THEN 920 # 921 # process the DateTime object for number of days 922 # since Jan,1 2000 (counted in days) 923 # Day 0.0 is at Jan 1 2000 0.0 UT 924 # 925 # _RETURN 926 # 927 # day number 928 # 929sub days_since_2000_Jan_0 { 930 my ($dt) = @_; 931 return int($dt->jd - $jd_2000_Jan_0); 932} 933 934sub sind { 935 sin( ( $_[0] ) * $DEGRAD ); 936} 937 938sub cosd { 939 cos( ( $_[0] ) * $DEGRAD ); 940} 941 942sub tand { 943 tan( ( $_[0] ) * $DEGRAD ); 944} 945 946sub atand { 947 ( $RADEG * atan( $_[0] ) ); 948} 949 950sub asind { 951 ( $RADEG * asin( $_[0] ) ); 952} 953 954sub acosd { 955 ( $RADEG * acos( $_[0] ) ); 956} 957 958sub atan2d { 959 ( $RADEG * atan2( $_[0], $_[1] ) ); 960} 961 962 # 963 # 964 # FUNCTIONAL SEQUENCE for revolution 965 # 966 # _GIVEN 967 # any angle in degrees 968 # 969 # _THEN 970 # 971 # reduces any angle to within the first revolution 972 # by subtracting or adding even multiples of 360.0 973 # 974 # 975 # _RETURN 976 # 977 # the value of the input is >= 0.0 and < 360.0 978 # 979sub revolution { 980 981 my $x = $_[0]; 982 return ( $x - 360.0 * floor( $x * $INV360 ) ); 983} 984 985# 986# 987# FUNCTIONAL SEQUENCE for _rev_lon 988# 989# _GIVEN 990# 991# two angles in degrees, the variable angle and the reference angle (longitude) 992# 993# _THEN 994# 995# Reduce input variable angle to within reference-180 .. reference+180 degrees 996# 997# 998# _RETURN 999# 1000# angle that was reduced 1001# 1002sub _rev_lon { 1003 my ($x, $lon) = @_; 1004 return $lon + rev180($x - $lon); 1005} 1006 1007# 1008# 1009# FUNCTIONAL SEQUENCE for rev180 1010# 1011# _GIVEN 1012# 1013# any angle in degrees 1014# 1015# _THEN 1016# 1017# Reduce input to within -180..+180 degrees 1018# 1019# 1020# _RETURN 1021# 1022# angle that was reduced 1023# 1024sub rev180 { 1025 1026 my ($x) = @_; 1027 1028 return ( $x - 360.0 * floor( $x * $INV360 + 0.5 ) ); 1029} 1030 1031 # 1032 # 1033 # FUNCTIONAL SEQUENCE for equal 1034 # 1035 # _GIVEN 1036 # 1037 # Two floating point numbers and Accuracy 1038 # 1039 # _THEN 1040 # 1041 # Use sprintf to format the numbers to Accuracy 1042 # number of decimal places 1043 # 1044 # _RETURN 1045 # 1046 # True if the numbers are equal 1047 # 1048sub equal { 1049 1050 my ( $A, $B, $dp ) = @_; 1051 1052 return sprintf( "%.${dp}g", $A ) eq sprintf( "%.${dp}g", $B ); 1053} 1054 1055 # 1056 # 1057 # FUNCTIONAL SEQUENCE for convert_hour 1058 # 1059 # _GIVEN 1060 # Hour_rise, Hour_set 1061 # hours are in UT 1062 # 1063 # _THEN 1064 # 1065 # split out the hours and minutes 1066 # Oct 20 2003 1067 # will convert hours to seconds and return this 1068 # let DateTime handle the conversion 1069 # 1070 # _RETURN 1071 # 1072 # number of seconds 1073sub convert_hour { 1074 1075 my ( $hour_rise_ut, $hour_set_ut ) = @_; 1076 my $seconds_rise = floor( $hour_rise_ut * 60 * 60 ); 1077 my $seconds_set = floor( $hour_set_ut * 60 * 60 ); 1078 1079 return ( $seconds_rise, $seconds_set ); 1080} 1081# 1082# 1083# FUNCTIONAL SEQUENCE for _convert_1_hour 1084# 1085# _GIVEN 1086# A UT time in hours 1087# 1088# _THEN 1089# 1090# split out the hours and minutes 1091# Oct 20 2003 1092# will convert hours to seconds and return this 1093# let DateTime handle the conversion 1094# 1095# _RETURN 1096# 1097# number of seconds 1098sub _convert_1_hour { 1099 1100 my ( $hour ) = @_; 1101 my $seconds = floor( $hour * 3600 ); 1102 1103 return $seconds; 1104} 1105 1106sub _fmt_hr { 1107 my ($hr, $lon, $is_lmt) = @_; 1108 my ($lmt, $utc); 1109 if ($is_lmt) { 1110 $lmt = $hr; 1111 $utc = $lmt - $lon / 15; 1112 } 1113 else { 1114 $utc = $hr; 1115 $lmt = $utc + $lon / 15; 1116 } 1117 my $hr_h_utc = $utc; my $hr_h_lmt = $lmt; 1118 my $hr_d_utc = $utc / 24; my $hr_d_lmt = $lmt / 24; 1119 my $hr_utc = floor($utc); my $hr_lmt = floor($lmt); 1120 $utc -= $hr_utc; $lmt -= $hr_lmt; 1121 $utc *= 60; $lmt *= 60; 1122 my $mn_utc = floor($utc); my $mn_lmt = floor($lmt); 1123 $utc -= $mn_utc; $lmt -= $mn_lmt; 1124 $utc *= 60; $lmt *= 60; 1125 my $sc_utc = floor($utc); my $sc_lmt = floor($lmt); 1126 return sprintf("UTC: %02d:%02d:%02d %f h %f d, LMT: %02d:%02d:%02d %f h %f d", $hr_utc, $mn_utc, $sc_utc, $hr_h_utc, $hr_d_utc 1127 , $hr_lmt, $mn_lmt, $sc_lmt, $hr_h_lmt, $hr_d_lmt); 1128} 1129 1130# 1131# Formatting a duration in hours, minutes, seconds. 1132# Also used for angles with a 24-hour scale instead of the usual 360-degree scale 1133# (right ascension, sidereal time...) 1134# 1135sub _fmt_dur { 1136 my ($dur) = @_; 1137 my $sign = ''; 1138 if ($dur < 0) { 1139 $sign = '-'; 1140 $dur *= -1; 1141 } 1142 my $hr = floor($dur); 1143 $dur -= $hr; 1144 $dur *= 60; 1145 my $mn = floor($dur); 1146 $dur -= $mn; 1147 $dur *= 60; 1148 my $sc = floor($dur); 1149 return sprintf("%s%02d h %02d mn %02d s", $sign, $hr, $mn, $sc); 1150} 1151 1152sub _fmt_angle { 1153 my ($angle) = @_; 1154 my $sign = ''; 1155 if ($angle < 0) { 1156 $sign = '-'; 1157 $angle *= -1; 1158 } 1159 my $hr = floor($angle); 1160 $angle -= $hr; 1161 $angle *= 60; 1162 my $mn = floor($angle); 1163 $angle -= $mn; 1164 $angle *= 60; 1165 my $sc = floor($angle); 1166 return sprintf(q<%s%02d°%02d'%02d">, $sign, $hr, $mn, $sc); 1167} 1168 11691962; # Hint: sung by RZ, better known as BD 1170 1171=encoding utf8 1172 1173=head1 NAME 1174 1175DateTime::Event::Sunrise - Perl DateTime extension for computing the sunrise/sunset on a given day 1176 1177=head1 SYNOPSIS 1178 1179 use DateTime; 1180 use DateTime::Event::Sunrise; 1181 1182 # generating DateTime objects from a DateTime::Event::Sunrise object 1183 my $sun_Kyiv = DateTime::Event::Sunrise->new(longitude => +30.85, # 30°51'E 1184 latitude => +50.45); # 50°27'N 1185 for (12, 13, 14) { 1186 my $dt_yapc_eu = DateTime->new(year => 2013, 1187 month => 8, 1188 day => $_, 1189 time_zone => 'Europe/Kiev'); 1190 say "In Kyiv (50°27'N, 30°51'E) on ", $dt_yapc_eu->ymd, " sunrise occurs at ", $sun_Kyiv->sunrise_datetime($dt_yapc_eu)->hms, 1191 " and sunset occurs at ", $sun_Kyiv->sunset_datetime ($dt_yapc_eu)->hms; 1192 } 1193 1194 # generating DateTime objects from DateTime::Set objects 1195 my $sunrise_Austin = DateTime::Event::Sunrise->sunrise(longitude => -94.73, # 97°44'W 1196 latitude => +30.3); # 30°18'N 1197 my $sunset_Austin = DateTime::Event::Sunrise->sunset (longitude => -94.73, 1198 latitude => +30.3); 1199 my $dt_yapc_na_rise = DateTime->new(year => 2013, 1200 month => 6, 1201 day => 3, 1202 time_zone => 'America/Chicago'); 1203 my $dt_yapc_na_set = $dt_yapc_na_rise->clone; 1204 say "In Austin (30°18'N, 97°44'W), sunrises and sunsets are"; 1205 for (1..3) { 1206 $dt_yapc_na_rise = $sunrise_Austin->next($dt_yapc_na_rise); 1207 $dt_yapc_na_set = $sunset_Austin ->next($dt_yapc_na_set); 1208 say $dt_yapc_na_rise, ' ', $dt_yapc_na_set; 1209 } 1210 1211 # If you deal with a polar location 1212 my $sun_in_Halley = DateTime::Event::Sunrise->new( 1213 longitude => -26.65, # 26°39'W 1214 latitude => -75.58, # 75°35'S 1215 precise => 1, 1216 ); 1217 my $Alex_arrival = DateTime->new(year => 2006, # approximate date, not necessarily the exact one 1218 month => 1, 1219 day => 15, 1220 time_zone => 'Antarctica/Rothera'); 1221 say $Alex_arrival->strftime("Alex Gough (a Perl programmer) arrived at Halley Base on %Y-%m-%d."); 1222 if ($sun_in_Halley->is_polar_day($Alex_arrival)) { 1223 say "It would be days, maybe weeks, before the sun would set."; 1224 } 1225 elsif ($sun_in_Halley->is_polar_night($Alex_arrival)) { 1226 say "It would be days, maybe weeks, before the sun would rise."; 1227 } 1228 else { 1229 my $sunset = $sun_in_Halley->sunset_datetime($Alex_arrival); 1230 say $sunset->strftime("And he saw his first antarctic sunset at %H:%M:%S."); 1231 } 1232 1233=head1 DESCRIPTION 1234 1235This module will computes the time of sunrise and sunset for a given date 1236and a given location. The computation uses Paul Schlyter's algorithm. 1237 1238Actually, the module creates a DateTime::Event::Sunrise object or a 1239DateTime::Set object, which are used to generate the sunrise or the sunset 1240times for a given location and for any date. 1241 1242=head1 METHODS 1243 1244=head2 new 1245 1246This is the DateTime::Event::Sunrise constructor. It takes keyword 1247parameters, which are: 1248 1249=over 4 1250 1251=item longitude 1252 1253This is the longitude of the location where the sunrises and sunsets are observed. 1254It is given as decimal degrees: no minutes, no seconds, but tenths and hundredths of degrees. 1255Another break with the normal usage is that Eastern longitude are positive, Western longitudes 1256are negative. 1257 1258Default value is 0, that is Greenwich or any location on the eponymous meridian. 1259 1260=item latitude 1261 1262This is the latitude of the location where the sunrises and sunsets are observed. 1263As for the longitude, it is given as decimal degrees. Northern latitudes are positive 1264numbers, Southern latitudes are negative numbers. 1265 1266Default value is 0, that is any location on the equator. 1267 1268=item altitude 1269 1270This is the height of the Sun at sunrise or sunset. In astronomical context, the altitude or 1271height is the angle between the Sun and the local horizon. It is expressed as degrees, usually 1272with a negative number, since the Sun is I<below> the horizon. 1273 1274Default value is -0.833, that is when the sun's upper limb touches the horizon, while 1275taking in account the light refraction. 1276 1277Positive altitude are allowed, in case the location is near a mountain range 1278behind which the sun rises or sets. 1279 1280=item precise 1281 1282Boolean to control which algorithm is used. A false value gives a simple algorithm, but 1283which can lead to inaccurate sunrise times and sunset times. A true value gives 1284a more elaborate algorithm, with a loop to refine the sunrise and sunset times 1285and obtain a better precision. 1286 1287Default value is 0, to choose the simple algorithm. 1288 1289This parameter replaces the C<iteration> deprecated parameter. 1290 1291=item upper_limb 1292 1293Boolean to choose between checking the Sun's upper limb or its center. 1294A true value selects the upper limb, a false value selects the center. 1295 1296This parameter is significant only when the altitude does not already deal with the sun radius. 1297When the altitude takes into account the sun radius, this parameter should be false. 1298 1299Default value is 0, since the upper limb correction is already 1300taken in account with the default -0.833 altitude. 1301 1302=item silent 1303 1304Boolean to control the output of some warning messages. 1305With polar locations and dates near the winter solstice or the summer solstice, 1306it may happen that the sun never rises above the horizon or never sets below. 1307If this parameter is set to false, the module will send warnings for these 1308conditions. If this parameter is set to true, the module will not pollute 1309your F<STDERR> stream. 1310 1311Default value is 0, for backward compatibility. 1312 1313=item trace 1314 1315This parameter should either be a false value or a filehandle opened 1316for output. In the latter case, a few messages are printed to the 1317filehandle, which allows the programmer to see step by step how the 1318sunrise and the sunset are computed. 1319 1320Used for analysis and debugging purposes. You need to read the text 1321F<doc/astronomical-notes.pod> in the sister module L<Astro::Sunrise> 1322to understand what the traced values represent. 1323 1324Default value is C<0>, which does not produce trace messages. 1325 1326=back 1327 1328=head2 sunrise, sunset 1329 1330Although they come from the DateTime::Event::Sunrise module, these methods 1331are C<DateTime::Set> constructors. They use the same parameters as the C<new> 1332constructor, but they give objects from a different class. 1333 1334=head2 sunrise_datetime, sunset_datetime 1335 1336These two methods apply to C<DateTime::Event::Sunrise> objects (that is, created 1337with C<new>, not C<sunrise> or C<sunset>). They receive one parameter in addition 1338to C<$self>, a C<DateTime> object. They return another C<DateTime> object, 1339for the same day, but with the time of the sunrise or sunset, respectively. 1340 1341=head2 sunrise_sunset_span 1342 1343This method applies to C<DateTime::Event::Sunrise> objects. It accepts a 1344C<DateTime> object as the second parameter. It returns a C<DateTime::Span> 1345object, beginning at sunrise and ending at sunset. 1346 1347=head2 is_polar_night, is_polar_day, is_day_and_night 1348 1349These methods apply to C<DateTime::Event::Sunrise> objects. They accept a 1350C<DateTime> object as the second parameter. They return a boolean indicating 1351the following condutions: 1352 1353=over 4 1354 1355=item * is_polar_night is true when the sun stays under the horizon. Or rather 1356under the altitude parameter used when the C<DateTime::Event::Sunrise> object was created. 1357 1358=item * is_polar_day is true when the sun stays above the horizon, 1359resulting in a "Midnight sun". Or rather when it stays above the 1360altitude parameter used when the C<DateTime::Event::Sunrise> object was created. 1361 1362=item * is_day_and_night is true when neither is_polar_day, nor is_polar_night 1363are true. 1364 1365=back 1366 1367=head2 next current previous contains as_list iterator 1368 1369See DateTime::Set. 1370 1371=head1 EXTENDED EXAMPLES 1372 1373 my $dt = DateTime->new( year => 2000, 1374 month => 6, 1375 day => 20, 1376 ); 1377 1378 my $sunrise = DateTime::Event::Sunrise ->sunrise ( 1379 longitude =>'-118', 1380 latitude =>'33', 1381 altitude => '-0.833', 1382 precise => '1' 1383 ); 1384 1385 my $sunset = DateTime::Event::Sunrise ->sunset ( 1386 longitude =>'-118', 1387 latitude =>'33', 1388 altitude => '-0.833', 1389 precise => '1' 1390 ); 1391 1392 my $tmp_rise = $sunrise->next( $dt ); 1393 1394 my $dt2 = DateTime->new( year => 2000, 1395 month => 12, 1396 day => 31, 1397 ); 1398 1399 # iterator 1400 my $dt_span = DateTime::Span->new( start =>$dt, end=>$dt2 ); 1401 my $set = $sunrise->intersection($dt_span); 1402 my $iter = $set->iterator; 1403 while ( my $dt = $iter->next ) { 1404 print ' ',$dt->datetime; 1405 } 1406 1407 # is it day or night? 1408 my $day_set = DateTime::SpanSet->from_sets( 1409 start_set => $sunrise, end_set => $sunset ); 1410 print $day_set->contains( $dt ) ? 'day' : 'night'; 1411 1412 my $dt = DateTime->new( year => 2000, 1413 month => 6, 1414 day => 20, 1415 time_zone => 'America/Los_Angeles', 1416 ); 1417 1418 my $sunrise = DateTime::Event::Sunrise ->new( 1419 longitude =>'-118' , 1420 latitude => '33', 1421 altitude => '-0.833', 1422 precise => '1' 1423 1424 ); 1425 1426 my $tmp = $sunrise->sunrise_sunset_span($dt); 1427 print "Sunrise is:" , $tmp->start->datetime , "\n"; 1428 print "Sunset is:" , $tmp->end->datetime; 1429 1430=head1 NOTES 1431 1432=head2 Longitude Signs 1433 1434Remember, contrary to the usual convention, 1435 1436EASTERN longitudes are POSITIVE, 1437 1438WESTERN longitudes are NEGATIVE. 1439 1440On the other hand, the latitude signs follow the usual convention: 1441 1442Northen latitudes are positive, 1443 1444Southern latitudes are negative. 1445 1446=head2 Sun Height 1447 1448There are a number of sun heights to choose from. The default is 1449-0.833 because this is what most countries use. Feel free to 1450specify it if you need to. Here is the list of values to specify 1451the sun height with: 1452 1453=over 4 1454 1455=item * B<0> degrees 1456 1457Center of Sun's disk touches a mathematical horizon 1458 1459=item * B<-0.25> degrees 1460 1461Sun's upper limb touches a mathematical horizon 1462 1463=item * B<-0.583> degrees 1464 1465Center of Sun's disk touches the horizon; atmospheric refraction accounted for 1466 1467=item * B<-0.833> degrees 1468 1469Sun's supper limb touches the horizon; atmospheric refraction accounted for 1470 1471=item * B<-6> degrees 1472 1473Civil twilight (one can no longer read outside without artificial illumination) 1474 1475=item * B<-12> degrees 1476 1477Nautical twilight (navigation using a sea horizon no longer possible) 1478 1479=item * B<-15> degrees 1480 1481Amateur astronomical twilight (the sky is dark enough for most astronomical observations) 1482 1483=item * B<-18> degrees 1484 1485Astronomical twilight (the sky is completely dark) 1486 1487=back 1488 1489=head2 Notes on the Precise Algorithm 1490 1491The original method only gives an approximate value of the Sun's rise/set times. 1492The error rarely exceeds one or two minutes, but at high latitudes, when the Midnight Sun 1493soon will start or just has ended, the errors may be much larger. If you want higher accuracy, 1494you must then select the precise variant of the algorithm. This feature is new as of version 0.7. Here is 1495what I (module creator) have tried to accomplish with this. 1496 1497 1498=over 4 1499 1500=item a) 1501 1502Compute sunrise or sunset as always, with one exception: to convert LHA from degrees to hours, 1503divide by 15.04107 instead of 15.0 (this accounts for the difference between the solar day 1504and the sidereal day. 1505 1506=item b) 1507 1508Re-do the computation but compute the Sun's RA and Decl, and also GMST0, for the moment 1509of sunrise or sunset last computed. 1510 1511=item c) 1512 1513Iterate b) until the computed sunrise or sunset no longer changes significantly. 1514Usually 2 iterations are enough, in rare cases 3 or 4 iterations may be needed. 1515 1516=back 1517 1518However, I (second module maintainer) have checked with a few external 1519sources, to obtain test data. And actually, using the value 15.0 gives 1520results closer to what Stellarium and the NOAA solar calculator give. 1521So I will use value 15.0, unless I find a bug in the precise algorithm 1522as presently implemented. 1523 1524=head2 Notes on polar locations 1525 1526If the location is beyond either polar circle, and if the date is near 1527either solstice, there can be midnight sun or polar night. In this 1528case, there is neither sunrise nor sunset, and the module C<carp>s 1529that the sun never rises or never sets. Then, it returns the time at 1530which the sun is at its highest or lowest point. 1531 1532When computing twilights instead of sunrises / sunsets, the limit for 1533polar locations extends a little beyond the polar circle. For example, 1534for nautical twilights (12 degrees below the horizon), the limits 1535where midnight sun might happen is 12 degrees southward of the Arctic 1536Circle and 12 degrees northward of the Antarctic Circle, that is, 1537about 54° latitude instead of 66°33′. 1538 1539 1540=head1 DEPENDENCIES 1541 1542This module requires: 1543 1544=over 4 1545 1546=item * 1547 1548DateTime 1549 1550=item * 1551 1552DateTime::Set 1553 1554=item * 1555 1556DateTime::Span 1557 1558=item * 1559 1560Params::Validate 1561 1562=item * 1563 1564Set::Infinite 1565 1566=item * 1567 1568POSIX 1569 1570=item * 1571 1572Math::Trig 1573 1574=back 1575 1576=head1 BUGS AND CAVEATS 1577 1578Using a latitude of 90 degrees (North Pole or South Pole) gives curious results. 1579I guess that it is linked with a ambiguous value resulting from a 0/0 computation. 1580 1581Using a longitude of 177 degrees, or any longitude near the 180 meridian, may also give 1582curious results, especially with the precise algorithm. 1583 1584The precise algorithm should be thoroughly analysed, to understand why 1585the value 15.04107 advised by Paul Schlyter does not give the expected 1586results. 1587 1588The precise algorithm is not tested with polar locations. At least, it 1589is tested with a near-polar location, Fairbanks, at the time when the 1590night is at its shortest, that is, in June. 1591 1592=head1 AUTHORS 1593 1594Original author: Ron Hill <rkhill@firstlight.net> 1595 1596Co-maintainer: Jean Forget <JFORGET@cpan.org> 1597 1598=head1 SPECIAL THANKS 1599 1600=over 4 1601 1602=item Robert Creager [Astro-Sunrise@LogicalChaos.org] 1603 1604for providing help with converting Paul's C code to perl. 1605 1606=item Flávio S. Glock [fglock@pucrs.br] 1607 1608for providing the the interface to the DateTime::Set 1609module. 1610 1611=item Eric Jensen 1612 1613for positive and interesting advices about the new version of the 1614module 1615 1616=back 1617 1618=head1 CREDITS 1619 1620=over 4 1621 1622=item Paul Schlyter, Stockholm, Sweden 1623 1624for his excellent web page on the subject. 1625 1626=item Rich Bowen (rbowen@rbowen.com) 1627 1628for suggestions. 1629 1630=item People at L<https://geocoder.opencagedata.com/> 1631 1632for noticing an endless loop condition in L<Astro::Sunrise> and for fixing it. 1633 1634=back 1635 1636=head1 COPYRIGHT and LICENSE 1637 1638=head2 Perl Module 1639 1640This program is distributed under the same terms as Perl 5.16.3: 1641GNU Public License version 1 or later and Perl Artistic License 1642 1643You can find the text of the licenses in the F<LICENSE> file or at 1644L<https://dev.perl.org/licenses/artistic.html> 1645and L<https://www.gnu.org/licenses/gpl-1.0.html>. 1646 1647Here is the summary of GPL: 1648 1649This program is free software; you can redistribute it and/or modify 1650it under the terms of the GNU General Public License as published by 1651the Free Software Foundation; either version 1, or (at your option) 1652any later version. 1653 1654This program is distributed in the hope that it will be useful, 1655but WITHOUT ANY WARRANTY; without even the implied warranty of 1656MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1657GNU General Public License for more details. 1658 1659You should have received a copy of the GNU General Public License 1660along with this program; if not, write to the Free Software Foundation, 1661Inc., <https://www.fsf.org/>. 1662 1663=head2 Original C program 1664 1665Here is the copyright information provided by Paul Schlyter 1666for the original C program: 1667 1668Written as DAYLEN.C, 1989-08-16 1669 1670Modified to SUNRISET.C, 1992-12-01 1671 1672(c) Paul Schlyter, 1989, 1992 1673 1674Released to the public domain by Paul Schlyter, December 1992 1675 1676Permission is hereby granted, free of charge, to any person obtaining a 1677copy of this software and associated documentation files (the "Software"), 1678to deal in the Software without restriction, including without limitation 1679the rights to use, copy, modify, merge, publish, distribute, sublicense, 1680and/or sell copies of the Software, and to permit persons to whom the 1681Software is furnished to do so, subject to the following conditions: 1682 1683The above copyright notice and this permission notice shall be included 1684in all copies or substantial portions of the Software. 1685 1686THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 1687IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 1688FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 1689THE AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, 1690WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT 1691OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 1692THE SOFTWARE. 1693 1694=head1 SEE ALSO 1695 1696perl(1). 1697 1698DateTime Web page at L<http://datetime.perl.org/> 1699 1700L<DateTime::Set> 1701 1702L<DateTime::SpanSet> 1703 1704L<Astro::Sunrise> 1705 1706L<DateTime::Event::Jewish::Sunrise> 1707 1708L<Astro::Coords> 1709 1710L<Astro::PAL> 1711 1712Paul Schlyter's homepage at L<https://stjarnhimlen.se/english.html> 1713 1714The NOAA solar calculator at L<https://www.esrl.noaa.gov/gmd/grad/solcalc/> 1715 1716=cut 1717 1718