Bestimmung der dynamischen Bahnelemente Periastronzeit (T) u. Umlaufzeit (U) mittels Least-Square-Fits.

M=E-e sin(E).
M,E,v = mittlere, exentrische u. wahre Anomalie, P=Positionswinkel (Ausgangselemente e,i,
z,W); z. B. nach der Zwiers-Methode).

v=ATN(TAN(P-W)/cos(i))-z
E=ATN(SQR((1-e)/(1+e))*TAN(v/2))*2.
DEFFN r(x)=x-INT(x/(2*PI))*(2*PI)
y=TAN(P-knoten)
x=COS(i)
z=SQR(x^2+y^2)
x=x/z
y=y/z
v=FN r(ATN(y/(1+x))*2+PI-
z) //wahre Anomalie
ex=FN r(ATN(SQR((1-e)/(1+e))*TAN((v/2))*2) //exentrische Anomalie

Mit der Methode der kleinsten Quadrate können die Thiele-Innes Elemente (A,B,F,G) aus den Ausgangselementen bestimmt werden: (i=1,2,3...,n)

X(i)=cos(E(i))-e
Y(i)=sin(E(i)) cos(
v)
v = ARCSIN(e).
x(i)=A X(i)+F Y(i)
y(i)=B X(i)+G Y(i)i
x(i)=d cos(P)
y(i)=d sin(P)

Bei bekannten Koeffizienten A,B,F,G wird umgekehrt die Wertpaare X(i),Y(i) durch Ausgleichsrechnung ermittelt. Mit den Konstanten A,B,F,G,X(i),Y(i) erhält man Px,Py,Qx,Qy usw., mittels Ausgleichsrechnung dA,dB,dF,dG, de, dT, dmy, und schließlich die dynamischen Elemente my (P Umlaufzeit = 360 Grad/my), T und e.

Vorgang mit den verbesserten Ausgangselementen wiederholen. Diese Methode wird vorzugsweise angewendet, wenn die Messungen in Rechteckkoordinaten [x(i),y(i)] vorliegen.

Alternative: Für jede Beobachtung (t) rechnet man t=T+(1/my)*(E-e sin(E)) oder (2PI/U)(t-T)=M=E-e sin(E). Doppelsternsystem Sirius: T=1894.13, U=50.09 Jahre, große Halbachse 7.499''.

REM GFA24
DEFFN r(x) = x - INT(x / (2 * PI)) * (2 * PI)
DIM t(50),m(50),g(50),e(50),r(50) //BEISPIEL DOPPELSTERN SIRIUS
kn = RAD(44.57) //EINTRAG KNOTEN
w = RAD(147.27) //EINTRAG PERIASTRONARGUMENT
i = RAD(136.53) //EINTRAG INCLINATION
e = 0.592       //EINTARG EXENTRIZITÄT
anz = 31 //EINTRAG ANZAHL DATEN
FOR ii = 1 TO anz
  READ t,p,d,g //t=Periastronzeit,p=Positionswinkel,d=Distanz,g=Gewicht =1
  t(ii) = t
  g(ii) = g
  p = RAD(p) //POSITIONSWINKEL ABNEHMEND (RETROGRADE BEWEGUNG, i=IM ZWEITEN QUADRANTEN)
  y = TAN(p - kn)
  x = COS(i)
  z = SQR(x ^ 2 + y ^ 2)
  x = x / z
  y = y / z
  v = FN r(ATN(y / (1 + x)) * 2 + PI - w) //wahre Anomalie
  ex = FN r(ATN(SQR((1 - e) / (1 + e)) * TAN(v / 2)) * 2) //exentr. Anomalie
  m(ii) = ex - e * SIN(ex) //mittl. Anomalie
  e(ii) = (1 - e ^ 2) / (1 + e * COS(v))
  u = ATN(TAN(p - kn) / COS(i)) //Argument der Breite
  r(ii) = d * (COS(p - kn) / COS(u)) //Radiusvektor
NEXT ii
g = 0
FOR i = 1 TO anz
  g = g + g(i)
  xo = xo + t(i) * g(i)
  yo = yo + m(i) * g(i)
  xx = xx + t(i) * t(i) * g(i)
  yy = yy + m(i) * m(i) * g(i)
  xy = xy + t(i) * m(i) * g(i)
NEXT i
GOSUB regre
PRINT "Mittl. Anomalie........: M=a+b*T Periastronzeit"
PRINT "Korrelationskoeffizient: ";ko
PRINT "Koeffizient a: ";a;" m.F. ";mfa
PRINT "Koeffizient b: ";b;" m.F. ";mfb
PRINT "Jährl. Eigenbewegung.....: ";b;" RAD"
bg = DEG(ABS(b))
u = 2 * PI / ABS(b)
PRINT "Jährl. Eigenbewegung.....: ";bg;" GRAD"
PRINT "Umlaufzeit...............: ";u;" Jahre"
t = ABS(a / b)
PRINT "Periastronzeit...........: ";t
PRINT
g = 0
xo = 0
yo = 0
xx = 0
yy = 0
xy = 0
FOR i = 1 TO anz
  g = g + g(i)
  xo = xo + e(i) * g(i)
  yo = yo + r(i) * g(i)
  xx = xx + e(i) * e(i) * g(i)
  yy = yy + r(i) * r(i) * g(i)
  xy = xy + e(i) * r(i) * g(i)
NEXT i
GOSUB regre
PRINT "Radiusvektor...........: r=a+b*e Periastronzeit" //r=a*(1-e^2)/(1-e*COS(v))
PRINT "Korrelationskoeffizient: ";ko
PRINT "Koeffizient a: ";a;" m.F. ";mfa
PRINT "Koeffizient b: ";b;" m.F. ";mfb
PRINT "Große halbe Bahnachse a: ";b;"´´"
KEYGET HALT%
PROCEDURE regre
  REM KORRELATIONSKOEFFIZIENT (ko = ;(&HF1);1 perfekte Anpassung, ko=0 kein Zusammenhang zwischen x u. y)
  ko = (g * xy - xo * yo) / (SQR(g * xx - xo * xo) * SQR(g * yy - yo * yo))
  REM NORMALGLEICHUNGEN N*A+XO*B=YO; XO*A+XX*B=XY
  a = (yo * xx - xo * xy) / (g * xx - xo ^ 2)
  b = (g * xy - yo * xo) / (g * xx - xo ^ 2)
  REM FEHLERQUADRATSUMME (F.Q.S.)
  vv = yy - yo * a - xy * b
  s = SQR(vv / (anz - 2))
  REM MITTLERER FEHLER (M.F.) a
  mfa = s * SQR(xx / (g * xx - xo * xo))
  REM MITTLERER FEHLER b
  mfb = s * SQR(g / (g * xx - xo * xo))
RETURN
REM DOPPELSTERNSYSTEM SIRIUS - AUS DEN BAHNELEMENTEN BERECHNETE >IDEALE< DATEN 1910-1940
REM BEOBACHTUNGSZEIT t, POSITIONSWINKEL p, DISTANZ d, GEWICHT g=1
mess:
DATA 1910,90.82,8.87,1,1911,87.93,9.22,1,1912,85.25,9.55,1,1913,82.73,9.85,1,1914,80.36,
10.12,1,1915,78.10,10.36,1,1916,75.95,10.58,1,1917,73.88,10.77,1,1918,71.87,10.93,
1,1919,69.91,11.06,1,1920,68.00,11.16,1,1921,66.12,11.23
DATA 1,1922,64.25,11.27,1,1923,62.29,11.27,1,1924,60.52,11.24,1,1925,58.64,11.18,
1,1926,56.73,11.07,1,1927,54.78,10.93,1,1928,52.77,10.75,1,1929,50.68,10.52,
1,1930,48.49,10.25,1,1931,46.17,9.94,1,1932,43.69,9.57,1,1933,40.99,9.15
DATA 1,1934,38.01,8.68,1,1935,34.67,8.14,1,1936,30.82,7.54,1,1937,26.27,6.87,
1,1938,20.66,6.13,1,1939,13.43,5.32,1,1940,3.47,4.46,0,0,0,0 

Bahnverbesserung durch Differenzierung des Elementesystems

Sind die Bahnelemente eines Doppelsterns mittels einer der zuvor beschriebenen Methoden bestimmt, werden diese als Ausgangselemente für eine Bahnverbesserung herangezogen, soweit kleine Bahnteile, die nur äußerst unsichere, provisorische Elemente ergeben, eine Ausgleichung bzw. Verbesserung des Elementensystems nicht sinnlos machen; denn die Bahn ist nur dann verbesserungsfähig, wenn um den ganzen Bahnbogen einigermaßen genaue Beobachtungen vorliegen, die den strengen Anforderungen einer Ausgleichung genügen.

Bei sehr unpräzisen Messungen sind einige Elemente festzuhalten und erst in der zweiten Iteration zu korrigieren. Notfalls vernachlässigt man die Distanzen und paßt nur die Positionswinkel (O-C) an, wenn die Ellipse nicht zu schmal ist; das Beobachtungsmaterial wird dadurch zwar nicht vollkommen ausgenutzt, ergibt aber meistens eine gute Annäherung. Da die Positionswinkel meistens sicherer als kleine Distanzen sind werden oftmals nur die Winkel angepaßt (s. Progr. STARFIT).

In dem die Polarkoordinaten Positionswinkel u. Distanz für jedes Beobachtungsdatum (t) aus dem Elementensystem berechnet wird, ist die Differenz des Ephemeridenorts (C) zum beobachteten Normalort (O) zu bilden:  O-C (Observation minus Calculation = B-R Beobachtung minus Rechnung).

dp=(O-C) gemessener Positionswinkel (p) minus berechneter Ephemeriden-Positionswinkel (aus den zuvor - Thiele-Innes, Hopmann, Rabe, Zwiers bestimmten vorläufigen Elementen).

57.29578 dd/d = d,dd = Distanz u. Distanzdifferenz O-C aus den vorläufigen Elementen.

Man kann mit den kartesischen Rechteckkoordinaten (x,y) rechnen, wenn mehr photographische oder Speckle-Beobachtungen vorliegen, oder mit den Polarkoordinaten (d,p), wenn hauptsächl. visuelle Beobachtungen vorliegen.

Für jede Beobachtung erfolgt die Aufstellung der Bestimmungsgleichungen (Anzahl i=1,2,3,...,n):

A(i) dW+B(i) di+C(i) dw+F(i) dn+G(i) dm+H(i) dv=R(i)=dp(i)
a(i) da° +b(i) di+c(i) dw+d(i) dn+e(i) dm+f(i) d
v=r(i)=57.29578 dd(i)/d(i)

Das Zeitmaß »N« (Jahre) liegt stets innerhalb des auszugleichenden Bahnbogens (N=(t1-t2)/2; t1 u. t2=Anfangs- u. Endjahr des auszugleichenen Bahnstücks, Epoche ep=t1+N). c = (t-to)/N <1; to=gewählte Epoche, t=Beobachtungszeit.
Die gewählte Oskulationsepoche (to) darf nicht mit der Epoche des auf- oder absteigenden Bahnknotens zusammenfallen, da das Verfahren sonst versagt.

Bahnelemente: e=Exentrizität, a=Große Bahnhalbachse, i=Inclination, my=l=jährl. Eigenbewegung, z=w=Periastronargument, W=Konten.

e = sin(e)
sec(
v)=1/cos(arcsin(e))
d
v = 57.29578 sec( v) de
dn=my sec(
v)2 dT
my=jährl. Eigenbewegung in GRAD
da° =57.29578 da/a
da=(da°*a)/57.29578
dm = N sec(
v)2 dmy
N=Zeitmaß: 
c =(t-T)/N.
de=d
v / (57.29578*(1/cos(arcsin(e))))
dT=dn/(my sec(
v)2)
dmy = dm/(N sec(
v)2

Koeffizienten A,B,C,D,E,F (P=pw=aus den zu verbessernden Elementen berechneter Positionswinkel:

A(i)=1,
B(i)=-cos(P-
W)2 tan(v+w) sin(i)
C(i)= (cos(P-
W)2 /cos(v+w)2) cos(i)
D(i)=-n C
E(i)=n C(i) ((t-T)/N)
F(i)=m C(i)
a(i)=1
b(i)=-sin(P-
W) tan(i)
c(i)=-sin(P-
W) cos(P-W) sin(i) tan(i)
d(i)=-(n c(i) + o)
e(i)=(n C(i) + o) ((t-T)/N)
f(i)=m c(i) + s

Hilfsfunktionen:

n=sec(v) (1+e cos(v))2
m=sec(
v) sin(v) (2+e*cos(v)); v=wahre Anomalie.
s=-sec(
v) cos(v) (1+e*cos(v))
o=tan(arcsin(e)) sin(v) (1+e*cos(v))

Bei z. B. n=20 Beobachtungen erhält man auf der linken Seite die 20 x 6 Matrix A (20 Beobachtungen x 6 Unbekannte) und auf rechten Seite der 20 x 1 Residuenvektor Matrix R = R(i) = r(i).

 Aufstellung einer Normalgleichung mit 6 Unbekannten (allgemein):

A(1) + B(1)v + C(1)w + D(1)x + E(1)y + F(1)z = X(1)
A(2) + B(2)v + C(2)w + D(2)x + E(2)y + F(2)z = X(2)
usw.  .....   .....   .....   .....   .....   .....   .....

1) A(n)=1; 1,2,...,n. Transponiere Matrix A, 2) die Multiplikation der transponierten Matrix AT mit der 20 X 6 Matrix A, ergibt die linke Seite der Normalgleichung Matrix 6 X 6, 3) die Multiplikation der transponierten Matrix A T mit dem 20 X 6 Residuenvektor Matrix R, bildet die rechte Seite der Normalgleichung Matrix 6 X 6.

