Autor Nachricht
TomS
BeitragVerfasst am: 30. Nov 2023 14:58    Titel:

sembel hat Folgendes geschrieben:
@TomS, du bist absolut genial, um nicht zu sagen ein Genie, zumindest für mich.

Danke.

Aber das ist so ein bisschen das kleine Einmaleins für Physiker.

sembel hat Folgendes geschrieben:
Mittels dem theta(t) Skript (in JavaScript) hatte ich bereits die t-Werte für volle Gradzahlen zwischen 0 und 360 ermittelt. Mit diesen Werten habe ich die Ergebnisse deines t(theta) Skripts gegengeprüft und sie stimmen überein. (also sind entweder beide falsch oder richtig ;) )

Ich habe ja drei unterschiedliche, voneinander unabhängige Methoden am Start, d.h. es sollte schon stimmen.

sembel hat Folgendes geschrieben:
Ich würde es gerne verstehen, und bis zu einem gewissen Grad gelingt mir das auch, aber dann kommen nur noch Fragezeichen und es verläuft sich in der Ahnungslosigkeit.

Du denkst m.E. zu sehr von der Lösung her.

Der Startpunkt der Überlegung ist eben nicht "Ellipse, dazu noch eine Gleichung", das ganze fängt viel früher an. Und wenn man es verstehen will, muss man da auch durch, sonst findet man sich in einem Wald voller Formeln wieder.

sembel hat Folgendes geschrieben:
Du bist der Erfinder, der Werkzeuge erfindet …

Nee, ich habe nur ein paar Formeln zusammengesucht.

sembel hat Folgendes geschrieben:
Vielen Dank auch für dieses Skript!

Gerne.
sembel
BeitragVerfasst am: 30. Nov 2023 12:31    Titel:

@TomS, du bist absolut genial, um nicht zu sagen ein Genie, zumindest für mich.

Mittels dem theta(t) Skript (in JavaScript) hatte ich bereits die t-Werte für volle Gradzahlen zwischen 0 und 360 ermittelt. Mit diesen Werten habe ich die Ergebnisse deines t(theta) Skripts gegengeprüft und sie stimmen überein. (also sind entweder beide falsch oder richtig Augenzwinkern )

Ich würde es gerne verstehen, und bis zu einem gewissen Grad gelingt mir das auch, aber dann kommen nur noch Fragezeichen und es verläuft sich in der Ahnungslosigkeit. Es wäre wohl auch sehr aufwendig da einzusteigen, zumindest für mich. Von daher lasse ich es gut sein damit. Du bist der Erfinder, der Werkzeuge erfindet, und ich der Handwerker der sie benutzt, bzw. derjenige, der Ideen hat und Werkzeuge braucht um sie zu verwirklichen.

Vielen Dank auch für dieses Skript!
TomS
BeitragVerfasst am: 30. Nov 2023 10:42    Titel:

Meine direkte Lösung anbei.

Getestet mittels Berechnung "1000 Zeitpunkte für einen vollständigen Orbit (T1) => daraus Winkel theta => daraus psi => daraus wieder die Zeiten (T2)"
das liefert einen relativen Fehler von T2 bezogen auf T1 < 2e-7, meist eher 1e-8

Code:
# geometry --------------------

# Verschieben des Winkels in eigener Funktion
def shift_angle(alpha):
    ''' set angle to correct interval '''   
    if np.isscalar(alpha):
        if alpha < 0.0:
            alpha += 2*np.pi
        return alpha
    else:
        alpha = np.where(alpha < 0.0, alpha+2*np.pi, alpha)
        alpha[0] = 0.0
        return alpha
   
def calc_theta(psi, ecc):
    ''' returns the angle as a function of te eccentric anomaly psi; currently restricted to [0, 2*pi] '''
    theta = 2 * np.arctan( np.sqrt((1+ecc)/(1-ecc)) * np.tan(psi/2))
    return shift_angle(theta)

# jetzt psi aus theta
def calc_psi(theta, ecc):
    ''' returns the eccentric anomaly psi as a function of the angle theta; currently restricted to [0, 2*pi] '''
    psi = 2 * np.arctan( np.sqrt((1-ecc)/(1+ecc)) * np.tan(theta/2))
    return shift_angle(psi)

# Zeit aus psi
def calc_time(psi, Om, ecc):
    return (psi - ecc * np.sin(psi)) / Om

def calc_radius(theta, a, ecc):
    return a / (1 + ecc * np.cos(theta))

# Test
theta_arr = np.linspace(0, 2*np.pi, 2000)

# so sieht der Aufruf aus
t_arr = calc_time(calc_psi(theta_arr, the_orb.ecc), the_orb.om_rev, the_orb.ecc)

plt.plot(theta_arr, (t_arr - theta_arr/the_orb.om_rev)/(24*3600))
plt.show()
TomS
BeitragVerfasst am: 30. Nov 2023 06:31    Titel:

Das ist auch nicht der genannte Beitrag.

Bisher ging es um ; jetzt fragst du aber nach .

Wenn ich mich nicht täusche, steht die Lösung in dem genannten Beitrag schon da:



Eine Formel zur Berechnung von psi aus theta u.u. ist ebenfalls genannt.


Weitere Alternativen wären:



Für das Integral findet Wolfram Alpha eine geschlossene Form, der ich aber noch nicht ganz traue.

Alternativ betrachtet man eine Taylorentwicklung des Integranden in e= 0.0167; der nullte Term liefert einfach Eins, alle weiteren Terme sind mit 0.0167, 0.0167^2= 0.00027889 … unterdrückt.

Die dritte Variante ist eine numerische Berechnung des Integrals; da weiß ich nicht, was JavaScript bietet. Das gleiche gilt, wenn man dies als Differentialgleichung auffasst und diese numerisch löst.
sembel
BeitragVerfasst am: 29. Nov 2023 23:50    Titel:

TomS hat Folgendes geschrieben:
Ich denke, die Genauigkeit hängt nicht von der Programmiersprache ab, höchstens die Ausgabe.

Wie "genau" JavaScript ist kann jeder selbst in der Browser Konsole testen:
Code:
1 - 0.9 // 0.09999999999999998

Das hatte ich weiter vorne schon erklärt (binär versus sichtbar).
Code:
(1 - 0.9) + 0.9 // 1


TomS hat Folgendes geschrieben:
Das kannst du selbst Prost

Die notwendigen Formeln stehen im Beitrag vom 21. Nov 2023 13:08, einiges ist schon programmiert.

Danke, dass du mir so viel zutraust, aber eher nein. Ich blicke da nicht durch. Es ist jedoch nicht so, dass ich mir keine Mühe geben würde. Meine Fähigkeiten dahingehend sind stark begrenzt.

Das einzige was ich gefunden habe, bei dem ein Bogenmaß eingegeben wird ist das hier aus #1.
JavaScript
Code:
var ecc = 0.0167086;
var per_rev = 31558432.539;
var om = ((2 * Math.PI) / per_rev);
var the_orb = {"ecc" : ecc, "om" : om};

// returns the derivative of theta w.r.t. time
function dot_theta(theta, Om, ecc) {
   return (Om / (1 - ecc**2)**(3/2)) * (1 + ecc * Math.cos(theta))**2;
}

var theta = (Math.PI / 2);
console.log(dot_theta(theta, the_orb.om, the_orb.ecc));
// 1.9918027937905798e-7

var theta = (Math.PI);
console.log(dot_theta(theta, the_orb.om, the_orb.ecc));
// 1.9257983876238298e-7

var theta = (Math.PI * 2);
console.log(dot_theta(theta, the_orb.om, the_orb.ecc));
// 2.0589193322651474e-7

