Doppelsterne7




		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/(2PI))(2PI) 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+bT 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+be Periastronzeit" //r=a(1-e^2)/(1-eCOS(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 NA+XOB=YO; XOA+XXB=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) dv=r(i)=57.29578 dd(i)/d(i) Das Zeitmaß »N« (Jahre) liegt stets innerhalb des auszugleichenden Bahnbogens (N=(t₁-t₂)/2; t₁ u. t₂=Anfangs- u. Endjahr des auszugleichenen Bahnstücks, Epoche ep=t₁+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)) dv = 57.29578 sec( v) de dn=my sec(v)² dT my=jährl. Eigenbewegung in GRAD da° =57.29578 da/a da=(da°a)/57.29578 dm = N sec(v)² dmy N=Zeitmaß: c =(t-T)/N. de=dv / (57.29578(1/cos(arcsin(e)))) dT=dn/(my sec(v)²) dmy = dm/(N sec(v)² Koeffizienten A,B,C,D,E,F (P=pw=aus den zu verbessernden Elementen berechneter Positionswinkel: A(i)=1, B(i)=-cos(P-W)² tan(v+w) sin(i) C(i)= (cos(P-W)² /cos(v+w)²) 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))² m=sec(v) sin(v) (2+ecos(v)); v=wahre Anomalie. s=-sec(v) cos(v) (1+ecos(v)) o=tan(arcsin(e)) sin(v) (1+ecos(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 A^T 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 A^T): 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 A^T = 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 A^T mit Matrix A (Zeile Matrix A^T 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)+... [B²] = 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)+... [C² ] = 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 A^T mal Zeile R) der transponierten Matrix A^T 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] dW + [B²] di + [BC] dw + [BD] dn + [BE] dm + [BF] dv = [BR] [C] dW + [CB] di + [C² ] dw + [CD] dn + [CE] dm + [CF] dv = [CR] [D] dW + [DB] di + [DC] dw + [D²] dn + [DE] dm + [DF] dv = [DR] [E] dW + [EB] di + [EC] dw + [ED] dn + [E²] dm + [EF] dv = [ER] [F] d W + [FB] di + [FC] dw + [FE] dn + [FE] dm + [F²] 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''Å(44)=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] dW+[pB² ] 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 A^T: Åp(1) A(1)+ Åp(2) A(2)+ Åp(3) A(3)+ Åp(4) A(4) A^T = Å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: [pA²] =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) [pB² ]=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 xx D xxx 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+Bx+Cxx+Dxxx. 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