n-Beobachtungen bzw. Bedingungsgleichungen des überbestimmten linearen Gleichungs- systems mit 6 Unbekannten o,v,w,x,y,z:

                           A(1) + B(1)v + C(1)w + D(1)x + E(1)y + F(1)z
                           A(2) + B(2)v + C(2)w + D(2)x + E(2)y + F(2)z
                           A(3) + B(3)v + C(3)w + D(3)x + E(3)y + F(3)z
   Matrix A =        A(4) + B(4)v + C(4)w + D(4)x + E(4)y + F(4)z
                           A(5) + B(5)v + C(5)w + D(5)x + E(5)y + F(5)z
                           A(6) + B(6)v + C(6)w + D(6)x + E(6)y + F(6)z
                           usw. .....   .....   .....   .....   .....   .....
                           A(n) + B(n)v + C(n)w + D(n)x + E(n)y + F(n)z

Residuenvektor R: 

                  X(1)
                  X(2)
                  X(3)
Matrix R  = X(4)
                  X(5)
                  X(6)
                  ....
                  X(n)

1) transponiere Matrix A (Spalte Marix A in Zeile Matrix AT):

                        A(1)   + A(2)   +  A(3)   + A(4)    + A(5).....A(n)
                        B(1)v + B(2)v  +  B(3)v + B(4)v  + B(5)v ...B(n)v
