Astronomische Berechnungen für Amateure/ Himmelsmechanik/ Planeten VSOP87/ Himmelsmechanik/ Astronomisches Rechenprogramm

Aus Wikibooks

Auch wenn es im Internet und Fachhandel zahlreiche Programme gibt, mit deren Hilfe die Position von Planeten berechnet werden kann, so besteht mitunter doch der Wunsch ein eigenes derartiges Programm zu schreiben. Leider gibt es wenige gute Anleitungen auf dieser Art. Auf http://neoprogrammics.com/vsop87/ kann man entsprechende Funktionen für diverse Programmiersprachen zur Berechnung von Planetenpositionen mit einer Genauigkeit besser als 1" herunterladen, doch ist eine Eingabe des Zeitpunkts für den die Planetenposition berechnet werden soll, als Julianische Tageszahl sehr unbequem und mit der Ausgabe als Bogenmaßwinkel im ekliptikalen System läßt sich das Ergebnis weder mit anderen Ephemeriden vergleichen, noch in übliche Sternkarten eintragen. Daher sind einige weitere Funktionen sinnvoll, welche die Umrechnung von normalen Datum in die Julianische Tageszahl und vom ekliptikalen System in die üblichen Rektaszensions- und Deklinationswerte durchführen.

Grundsätzlich müssen, wenn man die Position des Planeten von der Erde aus bestimmen möchte, die Position der Erde und des gewünschten Planeten ( im Beispiel Merkur) berechnet werden und für beide Objekte die nötigen Funktionen kopiert werden. Es müssen die Funktionen "Heliocentic XYZ" verwendet werden. Ob das Äquinoktium 2000 oder das aktuelle Äquinoktium verwendet wird, kann dem Programmierer überlassen werden.

Das Beispielprogramm ist in php geschrieben. Es dürfte nicht allzu schwierig sein, es in andere Programmiersprachen wie C++ umzuwandeln.

PHP-Programm[Bearbeiten]

<?PHP

$Jahr = 2017;
$Monat = 1;
$Tageswert = 1;
$Stunde = 0; //Weltzeit
$pi = 16*atan(1.0/5.0) - 4*atan(1.0/239.0);

$JD = JD_Wert ($Jahr, $Monat, $Tageswert, $Stunde);
$t = ($JD - 2451545) / 365250;


//Position der Erde

$X1 = Earth_X0($t) + Earth_X1($t) + Earth_X2($t) + Earth_X3($t) + Earth_X4($t) + Earth_X5($t);
$Y1 = Earth_Y0($t) + Earth_Y1($t) + Earth_Y2($t) + Earth_Y3($t) + Earth_Y4($t) + Earth_Y5($t);
$Z1 = Earth_Z0($t) + Earth_Z1($t) + Earth_Z2($t) + Earth_Z3($t) + Earth_Z4($t) + Earth_Z5($t);

//Position des Planeten, hier Merkur

$X2 = Mercury_X0($t) + Mercury_X1($t) + Mercury_X2($t) + Mercury_X3($t) + Mercury_X4($t) + Mercury_X5($t);
$Y2 = Mercury_Y0($t) + Mercury_Y1($t) + Mercury_Y2($t) + Mercury_Y3($t) + Mercury_Y4($t) + Mercury_Y5($t);
$Z2 = Mercury_Z0($t) + Mercury_Z1($t) + Mercury_Z2($t) + Mercury_Z3($t) + Mercury_Z4($t) + Mercury_Z5($t);


//Bestimmung der ekliptikalen Länge ( im Bogenmass)

if (($X2 - $X1) > 0)
{

 $L = atan(($Y2 -$Y1)/($X2 - $X1));

}
elseif (($X2 - $X1) == 0)
{

 if (($Y2 - $Y1) > 0)
{
$L = $pi/2.0;
}
else
{
$L = (3.0*$pi)/2.0;
}

}
else
{

 $L = $pi + atan(($Y2 -$Y1)/($X2 - $X1));

}


//Bestimmung der ekliptikalen Breite ( im Bogenmass)

$B = atan(($Z2 - $Z1)/sqrt((($Y2 - $Y1)*($Y2 - $Y1)) + (($X2 - $X1)*($X2 - $X1))));

