1 2# Copyright (c) 1999-2017 Rob Fugina <robf@fugina.com> 3# Distributed under the terms of the GNU Public License, Version 3.0 4 5package Astro::SunTime; 6use vars qw(@ISA @EXPORT $VERSION); 7$VERSION = 0.06; 8@ISA = qw(Exporter); 9@EXPORT = qw(sun_time); 10 11# Results can be checked with: http://aa.usno.navy.mil/data/docs/RS_OneYear.php 12 13# 09/03/00 :: winter Make ParseDate optional. It is overkill and I could not get it to 14# compile in perl2exe. It gave runaway comment errors :( 15# 10/12/00 :: winter Change time_zone check to defined, to allow for time_zone 0 16 17use POSIX; 18 19use strict; 20 21# sun_time takes: 22# type => 'rise' | 'set' 23# latitude 24# longitude 25# time_zone => hours from GMT 26# date => date parsable by Time::ParseDate::parsedate() 27# time => to feed to localtime 28 29sub sun_time 30{ 31 my %params = @_; 32 33 my $type = $params{type} || 'rise'; 34 my $latitude = (defined $params{latitude}) ? $params{latitude} : 38.74274; 35 my $longitude = (defined $params{longitude}) ? $params{longitude} : -90.560143; 36 my $time_zone = (defined $params{time_zone}) ? $params{time_zone} : -6; 37 38 my $time; 39 if ($params{date}) { 40 require Time::ParseDate; 41 $time = Time::ParseDate::parsedate($params{date}); 42 } 43 elsif ($params{time}) { 44 $time = $params{time}; 45 } 46 else { 47 $time = time; 48 } 49 my @suntime = localtime($time); 50 51 my $yday = $suntime[7] + 1; 52 53 my $A = 1.5708; 54 my $B = 3.14159; 55 my $C = 4.71239; 56 my $D = 6.28319; 57 my $E = 0.0174533 * $latitude; 58 my $F = 0.0174533 * $longitude; 59 my $G = 0.261799 * $time_zone; 60 61 # For astronomical twilight, use R = -.309017 62 # For nautical twilight, use R = -.207912 63 # For civil twilight, use R = -.104528 64 # For sunrise or sunset, use R = -.0145439 65 66 my $R = -.0145439; 67 if ($params{twilight}) { 68 if($params{twilight} eq 'astronomical') { 69 $R = -.309017; 70 } 71 elsif($params{twilight} eq 'nautical') { 72 $R = -.207912; 73 } 74 elsif($params{twilight} eq 'civil') { 75 $R = -.104528; 76 } 77 } 78 79 80 my $J = ($type eq 'rise') ? $A : $C; 81 my $K = $yday + (($J - $F) / $D); 82 my $L = ($K * .017202) - .0574039; # Solar Mean Anomoly 83 my $M = $L + .0334405 * sin($L); # Solar True Longitude 84 $M += 4.93289 + (3.49066E-04) * sin(2 * $L); 85 $M = &normalize($M, $D); # Quadrant Determination 86 $M += 4.84814E-06 if ($M / $A) - int($M / $A) == 0; 87 my $P = sin($M) / cos($M); # Solar Right Ascension 88 $P = atan2(.91746 * $P, 1); 89 90 # Quadrant Adjustment 91 if ($M > $C) 92 { 93 $P += $D; 94 } 95 elsif ($M > $A) 96 { 97 $P += $B; 98 } 99 100 my $Q = .39782 * sin($M); # Solar Declination 101 $Q = $Q / sqrt(-$Q * $Q + 1); # This is how the original author wrote it! 102 $Q = atan2($Q, 1); 103 104 my $S = $R - (sin($Q) * sin($E)); 105 $S = $S / (cos($Q) * cos($E)); 106 107 return 'none' if abs($S) > 1; # Null phenomenon 108 109 $S = $S / sqrt(-$S * $S + 1); 110 $S = $A - atan2($S, 1); 111 $S = $D - $S if $type eq 'rise'; 112 113 my $T = $S + $P - 0.0172028 * $K - 1.73364; # Local apparent time 114 my $U = $T - $F; # Universal timer 115 my $V = $U + $G; # Wall clock time 116 $V = &normalize($V, $D); 117 $V = $V * 3.81972; 118 119 my $hour = int($V); 120 my $min = int(($V - $hour) * 60 + 0.5); 121 122 @suntime[2,1,0] = ($hour, $min, 0); 123 124 @suntime = localtime(mktime(@suntime)); # normalize date structure 125 126 return sprintf("%d:%02d", @suntime[2,1]); 127} 128 129sub normalize 130{ 131 my $Z = shift; 132 my $D = shift; 133 134 die "Trying to normalize with zero offset..." if ($D == 0); 135 136 while ($Z < 0) {$Z = $Z + $D} 137 while ($Z >= $D) {$Z = $Z - $D} 138 139 return $Z; 140} 141 142 143 1441; 145 146__END__ 147 148 149