Matrix AT =      C(1)w + C(2)w + C(3)w + C(4)w + C(5)w...C(n)w
                        D(1)x + D(2)x  + D(3)x + D(4)x  + D(5)x....D(n)x
                        E(1)y  + E(2)y  + E(3)y  + E(4)y  + E(5)y.....E(n)y
                        F(1)z  + F(2)z   + F(3)z  + F(4)z  + F(5)z.....F(n)z

2) Multiplikation der transponierten Matrix AT mit Matrix A (Zeile Matrix AT mal Spalte Marix A):

 N   = A(1)*A(1)+A(2)*A(2)+A(3)*A(3)+A(4)*A(4)+...
[B] = B(1)*A(1)+B(2)*A(2)+B(3)*A(3)+B(4)*A(4)+...
[C] = C(1)*A(1)+C(2)*A(2)+C(3)*A(3)+...
[D] = D(1)*A(1)+D(2)*A(2)+D(3)*A(3)+...
[E] = E(1)*A(1)+E(2)*A(2)+E(3)*A(3)+...
[F] = F(1)*A(1)+F(2)*A(2)+F(3)*A(3)+...

Der 1. Block ergibt die erste Spalte der linken Seite der Normalgleichung Matrix 6 X 6. Da hier A(n)=1 gesetzt ist, ergibt die Ausrechnung der ersten Zeile die Anzahl der Bedingungsgleichungen N.

Der 2. Block der 6 X 6 Matrix ergibt die zweite Spalte der Normalgleichung.

[B]  = A(1)*B(1)+A(2)*B(2)+A(3)*B(3)+A(4)*B(4)+...
[B2] = B(1)*B(1)+B(2)*B(2)+B(3)*B(3)+B(4)*B(4)+...
[CB] = C(1)*B(1)+C(2)*B(2)+C(3)*B(3)+...
[DB] = D(1)*B(1)+D(2)*B(2)+D(3)*B(3)+...
[EB] = E(1)*B(1)+E(2)*B(2)+E(3)*B(3)+...
[FB] = F(1)*B(1)+F(2)*B(2)+F(3)*B(3)+...

Der 3. Block der 6 X 6 Matrix ergibt die dritte Spalte der Normalgleichung usw.

