![Welcome to the NavList Message Boards.](../sx1.jpg)
NavList:
A Community Devoted to the Preservation and Practice of Celestial Navigation and Other Methods of Traditional Wayfinding
Long-Term Sun Almanac - Perl version
From: Dan Allen
Date: 2000 Dec 09, 11:51 AM
From: Dan Allen
Date: 2000 Dec 09, 11:51 AM
Here is a program that I wrote using Perl (www.perl.org, www.perl.com) that will calculate the position of the sun given a date, time, latitude, longitude, and timezone correction. If you are not a Perl programmer, don't try and understand it, as Perl has a terrible syntax. --- # s2k.pl - Sun 2000 from Meeus Astronomical Algorithms # s2k is Copyright Daniel K. Allen, 2000. # All rights reserved. # # 6 Dec 2000 - Created by Dan Allen. # # Usage: perl s2k.pl mm.ddyyyy hh:mm:ss # # all arguments and results are in decimal degrees $, = " "; $\ = "\n"; $PI = 4*atan2(1,1); $D2R = $PI/180; $R2D = 180/$PI; $lat = 47.4807666667; $lon = 121.797433333; $tz = 8; $t = J2000($ARGV[0],$ARGV[1] ? $ARGV[1] : 12); $l = Mod(280.46645 + 36000.76983*$t + 0.0003032*$t*$t,360); $m = Mod(357.5291 + 35999.0503*$t + -0.0001559*$t*$t + -0.00000048*$t*$t*$t,360); $e = (0.016708617 + -0.000042037*$t + -0.0000001236*$t*$t); $c = (1.9146 + -0.004817*$t + -0.000014*$t*$t) * Sin($m); $c += (0.019993 + -0.000101*$t) * Sin($m*2) + (0.00029) * Sin($m*3); $sunLong = $l + $c; $ecc = (84381.448 + -46.815*$t - 0.00059*$t*$t + 0.001813*$t*$t*$t) / 3600; $ra = Mod(ATan2(Cos($ecc)*Sin($sunLong),Cos($sunLong)),360); $dec = ASin(Sin($ecc)*Sin($sunLong)); $gst = (280.46061837 + 0.000387933*$t*$t + -2.58331180573E-8*$t*$t*$t); $gst += $t*36525*360.985647366; $gst = Mod($gst,360); $lha = $gst - $lon - $ra; $alt = ASin(Sin($lat) * Sin($dec) + Cos($lat) * Cos($dec) * Cos($lha)); $az = 180 + ATan(Sin($lha) / (Cos($lha) * Sin($lat) - Cos($lat)*Tan($dec))); $# = "%10.6f"; print " Lon:",$lon," Lat:",$lat; print " RA:",$ra, " Dec:",$dec; print " AZ:",$az, " Alt:",$alt; sub Floor { my $x = shift; return $x < 0 ? int($x) - 1 : int($x) } sub Mod { my ($x,$y) = @_; return $x - $y * Floor($x/$y) } sub Sin { my $x = shift; return sin($x*$D2R) } sub Cos { my $x = shift; return cos($x*$D2R) } sub Tan { my $x = shift; return Sin($x)/Cos($x) } sub ASin { my $x = shift; return atan2($x,sqrt(1 - $x * $x))*$R2D } sub ATan { my $x = shift; return atan2($x,1)*$R2D } sub ATan2 { my ($y,$x) = @_; return atan2($y,$x)*$R2D } sub J2000 { my ($date,$t) = @_; my ($m,$dy,$y,@a); if (index($date,"/") > 0) { # mm/dd/yy or mm/dd/yyyy - Y2K compatible @a = split /\//,$date; $m = $a[0]; $d = $a[1]; $y = $a[2] < 50 ? 2000 + $a[2] : $a[2] < 100 ? 1900 + $a[2] : $a[2]; @a = (); } else { # mm.ddyyyy - HP calculator compatible $m = int($date); $d = int(($date - $m)*100); $y = int(1000000*($date - int($date) - $d/100)+0.5); } @a = split /\:/,$t; print "Time:",$y,$m,$d,$a[0],$a[1],$a[2]," TZ:",$tz; return (Julian($y,$m,$d,$a[0]+$tz,$a[1],$a[2]) - Julian(2000,1,1,12,0,0))/36525; } sub Julian { my ($year,$month,$day,$hr,$min,$sec) = @_; my ($m,$y,$t,$jd); if ($month <= 2) { $y = $year - 1; $m = $month + 12; } else { $y = $year; $m = $month; } if ($year < 1582) { $t = -2; } elsif ($year == 1582 && ($month < 10 || $month == 10 && $day <= 4)) { $t = -2; } else { $t = int($y/400) - int($y/100); } $jd = int(365.25*$y) + int(30.6001*($m+1)) + $t + 1720996.5 + $day; $jd += $hr/24.0 + $min/(24.0*60) + $sec/(24.0*3600); return $jd; }