NavList:

   

Reply

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;
}

   

Reply