//Bestimmung der Rektaszension ( in Stunden) aus den ekliptikalen Koordinaten ( im Bogenmass)

$R = Alpha ($L,$B,$Jahr + Jahresbruch($Tageswert,$Monat,$Jahr));

//Bestimmung der Deklination ( in Grad) aus den ekliptikalen Koordinaten ( im Bogenmass)

$D = Delta ($L,$B,$Jahr + Jahresbruch($Tageswert,$Monat,$Jahr));


// Ausgabe ( Rektaszension in Stunden und Minuten)

echo "R = ".intval($R)." h".(60*($R - intval($R)))." min D = ".$D." ";

// Ausgabe ( Rektaszension in Stunden, Minuten und Sekunden)

echo "R = ".intval($R)." h".intval(60*($R - intval($R)))." min ". 60*((60*($R - intval($R))) - intval(60*($R - intval($R)))) ." s D = ".$D." ";


function Jahresbruch ($Tageswert, $Monatswert, $Jahr) // Umrechnung des Datums in den Quotienten aus der Zahl des Tages im Jahreslauf und der Länge des Jahres.
{
if($Jahr % 4 == 0)
{

   if ($Jahr % 100 == 0)
{
if ($Jahr % 400 == 0)
{
$Schalttag = 1;
}
else
{
$Schalttag = 0;
}
}
else
{
$Schalttag = 1;
}

}
else
{

   $Schalttag = 0;

}

if ($Monatswert == 1)
{
}
elseif ($Monatswert == 2)
{

 $Tageswert = $Tageswert + 31.0;

}
elseif($Monatswert == 3)
{

 $Tageswert = $Tageswert + $Schalttag + 59.0;

}
elseif($Monatswert == 4)
{

 $Tageswert = $Tageswert + $Schalttag + 90.0;

}
elseif($Monatswert == 5)
{

 $Tageswert = $Tageswert + $Schalttag + 120.0;

}
elseif($Monatswert == 6)
{

 $Tageswert = $Tageswert + $Schalttag + 151.0;

}
elseif($Monatswert == 7)
{

 $Tageswert = $Tageswert + $Schalttag + 181.0;

}
elseif($Monatswert == 8)
{

 $Tageswert = $Tageswert + $Schalttag + 212.0;

}
elseif($Monatswert == 9) {

 $Tageswert = $Tageswert + $Schalttag + 243.0;

}
elseif($Monatswert == 10)
{

 $Tageswert = $Tageswert + $Schalttag + 273.0;

}
elseif($Monatswert == 11)
{

 $Tageswert = $Tageswert + $Schalttag + 304.0;

}
elseif($Monatswert == 12)
{

 $Tageswert = $Tageswert + $Schalttag + 334.0;

}
else
{

  echo "Fehler".$Monatswert;

}
return ($Tageswert/(365.0 + $Schalttag));

}


function Alpha ($l2, $b1, $epoche2) //Bestimmung der Rektaszension ( in Stunden) aus den ekliptikalen Koordinaten ( im Bogenmass)
{

   $pi=16.0*(atan(1.0/5.0))-4.0*(atan(1.0/239.0));
$s2=0.409319755-0.0002270867282*(($epoche2/100.0)-19.0)-0.00000002908882087*(($epoche2/100.0)-19.0)*(($epoche2/100.0)-19.0)+0.00000000872664626*(($epoche2/100.0)-19.0)*(($epoche2/100.0)-19.0)*(($epoche2/100.0)-19.0);
/*Schiefe der Ekliptik zum Zeitpunkt $epoche2*/
   /*Zurückumwandlung der Länge und Breite in Rektaszension und Deklination*/
   $delta2=asin(((cos($b1))*(sin($l2))*(sin($s2)))+((cos($s2))*(sin($b1))));
if (((cos($b1))*(cos($l2)))>0)
{
if ((atan((cos($s2)*tan($l2))-((sin($s2)*tan($b1))/cos($l2))))<0)/*$alpha2 darf nicht negativ werden!*/
{
$alpha2=2.0*$pi+atan((cos($s2)*tan($l2))-((sin($s2)*tan($b1))/cos($l2)));
}
else
{
$alpha2=atan((cos($s2)*tan($l2))-((sin($s2)*tan($b1))/cos($l2)));
}
}
if (((cos($b1))*(cos($l2)))==0)
{
if (((cos($s2))*(cos($b1))*(sin($l2)))-((sin($s2))*(sin($b1)))>0)
{
$alpha2=((3.0*$pi)/2.0);
}
else
{
$alpha2=$pi/2.0;
}
}
if (((cos($b1))*(cos($l2)))<0)
{
$alpha2=$pi+atan((cos($s2)*tan($l2))-((sin($s2)*tan($b1))/cos($l2)));
}
/*Umrechnung in Gradmass*/
$alpha2 = (($alpha2/$pi)*12);
$delta2 = (($delta2/$pi)*180);
return $alpha2;

}