[B]  = A(1)*C(1)+A(2)*C(2)+A(3)*C(3)+A(4)*C(4)+...
[BC] = B(1)*C(1)+B(2)*C(2)+B(3)*C(3)+B(4)*C(4)+...
[C2 ] = C(1)*C(1)+C(2)*C(2)+C(3)*C(3)+...
[DC] = D(1)*C(1)+D(2)*C(2)+D(3)*C(3)+...
[EC] = E(1)*C(1)+E(2)*C(2)+E(3)*C(3)+...
[FC] = F(1)*C(1)+F(2)*C(2)+F(3)*C(3)+...

Der 4,5 und 6 Block erfolgt anlaog.

3) Multiplikation (Spalte Matrix AT mal Zeile R) der transponierten Matrix AT mit der Spalte Residuenvektor Matrix R, ergit die rechte Seite der Normalgleichung 6 X 6 Matrix.

[R]  = A(1)*X(1)+A(2)*X(2)+A(3)*X(3)+A(4)*X(4)+...
[BR] = B(1)*X(1)+B(2)*X(2)+B(3)*X(3)+B(4)*X(4)+...
[CR] = C(1)*X(1)+C(2)*X(2)+C(3)*X(3)+...
[DR] = D(1)*X(1)+D(2)*X(2)+D(3)*X(3)+...
[ER] = E(1)*X(1)+E(2)*X(2)+E(3)*X(3)+...
[FR] = F(1)*X(1)+F(2)*X(2)+F(3)*X(3)+...

Daraus folgt die  Normalgleichung:

N   dW +    [B]  di + [C]   dw + [D]    dn + [E]   dm +    [F]  dv = [R]
[B] d
W +  [B2]  di + [BC] dw + [BD] dn + [BE] dm + [BF] dv = [BR]
[C] d
W + [CB] di +  [C2 ]  dw + [CD] dn + [CE] dm + [CF] dv = [CR]
[D] d
W + [DB] di + [DC] dw + [D2]  dn + [DE] dm + [DF] dv = [DR]
[E] d
W +  [EB] di + [EC] dw + [ED] dn + [E2]   dm + [EF] dv = [ER]
[F] d
W + [FB] di + [FC]  dw + [FE] dn + [FE]  dm +    [F2] dv = [FR]



Gewichtsfaktoren (p=g). Da die Beobachtungen von der Anzahl Beobachtungsnächte, Instrumentengröße u. -fehler, Beobachterkondition usw. abhängen, sind nicht alle Meßwerte von gleicher Güte. Jede Beobachtung erhält daher ein um so größeres »Gewicht«, je genauer und zuverlässiger eine Messung veranschlagt werden kann. Alle Messungen sind nach Möglichkeit zu berücksichtigen, allerdings existiert unter den Beobachtern kein einheitliches Wichtungs-System.

Gewichtsermittlung nach W.D. Heintz: 3 Messungen erhalten das Gewicht 1. Liegen z. B. 12 Messungen von 4 Beobachtern bei einer geglätteten Distanz von 1.2'' vor, erhält der Normalort (Äquinoktium B1950, J2000) das Gewicht 1.2''*Å(4*4)=4.8.

Distanzen und Winkel erhalten bei Speckle-Interferometrie oder photographischen Messungen gleiche Gewichte, die jedoch viel höher bewerten werden als vergleichbare visuelle Beobachtungen.



Normalgleichung oben mit Gewichtsfaktoren »p(n)«:

[p]    dW +[pB]  di+[pC]   dw+[pD]    dn+[pE]    dm+[pF]   dv = [pR]
[pB] d
W+[pB2 ] di+[pBC] dw+[pBD] dn+[pBE] dm+[pBF] dv = [pBR]

usw.

In der Praxis multipliziert man jede Bedingungsgleichung mit der Quadratwurzel des ihr zuerteilten Gewichts (Åp).

Beispiel: Bedingungsgleichung 4x3 Matrix A: Residuenvektor R:

          Åp(1) A(1)+Åp(1) B(1)+Åp(1) C(1)             Åp(1) X(1)
         
Åp(2) A(2)+Åp(2) B(2)+Åp(2) C(2)             Åp(2) X(2)
A =    
Åp(3) A(3)+Åp(3) B(3)+Åp(3) C(3)    R =   Åp(3) X(3)
         
Åp(4) A(4)+Åp(4) B(4)+Åp(4) C(4)             Åp(4) X(4)

Transponierte Matrix AT:

            Åp(1) A(1)+ Åp(2) A(2)+ Åp(3) A(3)+ Åp(4) A(4)
AT =    
Åp(1) B(1)+  Åp(2) B(2)+ Åp(3) B(3)+ Åp(4) B(4)
           
Åp(1) C(1)+ Åp(2) C(2)+ Åp(3) C(3)+ Åp(4) C(4)

Transponierte Matrix AT mal Residuenvektor R:

[pAX]=p(1) A(1) X(1)+p(2) A(2) X(2)+p(3) A(3) X(3)+p(4) A(4) X(4)
[pBX]=p(1) B(1) X(1)+p(2) B(2) X(2)+p(3) B(3) X(3)+p(4) B(4) X(4)
[pCX]=p(1) C(1) X(1)+p(2) C(2) X(2)+p(3) C(3) X(3)+p(4) C(4) X(4)

Normalgleichung:
[pA2] =p(1) A(1) A(1)+p(2) A(2) A(2)+p(3) A(3) A(3)+p(4) A(4) A(4)
[pBA]=p(1) B(1) A(1)+p(2) B(2) A(2)+p(3) B(3) A(3)+p(4) B(4) A(4)
[pCA]=p(1) C(1) A(1)+p(2) C(2) A(2)+p(3) C(3) A(3)+p(4) C(4) A(4)
[pAB]=p(1) A(1) B(1)+p(2) A(2) B(2)+p(3) A(3) B(3)+p(4) A(4) B(4)
[pB2 ]=p(1) B(1) B(1)+p(2) B(2) B(2)+p(3) B(3) B(3)+p(4) B(4) B(4)
[pCB]=p(1) C(1) B(1)+p(2) C(2) B(2)+p(3) C(3) B(3)+p(4) C(4) B(4)
usw. ...............................................................

