1# -*- encoding: utf-8; indent-tabs-mode: nil -*- 2# 3# Perl extension for computing the sunrise/sunset on a given day 4# Copyright (C) 1999-2003, 2013, 2015, 2017, 2019, 2021 Ron Hill and Jean Forget 5# 6# See the license in the embedded documentation below. 7# 8package Astro::Sunrise; 9 10use strict; 11use warnings; 12use POSIX qw(floor); 13use Math::Trig; 14use Carp; 15use vars qw( $VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $RADEG $DEGRAD ); 16 17require Exporter; 18 19@ISA = qw( Exporter ); 20@EXPORT = qw( sunrise sun_rise sun_set ); 21@EXPORT_OK = qw( DEFAULT CIVIL NAUTICAL AMATEUR ASTRONOMICAL sind cosd tand asind acosd atand atan2d equal ); 22%EXPORT_TAGS = ( 23 constants => [ qw/DEFAULT CIVIL NAUTICAL AMATEUR ASTRONOMICAL/ ], 24 trig => [ qw/sind cosd tand asind acosd atand atan2d equal/ ], 25 ); 26 27$VERSION = '0.99'; 28$RADEG = ( 180 / pi ); 29$DEGRAD = ( pi / 180 ); 30my $INV360 = ( 1.0 / 360.0 ); 31 32sub sun_rise { 33 my ($sun_rise, undef) = sun_rise_sun_set(@_); 34 return $sun_rise; 35} 36sub sun_set { 37 my (undef, $sun_set) = sun_rise_sun_set(@_); 38 return $sun_set; 39} 40 41sub sun_rise_sun_set { 42 my %arg; 43 if (ref($_[0]) eq 'HASH') { 44 %arg = %{$_[0]}; 45 } 46 else { 47 @arg{ qw/lon lat alt offset/ } = @_; 48 } 49 50 # This trick aims to fulfill two antagonistic purposes: 51 # -- do not load DateTime if the only function called is "sunrise" 52 # -- load DateTime implicitly if the user calls "sun_rise" or "sun_set". This is to be backward 53 # compatible with 0.92 or earlier, when Astro::Sunrise would load DateTime and thus, allow 54 # the user to remove this line from his script. 55 unless ($INC{DateTime}) { 56 eval "use DateTime"; 57 croak $@ 58 if $@; 59 } 60 61 my ($longitude, $latitude) = @arg{ qw/lon lat/ }; 62 my $alt = defined($arg{alt} ) ? $arg{alt} : -0.833; 63 my $offset = defined($arg{offset} ) ? int($arg{offset}) : 0 ; 64 my $tz = defined($arg{time_zone}) ? $arg{time_zone} : 'local'; 65 $arg{precise} ||= 0; 66 $arg{upper_limb} ||= 0; 67 $arg{polar} ||= 'warn'; 68 $arg{trace} ||= 0; 69 croak "Longitude parameter (keyword: 'lon') is mandatory" 70 unless defined $longitude; 71 croak "Latitude parameter (keyword: 'lat') is mandatory" 72 unless defined $latitude; 73 croak "Wrong value of the 'polar' argument: should be either 'warn' or 'retval'" 74 if $arg{polar} ne 'warn' and $arg{polar} ne 'retval'; 75 76 my $today = DateTime->today(time_zone => $tz); 77 $today->set( hour => 12 ); 78 $today->add( days => $offset ); 79 80 my( $sun_rise, $sun_set ) = sunrise( { year => $today->year, 81 month => $today->mon, 82 day => $today->mday, 83 lon => $longitude, 84 lat => $latitude, 85 tz => ( $today->offset / 3600 ), 86 # 87 # DST is always 0 because DateTime 88 # currently (v 0.16) adds one to the 89 # offset during DST hours 90 isdst => 0, 91 alt => $alt, 92 precise => $arg{precise}, 93 upper_limb => $arg{upper_limb}, 94 polar => $arg{polar}, 95 trace => $arg{trace}, 96 } ); 97 return ($sun_rise, $sun_set); 98} 99 100sub sunrise { 101 my %arg; 102 if (ref($_[0]) eq 'HASH') { 103 %arg = %{$_[0]}; 104 } 105 else { 106 @arg{ qw/year month day lon lat tz isdst alt precise/ } = @_; 107 } 108 my ( $year, $month, $day, $lon, $lat, $TZ, $isdst, $trace) 109 = @arg{ qw/year month day lon lat tz isdst trace/ }; 110 my $altit = defined($arg{alt} ) ? $arg{alt} : -0.833; 111 $arg{precise} ||= 0; 112 $arg{upper_limb} ||= 0; 113 $arg{polar} ||= 'warn'; 114 $trace ||= 0; 115 croak "Year parameter is mandatory" 116 unless defined $year; 117 croak "Month parameter is mandatory" 118 unless defined $month; 119 croak "Day parameter is mandatory" 120 unless defined $day; 121 croak "Longitude parameter (keyword: 'lon') is mandatory" 122 unless defined $lon; 123 croak "Latitude parameter (keyword: 'lat') is mandatory" 124 unless defined $lat; 125 croak "Wrong value of the 'polar' argument: should be either 'warn' or 'retval'" 126 if $arg{polar} ne 'warn' and $arg{polar} ne 'retval'; 127 128 if (! $arg{precise}) { 129 my $revsub = \&rev180; # normalizing angles around 0 degrees 130 if ($trace) { 131 printf $trace "\nBasic computation for %04d-%02d-%02d, lon %.3f, lat %.3f, altitude %.3f, upper limb %d\n" 132 , $year, $month, $day 133 , $lon 134 , $lat 135 , $altit 136 , $arg{upper_limb}; 137 } 138 my $d = days_since_2000_Jan_0( $year, $month, $day ) + 0.5 - $lon / 360.0; 139 my ($h1, $h2) = sun_rise_set($d, $lon, $lat, $altit, 15.0, $arg{upper_limb}, $arg{polar}, $trace, $revsub); 140 if ($h1 eq 'day' or $h1 eq 'night' or $h2 eq 'day' or $h2 eq 'night') { 141 return ($h1, $h2); 142 } 143 return convert_hour($h1, $h2, $TZ, $isdst); 144 } 145 else { 146 my $revsub = sub { _rev_lon($_[0], $lon) }; # normalizing angles around the local longitude 147 my $ang_speed = 15.0; 148 149 # This is the initial start 150 my $d = days_since_2000_Jan_0($year, $month, $day) - $lon / 360.0; 151 152 if ($trace) { 153 printf $trace "\nPrecise sunrise computation for %04d-%02d-%02d, lon %.3f, lat %.3f, altitude %.3f, upper limb %d angular speed %.5f\n" 154 , $year, $month, $day 155 , $lon, 156 , $lat, 157 , $altit 158 , $arg{upper_limb} 159 , $ang_speed; 160 } 161 my $h1_lmt = 12; # LMT decimal hours, noon then the successive values of sunrise 162 my $h1_utc; # UTC decimal hours, noon LMT then the successive values of sunrise 163 for my $counter (1..9) { 164 # 9 is a failsafe precaution against a possibly runaway loop 165 # but hopefully, we will leave the loop long before, with "last" 166 my $h2_utc; 167 ($h2_utc, undef) = sun_rise_set($d + $h1_lmt / 24, $lon, $lat, $altit, $ang_speed, $arg{upper_limb}, $arg{polar}, $trace, $revsub); 168 if ($h2_utc eq 'day' or $h2_utc eq 'night') { 169 $h1_utc = $h2_utc; 170 last; 171 } 172 $h1_utc = $h1_lmt - $lon / 15; 173 if (equal($h1_utc, $h2_utc, 5)) { 174 # equal within 1e-5 hour, a little less than a second 175 $h1_utc = $h2_utc; 176 last; 177 } 178 $h1_utc = $h2_utc; 179 $h1_lmt = $h1_utc + $lon / 15; 180 } 181 182 if ($trace) { 183 printf $trace "\nPrecise sunset computation for %04d-%02d-%02d, lon %.3f, lat %.3f, altitude %.3f, upper limb %d angular speed %.5f\n" 184 , $year, $month, $day 185 , $lon, 186 , $lat, 187 , $altit 188 , $arg{upper_limb} 189 , $ang_speed; 190 } 191 my $h3_lmt = 12; # LMT decimal hours, noon then the successive values of sunset 192 my $h3_utc; # UTC decimal hours, noon LMT then the successive values of sunset 193 for my $counter (1..9) { 194 # 9 is a failsafe precaution against a possibly runaway loop 195 # but hopefully, we will leave the loop long before, with "last" 196 my $h4_utc; 197 (undef, $h4_utc) = sun_rise_set($d + $h3_lmt / 24, $lon, $lat, $altit, $ang_speed, $arg{upper_limb}, $arg{polar}, $trace, $revsub); 198 if ($h4_utc eq 'day' or $h4_utc eq 'night') { 199 $h3_utc = $h4_utc; 200 last; 201 } 202 $h3_utc = $h3_lmt - $lon / 15; 203 if (equal($h3_utc, $h4_utc, 5)) { 204 # equal within 1e-5 hour, a little less than a second 205 $h3_utc = $h4_utc; 206 last; 207 } 208 $h3_utc = $h4_utc; 209 $h3_lmt = $h3_utc + $lon / 15; 210 } 211 212 return convert_hour($h1_utc, $h3_utc, $TZ, $isdst); 213 214 } 215} 216####################################################################################### 217# end sunrise 218################################################################################### 219 220# 221# 222# FUNCTIONAL SEQUENCE for days_since_2000_Jan_0 223# 224# _GIVEN 225# year, month, day 226# 227# _THEN 228# 229# process the year month and day (counted in days) 230# Day 0.0 is at Jan 1 2000 0.0 UT 231# Note that ALL divisions here should be INTEGER divisions 232# 233# _RETURN 234# 235# day number 236# 237sub days_since_2000_Jan_0 { 238 use integer; 239 my ( $year, $month, $day ) = @_; 240 241 my $d = 367 * $year 242 - int( ( 7 * ( $year + ( ($month + 9) / 12 ) ) ) / 4 ) 243 + int( (275 * $month) / 9 ) 244 + $day - 730530; 245 246 return $d; 247 248} 249 250# 251# 252# FUNCTIONAL SEQUENCE for convert_hour 253# 254# _GIVEN 255# Hour_rise, Hour_set, Time zone offset, DST setting 256# hours are in UT 257# 258# _THEN 259# 260# convert to local time 261# 262# 263# _RETURN 264# 265# hour:min rise and set 266# 267 268sub convert_hour { 269 my ($hour_rise_ut, $hour_set_ut, $TZ, $isdst) = @_; 270 return (convert_1_hour($hour_rise_ut, $TZ, $isdst), 271 convert_1_hour($hour_set_ut, $TZ, $isdst)); 272} 273# 274# 275# FUNCTIONAL SEQUENCE for convert_1_hour 276# 277# _GIVEN 278# Hour, Time zone offset, DST setting 279# hours are in UT 280# 281# _THEN 282# 283# convert to local time 284# 285# 286# _RETURN 287# 288# hour:min 289# 290 291sub convert_1_hour { 292 my ($hour_ut, $TZ, $isdst) = @_; 293 294 if ($hour_ut eq 'day' or $hour_ut eq 'night') { 295 return $hour_ut; 296 } 297 298 my $hour_local = $hour_ut + $TZ; 299 if ($isdst) { 300 $hour_local ++; 301 } 302 303 # The hour should be between 0 and 24; 304 if ($hour_local < 0) { 305 $hour_local += 24; 306 } 307 elsif ($hour_local > 24) { 308 $hour_local -= 24; 309 } 310 311 my $hour = int ($hour_local); 312 313 my $min = floor(($hour_local - $hour) * 60 + 0.5); 314 315 if ($min >= 60) { 316 $min -= 60; 317 $hour++; 318 $hour -= 24 if $hour >= 24; 319 } 320 321 return sprintf("%02d:%02d", $hour, $min); 322} 323 324 325sub sun_rise_set { 326 my ($d, $lon, $lat, $altit, $ang_spd, $upper_limb, $polar, $trace, $revsub) = @_; 327 328 if ($trace) { 329 printf $trace "\n"; 330 } 331 332 # Compute local sidereal time of this moment 333 my $gmst0 = GMST0($d); 334 my $sidtime = revolution( $gmst0 + 180.0 ); 335 336 # Compute Sun's RA + Decl + distance at this moment 337 my ( $sRA, $sdec, $sr ) = sun_RA_dec($d, $lon, $trace); 338 339 # Compute time when Sun is at south - in hours (LMT then UTC) 340 my $tsouth_lmt = 12.0 - $revsub->( $sidtime - $sRA ) / 15; 341 my $tsouth = $tsouth_lmt - $lon / 15; 342 343 if ($trace) { 344 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); 345 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); 346 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); 347 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); 348 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); 349 } 350 351 if ($upper_limb) { 352 # Compute the Sun's apparent radius, degrees 353 my $sradius = 0.2666 / $sr; 354 $altit -= $sradius; 355 } 356 357 # Compute the diurnal arc that the Sun traverses to reach 358 # the specified altitude altit: 359 my $cost = ( sind($altit) - sind($lat) * sind($sdec) ) 360 / ( cosd($lat) * cosd($sdec) ); 361 362 if ($trace) { 363 print $trace "altit = $altit, sind(altit) = ", sind($altit), ", lat = $lat, sind(lat) = ", sind($lat), "\n"; 364 print $trace "sdec = $sdec, sind(sdec) = ", sind($sdec), ", lat = $lat, cosd(lat) = ", cosd($lat), "\n"; 365 print $trace "sdec = $sdec, cosd(sdec) = ", cosd($sdec), ", cost = $cost\n"; 366 } 367 368 my $t; 369 if ( $cost >= 1.0 ) { 370 if ($polar eq 'retval') { 371 return ('night', 'night'); 372 } 373 carp "Sun never rises!!\n"; 374 $t = 0.0; # Sun always below altit 375 } 376 elsif ( $cost <= -1.0 ) { 377 if ($polar eq 'retval') { 378 return ('day', 'day'); 379 } 380 carp "Sun never sets!!\n"; 381 $t = 12.0; # Sun always above altit 382 } 383 else { 384 my $arc = acosd($cost); # The diurnal arc 385 $t = $arc / $ang_spd; # Time to traverse the diurnal arc, hours 386 if ($trace) { 387 printf $trace "Diurnal arc $arc -> $t hours (%s)\n", _fmt_dur($t); 388 } 389 } 390 391 # Store rise and set times - in hours UT 392 393 my $hour_rise_ut = $tsouth - $t; 394 my $hour_set_ut = $tsouth + $t; 395 if ($trace) { 396 printf $trace "For day $d (%s), sunrise at $hour_rise_ut (%s)\n", _fmt_hr(24 * ($d - int($d)), $lon), 397 _fmt_hr($hour_rise_ut, $lon); 398 printf $trace "For day $d (%s), sunset at $hour_set_ut (%s)\n", _fmt_hr(24 * ($d - int($d)), $lon), 399 _fmt_hr($hour_set_ut , $lon); 400 } 401 return($hour_rise_ut, $hour_set_ut); 402} 403 404######################################################################################################### 405# 406# 407# FUNCTIONAL SEQUENCE for GMST0 408# 409# _GIVEN 410# Day number 411# 412# _THEN 413# 414# computes GMST0, the Greenwich Mean Sidereal Time 415# at 0h UT (i.e. the sidereal time at the Greenwich meridian at 416# 0h UT). GMST is then the sidereal time at Greenwich at any 417# time of the day.. 418# 419# 420# _RETURN 421# 422# Sidtime 423# 424sub GMST0 { 425 my ($d) = @_; 426 427 my $sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) 428 + ( 0.9856002585 + 4.70935E-5 ) * $d ); 429 return $sidtim0; 430 431} 432 433# 434# 435# FUNCTIONAL SEQUENCE for sun_RA_dec 436# 437# _GIVEN 438# day number, $r and $lon (from sunpos) 439# 440# _THEN 441# 442# compute RA and dec 443# 444# 445# _RETURN 446# 447# Sun's Right Ascension (RA), Declination (dec) and distance (r) 448# 449# 450sub sun_RA_dec { 451 my ($d, $lon_noon, $trace) = @_; 452 453 # Compute Sun's ecliptical coordinates 454 my ( $r, $lon ) = sunpos($d); 455 if ($trace) { 456 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); 457 } 458 459 # Compute ecliptic rectangular coordinates (z=0) 460 my $x = $r * cosd($lon); 461 my $y = $r * sind($lon); 462 463 # Compute obliquity of ecliptic (inclination of Earth's axis) 464 my $obl_ecl = 23.4393 - 3.563E-7 * $d; 465 466 # Convert to equatorial rectangular coordinates - x is unchanged 467 my $z = $y * sind($obl_ecl); 468 $y = $y * cosd($obl_ecl); 469 470 # Convert to spherical coordinates 471 my $RA = atan2d( $y, $x ); 472 my $dec = atan2d( $z, sqrt( $x * $x + $y * $y ) ); 473 474 return ( $RA, $dec, $r ); 475 476} # sun_RA_dec 477 478 479# 480# 481# FUNCTIONAL SEQUENCE for sunpos 482# 483# _GIVEN 484# day number 485# 486# _THEN 487# 488# Computes the Sun's ecliptic longitude and distance 489# at an instant given in d, number of days since 490# 2000 Jan 0.0. 491# 492# 493# _RETURN 494# 495# ecliptic longitude and distance 496# ie. $True_solar_longitude, $Solar_distance 497# 498sub sunpos { 499 my ($d) = @_; 500 501 # Mean anomaly of the Sun 502 # Mean longitude of perihelion 503 # Note: Sun's mean longitude = M + w 504 # Eccentricity of Earth's orbit 505 # Eccentric anomaly 506 # x, y coordinates in orbit 507 # True anomaly 508 509 # Compute mean elements 510 my $Mean_anomaly_of_sun = revolution( 356.0470 + 0.9856002585 * $d ); 511 my $Mean_longitude_of_perihelion = 282.9404 + 4.70935E-5 * $d; 512 my $Eccentricity_of_Earth_orbit = 0.016709 - 1.151E-9 * $d; 513 514 # Compute true longitude and radius vector 515 my $Eccentric_anomaly = $Mean_anomaly_of_sun 516 + $Eccentricity_of_Earth_orbit * $RADEG 517 * sind($Mean_anomaly_of_sun) 518 * ( 1.0 + $Eccentricity_of_Earth_orbit * cosd($Mean_anomaly_of_sun) ); 519 520 my $x = cosd($Eccentric_anomaly) - $Eccentricity_of_Earth_orbit; 521 522 my $y = sqrt( 1.0 - $Eccentricity_of_Earth_orbit * $Eccentricity_of_Earth_orbit ) 523 * sind($Eccentric_anomaly); 524 525 my $Solar_distance = sqrt( $x * $x + $y * $y ); # Solar distance 526 my $True_anomaly = atan2d( $y, $x ); # True anomaly 527 528 my $True_solar_longitude = $True_anomaly + $Mean_longitude_of_perihelion; # True solar longitude 529 530 if ( $True_solar_longitude >= 360.0 ) { 531 $True_solar_longitude -= 360.0; # Make it 0..360 degrees 532 } 533 534 return ( $Solar_distance, $True_solar_longitude ); 535} 536 537 538 539sub sind { 540 sin( ( $_[0] ) * $DEGRAD ); 541} 542 543sub cosd { 544 cos( ( $_[0] ) * $DEGRAD ); 545} 546 547sub tand { 548 tan( ( $_[0] ) * $DEGRAD ); 549} 550 551sub atand { 552 ( $RADEG * atan( $_[0] ) ); 553} 554 555sub asind { 556 ( $RADEG * asin( $_[0] ) ); 557} 558 559sub acosd { 560 ( $RADEG * acos( $_[0] ) ); 561} 562 563sub atan2d { 564 ( $RADEG * atan2( $_[0], $_[1] ) ); 565} 566 567# 568# 569# FUNCTIONAL SEQUENCE for revolution 570# 571# _GIVEN 572# any angle 573# 574# _THEN 575# 576# reduces any angle to within the first revolution 577# by subtracting or adding even multiples of 360.0 578# 579# 580# _RETURN 581# 582# the value of the input is >= 0.0 and < 360.0 583# 584 585sub revolution { 586 my $x = $_[0]; 587 return ( $x - 360.0 * floor( $x * $INV360 ) ); 588} 589# 590# 591# FUNCTIONAL SEQUENCE for _rev_lon 592# 593# _GIVEN 594# 595# two angles in degrees, the variable angle and the reference angle (longitude) 596# 597# _THEN 598# 599# Reduce input variable angle to within reference-180 .. reference+180 degrees 600# 601# 602# _RETURN 603# 604# angle that was reduced 605# 606sub _rev_lon { 607 my ($x, $lon) = @_; 608 return $lon + rev180($x - $lon); 609} 610 611# 612# 613# FUNCTIONAL SEQUENCE for rev180 614# 615# _GIVEN 616# 617# any angle 618# 619# _THEN 620# 621# Reduce input to within -180..+180 degrees 622# 623# 624# _RETURN 625# 626# angle that was reduced 627# 628sub rev180 { 629 my ($x) = @_; 630 631 return ( $x - 360.0 * floor( $x * $INV360 + 0.5 ) ); 632} 633 634sub equal { 635 my ($A, $B, $dp) = @_; 636 637 return sprintf("%.${dp}g", $A) eq sprintf("%.${dp}g", $B); 638} 639 640sub _fmt_hr { 641 my ($hr, $lon, $is_lmt) = @_; 642 my ($lmt, $utc); 643 if ($is_lmt) { 644 $lmt = $hr; 645 $utc = $lmt - $lon / 15; 646 } 647 else { 648 $utc = $hr; 649 $lmt = $utc + $lon / 15; 650 } 651 my $hr_h_utc = $utc; my $hr_h_lmt = $lmt; 652 my $hr_d_utc = $utc / 24; my $hr_d_lmt = $lmt / 24; 653 my $hr_utc = floor($utc); my $hr_lmt = floor($lmt); 654 $utc -= $hr_utc; $lmt -= $hr_lmt; 655 $utc *= 60; $lmt *= 60; 656 my $mn_utc = floor($utc); my $mn_lmt = floor($lmt); 657 $utc -= $mn_utc; $lmt -= $mn_lmt; 658 $utc *= 60; $lmt *= 60; 659 my $sc_utc = floor($utc); my $sc_lmt = floor($lmt); 660 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 661 , $hr_lmt, $mn_lmt, $sc_lmt, $hr_h_lmt, $hr_d_lmt); 662} 663 664sub _fmt_dur { 665 my ($dur) = @_; 666 my $sign = ''; 667 if ($dur < 0) { 668 $sign = '-'; 669 $dur *= -1; 670 } 671 my $hr = floor($dur); 672 $dur -= $hr; 673 $dur *= 60; 674 my $mn = floor($dur); 675 $dur -= $mn; 676 $dur *= 60; 677 my $sc = floor($dur); 678 return sprintf("%s%02d h %02d mn %02d s", $sign, $hr, $mn, $sc); 679} 680 681sub _fmt_angle { 682 my ($angle) = @_; 683 my $sign = ''; 684 if ($angle < 0) { 685 $sign = '-'; 686 $angle *= -1; 687 } 688 my $hr = floor($angle); 689 $angle -= $hr; 690 $angle *= 60; 691 my $mn = floor($angle); 692 $angle -= $mn; 693 $angle *= 60; 694 my $sc = floor($angle); 695 return sprintf(q<%s%02d°%02d'%02d">, $sign, $hr, $mn, $sc); 696} 697 698sub DEFAULT () { -0.833 } 699sub CIVIL () { - 6 } 700sub NAUTICAL () { -12 } 701sub AMATEUR () { -15 } 702sub ASTRONOMICAL () { -18 } 703 704# Ending a module with whatever, which risks to be zero, is wrong. 705# Ending a module with 1 is boring. So, let us end it with: 7061950; 707# Hint: directed by BW, with GS, WH and EVS 708 709__END__ 710 711=encoding utf8 712 713=head1 NAME 714 715Astro::Sunrise - Perl extension for computing the sunrise/sunset on a given day 716 717=head1 VERSION 718 719This documentation refers to C<Astro::Sunrise> version 0.99. 720 721=head1 SYNOPSIS 722 723 # When did the sun rise on YAPC::Europe 2015? 724 use Astro::Sunrise; 725 my ($sunrise, $sunset) = sunrise( { year => 2015, month => 9, day => 2, # YAPC::EU starts on 2nd September 2015 726 lon => -3.6, lat => 37.17, # Granada is 37°10'N, 3°36'W 727 tz => 1, isdst => 1 } ); # This is still summer, therefore DST 728 729 # When does the sun rise today in Salt Lake City (home to YAPC::NA 2015)? 730 use Astro::Sunrise; 731 use DateTime; 732 $sunrise_today = sun_rise( { lon => -111.88, lat => 40.75 } ); # 40°45'N, 111°53'W 733 734 # And when does it set tomorrow at Salt Lake City? 735 use Astro::Sunrise; 736 use DateTime; 737 $sunset_tomorrow = sun_set( { lat => 40.75, # 40°45'N, 738 lon => -111.88, # 111°53'W 739 alt => -0.833, # standard value for the sun altitude at sunset 740 offset => 1 } ); # day offset up to tomorrow 741 742=head1 DESCRIPTION 743 744This module will return the sunrise and sunset for a given day. 745 746Months are numbered 1 to 12, in the usual way, not 0 to 11 as in 747C and in Perl's localtime. 748 749 Eastern longitude is entered as a positive number 750 Western longitude is entered as a negative number 751 Northern latitude is entered as a positive number 752 Southern latitude is entered as a negative number 753 754Please note that, when given as positional parameters, the longitude is specified before the 755latitude. 756 757The time zone is given as the numeric value of the offset from UTC. 758 759The C<precise> parameter is set to either 0 or 1. 760If set to 0 no Iteration will occur. 761If set to 1 Iteration will occur, which will give a more precise result. 762Default is 0. 763 764There are a number of sun altitudes to chose from. The default is 765-0.833 because this is what most countries use. Feel free to 766specify it if you need to. Here is the list of values to specify 767altitude (ALT) with, including symbolic constants for each. 768 769=over 770 771=item B<0> degrees 772 773Center of Sun's disk touches a mathematical horizon 774 775=item B<-0.25> degrees 776 777Sun's upper limb touches a mathematical horizon 778 779=item B<-0.583> degrees 780 781Center of Sun's disk touches the horizon; atmospheric refraction accounted for 782 783=item B<-0.833> degrees, DEFAULT 784 785Sun's upper limb touches the horizon; atmospheric refraction accounted for 786 787=item B<-6> degrees, CIVIL 788 789Civil twilight (one can no longer read outside without artificial illumination) 790 791=item B<-12> degrees, NAUTICAL 792 793Nautical twilight (navigation using a sea horizon no longer possible) 794 795=item B<-15> degrees, AMATEUR 796 797Amateur astronomical twilight (the sky is dark enough for most astronomical observations) 798 799=item B<-18> degrees, ASTRONOMICAL 800 801Astronomical twilight (the sky is completely dark) 802 803=back 804 805=head1 USAGE 806 807=head2 B<sunrise> 808 809 ($sunrise, $sunset) = sunrise( { year => $year, month => $month, 810 day => $day, 811 lon => $longitude, lat => $latitude, 812 tz => $tz_offset, isdst => $is_dst, 813 alt => $altitude, upper_limb => $upper_limb, 814 precise => $precise, polar => $action, 815 trace => $fh } ); 816 817 ($sunrise, $sunset) = sunrise(YYYY,MM,DD,longitude,latitude,Time Zone,DST); 818 819 ($sunrise, $sunset) = sunrise(YYYY,MM,DD,longitude,latitude,Time Zone,DST,ALT); 820 821 ($sunrise, $sunset) = sunrise(YYYY,MM,DD,longitude,latitude,Time Zone,DST,ALT,precise); 822 823Returns the sunrise and sunset times, in HH:MM format. 824 825The first form uses a hash reference to pass arguments by name. 826The other forms are kept for backward compatibility. The arguments are: 827 828=over 4 829 830=item year, month, day 831 832The three elements of the date for which you want to compute the sunrise and sunset. 833Months are numbered 1 to 12, in the usual way, not 0 to 11 as in C and in Perl's localtime. 834 835Mandatory, can be positional (#1, #2 and #3). 836 837=item lon, lat 838 839The longitude and latitude of the place for which you want to compute the sunrise and sunset. 840They are given in decimal degrees. For example: 841 842 lon => -3.6, # 3° 36' W 843 lat => 37.17, # 37° 10' N 844 845 Eastern longitude is entered as a positive number 846 Western longitude is entered as a negative number 847 Northern latitude is entered as a positive number 848 Southern latitude is entered as a negative number 849 850Mandatory, can be positional (#4 and #5). 851 852=item tz 853 854Time Zone is the offset from GMT 855 856Mandatory, can be positional (#6). 857 858=item isdst 859 8601 if daylight saving time is in effect, 0 if not. 861 862Mandatory, can be positional (#7). 863 864=item alt 865 866Altitude of the sun, in decimal degrees. Usually a negative number, 867because the sun should be I<under> the mathematical horizon. 868But if there is a high mountain range sunward (that is, southward if you 869live in the Northern hemisphere), you may need to enter a positive altitude. 870 871This parameter is optional. Its default value is -0.833. It can be positional (#8). 872 873=item upper_limb 874 875If this parameter set to a true value (usually 1), the algorithm computes 876the sun apparent radius and takes it into account when computing the sun 877altitude. This parameter is useful only when the C<alt> parameter is set 878to C<0> or C<-0.583> degrees. When using C<-0.25> or C<-0.833> degrees, 879the sun radius is already taken into account. When computing twilights 880(C<-6> to C<-18>), the sun radius is irrelevant. 881 882Since the default value for the C<alt> parameter is -0.833, the 883default value for C<upper_limb> is 0. 884 885This parameter is optional and it can be specified only by keyword. 886 887=item polar 888 889When dealing with a polar location, there may be dates where there is 890a polar night (sun never rises) or a polar day. The default behaviour of 891the module is to emit a warning in these cases ("Sun never rises!!" 892or "Sun never sets!!"). But some programmers may find this inconvenient. 893An alternate behaviour is to return special values reflecting the 894situation. 895 896So, if the C<polar> parameter is set to C<'warn'>, the module emits 897a warning. If the C<polar> parameter is set to C<'retval'>, the 898module emits no warning, but it returns either C<'day'> or C<'night'>. 899 900Example: 901 902 # Loosely based on Alex Gough's activities: scientist and Perl programmer, who spent a year 903 # in Halley Base in 2006. Let us suppose he arrived there on 15th January 2006. 904 my ($sunrise, $sunset) = sunrise( { year => 2006, month => 1, day => 15, 905 lon => -26.65, lat => -75.58, # Halley Base: 75°35'S 26°39'W 906 tz => 3, polar => 'retval' } ); 907 if ($sunrise eq 'day') { 908 say "Alex Gough saw the midnight sun the first day he arrived at Halley Base"; 909 } 910 elsif ($sunrise eq 'night') { 911 say "It would be days, maybe weeks, before the sun would rise."; 912 } 913 else { 914 say "Alex saw his first antarctic sunset at $sunset"; 915 } 916 917This parameter is optional and it can be specified only by keyword. 918 919=item trace 920 921This parameter should either be a false value or a filehandle opened for output. 922In the latter case, a few messages are printed to the filehandle, which allows 923the programmer to see step by step how the sunrise and the sunset are computed. 924 925Used for analysis and debugging purposes. You need to read the text 926F<doc/astronomical-notes.pod> to understand what the traced values 927represent. 928 929This parameter is optional and it can be specified only by keyword. 930 931=item precise 932 933Choice between a precise algorithm and a simpler algorithm. 934The default value is 0, that is, the simpler algorithm. 935Any true value switches to the precise algorithm. 936 937The original method only gives an approximate value of the Sun's rise/set times. 938The error rarely exceeds one or two minutes, but at high latitudes, when the Midnight Sun 939soon will start or just has ended, the errors may be much larger. If you want higher accuracy, 940you must then use the precise algorithm. This feature is new as of version 0.7. Here is 941what I have tried to accomplish with this. 942 943a) Compute sunrise or sunset as always, with one exception: to convert LHA from degrees to hours, 944 divide by 15.04107 instead of 15.0 (this accounts for the difference between the solar day 945 and the sidereal day). 946 947b) Re-do the computation but compute the Sun's RA and Decl, and also GMST0, for the moment 948 of sunrise or sunset last computed. 949 950c) Iterate b) until the computed sunrise or sunset no longer changes significantly. 951 Usually 2 iterations are enough, in rare cases 3 or 4 iterations may be needed. 952 953This parameter is optional. It can be positional (#9). 954 955=back 956 957=head3 I<For Example> 958 959 ($sunrise, $sunset) = sunrise( 2001, 3, 10, 17.384, 98.625, -5, 0 ); 960 ($sunrise, $sunset) = sunrise( 2002, 10, 14, -105.181, 41.324, -7, 1, -18); 961 ($sunrise, $sunset) = sunrise( 2002, 10, 14, -105.181, 41.324, -7, 1, -18, 1); 962 963=head2 B<sun_rise>, B<sun_set> 964 965 $sun_rise = sun_rise( { lon => $longitude, lat => $latitude, 966 alt => $altitude, upper_limb => $bool, 967 offset => $day_offset, 968 precise => $bool_precise, polar => $action } ); 969 $sun_set = sun_set ( { lon => $longitude, lat => $latitude, 970 alt => $altitude, upper_limb => $bool, 971 offset => $day_offset, 972 precise => $bool_precise, polar => $action } ); 973 $sun_rise = sun_rise( $longitude, $latitude ); 974 $sun_rise = sun_rise( $longitude, $latitude, $alt ); 975 $sun_rise = sun_rise( $longitude, $latitude, $alt, $day_offset ); 976 977Returns the sun rise time (resp. the sun set time) for the given location 978and for today's date (as given by DateTime), plus or minus some offset in days. 979The first form use all parameters and transmit them by name. The second form 980uses today's date (from DateTime) and the default altitude. The third 981form adds specifying a custom altitude. The fourth form allows for specifying 982an integer day offset from today, either positive or negative. 983 984The parameters are the same as the parameters for C<sunrise>. There is an additional 985parameter, C<offset>, which allows using a date other than today: C<+1> for 986to-morrow, C<-7> for one week ago, etc. 987 988The arguments are: 989 990=over 4 991 992=item lon, lat 993 994The longitude and latitude of the place for which you want to compute the sunrise and sunset. 995They are given in decimal degrees. For example: 996 997 lon => -3.6, # 3° 36' W 998 lat => 37.17, # 37° 10' N 999 1000 Eastern longitude is entered as a positive number 1001 Western longitude is entered as a negative number 1002 Northern latitude is entered as a positive number 1003 Southern latitude is entered as a negative number 1004 1005Mandatory, can be positional (#1 and #2). 1006 1007=item alt 1008 1009Altitude of the sun, in decimal degrees. Usually a negative number, 1010because the sun should be I<under> the mathematical horizon. 1011But if there is a high mountain range sunward (that is, southward if you 1012live in the Northern hemisphere), you may need to enter a positive altitude. 1013 1014This parameter is optional. Its default value is -0.833. It can be positional (#3). 1015 1016=item offset 1017 1018By default, C<sun_rise> and C<sun_set> use the current day. If you need another 1019day, you give an offset relative to the current day. For example, C<+7> means 1020next week, while C<-365> means last year. 1021 1022This parameter has nothing to do with timezones. 1023 1024Optional, 0 by default, can be positional (#4). 1025 1026=item time_zone 1027 1028Time Zone is the Olson name for a timezone. By default, the functions 1029C<sun_rise> and C<sun_set> will try to use the C<local> timezone. 1030 1031This parameter is optional and it can be specified only by keyword. 1032 1033=item upper_limb 1034 1035If this parameter set to a true value (usually 1), the algorithm computes 1036the sun apparent radius and takes it into account when computing the sun 1037altitude. This parameter is useful only when the C<alt> parameter is set 1038to C<0> or C<-0.583> degrees. When using C<-0.25> or C<-0.833> degrees, 1039the sun radius is already taken into account. When computing twilights 1040(C<-6> to C<-18>), the sun radius is irrelevant. 1041 1042Since the default value for the C<alt> parameter is -0.833, the 1043default value for C<upper_limb> is 0. 1044 1045This parameter is optional and it can be specified only by keyword. 1046 1047=item polar 1048 1049When dealing with a polar location, there may be dates where there is 1050a polar night (sun never rises) or a polar day. The default behaviour of 1051the module is to emit a warning in these cases ("Sun never rises!!" 1052or "Sun never sets!!"). But some programmers may find this inconvenient. 1053An alternate behaviour is to return special values reflecting the 1054situation. 1055 1056So, if the C<polar> parameter is set to C<'warn'>, the module emits 1057a warning. If the C<polar> parameter is set to C<'retval'>, the 1058module emits no warning, but it returns either C<'day'> or C<'night'>. 1059 1060This parameter is optional and it can be specified only by keyword. 1061 1062=item trace 1063 1064This parameter should either be a false value or a filehandle opened for output. 1065In the latter case, a few messages are printed to the filehandle, which allows 1066the programmer to see step by step how the sunrise and the sunset are computed. 1067 1068Used for analysis and debugging purposes. 1069 1070This parameter is optional and it can be specified only by keyword. 1071 1072=item precise 1073 1074Choice between a precise algorithm and a simpler algorithm. 1075The default value is 0, that is, the simpler algorithm. 1076Any true value switches to the precise algorithm. 1077 1078For more documentation, see the corresponding parameter 1079for the C<sunrise> function. 1080 1081This parameter is optional and it can be specified only by keyword. 1082 1083=back 1084 1085=head3 For Example 1086 1087 $sunrise = sun_rise( -105.181, 41.324 ); 1088 $sunrise = sun_rise( -105.181, 41.324, -15 ); 1089 $sunrise = sun_rise( -105.181, 41.324, -12, +3 ); 1090 $sunrise = sun_rise( -105.181, 41.324, undef, -12); 1091 1092=head2 Trigonometric functions using degrees 1093 1094Since the module use trigonometry with degrees, the corresponding functions 1095are available to the module user, free of charge. Just mention the 1096tag C<:trig> in the C<use> statement. These functions are: 1097 1098=over 4 1099 1100=item sind, cosd, tand 1101 1102The direct functions, that is, sine, cosine and tangent functions, respectively. 1103Each one receives one parameter, in degrees, and returns the trigonometric value. 1104 1105=item asind, acosd, atand 1106 1107The reverse functions, that is, arc-sine, arc-cosine, and arc-tangent. 1108Each one receives one parameter, the trigonometric value, and returns the corresponding 1109angle in degrees. 1110 1111=item atan2d 1112 1113Arc-tangent. This function receives two parameters: the numerator and the denominator 1114of a fraction equal to the tangent. Use this function instead of C<atand> when you 1115are not sure the denominator is not zero. E.g.: 1116 1117 use Astro::Sunrise qw(:trig); 1118 say atan2d(1, 2) # prints 26,5 1119 say atan2d(1, 0) # prints 90, without triggering a "division by zero" error 1120 1121=item equal 1122 1123Not really a trigonometrical function, but still useful at some times. This function 1124receives two floating values and an integer value. It compares the floating numbers, 1125and returns "1" if their most significant digits are equal. The integer value 1126specifies how many digits are kept. E.g. 1127 1128 say equal(22/7, 355/113, 3) # prints 1, because : 22/7 = 3.14285715286 rounded to 3.14 1129 # 355/113 = 3.14159292035 rounded to 3.14 1130 say equal(22/7, 355/113, 4) # prints 0, because : 22/7 = 3.14285715286 rounded to 3.143 1131 # 355/113 = 3.14159292035 rounded to 3.142 1132 1133=back 1134 1135=head1 EXPORTS 1136 1137By default, the functions C<sunrise>, C<sun_rise> and C<sun_set> are exported. 1138 1139The constants C<DEFAULT>, C<CIVIL>, C<NAUTICAL>, C<AMATEUR> and C<ASTRONOMICAL> are 1140exported on request with the tag C<:constants>. 1141 1142The functions C<sind>, C<cosd>, C<tand>, C<asind>, C<acosd>, C<atand>, C<atan2d> and C<equal> 1143exported on request with the tag C<:trig>. 1144 1145=head1 DEPENDENCIES 1146 1147This module requires only core modules: L<POSIX>, L<Math::Trig> and L<Carp>. 1148 1149If you use the C<sun_rise> and C<sun_set> functions, you will need also L<DateTime>. 1150 1151=head1 BUGS AND ISSUES 1152 1153Before reporting a bug, please read the text 1154F<doc/astronomical-notes.pod> because the strange behavior you observed 1155may be a correct one, or it may be a corner case already known and 1156already mentioned in the text. 1157 1158Nevertheless, patches and (justified) bug reports are welcome. 1159 1160See L<https://github.com/jforget/Astro-Sunrise/issues> 1161and L<https://github.com/jforget/Astro-Sunrise/pulls>. 1162 1163=head2 Precise Algorithm 1164 1165The explanations for the precise algorithm give a pretty good reason 1166for using an angular speed of 15.04107 instead of 15. Yet, 1167computations with 15.04107 do not give results conforming to what the 1168NOAA website and Stellarium give, while computations with 15 give 1169conforming results. The implementation of the precise algorithm should 1170be analysed and checked to find the reason why 15.04107 does not give 1171the proper results. 1172 1173=head2 Kwalitee 1174 1175The CPANTS tools do not recognize the LICENSE POD paragraph. But any 1176human reader will admit that this LICENSE paragraph exists and is valid. 1177 1178=head2 Haiku-OS CPAN Tester 1179 1180The built-in test F<t/06datetime.t> fails on Haiku-OS because there is 1181no way to extract the timezone name from the system parameters. This 1182failure does not affect the core functions of L<Astro::Sunrise>. 1183 1184Also reported from a user working on a partially configured FreeBSD machine, see 1185L<https://github.com/jforget/Astro-Sunrise/issues/16>. 1186 1187Hopefully, this will be fixed in the current version. 1188 1189=head1 AUTHOR 1190 1191Ron Hill 1192rkhill@firstlight.net 1193 1194Co-maintainer: Jean Forget (JFORGET at cpan dot org) 1195 1196=head1 SPECIAL THANKS 1197 1198Robert Creager [Astro-Sunrise@LogicalChaos.org] 1199for providing help with converting Paul's C code to Perl, 1200for providing code for sun_rise, sun_set subs. 1201Also adding options for different altitudes. 1202 1203Joshua Hoblitt [jhoblitt@ifa.hawaii.edu] 1204for providing the patch to convert to DateTime. 1205 1206Chris Phillips for providing patch for conversion to 1207local time. 1208 1209Brian D Foy for providing patch for constants :) 1210 1211Gabor Szabo for pointing POD mistakes. 1212 1213People at L<https://geocoder.opencagedata.com/> for noticing an endless 1214loop condition and for fixing it. 1215 1216=head1 CREDITS 1217 1218=over 4 1219 1220=item Paul Schlyter, Stockholm, Sweden 1221 1222for his excellent web page on the subject. 1223 1224=item Rich Bowen (rbowen@rbowen.com) 1225 1226for suggestions. 1227 1228=item Adrian Blockley [adrian.blockley@environ.wa.gov.au] 1229 1230for finding a bug in the conversion to local time. 1231 1232=item Slaven Rezić 1233 1234for finding and fixing a bug with DST. 1235 1236=back 1237 1238Lightly verified against L<http://aa.usno.navy.mil/data/docs/RS_OneYear.html> 1239 1240In addition, checked to be compatible with a C implementation of Paul Schlyter's algorithm. 1241 1242=head1 COPYRIGHT and LICENSE 1243 1244=head2 Perl Module 1245 1246Copyright (C) 1999-2003, 2013, 2015, 2017, 2019, 2021 Ron Hill and 1247Jean Forget, all rights reserved. This program is distributed under 1248the same terms as Perl 5.16.3: GNU Public License version 1 or later 1249and Perl Artistic License 1250 1251You can find the text of the licenses in the F<LICENSE> file or at 1252L<https://dev.perl.org/licenses/artistic.html> 1253and L<https://www.gnu.org/licenses/gpl-1.0.html>. 1254 1255Here is the summary of GPL: 1256 1257This program is free software; you can redistribute it and/or modify 1258it under the terms of the GNU General Public License as published by 1259the Free Software Foundation; either version 1, or (at your option) 1260any later version. 1261 1262This program is distributed in the hope that it will be useful, 1263but WITHOUT ANY WARRANTY; without even the implied warranty of 1264MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1265GNU General Public License for more details. 1266 1267You should have received a copy of the GNU General Public License 1268along with this program; if not, write to the Free Software Foundation, 1269Inc., L<https://www.fsf.org/>. 1270 1271=head2 Original C program 1272 1273Here is the copyright information provided by Paul Schlyter: 1274 1275Written as DAYLEN.C, 1989-08-16 1276 1277Modified to SUNRISET.C, 1992-12-01 1278 1279(c) Paul Schlyter, 1989, 1992 1280 1281Released to the public domain by Paul Schlyter, December 1992 1282 1283Permission is hereby granted, free of charge, to any person obtaining a 1284copy of this software and associated documentation files (the "Software"), 1285to deal in the Software without restriction, including without limitation 1286the rights to use, copy, modify, merge, publish, distribute, sublicense, 1287and/or sell copies of the Software, and to permit persons to whom the 1288Software is furnished to do so, subject to the following conditions: 1289 1290The above copyright notice and this permission notice shall be included 1291in all copies or substantial portions of the Software. 1292 1293THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 1294IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 1295FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 1296THE AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, 1297WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT 1298OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 1299THE SOFTWARE. 1300 1301=head1 SEE ALSO 1302 1303perl(1). 1304 1305L<DateTime::Event::Sunrise> 1306 1307L<DateTime::Event::Jewish::Sunrise> 1308 1309The text F<doc/astronomical-notes.pod> (or its original French version 1310F<doc/notes-astronomiques>) in this distribution. 1311 1312L<https://stjarnhimlen.se/comp/riset.html> 1313 1314=cut 1315