Das Ergebnis passt nicht.
DrStupid
BeitragVerfasst am: 29. Nov 2023 22:18    Titel:

TomS hat Folgendes geschrieben:
IEEE 754 wird doch heute von fast jeder Hardware unterstützt. Und Programmiersprachen verwenden üblicherweise die Float-Darstellung der Plattform.


IEEE 754 allein sagt noch nicht viel. Das kann von 16 bis 128 Bit gehen. Ob und wenn ja wie das genutzt wird, hängt vom Compiler bzw. Interpreter ab. Mit dem Typ Extended bei Delphi habe ich ja schon ein Beispiel genannt. Alle x86 FPUs unterstützen das, aber aktuelle Compiler nutzen es nicht und machen stattdessen Double daraus. Bei Programmiersprachen mit Typendeklaration sollte es eigentlich allein vom Programmierer abhängen, welcher Fließkommatyp verwendet wird. Aber leider kann man sich selbst darauf nicht verlassen.
TomS
BeitragVerfasst am: 29. Nov 2023 21:40    Titel:

IEEE 754 wird doch heute von fast jeder Hardware unterstützt. Und Programmiersprachen verwenden üblicherweise die Float-Darstellung der Plattform. Das betrifft insbs. die Präzision (im Gegensatz zu Stringdarstellung, inf, nan etc.)

Aber du hast recht, viele Sprachen bieten wohl Konformität zu IEEE 754, garantieren sie jedoch nicht.
DrStupid
BeitragVerfasst am: 29. Nov 2023 20:31    Titel:

TomS hat Folgendes geschrieben:
Ich denke, die Genauigkeit hängt nicht von der Programmiersprache ab, höchstens die Ausgabe.


Doch, das hängt u.a. auch von der Sprache ab. Beispielsweise verwendet JavaScript grundsätzlich immer 64bit Fließkommazahlen gemäß IEEE 754 Standard (auch bekannt als double precision). Da stehen für die signifikanten Ziffern 52 Bit zur Verfügung, was ca. 16 Dezimalziffern entspricht. Andere Fließkommazahlen gibt es in JavaScript nicht.

Es kann aber auch vom Interpreter bzw. Compiler abhängen. Beispielsweise haben ältere Delphi-Versionen den Datentyp Extended noch korrekt als 80bit Fließkommazahl übersetzt. Aktuelle Versionen machen mit demselben Quelltext nur noch ein schnödes Double daraus.
TomS
BeitragVerfasst am: 29. Nov 2023 19:41    Titel:

sembel hat Folgendes geschrieben:
@TomS , die Ergebnisse in JavaScript stimmen mit denen von dir zuletzt geposteten Testdaten und Testergebnissen 1:1 überein. JavaScript hat nur ein paar Nachkommastellen mehr.

Ich denke, die Genauigkeit hängt nicht von der Programmiersprache ab, höchstens die Ausgabe.

sembel hat Folgendes geschrieben:
Ist eine Funktion möglich, die die Anzahl der vergangenen Sekunden seit dem letzten Perihel (in der Funktion die Starttime) als Ergebnis liefert, wenn der Input eine Gradzahl zwischen 0-360 (eventuell auch Bogenmaß 0-2*Pi) ist? Also im Prinzip in die entgegengesetzte Richtung wie bisher. Kann das hier jemand programmieren? Ist es eine sehr komplexe oder zu komplexe Formel?

Das kannst du selbst Prost

Die notwendigen Formeln stehen im Beitrag vom 21. Nov 2023 13:08, einiges ist schon programmiert.
sembel
BeitragVerfasst am: 29. Nov 2023 17:26    Titel:

@TomS , die Ergebnisse in JavaScript stimmen mit denen von dir zuletzt geposteten Testdaten und Testergebnissen 1:1 überein. JavaScript hat nur ein paar Nachkommastellen mehr.

Frage: Am Anfang dieses Themas wurde dies schon einmal erwähnt und jetzt wird explizit danach gefragt: Ist eine Funktion möglich, die die Anzahl der vergangenen Sekunden seit dem letzten Perihel (in der Funktion die Starttime) als Ergebnis liefert, wenn der Input eine Gradzahl zwischen 0-360 (eventuell auch Bogenmaß 0-2*Pi) ist? Also im Prinzip in die entgegengesetzte Richtung wie bisher. Kann das hier jemand programmieren? Ist es eine sehr komplexe oder zu komplexe Formel?
sembel
BeitragVerfasst am: 26. Nov 2023 17:36    Titel:

@TomS , wo hast du den Fehler her? Den habe ich schon gefunden. Dusseligkeitsfehler beim 1:1 umschreiben von Python nach JavaScript deiner letzten Version von contract_Kepler_eq(). Macht sich sofort auf das Ergebnis spürbar. Konsistenzchecks und vieles vieles mehr laufen und kommen noch ... viel Arbeit.