Praktisches numerisches Beispiel. Die Verfahrensweise testet man am besten mit einer falschen »Ausgangsbahn« und »idealen« Meßpunkten aus definitiven Bahnelementen. Die Bahnelemente werden im Idealfall definitiv nach mehreren Iterationen erreicht. Bei mehr oder weniger groben Messungen sind reell viele Schritte zu rechnen, evtl. einige Parameter festzuhalten (s. Progr.).

Computer erledigen den immensen Rechenaufwand heute in Sekundenbruchteile. Früher machte sich eigentlich nur van den Bos die Mühe Diff.-Korrekturen in Rechteckkoordinaten unter Anlegung spezieller Tabellen zu berechnen.

REM GFA25  PROGR. BAHNELEMENTEFIT (GFA-BASIC PROGR. FÜR WINDOWS 95)
DEFFN r(x) = x - INT(x / (2 * PI)) * 2 * PI
DIM t(300),p1(300),d(300),g(300),p(30,30),ko(20)
REM DEFINITIVE BAHNELEMENTE DES DOPPELSTERNSYSTEMS CASTOR (GEMINORUM), ADS 6175
to = 1950.65  //PERIASTRONDURCHGANG T
u = 511.3    //UMLAUFZEIT
my = (2 * PI) / u //JÄHRL. BEWEGUNG
a = 7.37     //GROSSE HALBE BAHNACHSE
e = 0.36     //EXENTRIZITÄT
i = RAD(112.9)  //INCLINATION
knot = RAD(41.7) //AUFST. KNOTEN
w = RAD(239.8)    //PERIASTRONARGUMENT
//---------------------------------------------
s = 0
REM SATZ >IDEALER< MESSPUNKTE POSITIONSWINKEL p1 U. DISTANZ d AUS DEFINITIVEN BAHNELEMENTEN
FOR t = 1694 TO 2205 STEP 10   //t = BEOBACHTUNGSZEIT JAHR 1870 - 1940
  PRINT AT(1,1);"JAHR: ";t
  REM BAHNBERECHNUNG DEFINITIVE BAHNELEMENTE ---
  s = s + 1
  g(s) = 1  //GEWICHTSFAKTOR DER MESSUNG = 1
  t(s) = t
  m = my * (t - to) //MITTL. ANOMALIE (DURCHSCHNITTSGESCHWINDIGKEIT)
  ex = m
  FOR s1 = 1 TO 15
    ex = ex - (ex - e * SIN(ex) - m) / (1 - e * COS(ex)) //EXENTR. ANOMALIE
  NEXT s1
  v = 2 * ATN(SQR((1 + e) / (1 - e)) * TAN(ex / 2))  //WAHRE ANOMALIE
  r = a * (1 - e * COS(ex))
  x = r * COS(v + w)
  y = r * SIN(v + w) * COS(i)
  d = SQR(x * x + y * y) //DISTANZ ('')
  x1 = x / d
  y1 = y / d
  p1(s) = DEG(FN r(ATN(y1 / (1 + x1)) * 2 + knot))   //POSITIONSWINKEL GRAD IDEAL
  d(s) = d / 3600  //DISTANZ GRAD IDEAL
NEXT t
anz = s

