Code:
#!/usr/bin/perl -w
#
# Berechnung der Sonnen-Position sowie des Sonnenaufgangs und Sonnenuntergangs
# Kay Zobel, 6/2002
use Math::Trig;
return 1;
sub sun_pos # (Datum, LocalTime, GMT-Offset, Longitude, Latitude)
# Übergabeparameter:
# Datum ............ : TT.MM.JJJJ
# LocalTime ........ : HH:MM:SS
# GMT-Offset ....... : East = +hours | West = -hours
# Longitude ........ : GG.GGG°
# Latitude ......... : GG.GGG°
# Rückgabeparameter:
# Azimuth .......... : GG.GGG°
# Altitude ......... : GG.GGG°
{
my ($Datum, $LocalTime, $GMTOffset, $LON, $LAT) = @_;
my ($day, $month, $year) = split(/\./, $Datum);
my ($hour, $min, $sek) = split(/:/, $LocalTime);
my $LT = $hour + ($min*60+$sek)/3600;
my $UT = $LT - $GMTOffset;
my $lat = °2rad($LAT);
my $d = &days_since_2000_Jan_0($day, $month, $year);
my $w = 282.9404 + 4.70935E-5 * $d; # longitude of perihelion
my $a = 1.000000; # mean distance
my $e = 0.016709 - 1.151E-9 * $d; # eccentricity
my $M = 356.0470 + 0.9856002585 * $d; # mean anomaly
$M = &rev360($M);
my $M_rad = °2rad($M);
my $oblecl = 23.4393 - 3.563E-7 * $d; # obliquity of the ecliptic
my $oblecl_rad = °2rad($oblecl);
my $L = $M + $w; # sun's mean longitude
$L = &rev360($L);
my $E = $M + (180/pi) * $e * sin($M_rad) * (1+$e*cos($M_rad)); # eccentric anomaly
my $E_rad = °2rad($E);
my $x = cos($E_rad) - $e; # sun's rectangular coordinates
my $y = sin($E_rad) * sqrt(1-$e**2); # in the plane of the ecliptic
my $r = sqrt($x**2 + $y**2); # distance
my $v = &rad2deg(atan2($y, $x)); # true anomaly
my $lon = $v + $w; # longitude of the sun
$lon = &rev360($lon);
my $lon_rad = °2rad($lon);
my $x1 = $r * cos($lon_rad);
my $y1 = $r * sin($lon_rad);
my $z1 = 0.0;
my $xequat = $x1;
my $yequat = $y1 * cos($oblecl_rad) + 0.0 * sin($oblecl_rad);
my $zequat = $y1 * sin($oblecl_rad) + 0.0 * cos($oblecl_rad);
my $RA = &rad2deg(atan2($yequat, $xequat))/15;
my $Decl_rad = atan2($zequat, sqrt($xequat**2+$yequat**2));
my $Decl = &rad2deg($Decl_rad);
my $GMST0 = &rev360($L + 180) / 15;
my $SIDTIME = $GMST0 + $UT + $LON/15; # sidereal time
my $HA = ($SIDTIME - $RA) * 15;
my $HA_rad = °2rad($HA);
my $x2 = cos($HA_rad) * cos($Decl_rad);
my $y2 = sin($HA_rad) * cos($Decl_rad);
my $z2 = sin($Decl_rad);
my $xhor = $x2 * sin($lat) - $z2 * cos($lat);
my $yhor = $y2;
my $zhor = $x2 * cos($lat) + $z2 * sin($lat);
my $azimuth = atan2($yhor, $xhor) + pi;
my $AZIMUTH = &rad2deg($azimuth);
my $altitude = atan2($zhor, sqrt($xhor**2 + $yhor**2) );
my $ALTITUDE = &rad2deg($altitude);
return ($AZIMUTH, $ALTITUDE);
}
sub days_since_2000_Jan_0 #(DD.MM.YYYY)
{
my ($day, $month, $year) = @_;
$d = 367*$year - int((7*($year+(int(($month+9)/12))))/4)
+ int(275*$month/9) + $day - 730530;
return($d);
}
# reduziert einen Winkel auf einen Bereich zwischen 0 und 360°
sub rev360
{
my $a = $_[0];
my $b = $a - int($a/360.0)*360.0;
$b +=360 if $b<0;
return($b);
}
sub rev24
{
my $a = $_[0];
my $b = $a - int($a/24)*24;
$b +=24 if $b<0;
return($b);
}
sub sun_rise_set # (Datum, Rise/Set, Zenith, Longitude, Latitude, GMT-Offset, DaylightSavingTime)
# Übergabeparameter:
# Datum ............ : TT.MM.JJJJ
# Rise/Set ......... : Rise=1 | Set=0
# Zenith ........... : official | civil | nautical | astronomical
# Longitude ........ : GG.GGG°
# Latitude ......... : GG.GGG°
# GMT-Offset ....... : East = +hours | West = -hours
# DaylightSaving ... : East = +hours | West = -hours
# Rückgabeparameter:
# LocalTime ........ : HH:MM:SS
{
my ($Datum, $rise, $zenith, $LON, $LAT, $GMTOffset, $DST) = @_;
my ($day, $month, $year) = split(/\./, $Datum);
my $ZENITH;
if ($zenith eq 'civil') { $ZENITH = 96.0; }
elsif ($zenith eq 'nautical') { $ZENITH = 102.0; }
elsif ($zenith eq 'astronomical') { $ZENITH = 108.0; }
else { $ZENITH = 90.83333; }
$zenith = deg2rad($ZENITH);
my $lat = °2rad($LAT);
# calculate the day of the year
my $N1 = int(275 * $month / 9);
my $N2 = int(($month + 9) / 12);
my $N3 = (1 + int(($year - 4 * int($year/4) + 2) / 3));
my $N = $N1 - ($N2*$N3) + $day - 30;
# convert the longitude to hour value and
# calculate an approximate time
my $lngHour = $LON / 15;
my $t;
if ($rise) { $t = $N + (( 6 - $lngHour) / 24); }
else { $t = $N + ((18 - $lngHour) / 24); }
# calculate the sun's mean anomaly
my $M = (0.9856 * $t) - 3.289;
$M = &rev360($M);
$m = °2rad($M);
# calculate the sun's true longitude
my $L = $M + (1.916 * sin($m)) + (0.020 * sin(2*$m)) + 282.634;
$L = &rev360($L);
$l = °2rad($L);
# calculate the sun's right ascension
my $ra = atan(0.91764 * tan($l));
my $RA = rad2deg($ra);
$RA += 360 if $RA<0;
# right ascension value needs to be in the same quadrant as L
my $LQuadrant = int($L/90) * 90;
my $RAQuadrant = int($RA/90) * 90;
$RA = $RA + ($LQuadrant - $RAQuadrant);
# right ascension value needs to be in the same quadrant as L
$RA = $RA / 15;
# calculate the sun's declination
my $sinDec = 0.39782 * sin($l);
my $cosDec = cos(asin($sinDec));
# calculate the sun's local hour angle
my $cosH = (cos($zenith) - ($sinDec * sin($lat))) / ($cosDec * cos($lat));
if ($cosH > 1) { return('sun never rises'); }
if ($cosH < -1) { return('sun never sets'); }
# finish calculation H and convert into hours
my $H;
if ($rise) { $H = 360 - &rad2deg(acos($cosH)); }
else { $H = &rad2deg(acos($cosH)); }
$H = $H / 15;
# calculate local mean time of rising/setting
my $T = $H + $RA - (0.06571 * $t) - 6.622;
# adjust back to UTC
my $UT = $T - $lngHour;
$UT = &rev24($UT);
# convert UT value to local time zone of latitude/longitude
my $localT = $UT + $GMTOffset + $DST;
my $hour = int($localT);
my $min = int(($localT - $hour) * 60);
my $sek = 60 * (($localT - $hour) * 60 - $min);
my $localTime = sprintf("%02d:%02d:%02d", $hour, $min, $sek);
return ($localTime);
}
Foren auf CAD.de (alle Foren)
| 