Zum Benchmark
Diese schneidet ein bisschen besser ab als die 1:1 von deinem Python umgeschriebene:
Code:
var prec = 9; // position after decimal point
prec = Math.pow(10, prec);
function contract_Kepler_eq(t) {
  ...
  if (Math.round(psi * prec) == Math.round(psi_check * prec)) {
  ...

Der Unterschied sind wenige Millisekunden (Skriptlaufzeit bei 15 Inputs in Schleife).

Sowie beide besser als:
Code:
  if (psi == psi_check) {

Der Unterschied sind ~10 Millisekunden (Skriptlaufzeit bei 15 Inputs in Schleife).

PS:
- Klammern in Kalkulationen sind hilfreich zum besseren Codelesen.
- Shortcode nur bei Skripten verwenden, die nicht für alle gedacht sind.
- Feste Werte außerhalb von Funktionen / Schleifen vorberechnen.

Edit: Nochmals ein bisschen besser im Benchmark als alles bisherige ist:
Code:
var prec = 9; // position after decimal point
prec = Math.pow(10, (prec * -1));
function contract_Kepler_eq(t) {
  ...
  if (Math.abs(psi - psi_check) < prec) {
  ...


Edit: Zum Rechen-"Fehler" (bei JavaScript)
Dieser tritt immer mit den gleichen 2 Differenz-Werten auf:
Code:
4.440892098500626e-16
0.0000000000000004440892098500626

8.881784197001252e-16
0.0000000000000008881784197001252

Also ist eine Toleranz bis zu 15 Nachkommastellen möglich, um diesen Rechenfehler auszuschließen. Das sollte präzise genug sein.

Die meisten Werte liegen bereits unter 10 Iterationen in der Toleranz. Also kann zusätzlich die Iteration auf 10 gesetzt werden.
TomS
BeitragVerfasst am: 26. Nov 2023 16:53    Titel:

Sorry, jetzt wird's leider ganz falsch:

Code:
var psi_new = (om_t + ecc) * Math.sin(psi);


Es gibt da keine Klammer. So ist's richtig:

Code:
var psi_new = om_t + ecc * Math.sin(psi);


U.a. wg. Tests hatte ich drei völlig verschiedene Lösungsmethoden für das selbe Problem implementiert. Konsistenzchecks sind Pflicht, solltest du auch machen, andernfalls findest du bestimmte Fehler nicht. Und da du keine gute Testdaten im Web findest oder evtl. bei der Übersetzung solcher Daten in Tests nochmal Fehler machst, sind verschiedene Ansätze sinnvoll.
sembel
BeitragVerfasst am: 26. Nov 2023 16:28    Titel:

Don't panic!

Natürlich steckt ein Fehler drin. Hab ich doch gesagt, und zwar in JavaScript. Das hat in den letzten Nachkommastellen (derzeit maximal 17) — wie soll mans sagen — Differenzen zwischen binären Werten und angezeigten Werten. Zum Beispiel kann es sein, dass eine binäre 0.00000000000000099 als 0.00000000000000097 angezeigt wird. Das ist komplex zu erklären und am besten macht man sich darüber keine Gedanken, weil die Kalkulationen passen. Nach der Devise: Dem Koch nicht in den Topf schauen.

JavaScript hat jedoch noch weitere winzige Problemchen. Das betrifft Vergleiche wie if(a == b). Weil ... schwierig zu erklären, also weil jede Variable eine eigene Instanz ist. Was angezeigt wird ist nicht (unbedingt) das was dahinter steckt. Von daher ist das Problemchen mit dem fortlaufenden Unterschied von 4.440892098500626e-16 nach 25 Iterationen stark zu vernachlässigen, wenn es keine größeren Probleme schafft. Denn bedenke: Der Vergleich wurde/wird ohne Rundung ausgeführt. Also in voller Präzission.

Beispiel der Ungenauigkeit von JavaScript:
Code:
0.25 - 0.244
// 0.006000000000000005

Da kann man jetzt einen Aufstand machen! Oder locker bleiben und abrunden (Toleranz) — so wie du es machst.

Edit dazu:
Dieser Input erzeugt solch ein Problem:
Code:
Sekunden: 20889003
psi 48: 4.1448533131870215
psi_check 48: 4.144853313187022
psi 49: 4.144853313187022
psi_check 49: 4.1448533131870215
Differenz: 8.881784197001252e-16
Result: 4.130833583320317 (rad)

Noch einer:
Code:
Sekunden: 15798914
Differenz: 4.440892098500626e-16
Result: 3.145386104382406 (rad)


Viele Wege führen nach Rom. Es gibt sicherlich noch einige Methoden mehr die Aufgabe zu lösen. Da kann man sich streiten was besser ist, bis hin zu Benchmarks.
Code:
// solution via contraction of psi = Omega * t + e * sin(psi)
var num_prec_psi = 1e-8; // corresponds to 0.05 seconds
function Kepler_eq_iter(t, Om, ecc) {
   var om_t = (Om * t);
   var psi = om_t;
   for (var i = 0; i < 5; i++) {
      var psi_new = (om_t + ecc) * Math.sin(psi);
      var psi_dev = Math.max(Math.abs(psi - psi_new));
      if (psi_dev < num_prec_psi) {
         return psi_new;
      }
      psi = psi_new;
   }
   return psi
}

TomS hat Folgendes geschrieben:
Für welche Parameter verwendest du Zufallszahlen?

Für den Input, also in Sekunden.
Code:
Math.floor((Math.random() * Math.floor(per_rev)) + 1)

TomS hat Folgendes geschrieben:
Für ecc ist das heikel;

Aber welcher ecc soll jetzt verwendet werden?
TomS
BeitragVerfasst am: 26. Nov 2023 14:25    Titel:

@sembel –

Dein Code macht nicht das selbe wie meiner.

Auch dass die Schleife nach 25 Iterationen noch nicht abbricht, halte ich für ein Indiz, dass ein Fehler drin steckt.

Für welche Parameter verwendest du Zufallszahlen?

Für ecc ist das heikel; die Methode funktioniert hervorragend für sehr kleine ecc; für größere Werte sind evtl. andere Ansätze besser geeignet.
sembel
BeitragVerfasst am: 26. Nov 2023 12:14    Titel:

@TomS , in JavaScript könnte das so aussehen:
Code:
function Kepler_eq_iter(t, Om, ecc) {

  var om_t = Om * t;
  var psi = om_t;
  var psi_check = 0;

  for (var i = 0; i < 10; i++) {
    psi = om_t + ecc * Math.sin(psi);

    if (psi == psi_check) {
      break;
    }
    psi_check = psi;
  }
  return psi;
}

Es ist eine simple Lösung, die zur Absicherung einen maximalen Wert (10) für den Schleifendurchlauf setzt. Ein paar Test-Inputs haben ergeben, dass schon vor 10 Schleifendurchgängen die Gleichheit erreicht wird. Nachteil dieser Methode ist, dass Werte, die schon beim ersten Durchlauf die Gleichheit erreichen, noch einen zweiten machen müssen. Aber das sollten nur sehr wenige sein (bei 0, ½, 1 Vollkreis).

Edit dazu: Habe es nochmal mit Zufallswerten versucht. Habe dabei einen Ausreißer erwischt, der nach 25 noch nicht gleich war. Die Werte wurden jedoch schon beim 8. Durchlauf erreicht. Das ist vermutlich auf Ursachen auf binärer Ebene zurückzuführen. Die Toleranz sollte vertretbar sein. Weil Zufallszahl, kenne ich den Input nicht.
Code:
psi 9: 2.7726086283833244
psi_check 9: 2.772608628383325

psi 25: 2.7726086283833244
psi_check 25: 2.772608628383325


Es bleibt die Frage nach dem ecc Wert. Laut
https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html

gibt es zwei Angaben:
Code:
Orbital parameters
Orbit eccentricity    0.0167

Earth Mean Orbital Elements (J2000)
Orbital eccentricity 0.01671022

Welcher ist der passendste?
Laut
https://nssdc.gsfc.nasa.gov/planetary/factsheet/fact_notes.html
Code:
Orbital parameters

Instantaneous values (semimajor axis, eccentricity, perihelion, aphelion, inclination)
for planets referenced to Julian Date 2459000.5 (11 June 2020).
[Updated from Julian Date 2451800.5 (13 September 2000) on 22 Dec 2021]

Orbit eccentricity               A measure of the circularity of the orbit, equal to
                                 (aphelion - perihelion distance)/(2 x semi-major axis).
                                 For the Galilean satellites, the forced eccentricity is given. 
                                 For a circular orbit eccentricity = 0. Dimensionless.

Code:
Mean orbital elements

250-year least squares fit elements referenced to J2000 (Global Earth Physics, p. 14)

// ohne "eccentricity"

Dabei ist zu beachten, dass vom gedachten geraden Erdorbit anstatt der Schlangenlinie ausgegangen werden soll.
TomS
BeitragVerfasst am: 26. Nov 2023 08:08    Titel:

Hier mal Test-Daten für die von mir verwendeten Parameterwerte:

Code:
 t_arr =
[       0.           318769.18953076   637538.37906153   956307.56859229
  1275076.75812305  1593845.94765382  1912615.13718458  2231384.32671535
  2550153.51624611  2868922.70577687  3187691.89530764  3506461.0848384
  3825230.27436916  4143999.46389993  4462768.65343069  4781537.84296145
  5100307.03249222  5419076.22202298  5737845.41155375  6056614.60108451
  6375383.79061527  6694152.98014604  7012922.1696768   7331691.35920756
  7650460.54873833  7969229.73826909  8287998.92779985  8606768.11733062
  8925537.30686138  9244306.49639214  9563075.68592291  9881844.87545367
 10200614.06498444 10519383.2545152  10838152.44404596 11156921.63357673
 11475690.82310749 11794460.01263825 12113229.20216902 12431998.39169978
 12750767.58123054 13069536.77076131 13388305.96029207 13707075.14982283
 14025844.3393536  14344613.52888436 14663382.71841513 14982151.90794589
 15300921.09747665 15619690.28700742 15938459.47653818 16257228.66606894
 16575997.85559971 16894767.04513047 17213536.23466123 17532305.424192
 17851074.61372276 18169843.80325352 18488612.99278429 18807382.18231505
 19126151.37184582 19444920.56137658 19763689.75090734 20082458.94043811
 20401228.12996887 20719997.31949963 21038766.5090304  21357535.69856116
 21676304.88809193 21995074.07762269 22313843.26715345 22632612.45668422
 22951381.64621498 23270150.83574574 23588920.02527651 23907689.21480727
 24226458.40433803 24545227.5938688  24863996.78339956 25182765.97293032
 25501535.16246109 25820304.35199185 26139073.54152261 26457842.73105338
 26776611.92058414 27095381.11011491 27414150.29964567 27732919.48917643
 28051688.6787072  28370457.86823796 28689227.05776872 29007996.24729949
 29326765.43683025 29645534.62636102 29964303.81589178 30283073.00542254
 30601842.19495331 30920611.38448407 31239380.57401483 31558149.7635456 ]
delta_theta_arr_3 =
[ 0.          0.00216453  0.00431977  0.0064565   0.00856558  0.010638
  0.01266496  0.01463786  0.01654839  0.01838854  0.02015065  0.02182745
  0.02341207  0.02489809  0.02627956  0.02755101  0.02870751  0.02974463
  0.03065849  0.03144578  0.03210371  0.03263009  0.03302328  0.0332822
  0.03340634  0.03339572  0.03325092  0.03297305  0.03256374  0.03202513
  0.03135983  0.03057094  0.02966201  0.02863702  0.02750035  0.0262568
  0.02491152  0.02347001  0.02193809  0.02032188  0.01862779  0.01686247
  0.0150328   0.01314588  0.01120896  0.00922948  0.00721499  0.00517317
  0.00311175  0.00103856 -0.00103856 -0.00311175 -0.00517317 -0.00721499
 -0.00922948 -0.01120896 -0.01314588 -0.0150328  -0.01686247 -0.01862779
 -0.02032188 -0.02193809 -0.02347001 -0.02491152 -0.0262568  -0.02750035
 -0.02863702 -0.02966201 -0.03057094 -0.03135983 -0.03202513 -0.03256374
 -0.03297305 -0.03325092 -0.03339572 -0.03340634 -0.0332822  -0.03302328
 -0.03263009 -0.03210371 -0.03144578 -0.03065849 -0.02974463 -0.02870751
 -0.02755101 -0.02627956 -0.02489809 -0.02341207 -0.02182745 -0.02015065
 -0.01838854 -0.01654839 -0.01463786 -0.01266496 -0.010638   -0.00856558
 -0.0064565  -0.00431977 -0.00216453  0.        ]
TomS
BeitragVerfasst am: 26. Nov 2023 08:02    Titel:

So in etwa:

Code:
 def Kepler_eq_iter(t, Om, ecc):
    ''' solution via contraction of psi = Omega * t + e * sin(psi)  '''     
    num_prec_psi = 1e-8 # corresponds to 0.05 seconds
    om_t = Om * t
    psi = om_t
    for n in range(0, 5):   
        psi_new = om_t + ecc * np.sin(psi)
        psi_dev = np.max(np.abs(psi - psi_new))
        if psi_dev < num_prec_psi:
            return psi_new                     
        psi = psi_new       
    print('convergence warning: deviation =', psi_dev)       
    return psi
TomS
BeitragVerfasst am: 26. Nov 2023 07:10    Titel:

DrStupid hat Folgendes geschrieben:
TomS hat Folgendes geschrieben:
Alles was ich sagen wollte ist, dass man sich noch anschauen muss, ob 5 immer passt.

Kann man da nicht einfach eine while-Schleife laufen lassen, bis psi im Rahmen der Rechengenauigkeit konstant bleibt?

Klar sollte man das so machen. Ich ändere das.

@sembel –

Es geht um die Lösung der Gleichung aufgrund dieser Eigenschaft: https://de.m.wikipedia.org/wiki/Kontraktion_(Mathematik)
sembel
BeitragVerfasst am: 26. Nov 2023 00:56    Titel:

DrStupid hat Folgendes geschrieben:
Kann man da nicht einfach eine while-Schleife laufen lassen, bis psi im Rahmen der Rechengenauigkeit konstant bleibt?

Kann man. Also wenn es darauf ankommt, dass psi "gleich" (ausgeglichen) ist, also wenn es zweimal gleich ist, dann kann die Schleife unterbrochen werden. Aber heute nicht mehr. Morgen.
DrStupid
BeitragVerfasst am: 26. Nov 2023 00:49    Titel:

TomS hat Folgendes geschrieben:
Alles was ich sagen wollte ist, dass man sich noch anschauen muss, ob 5 immer passt.


Kann man da nicht einfach eine while-Schleife laufen lassen, bis psi im Rahmen der Rechengenauigkeit konstant bleibt?
sembel
BeitragVerfasst am: 26. Nov 2023 00:48    Titel:

TomS hat Folgendes geschrieben:

Dass du
Code:
for n in range(0, 5):

auskommentierst, ist falsch! Ich glaube, du verstehst nicht, was die Schleife bedeutet; sie ist essentiell. Alles was ich sagen wollte ist, dass man sich noch anschauen muss, ob 5 immer passt.

Ich glaube ich habe jetzt verstanden was diese Schleife bewirkt, bzw. dass sie nichts mit den Arrays fürs Plotten zu tun hat.

Der Wert "psi" wird also in sich verschachtelt mehrmals berechnet.
Code:
var om_t = (Om * t);

var psi = om_t;

for (var i = 0; i <= 5; i++) {
   
   psi = (om_t + (ecc * Math.sin(psi)));
}

Wie meinst du aber das "ob 5 immer passt"? Meinst du das so, ob 5 manchmal zu wenig oder manchmal zu viel ist? Mir ist sowas von anderen Berechnungen bekannt. Doch dann müsste gewusst werden, wann die Schleife beendet werden kann, oder wann sie fortgeführt werden soll. Hast du da was parat?
TomS
BeitragVerfasst am: 26. Nov 2023 00:28    Titel:

sembel hat Folgendes geschrieben:
TomS hat Folgendes geschrieben:
In der Funktion zur Iteration sollte man sich noch überlegen, ob und wann fünf Schritte ausreichend sind. Ich habe einfach mit scipy.optimize.root verglichen.

Welche fünf Schritte?

Diese hier:
Code:
#solution of the Kepler orbit: radius as function of the angle, ii) angle as function of time   

# 3. numerical solution via contraction --------------------

def Kepler_eq_iter(t, Om, ecc):
    ''' solution via contraction of psi = Omega * t + e * sin(psi)  '''         
    om_t = Om * t
    psi = om_t
#    for n in range(0, 5):   
#        psi = om_t + ecc * np.sin(psi)
    psi = om_t + ecc * np.sin(psi)
    return psi

Dass du
Code:
for n in range(0, 5):

auskommentierst, ist falsch! Ich glaube, du verstehst nicht, was die Schleife bedeutet; sie ist essentiell. Alles was ich sagen wollte ist, dass man sich noch anschauen muss, ob 5 immer passt.
sembel
BeitragVerfasst am: 25. Nov 2023 20:07    Titel:

Habe das JavaScript Skript (oben) nochmal um die nicht benötigten Zeilen bereinigt.

TomS hat Folgendes geschrieben:
In der Funktion zur Iteration sollte man sich noch überlegen, ob und wann fünf Schritte ausreichend sind. Ich habe einfach mit scipy.optimize.root verglichen.

Welche fünf Schritte? Du meinst das fürs Plotten? OK, dann bin ich raus.

TomS hat Folgendes geschrieben:
Die Parameter der NASA solltest du für com Erde-Mond anpassen, d.h. große Halbachse etc., auch per_rev und ecc.

Das ist ein wichtiger Hinweis! Weil wie sieht es denn aus mit dem ecc Wert, weil ja jetzt nicht mehr vom Earth-Barycenter und auch nicht vom Earth-Moon-Barycenter, sondern vom Barycenter einer gedachten gerade Ellipse (ohne Schlangenlinie) ausgegangen wird. Bei der NASA gibt es ecc nur für Earth-Barycenter. Oder?

Weitere Paramater, außer der Periodenlänge, werden nicht benötigt. Da weiterhin die offene Frage danach, welche Periodenlänge die passenste ist.

TomS hat Folgendes geschrieben:
Damit erledigt sich auch die Frage nach dem Namen.

Naja, der Programmierer sollte schon erwähnt werden.
TomS
BeitragVerfasst am: 25. Nov 2023 18:48    Titel:

Hi,

Name "Tom" ist ok.

In der Funktion zur Iteration sollte man sich noch überlegen, ob und wann fünf Schritte ausreichend sind. Ich habe einfach mit scipy.optimize.root verglichen.

Die Parameter der NASA solltest du für com Erde-Mond anpassen, d.h. große Halbachse etc., auch per_rev und ecc.

Damit erledigt sich auch die Frage nach dem Namen.
sembel
BeitragVerfasst am: 25. Nov 2023 18:09    Titel:

@TomS , alles gut. Siehe mein Edit in meinem letzten Kommentar.

Es ist jetzt geschafft! Ich danke euch recht herzlich! Insbesondere TomS für das Skript! Aber auch allen anderen Beteiligten. Habe viel gerlernt über Astronomie.

Hier das Skript in JavaScript:
Code:
// solution of the Kepler orbit: radius as function of the angle, ii) angle as function of time   
// Author: Tom
// Date: 2023-11-21
// Short version 2023-11-25 16:50


// npm import
// https://github.com/oleics/node-is-scalar/blob/master/index.js
var withSymbol = typeof Symbol !== 'undefined';
function isScalar(value) {
   
   var type = typeof(value);
   
   if(type === 'string') {return true;}
   if(type === 'number') {return true;}
   if(type === 'boolean') {return true;}
   if(withSymbol === true && type === 'symbol') {return true;}
   
   if(value == null) {return true;}
   if(withSymbol === true && value instanceof Symbol) {return true;}
   if(value instanceof String) {return true;}
   if(value instanceof Number) {return true;}
   if(value instanceof Boolean) {return true;}
   
   return false;
}


// Converter
function radtodegree(rad) {
   
   var result = rad * 180/ Math.PI;
   
   return result;
}


// Parameters
// NSSDC @ NASA https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
var ecc = 0.01671022;         // Orbit eccentricity
var per_rev = 31558149.7635456; // sidereal year (Fixstar) 365.256363004 * 86400;
// var per_rev = 31558432.539; // anomalistic year (Perihel) 365.25963586805557 * 86400;

var om_rev = (2 * Math.PI / per_rev); // angular frequency earth revolution


function theta(psi, ecc) {
   
    var theta = 2 * Math.atan(Math.sqrt((1+ecc)/(1-ecc)) * Math.tan(psi/2));
   
    if (isScalar(theta) === true) {
      if (theta < 0) {
         theta += 2*Math.PI;
      }
      return theta;
   }

    return theta;
}


// solution via contraction of psi = Omega * t + e * sin(psi)
function Kepler_eq_iter(t, Om, ecc) {

  var om_t = Om * t;
  var psi = om_t;
  var psi_check = 0;

  for (var i = 0; i < 10; i++) {
    psi = om_t + ecc * Math.sin(psi);

    if (psi == psi_check) {
      break;
    }
    psi_check = psi;
  }
  return psi;
}


// solution via contraction for Kepler eq.
function contract_Kepler_eq(t, Om, ecc) {
   
    psi = Kepler_eq_iter(t, Om, ecc);
    return theta(psi, ecc);
}


var t_arr_0 = 0;
var t_arr_90 = (per_rev / 2);
var t_arr_360 = per_rev;

// 3. 0
theta_arr_3_0 = contract_Kepler_eq(t_arr_0, om_rev, ecc);
console.log(theta_arr_3_0);

// 3. 90
theta_arr_3_90 = contract_Kepler_eq(t_arr_90, om_rev, ecc);
console.log(theta_arr_3_90);

// 3. 360
theta_arr_3_360 = contract_Kepler_eq(t_arr_360, om_rev, ecc);
console.log(theta_arr_3_360);

Eine Frage an die Beteiligten:
TomS: Genügt dir die Angabe deines Namens so?
Alle anderen: Wollt ihr auch erwähnt werden?


Letzte Frage: Im Skript wird ein siderisches Jahr verwendet. Perihel - Perihel ist jedoch ein anomalistisches Jahr. Zudem gibt es noch das Jahr nach Gregorianischen Kalender. Meine Frage ist nun, ob ein anomalistisches Jahr besser dazu passen würde, weil Perihels die Fixpunkte sind? Das Skript berechnet das Ergebniss nur anhand der Schwankungen innherhalb von Perihel-Perioden. Da wäre es doch logisch ein anomalistisches Jahr (Perihel-Perihel) zu verwenden. Oder?

Was ich bisher nicht verstanden habe, wieso Perihels/Aphels und Solstices/Equinoxe im Gregorianischen Kalender nicht kontinuierlich wandern, sondern um ein Datum schwanken. Oder ist der Effekt zu gering, um ihn in 100 oder 1000 Jahren zu spüren? Mir würde die Verwendung eines anomalistischen Jahres zusagen. Welche Bedenken gibt es dabei?
Code:
Bürgerliches Jahr / Gregorianisches Jahr
365 Tage, 5 Stunden, 49 Minuten, 12 Sekunden
365.2425 Tage
31556952 Sekunden

Tropisches Jahr (Frühlings-Tagundnachtgleiche - Frühlings-Tagundnachtgleiche)
365 Tage, 5 Stunden, 48 Minuten, 45,261 Sekunden
365.24219052083333 Tage
31556925.261 Sekunden

Anomalistisches Jahr (Perihel - Perihel)
1. Januar 2000 (J2000.0)
365 Tage, 6 Stunden, 13 Minuten, 52,539 Sekunden
365.25963586805557 Tage
31558432.539 Sekunden

Siderisches Jahr (Fixstern - Fixstern)
365 Tage, 6 Stunden, 9 Minuten, 9,76 Sekunden
365.25636296296295 Tage
31558149.759999998 Sekunden
TomS
BeitragVerfasst am: 25. Nov 2023 17:30    Titel:

@sembel – wenn du mir nochmal erklärst, was du willst, kann ich das gerne für dich erledigen. M.E. muss man nur ein paar Parameter für den com-Orbit anpassen. Alles andere kann so bleiben.


Dann: sämtliche verwendete Funktionen akzeptieren als erstes Argument eine Zahl oder ein Array; wenn letzteres, dann wird die Formel für alle Elemente des Arrays angewendet. Das gilt auch für +, -, *, **, /, np.sin() usw.

Bsp.:

np.sin(t_arr) liefert den Sinus aller tausend Werte in t_arr

Aber alle Funktionen akzeptieren auch eine Zahl. Also muss man da schon mal nix ändern.


np.where() liefert ein Array von Indizes des abgefragten Arrays zurück; und dieses Array kann man wieder in das ursprüngliche hineinstecken.

Beispiel:

arr_1 = [1.0, 2.0, 3.0, 4.0]
idx_arr = [0, 3]

arr_1 = [idx_arr] liefert [1.0, 4.0]

Aber wenn du nur mit einer Zahl arbeitest, brauchst du diesen Teil der Funktion nicht, der wird nur für Arrays durchlaufen.


Zuletzt: richtig, #1 und #2 können raus.
TomS
BeitragVerfasst am: 25. Nov 2023 17:25    Titel:

@DrStupid - ich hab das eben auch mal schnell berechnen lassen und sehe ein ähnliches Ergebnis – überraschend. Danke.

Letztlich gibt es zwei typische Massenskalen, nämlich die Masse der Erde sowie die des Mondes. Das liefert ein Verhältnis von 1/100, aus das sollte bei Taylor-Nährungen der Bahnkorrekturen ebenfalls wieder auftauchen. 1/100 eines Jahres sind aber ca. 3 1/2 Tage, kommt also hin.
sembel
BeitragVerfasst am: 25. Nov 2023 16:52    Titel:

@DrStupid , OK, der Unterschied im Unterschied (1 Minute, 1 Stunde) beträgt 1 Meter, aber der Unterschied zu ephemeris beträgt 60000 km.

Ich hatte Gestern eine Anpassung meiner Idee hier mitgeteilt und den Kommentar ein paar mal editiert. Ich hoffe das ist nicht an euch vorbei gegangen.
https://www.physikerboard.de/ptopic,393898.html#393898

Gleich darunter hatte ich noch eine Frage zum Unterschied zwischen siderischen und anomalistischen Jahr auf den Gregorianischen Kalender.
https://www.physikerboard.de/ptopic,393900.html#393900

Und in meinem letzten Kommentar habe ich die Grafik nachgereicht und nochmals ausgetauscht:
https://www.physikerboard.de/ptopic,393908.html#393908

Python → JavaScript
Zurzeit bin ich daran die Python-Formel von @TomS nach JavaScript umzuschreiben. Neben den Python-Funktionen und numpy-Funktionen, die es auch in JavaScript gibt und neben denen für die eine in npm oder eine manuell geschriebene existert, hat die Sache einige Tücken. Es werden Arrays erzeugt, ohne das dies (für mich) explizit ersichtlich ist.

Die Sache hängt gerade in der Funktion "theta" bei Zeile:
Code:
theta = np.where(theta < 0.0, theta+2*np.pi, theta)

Wegen "numpy.where".
Die Suche nach "numpy.where" gestaltet sich wegen "where" extrem schwierig.

Der Direktlink zur Formel: https://www.physikerboard.de/ptopic,393841.html#393841

Das hängt doch sicherlich mit dem "plot" zusammen, für den per "linspace" Werte zwischen 0 und 1000 generiert werden. "plot" fliegt bei mir komplett raus. Es wird nur eine Zahl als Ergebnis benötigt.

Also das Skript von TomS muss erst einmal auf das Wesentliche gekürzt werden.
Code:
#solution of the Kepler orbit: radius as function of the angle, ii) angle as function of time   
#Author: Tom
#Date: 2023-11-21
#Short version 2023-11-25 16:50

# imports --------------------

import numpy as np

# constants --------------------

class OrbParams:
    # NSSDC @ NASA, also wikipedia: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
    def __init__(self, ecc=0.0167086, a=149.598022960e9, per_rev=365.256363004*86400, per_rot=86164.10052):                                                                       
        self.ecc = ecc                                      # eccentricity
        self.a = a                                          # semi major axis   
        self.per_rev = per_rev                              # sideral year
        self.per_rot = per_rot                              # sideral day
      
        self.om_rev = 2 * np.pi / self.per_rev              # angular frequency earth revolution
        self.om_rot = 2 * np.pi / self.per_rot              # angular frequency earth rotation

# Creating an instance of MyObject
the_orb = OrbParams()


# geometry --------------------

def theta(psi, ecc):
    ''' returns the angle as a function of te eccentric anomaly psi; psi currently restricted to [0, 2*pi] '''
    ''' THERE'S STILL AN ISSUE WITH THE INTERVAL! '''
    theta = 2 * np.arctan( np.sqrt((1+ecc)/(1-ecc)) * np.tan(psi/2))
    if np.isscalar(theta):
        if theta < 0.0:
            theta += 2*np.pi
        return theta           
#    theta = np.where(theta < 0.0, theta+2*np.pi, theta)
#    theta[0] = 0.0
    return theta


# 3. numerical solution via contraction --------------------

def Kepler_eq_iter(t, Om, ecc):
    ''' solution via contraction of psi = Omega * t + e * sin(psi)  '''         
    om_t = Om * t
    psi = om_t
    for n in range(0, 5):   
        psi = om_t + ecc * np.sin(psi)
    psi = om_t + ecc * np.sin(psi)
    return psi

def contract_Kepler_eq(t, Om, ecc):
    ''' solution via contraction for Kepler eq. '''                 
    psi = Kepler_eq_iter(t, Om, ecc)
    return theta(psi, ecc)

   
# driver code --------------------

#t_arr = np.linspace(0.0, the_orb.per_rev, 1000)
t_arr_0 = 0.0
t_arr_90 = 365.256363004*86400 / 2
t_arr_360 = 365.256363004*86400

# 3. 0
theta_arr_3_0 = contract_Kepler_eq(t_arr_0, the_orb.om_rev, the_orb.ecc)
delta_theta_arr_3_0 = theta_arr_3_0 - the_orb.om_rev * t_arr_0
# 3. 90
theta_arr_3_90 = contract_Kepler_eq(t_arr_90, the_orb.om_rev, the_orb.ecc)
delta_theta_arr_3_90 = theta_arr_3_90 - the_orb.om_rev * t_arr_90
# 3. 360
theta_arr_3_360 = contract_Kepler_eq(t_arr_360, the_orb.om_rev, the_orb.ecc)
delta_theta_arr_3_360 = theta_arr_3_360 - the_orb.om_rev * t_arr_360


print(theta_arr_3_0)
print(delta_theta_arr_3_0)

print(theta_arr_3_90)
print(delta_theta_arr_3_90)

print(theta_arr_3_360)
print(delta_theta_arr_3_360)

Damit fällt "numpy.linspace" und "numpy.where" schon mal weg.
Problematisch sind jetzt noch die Ergebnisse:
Code:
0.0
0.0
3.141592653589793
0.0
6.283185307179586
0.0

Da sollten Zahlen zwischen 0 und 360 herauskommen.
Ich bitte um Hilfe.

Edit: Sorry, meine Dusseligkeit. Die Ergebnisse sind in Bogenmaß.
DrStupid
BeitragVerfasst am: 25. Nov 2023 14:06    Titel:

Mit den Werte aus meiner Berechnung oben habe ich jetzt auch mal ermittelt, welchen Fehler man macht, wenn man Erde und Mond als einen Körper betrachtet, anstatt den gemeinsamen Masseschwerpunkt aus ihren separaten Bahnen zu berechnen. Die Abweichung ist mit ca. 20 km tatsächlich sehr gering. Die Dominanz der jährlichen Schwankung hat mich allerding überrascht. Die monatlichen Schwankungen sind kaum zu erkennen.
DrStupid
BeitragVerfasst am: 25. Nov 2023 13:32    Titel:

sembel hat Folgendes geschrieben:
@DrStupid , deine Angabe "60000 km" vestehe ich nicht.


Ich habe mir von https://ssp.imcce.fr/forms/ephemeris Positionen und Geschwindigkeiten von Sonne, Erde, Mond und allen anderen Planeten für 2023-01-01 00:00:00 ausgeben lassen. Mit diesen Startwerten habe ich dann die Positionen für 2024-01-01 00:00:00 berechnet und mit den entsprechenden Werten aus https://ssp.imcce.fr/forms/ephemeris verglichen. Bei der Erde kommt dabei ein Unterschied von 60000 km raus - sowohl für die absolute Position, als auch für den Abstand zur Sonne.

sembel hat Folgendes geschrieben:
Ebenso nicht deine Angabe "Der Rechenfehler liegt unter einem Meter."


Ich habe die Berechnung mit Schrittweiten von einer Minute und einer Stunde durchgeführt und dann beide Ergebnisse vergleichen. Der Unterschied liegt für die Erde deutlich unter einem Meter. Das ist der numerische Fehler der Simulation. Die Ursache für die oben genannte Abweichung von 60000 km liegt also nicht in der Berechnung, sondern irgendwo anders.

sembel hat Folgendes geschrieben:
Ein weiteres Problem sind die verwendeten Baryzentren.


Es ist nur wichtig, dass sich alle Werte auf das gleiche Baryzentrum beziehen. Am besten verwendet man den gemeinsame Masseschwerpunkt alle beteiligten Körper. Es ist aber auch nicht schlimm, wenn man ein anderes Baryzentrum wählt. Außer, dass sich der Masseschwepunkt bewegt, passiert dabei nichts Schlimmes. Ich habe Werte für den Masseschwerpunkt des Sonnensystems verwendet.

sembel hat Folgendes geschrieben:
Einerseits müssten alle Einflüsse (Planeten, etc.) mit inbegriffen werden.


Das ist bei meiner Rechnung der Fall.
DrStupid
BeitragVerfasst am: 25. Nov 2023 12:43    Titel:

DrStupid hat Folgendes geschrieben:
Bei der geringen Krümmung und Exzentrizität der Erdbahn kann auch eine Schwankung des Abstandes um weniger als einen Erdradius große Auswirkungen haben.


Ich hab das jetzt mal für die Erdbahn (https://tinyurl.com/mrx6vhvr) und für den Masseschwerpunkt von Erde und Mond durchgerechnet (https://tinyurl.com/yckam5yk). Dabei erreicht die Erde den maximalen Abstand zur Sonne rund 90000 Sekunden später als der Masseschwerpunkt. Die vergleichsweise winzige Schwankung des Abstandes um weniger als 5000 km führt also tatsächlich zu einer zeitlichen Differenz um mehr als einen Tag. Weil das in beide Richtungen gehen kann, sind Unterschieden von mehreren Tagen beim Vergleich zwischen verschiedenen Jahren völlig plausibel.
sembel
BeitragVerfasst am: 25. Nov 2023 11:55    Titel:

Ein paar Erklärungen von meiner Seite. Es folgt in den nächsten Stunden noch eine Grafik.

Die EMB-Perihel-Zeitpunkte schwanken um bis zu 3 Stunden (2023-2025). Das ist schon immens weniger als die EB-Perihel-Zeitpunkte mit bis zu mehreren Tagen. Das bedeutet also, dass der EMB auf seiner Bahn auch eine Schlangenlinie vollführt, nur weniger schlänglich. Das hat zur Ursache, dass die Perihel-Zeitpunkte früher/später eintreten als auf einer Bahn ohne Schlangenlinien.

"Schlangenlinie" bedeutet, es wird "getrickst" in dem sich aus dem Fenster gelehnt wird, um früher/später am sonnennahesten zu sein.

@DrStupid , deine Angabe "60000 km" vestehe ich nicht. Ebenso nicht deine Angabe "Der Rechenfehler liegt unter einem Meter.".

Perihel 30,29 km/s (Internet)
Aphel 29,29 km/s (Internet)
Ø 29.79 km/s (rechnerisch)
Ø 29,7859 km/s (Internet)

Perihel = 147.095.000 Mio km (Internet)
Aphel = 152.100.000 Mio km (Internet)
Länge Erdbahn (EMB Bahn) = 940.000.000 km (Internet)

Angenommene Abweichung = 3 Stunden
((3 * (60 * 60)) * 30.29) = 327132 km

Das wäre eine rein rechnerische Abweichung vom theroretischen Gerade-Bahn Perihel. Aber dem ist nicht so, weil auch dies eine gerade Bahn voraussetzt. Aber es gibt an, dass die Abweichung kleiner sein muss, und das ist sie.

Ein weiteres Problem sind die verwendeten Baryzentren. Einerseits müssten alle Einflüsse (Planeten, etc.) mit inbegriffen werden. Andererseits ist das zu aufwenidg. Es finden sich jedoch nur Daten zu SSB (Solar-System-Barycenter) und EMB (Earth-Moon-Barycenter). Diese Werte können also mit der Formel von @TomS, die (nur) ein Sun-Earth-Barycenter anstatt einem SSB verwendet, als (grober) Vergleich herangezogen werden. Es wird also Abweichungen geben. Das ist tolerierbar.
DrStupid
BeitragVerfasst am: 25. Nov 2023 11:09    Titel:

TomS hat Folgendes geschrieben:
Damit stellt sich nicht nur die Frage, wie diese Bahnkorrekturen aufgrund des Mondes zu einer Abweichung der Position der Erde von einigen Tagen führen können, es stellt sich die Frage, bzgl. was denn diese Abweichung überhaupt gelten soll!


Ich habe oben versucht, es zu erklären. Das Problem liegt nicht bei der Position der Erde, sondern bei den Positionen von Aphel und Perihel. Bei der geringen Krümmung und Exzentrizität der Erdbahn kann auch eine Schwankung des Abstandes um weniger als einen Erdradius große Auswirkungen haben.

TomS hat Folgendes geschrieben:
Wenn der Durchgang der Erde durch den Perihel um einige Tage schwankt: ggü. was denn?


Beispielsweise im Vergleich der Zeiten von Perihel zu Perihel. Die unterscheiden sich von Jahr zu Jahr deutlich. Eine andere Möglichkeit wäre der Vergleich des Perihels der Erdbahn mit dem der Kepplerellipse des Erde-Mond-Systems. Das liegt in der gleichen Größenordnung.
TomS
BeitragVerfasst am: 25. Nov 2023 09:07    Titel:

DrStupid hat Folgendes geschrieben:
TomS hat Folgendes geschrieben:
Hatten wir das hier schon diskutiert? Wie kommst darauf, dass dieses Objekt deutlich besser einer Kepler-Ellipse folgt?


Dass das deutlich besser ist war meine Vermutung, weil Dipol-Wechselwirkungen üblicherweise mit 1/r³ abfallen. Der Fehler sollte also viel kleiner sein als das Herumeiern der Erde um den Masseschwerpunkt. Genau weiß man das natürlich erst, wenn man es ausgerechnet hat.

Die Näherung sieht absolut sinnvoll aus, allerdings habe ich ein grundsätzliches Problem, das auch die wortreiche Wikipedia-Artikel ohne eine einzige Gleichung nicht mal ansatzweise klären können – obwohl sie behaupten, das zu tun.

Ich habe einen gedachten Keplerorbit für die Erde alleine, mit a'(E), T'(E) sowie den Ortsvektor r'(E). Dann habe ich einen in ziemlich guter Näherung zutreffenden Keplerorbit für das System Erde-Mond um die Sonne, mit a(E,M), T(E,M) und r(E,M). Außerdem weiß ich, das für die wahre Erdposition r(E) und den wahren Massenmittelpunkt r(E,M) gilt



Diese Abweichung ist winzig!

Damit stellt sich nicht nur die Frage, wie diese Bahnkorrekturen aufgrund des Mondes zu einer Abweichung der Position der Erde von einigen Tagen führen können, es stellt sich die Frage, bzgl. was denn diese Abweichung überhaupt gelten soll! Wenn der Durchgang der Erde durch den Perihel um einige Tage schwankt: ggü. was denn?

Ich habe schon Probleme, die Frage zu formulieren.
TomS
BeitragVerfasst am: 24. Nov 2023 21:57    Titel:

Ich schaue mir die Störungstheorie zu DrStupids Idee mal an.
DrStupid
BeitragVerfasst am: 24. Nov 2023 20:53    Titel:

Ich habe jetzt mal das Sonnensystem mit allen Planeten und dem Mond mit mit Startwerten für 2023-01-01 00:00:00 aus https://ssp.imcce.fr/forms/ephemeris in Stundenschritten für ein Jahr simuliert:

https://tinyurl.com/bddkdwj7

Der Vergleich mit den gemäß https://ssp.imcce.fr/forms/ephemeris korrekten Werten für 2024-01-01 00:00:00 offenbart einen Fehler von rund 60000 km. Ich dachte erst, dass die Schrittweite zu groß ist und habe es nochmal mit einer Minute probiert, aber das Ergebnis bleibt fast unverändert. Der Rechenfehler liegt unter einem Meter.

Wenn die Daten aus https://ssp.imcce.fr/forms/ephemeris stimmen, dann fällt mir außer der Gravitationskonstante und den Massen der Himmelskörper keine Fehlerquelle ein.
sembel
BeitragVerfasst am: 24. Nov 2023 20:41    Titel:

Zusatz:

Bisher wird in der Formel ein siderisches Jahr verwendet. Welche Auswirkung hätte es, wenn ein anomalistisches Jahr verwendet werden würde? Die Idee ein anomalistisches Jahr zu verwenden rührt daher, weil der Beginn ein Perihel ist.

Code:
Anomalistisches Jahr (Perihel - Perihel)
1. Januar 2000 (Epoche J2000.0)
31558432.539 Sekunden
365 Tage, 6 Stunden, 13 Minuten und 52,539 Sekunden
(365 + (((6 * 60 * 60) + (13 * 60) + 52.539) / 86400))
365.25963586805557 Tage


Was ich dabei nicht ganz verstanden habe, was den Unterschied zwischen siderischen und anomalistischen auf den Gregorianischen Kalender ausmacht? Sowohl Perihel als auch Frühlings-Tagundnachtgleiche schwanken um die selben Daten, also sie wandern nicht in eine Richtung. Oder täusche ich mich da, weil es nur über mehrere Jahrhunderte spürbar ist? Da wären die Schalttage zu berücksichtigen, die alle 400 Jahre ausgeschalten werden.
sembel
BeitragVerfasst am: 24. Nov 2023 19:20    Titel:

DrStupid hat Folgendes geschrieben:
Doch, der eiert auch - aber wahrscheinlich viel weniger.

Super vielen Dank für diese anschauliche Info. Das bringt mich dazu meine Idee neu zu überdenken.

Weil sehr wahrscheinlich, gehe ich jetzt davon aus: @TomS hatte vollkommen recht damit, bzw. lag vollkommen richtig, dass die Perihel-Jahre alle gleich lang sind. Wie im Skript angegeben:
Code:
per_rev=365.256363004*24*60*60

Denn die Unterschiede in den EMB-Perihel-Jahren ergeben sich vermutlich aus der Schlangenlinie des EMB auf seiner Bahn um das Solar System Barycenter (SSB). Würde also diese Schwankung (Schlangenlinie) per Formel ausgeglichen werden, so dass von einer theoretischen "geraden" EMB-Bahn ausgegangen werden würde, dann wären die EMB-Perihel-Jahre alle gleichlang, bzw. zumindest sehr gleich.

Diese Überlegung ergibt folgende Anpassung meiner Idee:
Als initialer Perihel-Zeitpunkt wird der EMB-Perihel-Timestamp aus dem Jahr 2000 verwendet. Von diesem aus wird ein Wrapping mit "365.256363004*24*60*60" auf die ab dem EMB-Perihel-Timestamp-Year2000 fortlaufenen Sekunden aufgeschalten.

Perihelion Earth-Moon barycenter:
2000 01 03 23:48:11
Timestamp (seconds): 946943291

Ein bisschen JavaScript:
Code:
// Get timestamp in seconds
new Date(Date.UTC(2000, 0, 3, 23, 48, 11, 0)).getTime() / 1000
// 946943291
// months begin with 0 = January
// JavaScript timestamp is in Milliseconds

// Check:
new Date(946943291 * 1000).toUTCString()
// "Mon, 03 Jan 2000 23:48:11 GMT"

365.256363004*24*60*60
// 31558149.7635456 seconds (wrap)

Das Wrapping:
Code:
// Get seconds since Perihelion year 2000
var timestamp = Math.round(new Date().getTime() / 1000);
var incremental = 946943291; // Perihelion year 2000
var current = timestamp - incremental;
console.log(current); // 753906224 seconds since Perihelion year 2000

// Get seconds since the wrap Perihelion
var year_seconds = 31558149.7635456;
var current_wrap = Math.round(current % year_seconds);
console.log(current_wrap); // 28068779 seconds since the wrap Perihelion

// Check:
var check = (((28068222 / 60) / 60) / 24);
console.log(check ); // 324.86368055555556 days since the wrap Perihelion


DrStupid hat Folgendes geschrieben:
Vielleicht kann man dafür noch eine weniger hässliche Nährung finden, aber man sieht auch so, dass es mit der gleichen Geschwindigkeit rotiert, wie Mond und Erde umeinander.

Braucht es für meine Idee jetzt nicht mehr.

@TomS , die letzte Bitte der Anpassung des Skripts an schwankende Perihel-Jahre ist jetzt nicht mehr notwendig. Das Wrapping muss nicht implementiert werden, weil es nur den Input betrifft.

PS: So ist das manchmal in der Prozess-Entwicklung. Manchmal werden ganze Projekte kurz vor Fertigstellung gelöscht und neu angefangen, weil sie nicht mehr zur Idee-Anpassung passen. So schlimm ist es hier aber nicht. Sorry dafür.
DrStupid
BeitragVerfasst am: 24. Nov 2023 16:48    Titel:

sembel hat Folgendes geschrieben:
Ist auch irgendwie logisch, weil EMB nicht eiert


Doch, der eiert auch - aber wahrscheinlich viel weniger. Wenn beide Massen im Schwerpunkt vereint wären, dann würde



gelten. Wenn man Mond und Ede unterscheidet, dann gilt stattdessen



Dabei ist die Störung, die das Eiern verursacht. Mit









und



ergibt das


Vielleicht kann man dafür noch eine weniger hässliche Nährung finden, aber man sieht auch so, dass es mit der gleichen Geschwindigkeit rotiert, wie Mond und Erde umeinander.
sembel
BeitragVerfasst am: 24. Nov 2023 15:46    Titel:

TomS hat Folgendes geschrieben:
Kannst du die URLs einstellen?

Welche URLs meinst du? Oder meinst du die Quellen?

Das hier ist die Quelle der Perihel-Daten mit Erdmittelpunkt:
https://aa.usno.navy.mil/api/seasons?year=2023

In diesem Buch sind die Perihel-Daten mit Erdmittelpunkt und EMB:
Da gibt es nur eine Vorschau und wenn du Glück hast, wird dir die Seite 116 mit den Daten der Jahre 1989-2059 angezeigt. Das ist aber eher unwahrscheinlich. Ich hatte Gestern etwas Glück und habe die Seite sehen können (und gespeichert). Die Daten von 2023-2026 habe ich in Unix Timestamps umgewandelt und hier veröffentlicht.
https://www.google.de/books/edition/SEASONS_EASTER_AND_ASTRONOMICAL_MOTIONS/L1sjBgAAQBAJ
PS: Deswegen auch meine Frage nach den URL Parametern zur NASA Horizons API.

Powered by phpBB © 2001, 2005 phpBB Group