REM PROVISORISCHE AUSGANGSELEMENTE DER GROBEN BAHN --
to1 = 1940  //AUSGANGSELEMENTE  Epoche
u1 = 480    //Umlaufzeit
my1 = (2 * PI) / u1
a1 = 7   // a
s = 0.3     //Exentrizität
i1 = RAD(100)  // i
knot1 = RAD(30)  //Knoten
w1 = RAD(200)   //Periastronargument
REM ------------------
to1 = 1947  //FESTGEHALTENE WERTE - 1. NEUSTART
u1 = 480
my1 = (2 * PI) / u1
a1 = 7.37
es = 0.3
i1 = RAD(113)
knot1 = RAD(30)
w1 = RAD(207)
REM ------------------
to1 = 1949  //FESTGEHALTENE WERTE - 2. NEUSTART
u1 = 480
my1 = (2 * PI) / u1
a1 = 7.37
es = 0.3
i1 = RAD(113)   //DEG(x)=RAD*(180/PI) IN GRAD
knot1 = RAD(32) //RAD(x)=GRAD*(PI/180) IN RAD
w1 = RAD(228)
REM ------------------
to1 = 1950  //FESTGEHALTENE WERTE - 3. NEUSTART
u1 = 523
my1 = (2 * PI) / u1
a1 = 7.37
es = 0.37
i1 = RAD(113)
knot1 = RAD(40)
w1 = RAD(238)
REM -------------------------
REM BAHNVERBESSERUNG DURCH AUSGLEICHSRECHNUNG -------
FOR s2 = 1 TO 50  //ITERATIONSSCHLEIFE
  REM NULLSETZUNG DER PARAMETER [B]=b1,[CB]=cb usw.
  CLS
  g = 0
  ar = 0
  br = 0
  cr = 0
  dr = 0
  er = 0
  fr = 0
  b1 = 0
  c1 = 0
  d1 = 0
  e1 = 0
  f1 = 0
  bb = 0
  cb = 0
  db = 0
  eb = 0
  fb = 0
  bc = 0
  cc = 0
  dc = 0
  ec = 0
  fc = 0
  bd = 0
  cd = 0
  dd = 0
  ed = 0
  fd = 0
  be = 0
  ce = 0
  de = 0
  ee = 0
  fe = 0
  bf = 0
  cf = 0
  df = 0
  ef = 0
  ff = 0
  FOR s = 1 TO anz
    PRINT AT(1,2);s
    g = g + g(s)  //GEWICHTSFAKTOR DER MESSUNG
    m = my1 * (t(s) - to1) //MITTL. ANOMALIE
    ex = m
    FOR s1 = 1 TO 15
      ex = ex - (ex - es * SIN(ex) - m) / (1 - es * COS(ex)) //EXENTR. ANOMALIE
    NEXT s1
    v = 2 * ATN(SQR((1 + es) / (1 - es)) * TAN(ex / 2)) //WAHRE ANOMALIE
    r = a1 * (1 - es * COS(ex))
    x = r * COS(v + w1)
    y = r * SIN(v + w1) * COS(i1)
    d = SQR(x * x + y * y) //DISTANZ ('')
    x1 = x / d         //EPHEMERIDENWERTE
    y1 = y / d
    pw = DEG(FN r(ATN(y1 / (1 + x1)) * 2 + knot1)) //POSITIONSWINKEL GRAD FESTGEHALTENE WERTE
    REM PARAMETER --------------------
    n = (1 + es * COS(v)) ^ 2 * (1 / COS(ASIN(es)))
    m = (SIN(v) * (2 + es * COS(v))) * (1 / COS(ASIN(es)))
    bo = (-COS(RAD(pw) - knot1)^2 * TAN(v + w1) * SIN(i1))
    co = COS(RAD(pw) - knot1) ^ 2 * (1 / COS(v + w1)) ^ 2 * COS(i1)
    do = (-n * co)
    nn = 10 //EINTRAG N
    e2 = n * co * ((t(s) - to1) / nn)
    fo = m * co
    REM AKKUMULATION DER PARAMETER [B]=b1,[CB]=cb usw.
    b1 = b1 + bo * g(i)
    c1 = c1 + co * g(i)
    d1 = d1 + do * g(i)
    e1 = e1 + e2 * g(i)
    f1 = f1 + fo * g(i)
    bb = bb + bo * bo * g(i)
    cb = cb + co * bo * g(i)
    db = db + do * bo * g(i)
    eb = eb + e2 * bo * g(i)
    fb = fb + fo * bo * g(i)
    bc = bc + bo * co * g(i)
    cc = cc + co * co * g(i)
    dc = dc + do * co * g(i)
    ec = ec + e2 * co * g(i)
    fc = fc + fo * co * g(i)
    bd = bd + bo * do * g(i)
    cd = cd + co * do * g(i)
    dd = dd + do * do * g(i)
    ed = ed + e2 * do * g(i)
    fd = fd + fo * do * g(i)
    be = be + bo * e2 * g(i)
    ce = ce + co * e2 * g(i)
    de = de + do * e2 * g(i)
    ee = ee + e2 * e2 * g(i)
    fe = fe + fo * e2 * g(i)
    bf = bf + bo * fo * g(i)
    cf = cf + co * fo * g(i)
    df = df + do * fo * g(i)
    ef = ef + e2 * fo * g(i)
    ff = ff + fo * fo * g(i)
    IF p1(s) >= 0 AND p1(s) < 90 AND pw > 270 AND pw <= 360 THEN
      p1(s) = p1(s) + 360
    ENDIF
    IF pw >= 0 AND pw < 90 AND p1(s) > 270 AND p1(s) <= 360 THEN
      pw = pw + 360
    ENDIF
    rv = p1(s) - pw //RESIDUENVEKTOR R = BEOBACHTUNG MINUS RECHNUNG B-R
    ar = ar + rv * g(i)    //POSITIONSWINKELANPASSUNG RESIDUEN
    br = br + bo * rv * g(i)
    cr = cr + co * rv * g(i)
    dr = dr + do * rv * g(i)
    er = er + e2 * rv * g(i)
    fr = fr + fo * rv * g(i)
  NEXT s
  m = 6 //EINTRAG ANZAHL GLEICHUNGEN  MATRIX 6 X 6
  n = 6 //EINTRAG ANZAHL UNBEKANNTE
  p(1,1) = g
  p(1,2) = b1
  p(1,3) = c1
  p(1,4) = d1
  p(1,5) = e1
  p(1,6) = f1
  p(1,7) = ar
  REM --------
  p(2,1) = b1
  p(2,2) = bb
  p(2,3) = cb
  p(2,4) = db
  p(2,5) = eb
  p(2,6) = fb
  p(2,7) = br
  REM ----------
  p(3,1) = c1
  p(3,2) = bc
  p(3,3) = cc
  p(3,4) = dc
  p(3,5) = ec
  p(3,6) = fc
  p(3,7) = cr
  REM ----------------
  p(4,1) = d1
  p(4,2) = bd
  p(4,3) = cd
  p(4,4) = dd
  p(4,5) = ed
  p(4,6) = fd
  p(4,7) = dr
  REM ----------
  p(5,1) = e1
  p(5,2) = be
  p(5,3) = ce
  p(5,4) = de
  p(5,5) = ee
  p(5,6) = fe
  p(5,7) = er
  REM ----------------
  p(6,1) = f1
  p(6,2) = bf
  p(6,3) = cf
  p(6,4) = df
  p(6,5) = ef
  p(6,6) = ff
  p(6,7) = fr
  GOSUB elim
  PRINT "NORMALGLEICHUNG:"
  PRINT TAB(2);g * ko(1) + b1 * ko(2) + c1 * ko(3) + d1 * ko(4) + e1 * ko(5) + f1 * ko(6);TAB(25);ar
  PRINT TAB(2);b1 * ko(1) + bb * ko(2) + bc * ko(3) + bd * ko(4) + be * ko(5) + bf * ko(6);TAB(25);br
  PRINT TAB(2);c1 * ko(1) + cb * ko(2) + cc * ko(3) + cd * ko(4) + ce * ko(5) + cf * ko(6);TAB(25);cr
  PRINT TAB(2);d1 * ko(1) + db * ko(2) + dc * ko(3) + dd * ko(4) + de * ko(5) + df * ko(6);TAB(25);dr
  PRINT TAB(2);e1 * ko(1) + eb * ko(2) + ec * ko(3) + ed * ko(4) + ee * ko(5) + ef * ko(6);TAB(25);er
  PRINT TAB(2);f1 * ko(1) + fb * ko(2) + fc * ko(3) + fd * ko(4) + fe * ko(5) + ff * ko(6);TAB(25);fr
  dt = ko(4) / (57.29577951308 * my1 * (1 / COS(ASIN(es))) ^ 2)
  dmy = ko(5) / (nn * (1 / COS(ASIN(es))) ^ 2)
  de = ko(6) / (57.29577951308 * (1 / COS(ASIN(es))))
  PRINT "dknot: ";ko(1)
  PRINT "di   : ";ko(2)
  PRINT "dw:  : ";ko(3)
  PRINT "dt   : ";dt
  PRINT "my   : ";dmy
  PRINT "de   : ";de
  knot1 = knot1 + RAD(ko(1)) //KNOTEN
  i1 = i1 + RAD(ko(2))     // Inclination i
  w1 = w1 + RAD(ko(3))
  to1 = to1 + dt        //Periastrondurchgang
  my1 = my1 + RAD(dmy)  //Jährl. Bewegung
  um = (2 * PI) / my1
  es = es + de        //Exentrizität
  PRINT
  PRINT "i.....: ";DEG(i1);" GRAD"
  PRINT "KNOTEN: ";DEG(knot1);" GRAD"
  PRINT "OMEGA.: ";DEG(w1);" GRAD"
  PRINT "e.....: ";es
  PRINT "T.....: ";to1
  PRINT "U.....: ";um;" Jahre"
  PRINT "my....: ";DEG(my1);" Grad/Jahr"
  PRINT "a.....: ";a1
  PRINT "TASTE >W<"
  REPEAT
  UNTIL UPPER$(INKEY$) = "W" //WARTEN AUF TASTENDRUCK >W<