function Delta ($l2, $b1, $epoche2) //Bestimmung der Deklination ( in Grad) aus den ekliptikalen Koordinaten ( im Bogenmass)
{

   $pi=16.0*(atan(1.0/5.0))-4.0*(atan(1.0/239.0));
$s2=0.409319755-0.0002270867282*(($epoche2/100.0)-19.0)-0.00000002908882087*(($epoche2/100.0)-19.0)*(($epoche2/100.0)-19.0)+0.00000000872664626*(($epoche2/100.0)-19.0)*(($epoche2/100.0)-19.0)*(($epoche2/100.0)-19.0);
/*Schiefe der Ekliptik zum Zeitpunkt $epoche2*/
$delta2=asin(((cos($b1))*(sin($l2))*(sin($s2)))+((cos($s2))*(sin($b1))));
/*Umrechnung in Gradmass*/
$delta2 = (($delta2/$pi)*180);
return $delta2;

}

function JD_Wert ($Jahr, $Monat, $Tageswert, $Stunde) // Bestimmung des Julianischen Tageswertes
{
$Jahr = intval($Jahr);
if($Jahr % 4 == 0)
{

   if ($Jahr % 100 == 0)
{
if ($Jahr % 400 == 0)
{
$Schalttag = 1;
}
else
{
$Schalttag = 0;
}
}
else
{
$Schalttag = 1;
}

}
else
{

   $Schalttag = 0;

}


$Monatswert = intval($Monat);

if ($Monatswert == 1)
{
}
elseif ($Monatswert == 2)
{

 $Tageswert = $Tageswert + 31.0;

}
elseif($Monatswert == 3)
{

 $Tageswert = $Tageswert + $Schalttag + 59.0;

}
elseif($Monatswert == 4)
{

 $Tageswert = $Tageswert + $Schalttag + 90.0;

}
elseif($Monatswert == 5)
{

 $Tageswert = $Tageswert + $Schalttag + 120.0;

}
elseif($Monatswert == 6)
{

 $Tageswert = $Tageswert + $Schalttag + 151.0;

}
elseif($Monatswert == 7)
{

 $Tageswert = $Tageswert + $Schalttag + 181.0;

}
elseif($Monatswert == 8)
{

 $Tageswert = $Tageswert + $Schalttag + 212.0;

}
elseif($Monatswert == 9)
{

 $Tageswert = $Tageswert + $Schalttag + 243.0;

}
elseif($Monatswert == 10)
{

 $Tageswert = $Tageswert + $Schalttag + 273.0;

}
elseif($Monatswert == 11)
{

 $Tageswert = $Tageswert + $Schalttag + 304.0;

}
elseif($Monatswert == 12)
{

 $Tageswert = $Tageswert + $Schalttag + 334.0;

}
else
{

  echo "Fehler".$Monatswert;

}


$JTage = ($Jahr - 2000)*365;

if ($Jahr > 2000)
{

   $JTage = $JTage + 1 + intval(($Jahr - 2001)/4);
$JTage = $JTage - intval(($Jahr - 2001)/100);
$JTage = $JTage + intval(($Jahr - 2001)/400);

}
else
{

  $JTage = $JTage + intval(($Jahr - 2000)/4);
$JTage = $JTage - intval(($Jahr - 2000)/100);
$JTage = $JTage + intval(($Jahr - 2000)/400);

}

$JTage = 2451543.5 + $JTage + $Tageswert + ($Stunde/24.0);

return $JTage;
}

// Hier Funktionen aus http://neoprogrammics.com/vsop87/ einsetzen

?>