Lösung des Quadrantenproblems mit tan(x/2)

[Specht]

Hallo,

Anfang der 1980er Jahre begann ich, Ephemeriden an meinem damals neuen Commodore C64 zu berechnen. Man hatte noch ein graues Telefon mit Wählscheibe und Hörer, es gab noch keine Smartphones und keine Apps. Es war faszinierend, die ersten selbstberechneten Daten mit den Jahrbüchern zu vergleichen.

Der Weg war steinig, es gab eine Menge Hürden, z.B. das Quadrantenproblem. Man stößt darauf beispielsweise bei der Umrechnung von Planetenkoordinaten in verschiedene Koordinatensysteme. Bei der Umrechnung von ekliptikalen Koordinaten in äquatoriale Koordinaten gelten die bekannten Formeln

sin(r)*cos(d) = cos(e)*cos(b)*sin(l) - sin(e)*sin(b) (I)
cos(r)*cos(d) = cos(b)*cos(l) (II)
sin(d) = cos(e)*sin(b) + sin(e)*cos(b)*sin(l) (III)

r: Rektaszension, d: Deklination, l: ekliptikale Länge, b: ekliptikale Breite, e: Schiefe der Ekliptik

Die Größen sind je nach Programmiersprache im Gradmaß oder im Bogenmaß anzugeben. Die Umrechnung eines Winkels vom Grad- ins Bogenmaß erfolgt durch Multiplikation mit pi/180. Im Folgenden gelte das Gradmaß.

Die Problematik des Quadrantenproblems liegt darin, dass durch die Symmetrie der Graphen die Winkelangaben nicht eindeutig sind. Ein Sinuswert von 0.5 ergibt sich beispielsweise für einen Winkel von 30°, aber auch für einen Winkel von 150°. Ferner liegt der Wertebereich der arcsin-Funktion bei Rechnern im Intervall [-90°, 90°], der Wertebereich des Arcuscosinus im Intervall [0°, 180°]. Somit kann man nur die Deklination d in Formel (III) eindeutig über den Arcussinus berechnen, da die Deklination ja stets im Intervall [-90°, 90°] liegt.