NEXT s2
END
PROCEDURE elim  //AUSGLEICHSRECHNUNG NACH DER METHODE DER KLEINSTEN QUADRATE
  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
  ko(n) = p(n,n + 1) / p(n,n )  //RÜCKSUBSTITUTION
  j = n
  REPEAT
    j = j - 1
    ko(j) = p(j,n + 1)
    FOR i = j + 1 TO n
      ko(j) = ko(j) - p(j,i) * ko(i)
    NEXT i
  UNTIL j < 2
RETURN

Die Koeffizienten der Normalgleichung sind explizit angeführt. Um Arbeit und Speicherplatz zu sparen, kann die Matrix selbstverständlich zeilen- und spaltenweise dimensionieret werden.

Beispiel Bedingungsgleichung Polynom 3. Grades: y = A + B x + C x*x D x*x*x
NORMALGLEICHUNG EXPLIZIT (4 UNBEKANNTE, 4 RESIDUEN) POLYNOM 3. GRADES:

N       A [x]       B [xx]      C [xxx]        D = [y],1,0,0,0
[x]     A [xx]     B [xxx]     C [xxxx]     D = [xy],0,1,0,0
[xx]   A [xxx]   B [xxxx]   C [xxxxx]   D = [xxy],0,0,1,0
[xxx] A [xxxx] B [xxxxx] C [xxxxxx] D = [xxxy],0,0,0,1

