Mars4




			Translation

				Horst Schumacher

Lage der auf einer festen (Laplaceschen) Ebene bezogenen Mondbahn

Die Bahnelemente J u. N (Äquinoktium des Datums) liegen nicht fest, sondern unterliegen lfd. Störungen. Die Satellitenbahn wird daher auf eine feste Referenzebene (in der Regel die Laplacesche Ebene) bezogen, wodurch die Schwingungen des gestörten Bahnpols auf den Pol der festliegenden Ebene bezogen werden können. Im vorliegenden Fall umbahnen die Marsmonde den Planeten sehr eng, wodurch hauptsächlich die Sonne und die Abweichung des Mars von der Kugelform säkulare Störungen erzeugen, deren Berücksichtigung relativ zur Laplaceschen Ebene im vorliegenden Fall bereits eine gute Bahnwiedergabe ergeben.
Da die Laplacesche Ebene stets durch beide Schnittpunkte (Äquinoktien) des Marsäquators auf der Marsbahn verläuft (die gegenseitige Anziehung der Satelliten verursacht geringfügige Fluktuationen), unterliegt jedoch auch die Laplacesche Ebene der Präzessionsbewegung der Äquinoktien mit -0.1975° pro Jahrhundert und der säkularen Störung der Marsbahn.

Bahnelemente L (Länge des Satelliten),P (Länge des Perihels), K (Polwinkel) zur Epoche JDo.

dN1,dJ1, c und ho folgen in erster Näherung aus der auf die Erdäquatorebene bezogenen Bedingungsgleichung:

A) (N-N1)*sin(J1)=sin(J1)*dN1+ c*sin(ho)*cos(K*t)- c*cos(h o)*sin(K*t)
B) J-J1=dJ1+ c*sin(ho)*sin(K*t)+c*cos(h o)*cos(K*t)

Bedingungsgleichung A:
a(n)=cos(K*t)
b(n)=-sin(K*t)
c(n)=(N-N1)*sin(J1)

[px]   = S g(i)*a(n), 1,2,3,...,n
[pxx] = S g(i)*a(n)*a(n), 1,2,3,...,n
[pxy] = S g(i)*a(n)*b(n), 1,2,3,...,n
[pyy] = S g(i)*b(n)*b(n), 1,2,3,...,n
[py]   = S g(i)*b(n), 1,2,3...,n
[pd]   = S g(i)*c(n),1,2,3,...,n
[pdx] = S g(i)*a(n)*c(n),1,2,3,...,n
[pdy] = S g(i)*b(n)*c(n),1,2,3,...,n

Normalgleichung:
   p c + [px] a + [py]   + b = [pd],1,0,0
[px] c + [pxx] a + [pxy] + b = [pdx],0,1,0
[py] c + [pxy] a + [pyy] + b = [pdy],0,0,1

c=dN1/sin(J1)
a = c sin(h)
b = c cos(h)
c = Å(a²+b²)
y=a/c
x=b/c
ho = ARCTAN(y/(1+x))*2.

Bedingungsgleichung B:
a(n)=sin(K*t)
b(n)=cos(K*t)
c(n)=J-J1

[px1]   = S g(i)*a(n), 1,2,3,...,n
[pxx1] = S g(i)*a(n)*a(n), 1,2,3,...,n
[pxy1] = S g(i)*a(n)*b(n), 1,2,3,...,n
[pyy1] = S g(i)*b(n)*b(n), 1,2,3,...,n
[py1]   = S g(i)*b(n),1,2,3,...,n
[pd1]   = S g(i)*c(n),1,2,3,...,n
[pdx1] = S g(i)*a(n)*c(n),1,2,3,...,n
[pdy1] = S g(i)*b(n)*c(n),1,2,3,...,n

Normalgleichung:
    p   c + [px1] a + [py1] + b = [pd1],1,0,0
[px1] c + [pxx1] a + [pxy1] + b = [pdx1],0,1,0
[py1] c + [pxy1] a + [pyy1] + b = [pdy1],0,0,1

c=dJ1, a = c sin(h), b = c cos(h).
c = Å(a²+b²); y=a/c , x=b/c , ho = ARCTAN(y/(1+x))*2

t in Jahren ab der Epoche JDo: t=(JD-JDo)/365.25.
J1 = (Jm) Winkel der Marsäquatorebene mit der Erdäquatorebene, N1 (Nm) = auf den Erdäquator gemessene Rektaszension des aufsteigenden Knotens der Marsäquatorebene auf der Erdäquatorebene.
c = Winkel der Trabantenbahn mit der festen Ebene, ho = auf die feste Ebene bezogene Länge zwischen dem Schnittpunkt der festen Ebene auf dem Erdäquator (Nf) und dem der festen Ebene auf der Mondbahn; zur Epoche JDo (t=0).
J,N=Winkel u. Rektaszension (N) des aufst. Knoten der Mondbahnebene auf der Erdäquatorebene (in Rad).

Bei bekanntem Wert ho bezieht man die Bahnelemente N,J,P,L der Epoche JDo auf die Laplacesche Ebene (Nf,Jf, ho).
Jf=Winkel der Laplaceschen Ebene gegen die Erdäquatorebene, Nf=Rektaszension des Knotens der Laplaceschen Ebene auf der Erdäquatorebene.

Numerisches Beispiel Phobos:
Epoche 10.1.1967, 21h TT = JDo 2439501.375. In erster Näherung für die jährl. Bewegung des Knotens K, wird der berechnete Wert K (Kj) 159,003° angenommen: K=jährl. Änderung von ho (genäherter Wert K s. Drehung der Apsidenachse und der festen Ebene der Mondbahn infolge der Störung durch die Abplattung des Planeten ).

JD = julian. Datum der beobachteten J,N.
JD=2440952.5, J=36.9196, N=49.332093, Gewicht g=1.
JD=2441865.5, J=37.61435, N=49.239938, g=1
K = 159.003°.