Um die Rektaszension r aus (I) und (II) zu errechnen, findet man in der Literatur zwei Wege. Erstens durch umständliche Fallunterscheidungen der rechten Seiten der Gleichungen (I) und (II). Zweitens durch Division der Gleichung (I) durch Gleichung (II). So gelangt man durch die Arcustangens-Funktion zu r, was aber nicht wirklich gut ist, da diese Funktion auch nur Werte im Intervall ]-90°, 90°[ liefert. Manche Programmiersprachen bieten eine zusätzliche spezielle arctan-Funktion mit zwei Argumenten für dieses Problem an.

Einen eigenen, optimalen Weg fand ich in der Halbwinkelfunktion

tan(x/2) = sin(x) / (1+cos(x)) (IV)

Diese Möglichkeit zur Lösung des Quadrantenproblems hatte ich noch nirgendwo sonst gelesen. Man teilt dazu die Gleichungen (I) und (II) jeweils durch cos(d) und setzt dann die erhaltenen Ausdrücke für sin(r) und cos(r) in Gleichung (IV) ein. Man erhält:

tan(r/2) = (cos(e)*cos(b)*sin(l) - sin(e)*sin(b)) / (cos(d) + cos(b)*cos(l))

und somit

r = 2*arctan[(cos(e)*cos(b)*sin(l) - sin(e)*sin(b)) / (cos(d) + cos(b)*cos(l))] (V)

Diese Funktion für die Rektaszension hat nun den Wertebereich im Intervall ]-180°, 180°[ und liefert stets ein eindeutiges Ergebnis für r. Bei einem negativen Wert addiert man einfach 360°. Zuvor muss man noch die Deklination d aus (III) berechnen.

Fast alle Gleichungen für Koordinatentransformationen haben eine Form von (I), (II) und (III). Eine der drei Gleichungen liefert sofort ein eindeutiges Ergebnis, den anderen Wert berechnet man dann über das Einsetzen der jeweiligen Ausdrücke von Sinus und Cosinus des Wertes in die tan(x/2)-Funktion.

Mittlerweile laufen meine kleinen Programme schon fast 40 Jahre mit diesem kleinen Trick. Vom C64 schrieb ich sie um für den späteren Atari ST, dann in Visual Basic unter Windows XP, später Vista und jetzt Windows 10...

 

[Auriga_HH]

Mit dem Problem musste ich mich bei der Excel-Tabelle für meine Einnordungs-Routine auch herumschlagen.

 

[h_c_greier]

Hallo Volker,

Tja das Herumschlagen mit diesen sublimen Dingen sind diese lästigen Basics in der Programmierung.
Gerade in der rechnenden Astronomie kommt das oft vor, wie du schreibst zB bei der Koordinatenumrechnung. Wenn man nach der einschlägigen Literatur astronomische Sachverhalte programmiert, stolpert man immer wieder über solche nervigen Sachen wie Zahlen sortieren, interpolieren oder eben die Haupt- und Nebenwerte trigonometrischer Funktionen, die einen nur aufhalten

Zu deiner Lösung über den halben Tangens: Das kannte ich, habe ich aber nie verwendet. Es gibt eine ähnliche Situation zB bei kleinen Winkeln nahe 0 (sinus) oder nahe 90°/270° (cosinus), wenn diese über arcsin/arccos berechnet werden sollen. Durch die langsame Konvergenz in diesem Bereich benutzt(e) man dort den "Semiversus", im Englischen "haversine", also die Hälfte des Sinus versus.

haversin(x) = (1 - cos(x)) / 2
bzw.
havercos(x) = (1 + cos(x)) / 2

Damit kann man dann ein paar Kommastellen retten. Ich weiß nicht wie sehr das heute noch gebraucht wird.

Bei arctan2, atan2 usw. - je nach Sprache heißt das ja immer irgendwie anders - habe ich mir immer einfach eine eigene Funktion geschrieben. Man muss ja auch den Error 0/0 abfangen. In JavaScript bekommt man den Quadranten-Tangens mitgeliefert, nämlich

Math.atan2(y, x)
y...Zähler
x...Nenner

Beispiel:

const DEGS = Math.PI / 180;

function getArctan2(y, x) {
  return Math.atan2(y, x) / DEGS;
}

console.log(getArctan2(5, 5)); // => 45
console.log(getArctan2(-5, 5)); // => 135
console.log(getArctan2(0, 10)); // => 90
console.log(getArctan2(0, -10)); // => -90

Und dann, wie in der Astronomie oft notwendig, den Winkel noch in das Intervall [0-360°] bringen.

Aber in Zeiten von denen du sprichst gab es noch kein JavaScript und da musste man sich das selbst zusammenbasteln. Von der Performanz her war das wahrscheinlich nicht so toll, denn in der Berechnung tan(r/2) = ... steckt 8(!) mal die Auswertung einer trigonometrischen Funktion. So gesehen würde ich heute zu den in den Sprachen verbauten Routinen greifen und ein Error handling draufsetzen. Aber mit voller Hose ist leicht stinken

In Excel hat man ARCTAN2(x;y), wobei hier die Parameter vertauscht sind! Super, Microsoft!
In JavaScript gibt es Math.atan2(y, x)
In Python ist das glaub ich numpy.arctan2(y, x)


Übrigens ein schöner Post, der wieder mal zeigt, wie gut ein vernünftiger Formeleditor in diesem Forum fehlt!
Es mag ja Menschen geben, die (La)TeX fließen lesen können. Ich kann das nicht, zumindest nicht fließend.
Also liebe dafür zuständige Menschen, wie wär's? Im Parallelforum funktioniert das ohne Probleme.

cs,
harald

--

 

[P_E_T_E_R]

Übrigens ein schöner Post, der wieder mal zeigt, wie gut ein vernünftiger Formeleditor in diesem Forum fehlt!

Das möchte ich mit Nachdruck sekundieren, wobei das dann auch eine brauchbare Sammlung von Sonderzeichen enthalten sollte. Das wurde ja bereits mehrfach angemahnt, aber nie zufriedenstellend erfüllt.

Gruß, Peter

 

[Specht]

Hallo,

Aber in Zeiten von denen du sprichst gab es noch kein JavaScript und da musste man sich das selbst zusammenbasteln. Von der Performanz her war das wahrscheinlich nicht so toll, denn in der Berechnung tan(r/2) = ... steckt 8(!) mal die Auswertung einer trigonometrischen Funktion. So gesehen würde ich heute zu den in den Sprachen verbauten Routinen greifen und ein Error handling draufsetzen.

ja, man musste wirklich oft das Rad immer wieder selbst erfinden. Und die Rechengeschwindigkeit, das dauerte alles wirklich sehr lange . Mein erstes Planetariumprogramm für den C64 umfasste ganze 300(!) Sterne, ein Durchlauf dauerte etwa 4 Minuten, also Zeit genug für einen ...

Es war trotzdem eine abenteuerliche Zeit. Mit zunehmender Rechnerleistung wurde dann wirklich jedes Zeichen und jede Sterngröße auf Millimeterpapier als kleine Matrix selbst entworfen. Ich entwickele mein Hauptprogramm (eine Planetariumsoftware mit Sternen bis 7.0mag) sogar immer noch weiter, benutze es auch noch, zuletzt am 23. Mai abends am Nordwesthimmel, als sich Venus, Mond, Castor und Pollux ein schönes Stelldichein gaben. Die Grafik wirkt ggb modernen Apps wirklich sehr altbacken und rudimentär, funktioniert aber...

 

[h_c_greier]

Das möchte ich mit Nachdruck sekundieren, wobei das dann auch eine brauchbare Sammlung von Sonderzeichen enthalten sollte. Das wurde ja bereits mehrfach angemahnt, aber nie zufriedenstellend erfüllt.

Gruß, Peter

"Ich sehe da Wolken schwarze...."
(Ganz so wie dem momentanen Wetter entsprechend.)

Die Grafik wirkt ggb modernen Apps wirklich sehr altbacken und rudimentär, funktioniert aber...

Wieso, die Grafik sieht eh sehr passabel aus. Hat ein bisschen GUIDE Feeling, und die von CdC ist auch nicht viel anders.
Hut auf!

cs,
harald

--

 

[Ehemaliges Mitglied 77903]

Hallo Sternfreunde,

zur Lösung des Problems "richtiger Quadrant" verwende ich mit großem Nutzen eine Funktion aus den Dunstkreisen Rechnen mit imaginären Zahlen bzw. Umwandlung polare in kartesische Koordinaten, die alle meine Rechenhilfen haben. Ich gebe hier keine genaue Betriebsanleitung, weil die gewünschte Funktion in buchstäbliche jeder Maschine anders heißt, anders untergebracht und anders anzusprechen ist. Stattdessen gebe ich einen Scan meiner internen Notiz, wo ich für mich und alle 6 Rechenhilfen das nötige eindampfte und zusammenstellte. Jeder interessierte Sternfreund muss die fragliche Funktion in seiner eigenen Rechenhilfe erst einmal finden (sie ist garantiert da - sie kann nur verschieden heißen....), sich mit ihrer Funktionsweise und ihren Konventionen befassen.

Meine Rechenhilfen stammen von Texas Instruments (TI), Hewlett Packard (HP), Casio und Libre Office, einem Tabellenkalkulationsprogramm auf dem PC. TI-typisch ist die Benennung R<>P für kartesisch in polar, HP-typisch ist der Stack (umgekehrte polnische Notation). Ersteres muss man wissen, um die Funktion überhaupt zu finden, letzteres hat massive Folgen für die Reihenfolge und Art, wie man die Argumente Z und N der Maschine beizubringen hat. Da es sich um eine zweiwertige Funktion handelt, führt eine verkehrte Reihenfolge der Eingabe von Z und N zu falschen Ergebnissen. Der Anwender MUSS sich diese Dinge für JEDE beliebige Rechenhilfe klarmachen mit Hilfe eines einfachen Zahlenversuchs. In Mathematik-Kompendien suche man unter den Stichworten "kartesisch in polar" und "Darstellungsarten imaginärer Zahlen". Tiefere Kenntnisse des Themas "imaginäre Zahlen" sind nicht notwendig.

Die mathematische Grundidee steht bei mir oben und man sieht sofort die Analogie zu den astronomischen Formeln:

TAN(irgendwas, z.B. Azimut) = Zähler / Nenner.

In meinem Beispiel wählte ich 30 / 18, was zu einem gesuchten Winkel von etwas über 59° führt. +30 und -18 führen zu etwas unter +121°. -30 und -18 gibt -121, -30 und +18 gibt -59. In den astronomischen Formeln hat man es oft mit numerischen Werten von Z und N zwischen -1 und 1 zu tun (cosA mal cosB kann nie größer als absolut 1 sein...), was hier aber ohne Belang ist. BELIEBIGE Werte von Z und N ergeben immer TAN(irgendwas) zwischen +180 über 0 bis -180. Das Problem "richtiger Quadrant" ist vom Acker, man muss nur mehr negative Werte um 360° vermehren, um den in der Astronomie gewohnten Wertevorrat 0 bis 360 zu erreichen.

Bis demnächst, R.M.

 

[astrowolf67]

Hallo,
interessante Lösung für das Quadrantenproblem
Ich habe selbst schon Mitte der 80er begonnen auf dem C64 in Pascal astronomische Berechnungen durchzuführen.
Zumindest in ANSI C gab es schon am Atari ST um 1988 die atan2() Funktion, mit der das Problem gelöst war...

In was habt ihr so programmiert?

CS
Wolfgang

 

[maximilian]

Hallo!

In was habt ihr so programmiert?

Programmierbare Taschenrechner (vor allem HP und Texas Instruments), BASIC, FORTRAN, C, C++, Pascal und ein bisschen Java. Richtig gelernt (an der Universität mit Lochkarten!) habe ich FORTRAN, beruflich verwendet vor allem FORTRAN und C. Von daher kenne ich eigentlich schon immer die Funktion ATAN2, mit der das hier besprochene Problem gar nicht erst auftaucht Diese Funktion habe ich mir in den anderen Programmiersprachen, sofern die nicht von sich aus schon etwas vergleichbares hatten, selbst nachprogrammiert.

Grüße
Maximilian

 

[K-H_49]

Hallo,
ich kenne eine einfache Lösung des Quadrantenproblems bei folgender Gleichungsart:
cos x sin y = A
cos x cos y = B
y erhält man mit tan y = A/B.
Das nicht eindeutige Ergebnis wird mit folgender Regel eindeutig: ist B kleiner Null, dann addiert man zum Ergebnis y einfach 180°.
Ich kenne auch die Halbwinkelberechnung. Aber weil die Berechnung etwas aufwendiger ist, habe ich diesen Weg nicht verwendet.

 

[h_c_greier]

cos x sin y = A
cos x cos y = B

In JS umgesetzt wäre das also sowas

function getArctan2(y, x) {
  let a = Math.cos(x) * Math.sin(y);
  let b = Math.cos(x) * Math.cos(y);
  return b < 0 ? (a / b) + Math.PI : (a / b);
}

JS ES6 hat ja auch eine arctan2 Funktion namens atan2(y,x) eingebaut, die das Ergebnis im richtigen Quadranten liefert.
Die Rückgabe ist natürlich in Radiant und liegt im Intervall [-PI ; PI]

Die Sonderfälle:
x=0, y>0 => Winkel = PI/2
x=0, y<0 => Winkel = -PI/2 (= 3*PI/2)
x=0, y=0 => nicht def., Rückgabe 0

War das auch in FORTRAN schon so, dass beim unbestimmten Ausdruck 0/0 eine 0 zurückgegeben wurde?

cs,
harald

--

 

[Specht]

Hallo,

was mir sehr gut gefällt, ist der Kniff von Karl-Heinz (@K-H_49) mit der Abfrage des Nenners:

ich kenne eine einfache Lösung des Quadrantenproblems bei folgender Gleichungsart:
cos x sin y = A
cos x cos y = B
y erhält man mit tan y = A/B.
Das nicht eindeutige Ergebnis wird mit folgender Regel eindeutig: ist B kleiner Null, dann addiert man zum Ergebnis y einfach 180°.

Astrein ...

Der Kniff mit der tan(x/2)-Funktion ist halt, dass diese Funktion direkt eindeutig ins Intervall ]-180°, 180°[ abbildet und gar keine Abfrage mehr nötig ist. So kann man Planetenkoordinaten ohne Abfragen mit nur wenigen Programmzeilen vom System der Bahnebene zu äquatorialen oder azimutalen Koordinaten transformieren.

 

[Ehemaliges Mitglied 77903]

Hallo Sternfreunde,

den

Kniff von Karl-Heinz

kannte schon der alte Meeus. Ich gebe einen Scan aus seinem Buch von 1991. Bis demnächst R.M.