Zeile i (1,2,3,...,n1 Anzahl Bedingungsgleichungen) und Salte (Koeffizienten) einer Bedingungsgeleichung mit z. B. 4 Unbekanten A,B,C,D (SQR(p1(i)=Gewicht der Messung).

FOR i=1 TO n1
REM DIMENSIONIERUNG (i-te Zeile, Spalte) DER BEDINGUNGSGLEICHUNG
kof(i,1)=SQR(p1(i))*1
kof(i,2)=SQR(p1(i))*x(i)  //SPALTE 2
kof(i,3)=SQR(p1(i))*x(i)*x(i)  //SPALTE 3
kof(i,4)=SQR(p1(i))*x(i)*x(i)*x(i) //SPALTE 4
kof(i,5)=SQR(p1(i))*y(i) // RESIDUUM SPALTE 5
NEXT i

DIMENSIONIERUNG NORMALGLEICHUNG:
REM i Zeile, j Spalte der Normalgleichung p(i,j)
n=4    //EINTRAG ANZAHL UNBEKANNTE DER BEDINGUNGSGLEICHUNG
n1=7    //EINTRAG ANZAHL x(i)
REM NORMALGLEICHUNG p(i,j) ------------
FOR i=1 TO n
FOR j=1 TO n+1
p(i,j)=0
FOR k=1 TO n1
p(i,j)=p(i,j)+kof(k,i)*kof(k,j) //AKKUMULATION
NEXT k
NEXT j
 NEXT i
REM -------------------------------------
FOR i=1 TO n
p(i,n+1+i)=1 //RESIDUEN DER GRÖSSE 1
NEXT i
REM ------------------------------------------------

Numerisches Beispiel. y = A+B*x+C*x*x+D*x*x*x. A=10, B=3, C=5, D=15.5. x(1)=3, x(2)=5, x(3)=10, x(4)=8, x(5)=3, x(6)=9, x(7)=2.2. y(1)=482.5, y(2)=2087.5, y(3)=16040, y(4)=8290, y(5)=482.5, y(6)=11741.5, y(7)=205.844.

REM GFA26  MULTIPLE REGRESSION
DIM x(10),y(10),kof(10,10),p1(10),p(10,10),ko(10,10)
A = 10
b = 3
c = 5
d = 15.5
//---------------------------------------------------
x(1) = 3
x(2) = 5
x(3) = 10
x(4) = 8
x(5) = 3
x(6) = 9
x(7) = 2.2
//-----------------------------------------------------
y(1) = 482.5
y(2) = 2087.5
y(3) = 16040
y(4) = 8290
y(5) = 482.5
y(6) = 11741.5
y(7) = 205.844
//------------------------------------------
p1(1) = 1 //GEWICHT
p1(2) = 1
p1(3) = 1
p1(4) = 1
p1(5) = 1
p1(6) = 1
p1(7) = 1
//-----------------------------------
REM EXPLIZITE AKKUMULATION POLYNOM 3. GRADES
n1 = 7
x = 0    //INITIALISIERUNG
xx = 0
xxx = 0
xxxx = 0
xxxxx = 0
xxxxxx = 0
p = 0
y = 0
yy = 0
xy = 0
xxy = 0
xxxy = 0
FOR i = 1 TO n1
  p = p + p1(i)
  x = x + x(i) * p1(i)
  xx = xx + x(i) ^ 2 * p1(i)
  xxx = xxx + x(i) ^ 3 * p1(i)
  xxxx = xxxx + x(i) ^ 4 * p1(i)
  xxxxx = xxxxx + x(i) ^ 5 * p1(i)
  xxxxxx = xxxxxx + x(i) ^ 6 * p1(i)
  y = y + y(i) * p1(i)
  yy = yy + y(i) ^ 2 * p1(i)            //RESIDUEN
  xy = xy + x(i) * y(i) * p1(i)
  xxy = xxy + x(i) ^ 2 * y(i) * p1(i)
  xxxy = xxxy + x(i) ^ 3 * y(i) * p1(i)
NEXT i
REM ---------------
resi = 5 //ANZAHL RESIDUEN
m = 4 //ANZAHL GLEICHUNGEN
n = 4 //ANZAHL UNBEKANNTE
m = m + resi
p(1,1) = p
p(1,2) = x
p(1,3) = xx
p(1,4) = xxx
p(1,5) = y //1. RESIDUENVEKTOR
p(1,6) = 1  //2. RESIDUENVEKTOR
p(1,7) = 0
p(1,8) = 0
p(1,9) = 0
p(2,1) = x
p(2,2) = xx
p(2,3) = xxx
p(2,4) = xxxx
p(2,5) = xy //RESIDUENVEKTOR
p(2,6) = 0
p(2,7) = 1
p(2,8) = 0
p(2,9) = 0
p(3,1) = xx
p(3,2) = xxx
p(3,3) = xxxx
p(3,4) = xxxxx
p(3,5) = xxy //RESIDUENVEKTOR
p(3,6) = 0
p(3,7) = 0
p(3,8) = 1
p(3,9) = 0
p(4,1) = xxx
p(4,2) = xxxx
p(4,3) = xxxxx
p(4,4) = xxxxxx
p(4,5) = xxxy //RESIDUENVEKTOR
p(4,6) = 0
p(4,7) = 0
p(4,8) = 0
p(4,9) = 1
GOSUB elim
PRINT "KOEFFIZIENT a: ";ko(1,1)
PRINT "KOEFFIZIENT b: ";ko(2,1)
PRINT "KOEFFIZIENT c: ";ko(3,1)
PRINT "KOEFFIZIENT d: ";ko(4,1)
REM EINSETZEN IN NORMALGLEICHUNG (KONTROLLE):
PRINT n1 * ko(1,1) + x * ko(2,1) + xx * ko(3,1) + xxx * ko(4,1),y
PRINT x * ko(1,1) + xx * ko(2,1) + xxx * ko(3,1) + xxxx * ko(4,1),xy
PRINT xx * ko(1,1) + xxx * ko(2,1) + xxxx * ko(3,1) + xxxxx * ko(4,1),xxy
PRINT xxx * ko(1,1) + xxxx * ko(2,1) + xxxxx * ko(3,1) + xxxxxx * ko(4,1),xxxy
REM MESSUNG
x = 3  //EINTRAG
y = ko(1,1) + ko(2,1) * x + ko(3,1) * x * x + ko(4,1) * x * x * x
PRINT y
PRINT
PRINT
// DIMENSIONIERUNG-----------------------------------------------------------------------
n1 = 7 //ANZAHL
FOR i = 1 TO n1
  REM DIMENSIONIERUNG (i-te Zeile, Spalte) DER BEDINGUNGSGLEICHUNG
  kof(i,1) = SQR(p1(i)) * 1
  kof(i,2) = SQR(p1(i)) * x(i)  //SPALTE 2
  kof(i,3) = SQR(p1(i)) * x(i) ^ 2  //SPALTE 3
  kof(i,4) = SQR(p1(i)) * x(i) ^ 3    //SPALTE 4
  kof(i,5) = SQR(p1(i)) * y(i) // RESIDUUM SPALTE 5
NEXT i
REM DIMENSIONIERUNG NORMALGLEICHUNG:
REM i Zeile, j Spalte der Normalgleichung p(i,j)
n = 4    //EINTRAG ANZAHL UNBEKANNTE DER BEDINGUNGSGLEICHUNG
REM NORMALGLEICHUNG p(i,j) ------------
FOR i = 1 TO n
  FOR j = 1 TO n + 1
    p(i,j) = 0
    FOR k = 1 TO n1
      p(i,j) = p(i,j) + kof(k,i) * kof(k,j) //AKKUMULATION
    NEXT k
  NEXT j
NEXT i
REM -------------------------------------
FOR i = 1 TO n
  p(i,n + 1 + i) = 1 //RESIDUEN DER GRÖSSE 1
NEXT i
//--------------------------------------------------------
resi = 5 //ANZAHL RESIDUEN
m = 4 //ANZAHL GLEICHUNGEN
n = 4 //ANZAHL UNBEKANNTE
m = m + resi
GOSUB elim
PRINT "KOEFFIZIENT a: ";ko(1,1)
PRINT "KOEFFIZIENT b: ";ko(2,1)
PRINT "KOEFFIZIENT c: ";ko(3,1)
PRINT "KOEFFIZIENT d: ";ko(4,1)
REM MESSUNG
x = 3  //EINTRAG
y = ko(1,1) + ko(2,1) * x + ko(3,1) * x * x + ko(4,1) * x * x * x
PRINT y
STOP
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

Sternbeobachter - Sterntagebuch - Produktinformation - www.spaceglobe.de