Nach dem folg. Progr. ist die Fehlerquadratsumme und der mittl. Fehler der Koeffizienten a,b,c der Bedingungsgleichung B kleiner als die Fehlergrenze der Bedingungsgleichung A.
Akzeptierte Werte der Bedingungsgleichung B in erster Näherung: Konstanter Winkel der Mondbahnebene Phobos gegen die Laplaceschen Ebene: c =1.03°! 0.005°.
Aufsteigender Knoten der Mondbahnebene des Phobos auf der Laplaceschen Ebene Epoche JDo=2439501.375: ho=16.37° ± 44.02°.
Damit erhält man den aufst. Knoten für einen beliebigen Zeitpunkt (JD) ab der Epoche JDo in erster Näherung: h = (16.37°-159.003° * ((JD-2439501.375)/365.25).
Deklin. Nordpol Mars J2000 in erster Näherung: +52.926° ± 0.0019°.
Rektaszension Nordpol Mars J2000 in erster Näherung: 317.782° ± 0.0043°.

Koordinaten des Marsnordpols lt. A.T. Sinclair, >The orbits of the satellites of Mars determined from Earth-based and spacecraft observations< [Astron. Astrophys. 220, 321-328 (1989)]
d 1 +52.7226° (Äquinoktium J1950), d1 +52.928° (J2000)
a 1 317.3966° (J1950), a1 317.789° (J2000).

REM ERDÄQUATOREBENE/LAPLACESCHE EBENE PHOBOS
DIM p(10,10),ko(10,10),a(20),b(20),c(20),d(20),a1(20),b1(20),c1(20),g(20)
DEFFN RAD(x) = x * (PI / 180)
DEFFN deg(x) = x * (180 / PI)
DEFFN r(x) = x - INT(x / (2 * PI)) * (2 * PI)
REM --------------------
n1 = 7 //EINTRAG ANZAHL WERTPAARE J,N
jm = PI / 2 - FN RAD(52) //EINTRAG UNGEFHRE DEKLINATION NORDPOL MARS QUINOKTIUM J2000
nm = FN r(RAD(317) + PI / 2) //EINTRAG UNGEFHRE REKTASZENSION NORDPOL MARS QUINOKTIUM J2000
k = FN RAD(159.003)
RESTORE dat
FOR i = 1 TO n1
REM EPOCHE ÄQUINOKTIUM J2000 JDe=J2000=JD 2451545------------------------
jde = 2451545
READ jd,j,n,g
g(i) = g
j = FN RAD(j) //J ÄQUINOKTIUM DES DATUMS
n = FN RAD(n) //N ÄQUINOKTIUM DES DATUMS
GOSUB reduk
j = jj   //J=QUINOKTIUM J2000
n = nn //N=QUINOKTIUM J2000
REM EPOCHE 10.1.1967, 21h TT, JDo=2439501.375
t = (jd - 2439501.375) / 365.25
REM BEDINGUNGSGLEICHUNG A
a(i) = COS(k * t)
b(i) = -SIN(k * t)
c(i) = (n - nm) * SIN(jm)
REM BEDINGUNGSGLEICHUNG B
a1(i) = SIN(k * t)
b1(i) = COS(k * t)
c1(i) = j - jm
NEXT i
dat: //DATA-DATEI JD,J,N,ÄQUINOKTIUM DES DATUMS JD
DATA 2440952.5,36.9196,49.332093,1
DATA 2441865.5,37.61435,49.239938,1
DATA 2442747.5,37.97487,48.820197,1
DATA 2443478.5,37.26688,49.423955,1
DATA 2444025.5,36.1684,46.77792,1
DATA 2444786.5,36.52682,46.109697,1
DATA 2445578.5,36.7814,45.9154,1
REM ----------------------------
REM METHODE DER KLEINSTEN QUADRATE
REM AKKUMULATION
x = 0
xx = 0
xy = 0
y = 0
yy = 0
yz = 0
d = 0
dx = 0
dy = 0
x1 = 0
xx1 = 0
xy1 = 0
yy1 = 0
y1 = 0
yz1 = 0
d1 = 0
dx1 = 0
dy1 = 0
dd = 0
dd1 = 0
p = 0
FOR i = 1 TO n1
REM BEDINGUNGSGLEICHUNG A
p = p + g(i) //GEWICHT
x = x + a(i) * g(i)
xx = xx + a(i) * a(i) * g(i)
xy = xy + a(i) * b(i) * g(i)
y = y + b(i) * g(i)
yy = yy + b(i) * b(i) * g(i)
d = d + c(i) * g(i)
dx = dx + a(i) * c(i) * g(i)
dy = dy + b(i) * c(i) * g(i)
dd = dd + c(i) * c(i) * g(i)
REM BEDINGUNGSGLEICHUNG B
x1 = x1 + a1(i) * g(i)
xx1 = xx1 + a1(i) * a1(i) * g(i)
xy1 = xy1 + a1(i) * b1(i) * g(i)
y1 = y1 + b1(i) * g(i)
yy1 = yy1 + b1(i) * b1(i) * g(i)
d1 = d1 + c1(i) * g(i)
dx1 = dx1 + a1(i) * c1(i) * g(i)
dy1 = dy1 + b1(i) * c1(i) * g(i)
dd1 = dd1 + c1(i) * c1(i) * g(i)
NEXT i
REM BEDINGUNGSGLEICHUNG A --------
resi = 4 //EINTRAG ANZAHL RESIDUEN
m = 3    //EINTRAG ANZAHL GLEICHUNGEN
n = 3    //EINTRAG ANZAHL UNBEKANNTE
REM ----------------
m = m + resi
p(1,1) = p
p(1,2) = x
p(1,3) = y
p(1,4) = d   //1. RESIDUENVEKTOR
p(1,5) = 1   //2. RESIDUENVEKTOR
p(1,6) = 0
p(1,7) = 0
p(2,1) = x
p(2,2) = xx
p(2,3) = xy
p(2,4) = dx   //1. RESIDUENVEKTOR
p(2,5) = 0    //2. RESIDUENVEKTOR
p(2,6) = 1
p(2,7) = 0
p(3,1) = y
p(3,2) = xy
p(3,3) = yy
p(3,4) = dy //1. RESIDUENVEKTOR
p(3,5) = 0   //2. RESIDUENVEKTOR
p(3,6) = 0
p(3,7) = 1
GOSUB elim
c = ko(1,1)
a = ko(2,1)
b = ko(3,1)
REM FEHLERQUADRATSUMME
vv = dd - d * ko(1,1) - dx * ko(2,1) - dy * ko(3,1)
REM FEHLERQUADRATSUMME BEOBACHTUNG MINUS RECHNUNG
vv1 = 0
FOR i = 1 TO n1
vv1 = vv1 + (c(i) - (ko(1,1) + a(i) * ko(2,1) + b(i) * ko(3,1))) ^ 2 * g(i)
NEXT i
s = ABS(vv / (n1 - 3))
s1 = SQR(s)
mfa = SQR(ABS(s * ko(1,2)))
mfb = SQR(ABS(s * ko(2,3)))
mfc = SQR(ABS(s * ko(3,4)))
IF n1 <> p THEN
s1 = s1 / SQR(p)
ENDIF
REM -------------
REM BEDINGUNGSGLEICHUNG B ----------
resi = 4 //EINTRAG ANZAHL RESIDUEN
m = 3    //EINTRAG ANZAHL GLEICHUNGEN
n = 3    //EINTRAG ANZAHL UNBEKANNTE
REM ----------------
m = m + resi
p(1,1) = p
p(1,2) = x1
p(1,3) = y1
p(1,4) = d1 //1. RESIDUENVEKTOR
p(1,5) = 1   //2. RESIDUENVEKTOR
p(1,6) = 0
p(1,7) = 0
p(2,1) = x1
p(2,2) = xx1
p(2,3) = xy1
p(2,4) = dx1 //1. RESIDUENVEKTOR
p(2,5) = 0    //2. RESIDUENVEKTOR
p(2,6) = 1
p(2,7) = 0
p(3,1) = y1
p(3,2) = xy1
p(3,3) = yy1
p(3,4) = dy1 //1. RESIDUENVEKTOR
p(3,5) = 0    //2. RESIDUENVEKTOR
p(3,6) = 0   //3. RESDIUENVEKTOR
p(3,7) = 1   //4. RESIDUENVEKTOR
GOSUB elim
c1 = ko(1,1)
a1 = ko(2,1)
b1 = ko(3,1)
REM FEHLERQUADRATSUMME
vv2 = dd1 - d1 * ko(1,1) - dx1 * ko(2,1) - dy1 * ko(3,1)
REM FEHLERQUADRATSUMME BEOBACHTUNG MINUS RECHNUNG
vv3 = 0
FOR i = 1 TO n1
vv3 = vv3 + (c1(i) - (ko(1,1) + a1(i) * ko(2,1) + b1(i) * ko(3,1))) ^ 2 * g(i)
NEXT i
s2 = ABS(vv2 / (n1 - 3))
s3 = SQR(s2)
mfa1 = SQR(ABS(s2 * ko(1,2)))
mfb1 = SQR(ABS(s2 * ko(2,3)))
mfc1 = SQR(ABS(s2 * ko(3,4)))
IF n1 <> p THEN
s3 = s3 / SQR(p)
ENDIF
REM ---------------------
j1 = jm + c1
n1 = nm + c / SIN(j1)
dp1 = PI / 2 - j1 //DEKLIN. MARSNORDPOL J2000 IN ERSTER NHERUNG
ap1 = FN r(n1 - PI / 2) //REKTASZENSION MARSNORDPOL J2000 IN ERSTER NHERUNG
REM ------------------
g = SQR(a * a + b * b)
ao = a / g
bo = b / g
ko = FN r(ATN(ao / (1 + bo)) * 2)
mfg = SQR(mfa * mfa + mfb * mfb)
fa = mfa / mfg
fb = mfb / mfg
fko = FN r(ATN(fa / (1 + fb)) * 2)
g1 = SQR(a1 * a1 + b1 * b1)
ao = a1 / g1
bo = b1 / g1
ko1 = FN r(ATN(ao / (1 + bo)) * 2)
mfg1 = SQR(mfa1 * mfa1 + mfb1 * mfb1)
fa = mfa1 / mfg1
fb = mfb1 / mfg1
fko1 = FN r(ATN(fa / (1 + fb)) * 2)
REM -------------------------------;;M
PRINT "BEDINGUNGSGLEICHUNG A:"
PRINT "DEKLINATION   NORDPOL MARS J2000: ";FN deg(dp1);" +/- ";FN deg(mfc1);" GRAD"
PRINT "REKTASZENSION NORDPOL MARS J2000: ";FN deg(ap1);" +/- ";FN deg(mfc) / SIN(j1);" GRAD"
PRINT "WINKEL LAPLAC. EBENE/MONDBAHNEBENE: ";FN deg(g);" +/- ";FN deg(mfg);" GRAD"
PRINT "KNOTEN LAPLAC. EBENE/MONDBAHNEBENE: ";FN deg(ko);" +/- ";FN deg(fko);" GRAD"
PRINT "KOEFFIZIENT;a ";FN deg(a);" GRAD"
PRINT "KOEFFIZIENT;b ";FN deg(b);" GRAD"
PRINT "KOEFFIZIENT;c ";FN deg(c);" GRAD"
PRINT "MITTL. FEHLER EINER EINZELMESSUNG: ";FN deg(s1)
PRINT "MITTL. FEHLER a....: ";FN deg(mfa)
PRINT "MITTL. FEHLER b....: ";FN deg(mfb)
PRINT "MITTL. FEHLER c....: ";FN deg(mfc) / SIN(j1)
PRINT "FEHLERQUADRATSUMME.: ";vv;" / ";vv1
PRINT "BEDGINGUNGSGLEICHUNG B:"
PRINT "WINKEL LAPLACE. EBENE/MONDBAHNEBENE: ";FN deg(g1);" +/- ";FN deg(mfg1);" GRAD"
PRINT "KNOTEN LAPLACE. EBENE/MONDBAHNEBENE: ";FN deg(ko1);" +/- ";FN deg(fko1);" GRAD"
PRINT "KOEFFIZIENT;a1 ";FN deg(a1);" GRAD"
PRINT "KOEFFIZIENT;b1 ";FN deg(b1);" GRAD"
PRINT "KOEFFIZIENT;c1 ";FN deg(c1);" GRAD"
PRINT "MITTL. FEHLER EINER EINZELMESSUNG: ";FN deg(s2)
PRINT "MITTL. FEHLER a1...: ";FN deg(mfa1)
PRINT "MITTL. FEHLER b1...: ";FN deg(mfb1)
PRINT "MITTL. FEHLER c1...: ";FN deg(mfc1)
PRINT "FEHLERQUADRATSUMME.:";vv2;"/ ";vv3
REPEAT
UNTIL UPPER$(INKEY$) = " "
END
PROCEDURE reduk
to = (jd - jde) / 36525
t1 = (jde - jd) / 36525
w1 = RAD(((2306.218 + 1.3966 * to - 0.00014 * to * to) * t1 + (0.3019 - 0.00035 * to) * t1 * t1 + 0.018 * t1 * t1 * t1) / 3600)
w2 = RAD(((2306.218 + 1.3966 * to - 0.00014 * to * to) * t1 + (1.095 + 0.00007 * to) * t1 * t1 + 0.0182 * t1 * t1 * t1) / 3600)
w3 = RAD(((2004.311 - 0.8533 * to - 0.00022 * to * to) * t1 + (-0.4267 - 0.00022 * to) * t1 * t1 - 0.04183 * t1 * t1 * t1) / 3600)
REM -----------------------------------------
jj = SIN(w3) * SIN(j) * SIN(n + w1) + COS(w3) * COS(j)
IF ABS(jj) >= 1 THEN
    jj = 0.999999999 * SGN(jj)
ENDIF
jj = ACOS(jj)
x = (SIN(j) * COS(n + w1)) / SIN(jj)
y = (COS(w3) * SIN(j) * SIN(n + w1) - SIN(w3) * COS(j)) / SIN(jj)
IF x <= -0.9999999999 THEN
    nn = FN r(PI + w2)
ELSE
    nn = FN r(ATN(y / (1 + x)) * 2 + w2)
ENDIF
RETURN
PROCEDURE elim
FOR j = 1 TO n - 1     //GAUSS ELIMINATION
    nr = j
    no = ABS(p(j,j))
    FOR i = j + 1 TO n    //ZEILENPIVOT
      noo = ABS(p(i,j))
      EXIT IF (noo - no) < 0
      no = noo
      nr = i
    NEXT i
    IF nr = j THEN
      GOTO jum2
    ENDIF
    FOR i = j TO m + 1
      no = p(nr,i)
      p(nr,i) = p(j,i)
      p(j,i) = no
    NEXT i
    jum2:
    FOR i = j + 1 TO m + 1    //ELIMINATION
      p(j,i) = p(j,i) / p(j,j)
    NEXT i
    FOR i = j + 1 TO n
      FOR k = j + 1 TO m + 1
        p(i,k) = p(i,k) - p(j,k) * p(i,j)
      NEXT k
    NEXT i
NEXT j
FOR i = 1 TO (m + 1) - n
    ko(n,i) = p(n,n + i) / p(n,n) //RÜCKSUBSTITUTION
    FOR k = 1 TO n - 1
      l = n - k
      ko(l,i) = p(l,n + i)
      FOR s = l + 1 TO n
        ko(l,i) = ko(l,i) - p(l,s) * ko(s,i)
      NEXT s
    NEXT k
NEXT i
RETURN

Die zuvor in erster Näherung erhaltenen Ausgangselemente der Laplaceschen Ebene (K, c ,ho) verbessert man durch folg. Bedingungsgleichung [Verbesserung der Ausgangselemente: dK, dNf, dJf, d(c sin ho), d(c cos ho)], woraus auch der Winkel (Jf) und der Schnittpunkt (Nf) der festen Ebene auf der Erdäquatorebene folgt.

Bedingungsgleichung:

A) sin(Jf)*dNf+d[c sin(ho)]*cos(K*t)-d[c cos(h o)]*sin(K*t)-t*sin(c)*cos(ho-K*t)*dK=(N-Nj)*sin(Jf)
B) dJf+d[ sin(ho)]*sin(K*t)+d[ cos(ho)]*cos(K*t)+t*sin(c*sin(ho-K*t)*dK=J-Jf

Bedingungsgleichung A:
a(n)=cos(K*t)
b(n)=-sin(K*t)
c(n)=-t*sin(c)*cos(ho-K*t)
d(n)=(N-Nf)*sin(Jf)

[px] = S g(n)*a(n), 1,2,3,...,n
[py] = S g(n)*b(n), 1,2,3,...,n
[pz] = S g(n)*c(n), 1,2,3,...,n
[pxx] = S g(n)*a(n)*a(n), 1,2,3,...,n
[pyy] = S g(n)*b(n)*b(n), 1,2,3,...,n
[pzz] = S g(n)*c(n)*c(n), 1,2,3,...,n
[pxy] = S g(n)*a(n)*b(n), 1,2,3,...,n
[pxz] = S g(n)*a(n)*c(n), 1,2,3,...,n
[pyz] = S g(n)*b(n)*c(n), 1,2,3,...,n
[pd] =   S g(n)*d(n), 1,2,3,...,n
[pxd] = S g(n)*a(n)*d(n), 1,2,3,...,n
[pyd] = S g(n)*b(n)*d(n), 1,2,3,...,n
[pzd] = S g(n)*c(n)*d(n), 1,2,3,...,n

Normalgleichung für vier Unbekannte der Bedingungsgleichung A:
[p] dNf+[px] d(c sin h)+[py] d(c cos h)+[pz ] dK=[pd]
[px] dNf+[pxx] d(c sin h)+[pxy] d(c cos h)+[pxz] dK=[pxd]
[py] dNF+[pxy] d( c sin h )+[pyy] d(c cos h)+[pyz] dK=[pyd]
[pz] dNf+[pxz] d(c sin h)+[pyz] d(c cos h)+[pzz] dK=[pzd]

Bedingungsgleichung B:
a1(n)=sin(K*t)
b1(n)=cos(K*t)
c1(n)=t*sin(c)*sin(ho-K*t)
d1(n)=J-Jf

[px1] = S g(n)*a1(n), 1,2,3,...,n
[py1] = S g(i)*b1(n), 1,2,3,...,n
[pz1] = S g(i)*c1(n), 1,2,3,...,n
[pxx1] = S g(n)*a1(n)*a1(n), 1,2,3,...,n
[pyy1] = S g(i)*b1(n)*b1(n), 1,2,3,...,n
[pzz1] = S g(i)*c1(n)*c1(n), 1,2,3,...,n
[pxy1] = S g(n)*a1(n)*b1(n), 1,2,3,...,n
[pxz1] = S g(i)*a1(n)*c1(i), 1,2,3,...,n
[pyz1] = S g(i)*b1(n)*c1(i), 1,2,3,...,n
[pd1] = S g(i)*d1(n), 1,2,3,...,n
[pxd1] = S g(n)*a1(n)*d1(n), 1,2,3,...,n
[pyd1] = S g(i)*b1(n)*d1(i), 1,2,3,...,n
[pzd1] = S g(i)*c1(n)*d1(i), 1,2,3,...,n

Normalgleichung für vier Unbekannte der Bedingungsgleichung A:
[p1] dJf+[px1]    d(c sin h)+[py1]   d( c cos h )+[pz1 ] dK=[pd1]
[px1] dJf+[pxx1] d(c sin h)+[pxy1] d(c cos h)+[pxz1] dK=[pxd1]
[py1] dJF+[pxy1] d(c sin h)+[pyy1] d(c cos h)+[pyz1] dK=[pyd1]
[pz1] dJf+[pxz1] d( c sin h)+[pyz1] d(c cos h )+[pzz1] dK=[pzd1]

t in Jahren ab Epoche JDo: t=(JD-JDo)/365.25.

Beispiel.
Angenommene Ausgangselemente: Jf = 37°, Nf = 48°.
c (go) = 1.103°, h o (kn) = 16.37°, K (ko) = 159.003°.

Bedingungsgleichung B verweist auf kleinere Fehlergrenzen.
Da c 1.103° ± 0.003°, ho = 16.40° ± 50.094°, K = 159.001° ± 0.0173° kaum abweichende Werte sind, können die Ausgangswerte beibehalten werden.

Für die Länge auf der Laplaceschen Ebene (h) zwischen dem Schnittpunkt der festen Ebene auf dem Erdäquator und der Mondbahn für einen beliebigen Zeitpunkt (JD) ab der Epoche JDo erhält man: h = 16.37°-159.003° * (JD-2439501.375)/365.25.

Winkel der Laplaceschen Ebene gegen die Erdäquatorebene: Jf=37.07° ± 0.001°. Aufsteigender Knoten der Laplaceschen Ebene auf der Erdäquatorebene: Nf=47.78° ± 0.004°.

Akzeptierte Werte der Epoche JDo=2439501.375:
Jf=37.07° ± 0.001°, Äquinoktium J2000.
Nf=47.78° ± 0.004°, Äquinoktium J2000.
c =1.103° ± 0.003°.
ho = 16.37° ± 50.094°.
K=159.003° ± 0.0173°.
h =h o 16.37°-159.003°*(JD-2439501.375)/365.25, Äquinoktium J2000.

REM Programmabbruch durch Leertaste
DEFFN rad(x) = x * (PI / 180)
DEFFN deg(x) = x * (180 / PI)
DEFFN r(x) = x - INT(x / (2 * PI)) * (2 * PI)
REM TRANSFORMATION -------------------
REM J=WINKEL DER MONDBAHN AUF DER ERDÄQUATOREBENE
REM N=REKTASZENSION DES AUFST. KNOTENS DER MONDBAHN AUF DER ERDÄQUATOREBENE
REM Jf=WINKEL DER LAPLACESCHEN EBENE AUF DER ERDÄQUATOREBENE
REM Nf=REKTASZENSION DES AUFST. KNOTENS DER LAPLACESCHEN EBENE AUF DER ERDQUATOREBENE
REM go, =WINKEL DER MONDBAHN AUF DER LAPLACESCHEN EBENE
REM ko, theta o=AUF DER LAPLAC. EBENE GEMESSENE LÄNGENDIFFERENZ ZWISCHEN DEM DEM KNOTEN DER LAPLACESCHEN EBENE AUF DER ERDÄQUATOR- UND MONDBAHNEBENE
REM w = AUF DER MONDBAHNEBENE GEMESSENE LÄNGENDIFFERENZ ZWISCHEN DEM KNOTEN DER MONDBAHNEBENE AUF DER ERDQUATOREBENE U. LAPLAC. EBENE
REM -----------------------------
REM TRANSFORMATION J,N,w aus Jf,Nf,ko;(&HE9);o,g=
REM j = ACOS(COS(g) * COS(jf) - SIN(g) * SIN(jf) * COS(ko))
REM x = (COS(g) * SIN(jf) + SIN(g) * COS(jf) * COS(ko)) / SIN(j)
REM y = (SIN(g) * SIN(ko)) / SIN(j)
REM n = FN r(ATN(y / (1 + x)) * 2 + nf)
REM x = (COS(jf) * SIN(g) + SIN(jf) * COS(g) * COS(ko)) / SIN(j)
REM y = (SIN(jf) * SIN(ko)) / SIN(j)
REM F x <= -0.9999999 THEN
REM w = PI
REM ELSE
REM w = FN r(ATN(y / (1 + x)) * 2)
REM ENDIF
REM UMKEHRUNG: Jf,Nf,ko (= theta o) folgt aus J,N,w,g
REM jf = ACOS(COS(g) * COS(j) + SIN(g) * SIN(j) * COS(w))
REM x = (COS(g) * SIN(j) - SIN(g) * COS(j) * COS(w)) / SIN(jf)
REM y = (SIN(g) * SIN(w)) / SIN(jf)
REM nf = FN r(n - ATN(y / (1 + x)) * 2)
REM x = (-COS(j) * SIN(g) + SIN(j) * COS(g) * COS(w)) / SIN(jf)
REM y = (SIN(j) * SIN(w)) / SIN(jf)
REM IF x <= -0.99999999 THEN
REM ko = PI
REM ELSE
REM ko = FN r(ATN(y / (1 + x)) * 2)
REM ENDIF
REM TRANSFORMATION g= ,ko=;(&HE9);o,w aus J,N,Jf,Nf
REM g = ACOS(COS(j) * COS(jf) + SIN(j) * SIN(jf) * COS(n - nf))
REM x = (-COS(j) * SIN(jf) + SIN(j) * COS(jf) * COS(n - nf)) / SIN(g)
REM y = (SIN(j) * SIN(n - nf)) / SIN(g)
REM ko = FN r(ATN(y / (1 + x)) * 2)
REM x = (SIN(j) * COS(jf) - COS(j) * SIN(jf) * COS(n - nf)) / SIN(g)
REM y = (SIN(jf) * SIN(n - nf)) / SIN(g)
REM w = FN r(ATN(y / (1 + x)) * 2)
REM LAPLACESCHE EBENE - VERBESSERUNG DER AUSGANGSWERTE PHOBOS
DIM p(30,30),ko(30,30),a(20),b(20),c(20),d(20),a1(20),b1(20),c1(20),d1(20),g(20)
REM ----------------------
n1 = 20 //EINTRAG ANZAHL WERTPAARE J,N
REM AUSGANSELEMENTE
jf = FN rad(37)      //EINTRAG Jf UNGEFHR
nf = FN rad(48)      //EINTRAG Nf UNGEFHR
go = FN rad(1.103)   //EINTRAG = 1.03 GRAD (s. S. 102)
kn = FN rad(16.37)   //EINTRAG =16.37 GRAD (s. S. 102)
ko = FN rad(159.003) //EINTRAG K
RESTORE dat
FOR i = 1 TO n1
REM EPOCHE JDo=J2000=JD 2451545------------------------
jde = 2451545
READ jd,j,n,g
g(i) = g
j = FN rad(j)
n = FN rad(n)
GOSUB reduk
j = jj
n = nn
REM EPOCHE JDo=2439501.375=10.1.1967, 21h TT.
t = (jd - 2439501.375) / 365.25
REM BEDINGUNGSGLEICHUNG A
a(i) = COS(ko * t)
b(i) = -SIN(ko * t)
c(i) = -t * SIN(go) * COS(kn - ko * t)
d(i) = (n - nf) * SIN(jf)
REM BEDINGUNGSGLEICHUNG B
a1(i) = SIN(ko * t)
b1(i) = COS(ko * t)
c1(i) = t * SIN(go) * SIN(kn - ko * t)
d1(i) = j - jf
NEXT i
dat: //DATA-DATEI JD,J,N, ÄQUINOKTIUM DES DATUMS JD
DATA 2439856.5,36.38022,46.29004,1
DATA 2439916.5,36.135227,47.02587,1
DATA 2439977.5,36.11353,47.88572,1
DATA 2440038.5,36.325607,48.66509,1
DATA 2440100.5,36.731237,49.197308,1
DATA 2440161.5,37.228109,49.35683,1
DATA 2440222.5,37.715103,49.13287,1
DATA 2440281.5,38.081532,48.605161,1
DATA 2440342.5,38.277769,47.851407,1
DATA 2440403.5,38.25006,47.03588,1
DATA 2440465.5,37.99597,46.3057832,1
DATA 2440526.5,37.57985,45.842756,1
DATA 2440587.5,37.077729,45.725376,1
DATA 2440952.5,36.9196,49.332093,1
DATA 2441865.5,37.61435,49.239938,1
DATA 2442747.5,37.97487,48.820197,1
DATA 2443478.5,37.26688,49.423955,1
DATA 2444025.5,36.1684,46.77792,1
DATA 2444786.5,36.52682,46.109697,1
DATA 2445578.5,36.7814,45.9154,1
REM ----------------------------
REM METHODE DER KLEINSTEN QUADRATE
REM AKKUMULATION
x = 0
y = 0
z = 0
xx = 0
yy = 0
zz = 0
xy = 0
xz = 0
yz = 0
d = 0
xd = 0
yd = 0
zd = 0
x1 = 0
y1 = 0
z1 = 0
xx1 = 0
yy1 = 0
zz1 = 0
xy1 = 0
xz1 = 0
yz1 = 0
d1 = 0
xd1 = 0
yd1 = 0
zd1 = 0
dd = 0
dd1 = 0
p = 0
FOR i = 1 TO n1
REM BEDINGUNGSGLEICHUNG A
p = p + g(i)
x = x + g(i) * a(i)
y = y + g(i) * b(i)
z = z + g(i) * c(i)
xx = xx + g(i) * a(i) * a(i)
yy = yy + g(i) * b(i) * b(i)
zz = zz + g(i) * c(i) * c(i)
xy = xy + g(i) * a(i) * b(i)
xz = xz + g(i) * a(i) * c(i)
yz = yz + g(i) * b(i) * c(i)
d = d + g(i) * d(i)
dd = dd + g(i) * d(i) * d(i)
xd = xd + g(i) * a(i) * d(i)
yd = yd + g(i) * b(i) * d(i)
zd = zd + g(i) * c(i) * d(i)
REM BEDINGUNGSGLEICHUNG B
x1 = x1 + g(i) * a1(i)
y1 = y1 + g(i) * b1(i)
z1 = z1 + g(i) * c1(i)
xx1 = xx1 + g(i) * a1(i) * a1(i)
yy1 = yy1 + g(i) * b1(i) * b1(i)
zz1 = zz1 + g(i) * c1(i) * c1(i)
xy1 = xy1 + g(i) * a1(i) * b1(i)
xz1 = xz1 + g(i) * a1(i) * c1(i)
yz1 = yz1 + g(i) * b1(i) * c1(i)
d1 = d1 + g(i) * d1(i)
dd1 = dd1 + g(i) * d1(i) * d1(i)
xd1 = xd1 + g(i) * a1(i) * d1(i)
yd1 = yd1 + g(i) * b1(i) * d1(i)
zd1 = zd1 + g(i) * c1(i) * d1(i)
NEXT i
REM BEDINGUNGSGLEICHUNG A --------
REM Normalgleichung fr vier Unbekannte Bedingungsgleichung A:
REM [p] dNf+[px] d( sin(;(&HE9);)+[py] d( cos(theta)+[pz] dK=[pd]
REM [px] dNf+[pxx] d( sin(;(&HE9);)+[pxy] d( cos(theta)+[pxz] dK=[pxd]
REM [py] dNF+[pxy] d( sin(;(&HE9);)+[pyy] d( cos(theta)+[pyz] dK=[pyd]
REM [pz] dNf+[pxz] d( sin(;(&HE9);)+[pyz] d( cos(theta)+[pzz] dK=[pzd]
resi = 5 //EINTRAG ANZAHL RESIDUEN
m = 4    //EINTRAG ANZAHL GLEICHUNGEN
n = 4    //EINTRAG ANZAHL UNBEKANNTE
REM ----------------
m = m + resi
p(1,1) = p
p(1,2) = x
p(1,3) = y
p(1,4) = z
p(1,5) = d   //1. RESIDUENVEKTOR
p(1,6) = 1
p(1,7) = 0
p(1,8) = 0
p(1,9) = 0
p(2,1) = x
p(2,2) = xx
p(2,3) = xy
p(2,4) = xz
p(2,5) = xd    //1. RESIDUENVEKTOR
p(2,6) = 0
p(2,7) = 1
p(2,8) = 0
p(2,9) = 0
p(3,1) = y
p(3,2) = xy
p(3,3) = yy
p(3,4) = yz
p(3,5) = yd   //1. RESIDUENVEKTOR
p(3,6) = 0
p(3,7) = 0
p(3,8) = 1
p(3,9) = 0
p(4,1) = z
p(4,2) = xz
p(4,3) = yz
p(4,4) = zz
p(4,5) = zd   //1. RESIDUENVEKTOR
p(4,6) = 0
p(4,7) = 0
p(4,8) = 0
p(4,9) = 1
GOSUB elim
REM FEHLERQUADRATSUMME
vv = dd - d * ko(1,1) - xd * ko(2,1) - yd * ko(3,1) - zd * ko(4,1)
REM FEHLERQUADRATSUMME BEOBACHTUNG MINUS RECHNUNG
vv1 = 0
FOR i = 1 TO n1
vv1 = vv1 + (d(i) - (ko(1,1) + a(i) * ko(2,1) + b(i) * ko(3,1) + c(i) * ko(4,1))) ^ 2 * g(i)
NEXT i
s = ABS(vv / (n1 - 4))
s1 = SQR(s)
mfa = SQR(ABS(s * ko(1,2)))
mfb = SQR(ABS(s * ko(2,3)))
mfc = SQR(ABS(s * ko(3,4)))
mfd = SQR(ABS(s * ko(4,5)))
IF n1 <> p THEN
s1 = s1 / SQR(p)
ENDIF
a = ko(1,1)
b = ko(2,1)
c = ko(3,1)
d = ko(4,1)
REM BEDINGUNGSGLEICHUNG B ----------
resi = 5 //EINTRAG ANZAHL RESIDUEN
m = 4    //EINTRAG ANZAHL GLEICHUNGEN
n = 4    //EINTRAG ANZAHL UNBEKANNTE
REM ---------------------
m = m + resi
p(1,1) = p
p(1,2) = x1
p(1,3) = y1
p(1,4) = z1
p(1,5) = d1   //1. RESIDUENVEKTOR
p(1,6) = 1
p(1,7) = 0
p(1,8) = 0
p(1,9) = 0
p(2,1) = x1
p(2,2) = xx1
p(2,3) = xy1
p(2,4) = xz1
p(2,5) = xd1   //1. RESIDUENVEKTOR
p(2,6) = 0
p(2,7) = 1
p(2,8) = 0
p(2,9) = 0
p(3,1) = y1
p(3,2) = xy1
p(3,3) = yy1
p(3,4) = yz1
p(3,5) = yd1 //1. RESIDUENVEKTOR
p(3,6) = 0
p(3,7) = 0
p(3,8) = 1
p(3,9) = 0
p(4,1) = z1
p(4,2) = xz1
p(4,3) = yz1
p(4,4) = zz1
p(4,5) = zd1 //1. RESIDUENVEKTOR
p(4,6) = 0
p(4,7) = 0
p(4,8) = 0
p(4,9) = 1
GOSUB elim
REM FEHLERQUADRATSUMME
vv2 = dd1 - d1 * ko(1,1) - xd1 * ko(2,1) - yd1 * ko(3,1) - zd1 * ko(4,1)
REM FEHLERQUADRATSUMME BEOBACHTUNG MINUS RECHNUNG
vv3 = 0
FOR i = 1 TO n1
vv3 = vv3 + (d1(i) - (ko(1,1) + a1(i) * ko(2,1) + b1(i) * ko(3,1) + c1(i) * ko(4,1))) ^ 2 * g(i)
NEXT i
s = ABS(vv2 / (n1 - 4))
s2 = SQR(s)
mfa1 = SQR(ABS(s * ko(1,2)))
mfb1 = SQR(ABS(s * ko(2,3)))
mfc1 = SQR(ABS(s * ko(3,4)))
mfd1 = SQR(ABS(s * ko(4,5)))
IF n1 <> p THEN
s2 = s2 / SQR(p)
ENDIF
a1 = ko(1,1)
b1 = ko(2,1)
c1 = ko(3,1)
d1 = ko(4,1)
REM ---------------------
jf = jf + a1
nf = FN r(nf + a / SIN(jf))
k = FN r(ko + d)
k1 = FN r(ko + d1)
REM -------------
y = b
x = c
g1 = SQR(b * b + c * c)
y1 = y / g1
x1 = x / g1
kno = FN r(ATN(y1 / (1 + x1)) * 2)
y = b1
x = c1
g2 = SQR(b1 * b1 + c1 * c1)
y1 = y / g2
x1 = x / g2
kn1 = FN r(ATN(y1 / (1 + x1)) * 2)
fg1 = SQR(mfb * mfb + mfc * mfc)
y1 = mfb / fg1
x1 = mfc / fg1
fkn = FN r(ATN(y1 / (1 + x1)) * 2)
fg2 = SQR(mfb1 * mfb1 + mfc1 * mfc1)
y1 = mfb1 / fg2
x1 = mfc1 / fg2
fkn2 = FN r(ATN(y1 / (1 + x1)) * 2)
REM -------------------------------
PRINT "WINKEL Jf J2000 ";FN deg(jf);" +/- ";FN deg(mfa1);" GRAD"
PRINT "KNOTEN Nf J2000 ";FN deg(nf);" +/- ";FN deg(mfa) / SIN(jf);" GRAD'"
PRINT "BEDINGUNGSGLEICHUNG A:"
PRINT "KNOTEN K ";FN deg(k);" +/- ";FN deg(mfd);" GRAD"
PRINT "WINKEL LAPLAC. EBENE/MONDBAHNEBENE: ";FN deg(g1);" +/- ";FN deg(fg1);" GRAD"
PRINT "LÄNGE: ";FN deg(kno);" +/- ";FN deg(fkn);" GRAD"
PRINT "MITTL. FEHLER EINER EINZELMESSUNG: ";FN deg(s1);" GRAD"
PRINT "FEHLERQUADRATSUMME.: ";vv;" / ";vv1
PRINT "BEDGINGUNGSGLEICHUNG B:"
PRINT "KNOTEN K ";FN deg(k1);" +/- ";FN deg(mfd1);" GRAD"
PRINT "WINKEL LAPLAC. EBENE/MONDBAHNEBENE: ";FN deg(g2);" +/- ";FN deg(fg2);" GRAD"
PRINT "LÄNGE: ";FN deg(kn1);" +/- ";FN deg(fkn2);" GRAD"
PRINT "MITTL. FEHLER EINER EINZELMESSUNG: ";FN deg(s2);" GRAD"
PRINT "FEHLERQUADRATSUMME.: ";vv2;" / ";vv3
REPEAT
UNTIL UPPER$(INKEY$) = " "
END
PROCEDURE reduk
to = (jd - jde) / 36525
t1 = (jde - jd) / 36525
w1 = RAD(((2306.218 + 1.3966 * to - 0.00014 * to * to) * t1 + (0.3019 - 0.00035 * to) * t1 * t1 + 0.018 * t1 * t1 * t1) / 3600)
w2 = RAD(((2306.218 + 1.3966 * to - 0.00014 * to * to) * t1 + (1.095 + 0.00007 * to) * t1 * t1 + 0.0182 * t1 * t1 * t1) / 3600)
w3 = RAD(((2004.311 - 0.8533 * to - 0.00022 * to * to) * t1 + (-0.4267 - 0.00022 * to) * t1 * t1 - 0.04183 * t1 * t1 * t1) / 3600)
REM -----------------------------------------
jj = SIN(w3) * SIN(j) * SIN(n + w1) + COS(w3) * COS(j)
IF ABS(jj) >= 1 THEN
    jj = 0.999999999 * SGN(jj)
ENDIF
jj = ACOS(jj)
x = (SIN(j) * COS(n + w1)) / SIN(jj)
y = (COS(w3) * SIN(j) * SIN(n + w1) - SIN(w3) * COS(j)) / SIN(jj)
IF x <= -0.9999999999 THEN
    nn = FN r(PI + w2)
ELSE
    nn = FN r(ATN(y / (1 + x)) * 2 + w2)
ENDIF
RETURN
PROCEDURE elim
FOR j = 1 TO n - 1     //GAUSS ELIMINATION
    nr = j
    no = ABS(p(j,j))
    FOR i = j + 1 TO n    //ZEILENPIVOT
      noo = ABS(p(i,j))
      EXIT IF (noo - no) < 0
      no = noo
      nr = i
    NEXT i
    IF nr = j THEN
      GOTO jum2
    ENDIF
    FOR i = j TO m + 1
      no = p(nr,i)
      p(nr,i) = p(j,i)
      p(j,i) = no
    NEXT i
    jum2:
    FOR i = j + 1 TO m + 1    //ELIMINATION
      p(j,i) = p(j,i) / p(j,j)
    NEXT i
    FOR i = j + 1 TO n
      FOR k = j + 1 TO m + 1
        p(i,k) = p(i,k) - p(j,k) * p(i,j)
      NEXT k
    NEXT i
NEXT j
FOR i = 1 TO (m + 1) - n
    ko(n,i) = p(n,n + i) / p(n,n) //RÜCKSUBSTITUTION
    FOR k = 1 TO n - 1
      l = n - k
      ko(l,i) = p(l,n + i)
      FOR s = l + 1 TO n
        ko(l,i) = ko(l,i) - p(l,s) * ko(s,i)
      NEXT s
    NEXT k
NEXT i
RETURN

REM DEIMOS ERDÄQUATORBENE
DIM p(10,10),ko(10,10),a(40),b(40),c(40),d(40),a1(40),b1(40),c1(40),g(40)
DEFFN rad(x) = x * (PI / 180)
DEFFN deg(x) = x * (180 / PI)
DEFFN r(x) = x - INT(x / (2 * PI)) * (2 * PI)
REM ----------------------
n1 = 27 //EINTRAG ANZAHL WERTPAARE J,N
jb = PI / 2 - FN rad(52) //EINTRAG UNGEFÄHRE DEKLINATION NORDPOL MARS ÄQUINOKTIUM J2000
nb = FN r(RAD(317) + PI / 2) //EINTRAG UNGEFÄHRE REKTASZENSION NORDPOL MARS ÄQUINOKTIUM J2000
k1 = RAD(0.2673304)
k = RAD(6.577422)
RESTORE dat
FOR i = 1 TO n1
REM EPOCHE JDo=J2000=JD 2451545------------------------
jdo = 2451545
READ jd,j,n,g
g(i) = g
j = FN rad(j)
n = FN rad(n)
GOSUB reduk
j3 = jj
n3 = nn
t = (jd - 2441266.5) / 365.25
REM BEDINGUNSGLEICHUNG A
a(i) = COS(k * t)
b(i) = -SIN(k * t)
c(i) = (n3 - nb) * SIN(jb)
REM BEDINGUNSGLEICHUNG B
a1(i) = SIN(k * t)
b1(i) = COS(k * t)
c1(i) = j3 - jb + (k1 / k) * SIN(jb) * COS(jb)
NEXT i
dat: //DATA-DATEI JD,J,N,¥QUINOKTIUM DES DATUMS JD
DATA 2440952.5,35.8433,43.73345,1
DATA 2441865.5,35.390465,44.19275,1
DATA 2442747.5,35.04149,44.834776,1
DATA 2443478.5,34.839039,45.480678,1
DATA 2444025.5,34.746705,46.00889,1
DATA 2444786.5,37.70838,46.7726,1
DATA 2445578.5,34.7812,47.555675,1
DATA 2446309.5,34.9457151,48.222398,1
DATA 2447039.5,35.19279327,48.79599,1
DATA 2447892.5,35.5663778,49.309379,1
DATA 2448622.5,35.9375174,49.5907038,1
DATA 2449718.5,36.535842,49.7169278,1
DATA 2450449.5,36.9308034,49.6046091,1
DATA 2451179.5,37.297719,49.346842,1
DATA 2452275.5,37.7566990388,48.7253933,1
DATA 2453005.5,37.978977,48.1889985,1
DATA 2453736.5,38.1214043,47.5852507,1
DATA 2454466.5,38.1765344,46.9452509,1
DATA 2455197.5,38.14157,46.2978717,1
DATA 2457023.5,37.67842,44.852446,1
DATA 2457754.5,37.36459,44.418031,1
DATA 2458119.5,37.186939,44.24478,1
DATA 2459580.5,36.3876677,43.899106,1
DATA 2460310.5,35.97062,43.957077,1
DATA 2461041.5,35.56976,44.174942,1
DATA 2461771.5,35.20751,44.54556,1
DATA 2462502.5,34.902707,45.053875,1
DATA 2463232.5,34.673164,45.673083,1
DATA 2463963.5,34.53122,46.371183,1
DATA 2464693.5,34.485194,47.10654,1
DATA 2465424.5,34.53722,48.838264,1
DATA 2466154.5,34.683587,48.521888,1
DATA 2466885.5,34.915533,49.120511,1
DATA 2467615.5,35.21846,49.600402,1
DATA 2468346.5,35.575445,49.939292,1
DATA 2469076.5,35.965231,50.1221802,1
DATA 2469807.5,36.36735,50.1449378,1
REM ----------------------------
REM METHODE DER KLEINSTEN QUADRATE
REM AKKUMULATION
x = 0
xx = 0
xy = 0
y = 0
yy = 0
yz = 0
d = 0
dx = 0
dy = 0
x1 = 0
xx1 = 0
xy1 = 0
yy1 = 0
y1 = 0
yz1 = 0
d1 = 0
dx1 = 0
dy1 = 0
dd = 0
dd1 = 0
p = 0
FOR i = 1 TO n1
REM BEDINGUNGSGLEICHUNG A
p = p + g(i)
x = x + a(i) * g(i)
xx = xx + a(i) * a(i) * g(i)
xy = xy + a(i) * b(i) * g(i)
y = y + b(i) * g(i)
yy = yy + b(i) * b(i) * g(i)
d = d + c(i) * g(i)
dx = dx + a(i) * c(i) * g(i)
dy = dy + b(i) * c(i) * g(i)
dd = dd + c(i) * c(i) * g(i)
REM BEDINGUNGSGLEICHUNG B
x1 = x1 + a1(i) * g(i)
xx1 = xx1 + a1(i) * a1(i) * g(i)
xy1 = xy1 + a1(i) * b1(i) * g(i)
y1 = y1 + b1(i) * g(i)
yy1 = yy1 + b1(i) * b1(i) * g(i)
d1 = d1 + c1(i) * g(i)
dx1 = dx1 + a1(i) * c1(i) * g(i)
dy1 = dy1 + b1(i) * c1(i) * g(i)
dd1 = dd1 + c1(i) * c1(i) * g(i)
NEXT i
REM BEDINGUNGSGLEICHUNG A --------
resi = 4 //EINTRAG ANZAHL RESIDUEN
m = 3    //EINTRAG ANZAHL GLEICHUNGEN
n = 3    //EINTRAG ANZAHL UNBEKANNTE
REM ----------------
m = m + resi
p(1,1) = p
p(1,2) = x
p(1,3) = y
p(1,4) = d   //1. RESIDUENVEKTOR
p(1,5) = 1   //2. RESIDUENVEKTOR
p(1,6) = 0
p(1,7) = 0
p(2,1) = x
p(2,2) = xx
p(2,3) = xy
p(2,4) = dx //1. RESIDUIENVEKTOR
p(2,5) = 0    //2. RESIDUENVEKTOR
p(2,6) = 1
p(2,7) = 0
p(3,1) = y
p(3,2) = xy
p(3,3) = yy
p(3,4) = dy //1. RESIDUIENVEKTOR
p(3,5) = 0   //2. RESIDUENVEKTOR
p(3,6) = 0
p(3,7) = 1
GOSUB elim
c = ko(1,1)
a = ko(2,1)
b = ko(3,1)
REM FEHLERQUADRATSUMME
vv = dd - d * ko(1,1) - dx * ko(2,1) - dy * ko(3,1)
REM FEHLERQUADRATSUMME BEOBACHTUNG MINUS RECHNUNG
vv1 = 0
FOR i = 1 TO n1
vv1 = vv1 + (c(i) - (ko(1,1) + a(i) * ko(2,1) + b(i) * ko(3,1))) ^ 2
NEXT i
s = ABS(vv / (n1 - 3))
s1 = SQR(s)
mfa = SQR(ABS(s * ko(2,3)))
mfb = SQR(ABS(s * ko(3,4)))
mfc = SQR(ABS(s * ko(1,2)))
IF n1 <> p THEN
s1 = s1 / SQR(p)
ENDIF
REM -------------
REM BEDINGUNGSGLEICHUNG B ----------
resi = 4 //EINTRAG ANZAHL RESIDUEN
m = 3    //EINTRAG ANZAHL GLEICHUNGEN
n = 3    //EINTRAG ANZAHL UNBEKANNTE
REM ----------------
m = m + resi
p(1,1) = p
p(1,2) = x1
p(1,3) = y1
p(1,4) = d1 //1. RESIDUENVEKTOR
p(1,5) = 1   //2. RESIDUENVEKTOR
p(1,6) = 0
p(1,7) = 0
p(2,1) = x1
p(2,2) = xx1
p(2,3) = xy1
p(2,4) = dx1 //1. RESIDUIENVEKTOR
p(2,5) = 0    //2. RESIDUENVEKTOR
p(2,6) = 1
p(2,7) = 0
p(3,1) = y1
p(3,2) = xy1
p(3,3) = yy1
p(3,4) = dy1 //1. RESIDUIENVEKTOR
p(3,5) = 0   //2. RESIDUENVEKTOR
p(3,6) = 0
p(3,7) = 1
GOSUB elim
c1 = ko(1,1)
a1 = ko(2,1)
b1 = ko(3,1)
REM FEHLERQUADRATSUMME
vv2 = dd1 - d1 * ko(1,1) - dx1 * ko(2,1) - dy1 * ko(3,1)
REM FEHLERQUADRATSUMME BEOBACHTUNG MINUS RECHNUNG
vv3 = 0
FOR i = 1 TO n1
vv3 = vv3 + (c1(i) - (ko(1,1) + a1(i) * ko(2,1) + b1(i) * ko(3,1))) ^ 2
NEXT i
s2 = ABS(vv2 / (n1 - 3))
s3 = SQR(s2)
mfa1 = SQR(ABS(s2 * ko(2,3)))
mfb1 = SQR(ABS(s2 * ko(3,4)))
mfc1 = SQR(ABS(s2 * ko(1,2)))
IF n1 <> p THEN
s3 = s3 / SQR(p)
ENDIF
REM ---------------------
j1 = jb + c1
n1 = nb + c / SIN(j1)
dp1 = PI / 2 - j1
ap1 = FN r(n1 - PI / 2)
REM ------------------
g = SQR(a * a + b * b)
ao = a / g
bo = b / g
ko = FN r(ATN(ao / (1 + bo)) * 2)
g1 = SQR(a1 * a1 + b1 * b1)
ao = a1 / g1
bo = b1 / g1
ko1 = FN r(ATN(ao / (1 + bo)) * 2)
mfg = SQR(mfa * mfa + mfb * mfb)
fa = mfa / mfg
fb = mfb / mfg
fko = FN r(ATN(fa / (1 + fb)) * 2)
mfg1 = SQR(mfa1 * mfa1 + mfb1 * mfb1)
fa = mfa1 / mfg1
fb = mfb1 / mfg1
fko1 = FN r(ATN(fa / (1 + fb)) * 2)
REM -------------------------------
PRINT "BEDINGUNGSGLEICHUNG A:"
PRINT "DEKLINATION   NORDPOL MARS J2000: ";FN deg(dp1);" +/-";FN deg(mfc1);" GRAD"
PRINT "REKTASZENSION NORDPOL MARS J2000: ";FN deg(ap1);" +/-";FN deg(mfc) / SIN(j1);" GRAD"
PRINT "WINKEL LAPLAC. EBENE/MONDBAHNEBENE: ";FN deg(g);" +/-";FN deg(mfg);" GRAD"
PRINT "KNOTEN LAPLAC. EBENE/MONDBAHNEBENE: ";FN deg(ko);" +/-";FN deg(fko);" GRAD"
PRINT "KOEFFIZIENT a ";FN deg(a);" GRAD"
PRINT "KOEFFIZIENT b ";FN deg(b);" GRAD"
PRINT "KOEFFIZIENT c ";FN deg(c);" GRAD"
PRINT "MITTL. FEHLER EINER EINZELMESSUNG: ";FN deg(s1)
PRINT "MITTL. FEHLER a....: ";FN deg(mfa)
PRINT "MITTL. FEHLER b....: ";FN deg(mfb)
PRINT "MITTL. FEHLER c....: ";FN deg(mfc) / SIN(j1)
PRINT "FEHLERQUADRATSUMME.: ";vv;" / ";vv1
PRINT "BEDGINGUNGSGLEICHUNG B:"
PRINT "WINKEL LAPLACE. EBENE/MONDBAHNEBENE: ";FN deg(g1);" +/-";FN deg(mfg1);" GRAD"
PRINT "KNOTEN LAPLACE. EBENE/MONDBAHNEBENE: ";FN deg(ko1);" +/-";FN deg(fko1);" GRAD"
PRINT "KOEFFIZIENT a1 ";FN deg(a1);" GRAD"
PRINT "KOEFFIZIENT b1 ";FN deg(b1);" GRAD"
PRINT "KOEFFIZIENT c1 ";FN deg(c1);" GRAD"
PRINT "MITTL. FEHLER EINER EINZELMESSUNG: ";FN deg(s2)
PRINT "MITTL. FEHLER a1...: ";FN deg(mfa1)
PRINT "MITTL. FEHLER b1...: ";FN deg(mfb1)
PRINT "MITTL. FEHLER c1...: ";FN deg(mfc1)
PRINT "FEHLERQUADRATSUMME.: ";vv2;" / ";vv3
REPEAT
UNTIL UPPER$(INKEY$) = " "
END
PROCEDURE reduk
to = (jd - 2451545) / 36525
t1 = (jdo - jd) / 36525
w1 = RAD(((2306.218 + 1.3966 * to - 0.00014 * to * to) * t1 + (0.3019 - 0.00035 * to) * t1 * t1 + 0.018 * t1 * t1 * t1) / 3600)
w2 = RAD(((2306.218 + 1.3966 * to - 0.00014 * to * to) * t1 + (1.095 + 0.00007 * to) * t1 * t1 + 0.0182 * t1 * t1 * t1) / 3600)
w3 = RAD(((2004.311 - 0.8533 * to - 0.00022 * to * to) * t1 + (-0.4267 - 0.00022 * to) * t1 * t1 - 0.04183 * t1 * t1 * t1) / 3600)
REM -----------------------------------------
jj = SIN(w3) * SIN(j) * SIN(n + w1) + COS(w3) * COS(j)
IF ABS(jj) >= 1 THEN
    jj = 0.999999999 * SGN(jj)
ENDIF
jj = ACOS(jj)
x = (SIN(j) * COS(n + w1)) / SIN(jj)
y = (COS(w3) * SIN(j) * SIN(n + w1) - SIN(w3) * COS(j)) / SIN(jj)
IF x <= -0.9999999999 THEN
    nn = FN r(PI + w2)
ELSE
    nn = FN r(ATN(y / (1 + x)) * 2 + w2)
ENDIF
RETURN
PROCEDURE elim
FOR j = 1 TO n - 1     //GAUSS ELIMINATION
    nr = j
    no = ABS(p(j,j))
    FOR i = j + 1 TO n    //ZEILENPIVOT
      noo = ABS(p(i,j))
      EXIT IF (noo - no) < 0
      no = noo
      nr = i
    NEXT i
    IF nr = j THEN
      GOTO jum1
    ENDIF
    FOR i = j TO m + 1
      no = p(nr,i)
      p(nr,i) = p(j,i)
      p(j,i) = no
    NEXT i
    jum1:
    FOR i = j + 1 TO m + 1   //ELIMINATION
      p(j,i) = p(j,i) / p(j,j)
    NEXT i
    FOR i = j + 1 TO n
      FOR k = j + 1 TO m + 1
        p(i,k) = p(i,k) - p(j,k) * p(i,j)
      NEXT k
    NEXT i
NEXT j
FOR i = 1 TO (m + 1) - n
    ko(n,i) = p(n,n + i) / p(n,n) //RÜCKSUBSTITUTION
    FOR k = 1 TO n - 1
      l = n - k
      ko(l,i) = p(l,n + i)
      FOR s = l + 1 TO n
        ko(l,i) = ko(l,i) - p(l,s) * ko(s,i)
      NEXT s
    NEXT k
NEXT i
RETURN

Alle Rechte vorbehalten (all rights reserved) , auch die der fotomechanischen Wiedergabe und der Speicherung in elektronischen Medien, Translation usw. Dasselbe gilt für das Recht der öffentlichen Wiedergabe. Copyright © by H. Schumacher, Spaceglobe.


	*<<< Auszug*


	*Auszug >>>*

Sternbeobachter - Sterntagebuch - Produktinformation - www.spaceglobe.de