Raumfahrt-Programmsammlung                       März 2019, Thomas Muetsch

Dokumentation zum Raumfahrtsammelprogramm RFP.EXE

Lauffähigkeit:

Das vorliegende Programm läuft seit 1988 auf IBM XT/AT und Kompatiblen
(unter MS-DOS 32-bit oder DOS-Emulation (z.B. DOSBox 0.74) auf allen
neueren Versionen) mit Festplatte und mindestens 640 kByte Arbeitsspeicher.
Ein Drucker wird nicht benötigt. Maus-Zeiger werden nicht unterstützt.
Für den fehlerlosen
Ablauf wird der Kommando-Interpreter (COMMAND.COM) im Root-Directory und
ein vollständiges Programmpaket mit den 6 EXE-Dateien RFP, RFP1, RFP2,
RFPGRAF1, RFPGRAF2 und RFPGRAF3 vorausgesetzt.
Die Raketenbeispiele "*.RAK" werden nicht benötigt.
Grafikteile laufen nur auf Hercules- und EGA(256k RAM)/VGA-Grafikkarte.
Sind diese nicht vorhanden, so wird z.T. eine Grobgrafik aus
ASCII-Zeichen verwendet. Die Programme und deren Namen dürfen i.a. nicht
verändert werden. Die in eckigen Klammern angegebenen DUMMY-Werte sind
vorgegeben und brauchen bei Nichtänderung nur bestätigt zu werden (mit ENTER).
Bei Zeichenketten lassen sich die DUMMY-Werte mit den Steuertasten bearbeiten.
Die TAB-Taste ermöglicht dabei ein Löschen des DUMMY-Inhalts. Das Programm
ist menügesteuert; die einzelnen Programme werden durch Betätigen der
davorstehenden Buchstaben (in Groß- oder Kleinschreibung) und Ziffern
gestartet. Dies gilt auch für die Auswahl der Himmelskörpers. Die Menüs
können durch die Eingabe einer nicht belegten Ziffer oder Taste wieder
verlassen werden. Es wurden teilweise Begriffe eingeführt, wie sie in der
einschlägigen Literatur nicht zu finden sind (z.B. Hohmannaufenthaltsdauer).
Diese Dokumentation soll dazu dienen, diese Begriffe durch Beschreibung des
formelmäßigen Zusammenhangs darzustellen. In den numerischen
Unterprogrammen kann jederzeit mit der PAUSE-Taste der Durchlauf gestoppt
und mit einer anderen Taste die Ausführung fortgesetzt werden. In dieser
Dokumentation erfolgt die Nennung der Unterprogramme immer mit dem jeweiligen
alphanumerischen Zeichen für dessen Aufruf. Einzelne Ergebnisse aus den
verschiedenen Unterprogrammen werden beim Aufruf eines anderen Programmteils
als Eingabe wiederverwendet; diese Werte (in eckigen Klammern angegeben)
brauchen lediglich mit ENTER bestätigt zu werden. Damit lassen sich
komplexere Zusammenhänge (z.B. für den Wechsel eines Koordinatensystems)
einfacher durchrechnen. Auch die Richtigkeit der verwendeten Berechnungen
kann z.T. durch Rückrechnung mit einem anderen Unterprogramm bestätigt
werden.

Definition der physikalischen Einheiten:

Für die Dimensionierung der physikalischen Größen wurden die SI-Einheiten
zugrunde gelegt. Diese sind in den Ein- und Ausgabeteilen nicht explizit
aufgeführt. Werden andere Dimensionionierungen verwendet, so sind diese in
den Eingabezeilen meist in runden Klammern vorgegeben. Die einzelnen
Dimensionen zur Erinnerung : kg (Masse), m (Länge), sec (Zeit), Ampere
(Stromstärke), mol (Stoffmenge, nicht in kmol !), Kelvin (Temperatur),
candela (Lichtstärke). Das Jahr als Zeiteinheit umfasst 365 volle Erdtage.
Um die Fehler- und damit Ausstiegshäufigkeit des Programms zu vermindern,
werden i.a. unsinnige Dateneingaben (z.B. negative Massen) nicht akzeptiert.
Ist der Definitionsbereich allerdings nicht von einfachem Zusammenhang, kann
u.U. eine unsinnige Berechnung erfolgen. Kommt es zu Laufzeitfehlern, wird
auf die nächsthöhere Programmebene gesprungen, bis zuletzt auch das
Hauptprogramm verlassen wird. Gleiches gilt analog für den stets möglichen
Abbruch mit CTRL-Break.

Copyright und Weitergabe des Programms:

*******************************************
** (c) 1988-2019 Copyright worldwide     **
**::::::::: by ::::::::::::::::::::::::::**
**          Thomas Muetsch, Windischbuch **
**:::::::::::::::::::::::::::::::::::::::**
** e-mail : thomas.muetsch @ freenet.de  **
*******************************************

Eine Weitergabe der EXE-Dateien ist uneingeschränkt erlaubt und ausdrücklich
erwünscht. Die Weitergabe der Quellcodes ist nur demjenigen gestattet, welchem
ich es persönlich ausgehändigt habe, allerdings unter Ausschluss der Vermarktung.
Diesen Anspruch behalte ich mir für einen späteren Zeitpunkt vor.
Die Programme wurden mit BORLAND Turbo-Pascal Version 6.0 SN P91121083
erzeugt.

Physikalische Annahmen und Vereinfachungen:

Die zur Verfügung stehenden Himmelskörper werden als kugelförmig angenommen.
Störeinflüsse anderer Himmelskörper werden vernachlässigt (Einkörperproblem).
Relativistische Einflüsse werden ebenfalls vernachlässigt.
Als Apsidenhöhe wird die Bahnordinate im Brennpunkt definiert.
Der Apsidenwinkel entspricht der wahren Anomalie der klassischen Bahnmechanik.
Auf Ausnahmen von diesen Grundsätzen wird besonders hingewiesen.

Beschreibung der Funktionen:

In dieser vorliegenden Dokumentation wurden z.T. Ausschnitte der
Quellprogramme verwendet und kommentiert. Der Funktionen-Syntax
besteht dabei aus dem Standardbefehlssatz von TURBO Pascal 6.0:
 abs       : Absolutbetrag-Funktion
 arctan   : Arcus-Tangens-Funktion in rad
 cos       : Cosinus-Funktion in rad
 exp       : Exponential-Funktion
 ln          : natürlicher Logarithmus
 max      : Maximum-Funktion
 min       : Minimum-Funktion
 round   : Gerundete Ganzzahl einer rationalen Zahl
 sin        : Sinus-Funktion in rad
 sqr        : Quadrat-Funktion
 sqrt       : Wurzel-Funktion
 trunc     : Ganzzahliger Anteil einer rationalen Zahl

und eigenen Erweiterungen:
 gatn           : Arcus-Tangens-Funktion in Grad
 gcos          : Cosinus-Funktion in Grad
 gsin           : Sinus-Funktion in Grad
 hoch          : Potenz-Funktion, z.B. hoch(a,b)=a^b
 racs,gacs : Arcus-Cosinus-Funktion rad,in Grad
 rasn,gasn : Arcus-Sinus-Funktion in rad,Grad
 rtan,gtan   : Tangens-Funktion in rad,Grad
 sgn            : Vorzeichen-Funktion (sgn={-1,0,+1})

Unter den Berechnungen werden lediglich die physikalischen Grundlösungen
aufgezeigt. Programmtechnisch müssen aber Ausnahmen von diesen Regeln
noch abgefangen werden, die hier jedoch nicht näher beschrieben werden.
Auch auf die Existenz von Doppellösungen bei bestimmten Gleichungen (häufig
bei Ellipsen- bzw. Hyperbelbahnen der klassischen Bahnmechanik) wird hier
nicht eingegangen.

Zur Verfügung stehende Unterprogramme:

 1  Ballistische Bahnen              7  Swingby-Manöver
 2  Kreisförmige Bahnen           8  Bodenspur zeichnen
 3  Elliptische Bahnen                9  Seile
 4  Allgemeine Bahnen
 5  Bahnbestimmung aus Ortsgrößen
 6  Bahnbestimmung aus Flugdauer

 a  Planetendaten                                  e  Raketengröße für Einstufenträger
 b  Erdabplattungseinflüsse                  f  Raketengeschwindigkeitsvermögen
 c  Manövergeschwindigkeitsbedarf   g  Raumschlepper
 d  Bahnelemente ermitteln                  h  Raumfahrtantriebe

 Numerische Programme : i  Freiflug im Erde-Mond-System
                                              j  Wiedereintrittsprogramm
                                              k  Raketenstartprogramm
                                              l  Planetensystem

Belegung der Konstanten :

Name des Himmelskörpers:
   pls: array[0..23] of string[10] =('Sonne','Merkur','Venus','Erde',
        'Mars','Jupiter','Saturn','Uranus','Neptun','Pluto','Ceres',
        'Eros','Hermes','Mond','Phobos','Deimos','Io','Europa',
        'Ganymed','Callisto','Titan','Oberon','Triton','Weltall');

Radius des Himmelskörpers:
   rv : array[0..23] of real =(696000e3,2439e3,6052e3,6371.2e3,
        3394e3,71398e3,60435e3,25900e3,24712e3,1142e3,380e3,
        7e3,300,1738e3,20e3,13e3,1818e3,1530e3,
        2610e3,2445e3,2575e3,775e3,1360e3,1e16);

Mittlere große Halbachse der Himmelskörperbahn:
   tv : array[0..23] of real =(0,5.791e10,1.0821e11,1.49598e11,
        2.2794e11,7.784e11,1.427e12,2.87e12,4.496e12,5.9135e12,4.138e11,
        2.181e11,1.93e11,384403e3,938e4,2346e4,422e6,671e6,
        1.07e9,1.883e9,1.22186e9,5834e5,3546e5,0);

Masse des Himmelskörpers:
   mv : array[0..23] of real =(1.9891e30,3.302e23,4.869e24,5.974e24,
        6.419e23,1.8988e27,5.684e26,8.698e25,1.028e26,1.3e22,1.17e21,
        5e15,4e11,7.3432e22,1.34e16,2.3e15,8.907e22,4.869e22,
        1.488e23,1.063e23,1.36e23,6e21,3e22,1.9891e30);

Rotationsdauer des Himmelskörpers:
   lv : array[0..23] of real =(2200000,5100000,21000000,86164,
        88643,35430,36840,38940,56400,550000,32688,
        18972,0,2360592,0,0,0,0,
        0,0,0,0,0,0);

Exzentrizität der Himmelskörperbahn:
   ev : array[0..23] of real =(0,0.206,0.006787,0.016721,
        0.093,0.048,0.056,0.047,0.009,0.25,0.079,
        0.223,0.475,0.0549,0.021,0.0028,0.0001,0.0001,
        0.0014,0.0074,0.0292,0.001,0,0);

Reziproker Wert der Abplattung des Himmelskörpers:
   abv: array[0..23] of real =(0,0,0,297,
        170,15,9.5,42,40,0,0,
        0,0,500,0,0,0,0,
        0,0,0,0,0,0);

Bahnneigung des Himmelskörpers:
   nwv: array[0..23] of real =(0,10,6,23.5,
        24,3,27,98,29.6,0,0,
        0,0,5.15,0,0,0,0,
        0,0,0,0,0,0);

Gravitationskonstante:
   gr:  6.672e-11

Zur Verfügung stehende Himmelskörper:

 0 Sonne
 1 Merkur
 2 Venus
 3 Erde      a Mond
 4 Mars      b Phobos   c Deimos
 d Ceres    e Eros         f Hermes
 5 Jupiter   g Io              h Europa   i Ganymed  j Callisto
 6 Saturn    k Titan
 7 Uranus   l Oberon
 8 Neptun   m Triton
 9 Pluto
In einer Aufrufroutine wird dem gewählten Himmelskörper die
INTEGER-Variable vc zugewiesen. Außerdem werden die folgenden häufig
benutzten Variablen für den jeweiligen Himmelskörper zugewiesen:
 Körperradius        re = rv [vc]
 Körpermasse       me = mv [vc]
 Rotationsdauer    te = lv [vc]
 Körperkonstante  mu = gr*me

Beschreibung der Unterprogramme:

1 Ballistische Bahnen:

Unterprogramm zur Ermittlung charakteristischer Daten einer ballistischen
Flugbahn im Schwerefeld eines Himmelskörpers. Dabei werden die Einflüsse der
Abplattung und des Luftwiderstands, sowie die Dauer der Antriebsphase
vernachlässigt.
Eingaben:
 Brennschlussgeschwindigkeit               ve
 Abschussneigungswinkel                      wi  (gegenüber der Horizontalen)
 Abschussrichtung (0=Osten,-90=Süd)  wb
 Abschussbreitengrad                            bg1
 Abschusslängengrad (Osten > 0)         lg1
Berechnungen:
  vpl  = 2*pi*re*gcos(bg1)/te
  vx   = vpl+ve*gcos(wb)*gcos(wi)
  vy   = ve*gsin(wb)*gcos(wi)
  vg   = sqrt(sqr(vx)+sqr(vy)+sqr(ve*gsin(wi)))
  wi1  = gasn(ve*gsin(wi)/vg)
  wb1  = gacs(vx/vg/gcos(wi1))
  ex   = sqrt(sqr(vg*vg/mu*re-1)*sqr(gcos(wi1))+sqr(gsin(wi1)))
  p    = vg*vg*sqr(gcos(wi1))/mu*re*re
  hg   = p/(1-ex)-re
  r9   = p/(1-ex)
  r8   = p/(1+ex)
  gh   = (r9+r8)/2
  t0   = 2*pi*sqrt(gh*gh/mu*gh)
  we   = gacs((1-p/re)/ex)
  in1  = gacs(gcos(abs(bg1))*gsin(90-wb1))
  y1   = 2*arctan(gtan((180-we)/2)/sqrt((1+ex)/(1-ex)))
  t1   = (y1-ex*sin(y1))/sqrt(mu/gh/gh/gh)
  t4   = (t0-2*t1)/lv[vc]*360
  om   = gasn(gsin(bg1)/gsin(in1))-(180-we)
  rek  = lg1-gacs(gcos(om+180-we)/gcos(bg1))
  bg2  = gasn(gsin(in1)*gsin(om+180+we))
  llg  = gacs(gcos(om+180+we)/gcos(bg2))
  lg2  = llg+rek-t4
  x1   = gcos(bg1)*gcos(lg1)
  y1   = gcos(bg1)*gsin(lg1)
  z1   = gsin(bg1)
  x2   = gcos(bg2)*gcos(lg2)
  y2   = gcos(bg2)*gsin(lg2)
  z2   = gsin(bg2)
  l5   = sqrt(sqr(x1-x2)+sqr(y1-y2)+sqr(z1-z2))
  rwe  = 2*re*arctan(l5/2/sqrt(1-l5*l5/4))
  kap1 = gatn(sqrt(1-vg*vg*re/mu))
Ausgaben:
 Bahnanstiegswinkel                                           wi1
 Flugrichtung für ruhenden Beobachter              wb1
 effektive Startgeschwindigkeit                            vg
 Startplatzdrehgeschwindigkeit                            vpl
 Abschusswinkel für max.Reichweite ca.             kap1
 Gipfelhöhe                                                           hg/1000   km
 Reichweite                                                           rwe/1000  km
 Landebreitengrad                                                  bg2
 Landelängengrad                                                  lg2
 Winkelsprung auf der Himmelskörperoberfläche  rwe/re*180/pi
 Winkelsprung der Bahn                                         2*we
 Flugdauer                                                              t0-2*t1
 Inklination                                                              in1
 Exzentrizität                                                          ex
 Argument des Perigäums                                     om
 Rektaszension                                                      rek

2 Kreisförmige Bahnen:

Unterprogramm zur Ermittlung charakteristischer Größen einer kreisförmigen
Umlaufbahn um einen Himmelskörper, sowie spezifischer Oberflächendaten.
Eingaben:
 Höhe (km)                               h1
 Inklination bzw.Breitengrad   in1
Berechnungen:
 r1   = re+1000*h1
 inl  = in1/180*pi
 t0   = sqrt(r1*r1/mu*r1)*2*pi
 rh   = r1*sin(inl)
 r3   = sqrt(re*re-rh*rh)
 t1   = 2*arctan((r3/r1)/sqrt(1-sqr(r3/r1)))/sqrt(mu/r1/r1/r1)
 t4   = 2*pi*sqrt(sqr(tv[vc])/mv[0]*tv[vc]/gr)
 s3   = round(lv[vc]*t4/(t4-lv[vc]))
 um   = 2*pi*re*cos(inl)
 rd   = re*cos(inl)
 bb   = rd*rtan(nwv[vc]*pi/180)/rtan((90-in1)*pi/180)
 gm   = 90-((pi/2-arctan((bb/rd)/sqrt(1-sqr(bb/rd))))*180/pi)
 tl   = 2*pi*rd*(1+gm/90)/um
 aa3  = rv[vc]*(1-1/abv[vc]*in1/90)
 um2  = 2*pi*aa3*cos(inl)
 ve   = um2/te
 ge   = 9.8144+(0.0178/90*in1)          (für Erde)
 ge   = gr*mv[vc]/aa3/aa3               (für andere Himmelskörper)
 gz   = abs(sqr(2*pi/te)*aa3*cos(inl))
 aef  = sqrt(sqr(cos(inl)*ge-gz)+sqr(sin(inl)*ge))
 b1   = pi/2-rasn(re/r1)
 s6   = abs(lv[vc]*t0/(lv[vc]-t0*cos(inl)))
 s4   = lv[vc]*t0/(lv[vc]-t0*cos(inl))*b1/pi
Ausgaben:
 Umlaufzeit                                              t0
 Dunkelzeit bei Tag- u.Nachtgleiche    t1 <=> t1/t0*100 %
 Kreisbahngeschwindigkeit                 sqrt(mu/r1)
 Fluchtkickgeschwindigkeit                 sqrt(mu/r1)*(sqrt(2)-1)
 Kinetische Energie (kWh/kg)             mu/r1/7200000
 Potentielle Energie (kWh/kg)             mu*(1/re-1/r1)/3.6e6 = (1-re/r1)*100 %
 Dauer eines Tages                             s3
 Tageslänge bei Sonnenwende           +/- abs(tl-1)*100 %
 Himmelskörperdrehgeschwindigkeit      ve
 Oberflächenanziehung                          ge - gz = aef
 Himmelskörperblickwinkel                 (pi-2*b1)*180/pi
 Himmelskörperbreitengrade              b1*360/pi (ohne Neigungsbreitengrade)
 Breitenkilometer                                  2*re*b1/1000
 Blickfläche                                            50*(1-cos(pi/2-b1)) %
 Himmelskörpersichtbarkeit                 (1-cos(b1))*50 %
 Synodische Oberflächenperiode         s6
 Satellitensichtbarkeit (max.)                -s4 <=> b1/pi*100 %
 Satellitenentfernung (max.)                  sqrt(sqr(re+1000*h1)-re*re)/1000 km

3 Elliptische Bahnen:

Unterprogramm zur Ermittlung charakteristischer Größen einer elliptischen
Umlaufbahn um einen Himmelskörper. Eine Hilfsroutine ermöglicht dabei die
Zuweisung von verschiedenen Monddaten zum jeweiligen Gravitationszentrum.
Eingaben:
 Peri-höhe (km)                           h1
 Apo-höhe  (km)                          h2
 Apsidenwinkel (0 - 180 Grad)  q
Berechungen:
 qw   = q/180*pi
 r1   = re+1000*h1
 r2   = re+1000*h2
 gh   = (r1+r2)/2
 ex   = (r2-r1)/(r2+r1)
 pp   = (1+ex)*r1
 rr   = pp/(1+ex*cos(qw))
 s2   = sqrt(mu*pp)
 vrad = sqrt(mu/pp)*ex*sin(qw)
 vnor = sqrt(mu*pp)/rr
 vt   = s2/pp*sqrt(1+2*ex*cos(qw)+ex*ex)
 v3   = s2/pp*sqrt(1+2*ex+ex*ex)
 t0   = 2*pi*sqrt(gh*gh/mu*gh)
 y1   = 2*arctan(rtan(qw/2)/sqrt((1+ex)/(1-ex)))
 t1   = (y1-ex*sin(y1))/sqrt(mu/gh/gh/gh)
 ps   = pi-arctan((1+ex*cos(qw))/ex/sin(qw))
 psi  = ps*180/pi
 hr   = (rr-re)/1000
 vk   = sqrt(mu/rr)
 kvb  = sqrt(vk*vk+vt*vt-2*vt*vk*sin(psi/180*pi))
 vk1  = sqrt(mu/r1)
 kvb1 = v3-vk1
 v2   = sqrt(pp*mu)/pp*(1-ex)
 v1   = sqrt(pp*mu)/pp*(1+ex)
 t3   = 2*pi*sqrt(r1*r1/mu*r1)
 t4   = 2*pi*sqrt(r2*r2/mu*r2)
 dv   = sqrt(mu/r1)*(sqrt(2*r2/(r1+r2))-1)
        +sqrt(mu/r2)*(1-sqrt(2*r1/(r1+r2)))
 dvo  = sqrt(mu/r1)*(sqrt(2*rr/(r1+rr))-1)
        +sqrt(mu/rr)*(1-sqrt(2*r1/(r1+rr)))
 dv1  = (sqrt(2)-1)*(sqrt(mu/r1)+sqrt(mu/r2))
 dv2  = (sqrt(2)-1)*(sqrt(mu/r1)+sqrt(mu/rr))
 t5   = t3*t4/(t4-t3)
 tj   = t0/(t4/t3-1)
 tb   = t0/(t3/t4-1)
 n8   = trunc(abs(tb)/t5)+1
 tb   = tb+t5*n8
Ausgaben:
 Exzentrizität                                             ex
 Ortstangentenwinkel                               psi
 Umlaufdauer                                            t0
 Ortszeit                                                     t1
 Höhe                                                         hr km
 Peri-Geschwindigkeit                             v1
 Apo-Geschwindigkeit                             v2
 Anziehung                                                mu/rr/rr
 Ortsgeschwindigkeit                               vt
 Radialgeschwindigkeit                           vrad
 Normalgeschwindigkeit                          nor
 Ortskreisbahngeschwindigkeit               vk
 Kreisbahngeschwindigkeitsbedarf        kvb
 Bahntransferbedarf für Ortshöhe            kvb+kvb1
 Hohmanngeschwindigkeitsbedarf          dv
  bzw. für Ortshöhe                                    dvo
 3-Impuls-Geschwindigkeitsbedarf         dv1
  bzw. für Ortshöhe                                    dv2
 Spiralgeschwindigkeitsbedarf               sqrt(mu/r1)*(1-sqrt(r1/r2))
  bzw. für Ortshöhe                                    sqrt(mu/r1)*(1-sqrt(r1/rr))
 Synodische Periode                                t5
 Hohmannaufenthalt innen mindestens   tj
                     außen mindestens               tb
 Kinetische Energie (kWh/kg)                 vt*vt/7.2e6
 Potentielle Energie (kWh/kg)                 mu*(1/re-1/rr)/3.6e6

4 Allgemeine Bahnen:

Unterprogramm zur Ermittlung charakteristischer Größen einer Umlaufbahn um
einen Himmelkörper unter Vorgabe des Flugzustands im Perizentrum.
Eingaben:
 Starthöhe (km)                 h1  (des Perizentrums)
 Startgeschwindigkeit       v1
 Ortshöhe (km)                  hr
Berechnungen:
 r2  = re+hr*1000
 r1  = re+h1*1000
 si  = r1*v1
 hh  = v1*v1-2*mu/r1
 ex  = sqrt(1+hh*si/mu*si/mu)
 pp  = v1*v1/mu*r1*r1
 f   = (pp-r2)/ex/r2
 vt  = sqrt(hh+2*mu/r2)
 Für ex=1:
  ph = pi/2-rasn(f)
  tu = pp*pp/sqrt(mu*pp)*(rtan(ph/2)/2+hoch((rtan(ph/2)),3)/6)
 Für ex<1:
  gh = pp/(1-ex*ex)
  ph = pi/2-rasn(f)
  t2 = 2/sqrt(1-ex*ex)*arctan((1-ex)*rtan(ph/2)/sqrt(1-ex*ex))
  tu = sqrt(pp/mu*pp*pp)
       *(ex*sin(ph)/(ex*ex-1)/(1+ex*cos(ph))-t2/(ex*ex-1))
  h2 = (gh*(1+ex)-re)/1000
  t4 = 2*pi*sqrt(gh/mu*gh*gh)
 Für ex>1:
  gh = pp/(ex*ex-1)
  ph = abs(pi/2-rasn(f))
  t2 = (ln(abs(((ex-1)*rtan(ph/2)+sqrt(ex*ex-1))/((ex-1)*rtan(ph/2)
        -sqrt(ex*ex-1)))))/sqrt(ex*ex-1)
  tu = sqrt(pp/mu*pp*pp)
       *(ex*sin(ph)/(ex*ex-1)/(1+ex*cos(ph))-t2/(ex*ex-1))
  vf = sqrt(2*mu/r2)
  vr = sqrt(vt*vt-vf*vf)
 x   = 1/(v1*v1/(mu/r1)-1)
 ps  = pi-arctan((1+ex*cos(ph))/ex/sin(ph))
 kvb = sqrt(vt*vt+mu/r2-2*vt*sqrt(mu/r2)*sin(ps))
 psi = ps*180/pi
Ausgaben:
 Abflugwinkel                                         rasn(x)*180/pi
 Exzentrizität                                          ex
 Parabolische Geschwindigkeit         v1-sqrt(2*mu/r1)
 Startkreisbahngeschwindigkeit        sqrt(mu/r1)
 Kickgeschwindigkeit                          v1-sqrt(mu/r1)
 Apo-höhe                                             h2 km
 Apsidenwinkel                                    ph*180/pi
 Ortstangentenwinkel                          psi
 Ortsgeschwindigkeit                          vt
 Ortskreisbahngeschwindigkeit         sqrt(mu/r2)
 Kreisbahngeschwindigkeitsbedarf  kvb
 Umlaufzeit                                           t4
 Ortszeit                                                tu
 Kinetische Energie  (kWh/kg)         vt*vt/7.2e6
 Potentielle Energie (kWh/kg)         mu*(1/re-1/r2)/3.6e6 = (1-re/r2)*100 %
 Anziehung                                        mu/r2/r2 = re*re/r2/r2*100 %
 Restgeschwindigkeit                      vr
 Geschwindigkeitsgewinnfaktor     vr/(v1-sqrt(2*mu/r1))

5 Bahnbestimmung aus Ortsgrößen:

Unterprogramm zur Ermittlung charakteristischer Größen einer Umlaufbahn um
einen Himmelskörper unter Vorgabe des aktuellen Flugzustands. Des weiteren
wird der Einfluss kleiner Störungen berechnet.
Eingaben:
 Ortshöhe (km)                     hr
 Ortsgeschwindigkeit           vt
 Ortstangentenwinkel (90-180)  psi (Winkel von 0 bis 90 Grad können aus
                                   Symmetriegründen transferiert werden)
Berechnungen:
 r  = hr*1000+re
 gh = r*mu/(2*mu-vt*vt*r)
 Numerische Lösung von fi(r,gh,psi):
  gh = gh/sqr(gcos(fi)+gtan(psi)*gsin(fi))-r
       +gcos(fi)*r/(gcos(fi)+gtan(psi)*gsin(fi))
 ex = abs(-1/(gcos(fi)+gtan(psi)*gsin(fi)))
 pp = r*(1+ex*gcos(fi))
 h1 = (pp/(1+ex)-re)/1000
 h2 = (pp/(1-ex)-re)/1000
 v1 = sqrt(mu*(2/(h1*1000+re)-1/gh))
Ausgaben:
 Exzentrizität                                     ex
 Apsidenwinkel                                 fi
 Peri-höhe                                         h1 km
 Apo-höhe                                          h2 km
 Peri-geschwindigkeit                      v1
 Störungseinfluss von +1 m/s (Näherung) tangential bzw. bahnnormal
 absolute Änderung der großen Halbachse    2/vt*gh
 relative Änderung der Umlaufzeit                      300/vt %
 absolute Änderung der Inklination                     180/vt/pi Grad
 absolute Änderung der Exzentrizität                  2/vt

6 Bahnbestimmung aus Flugdauer:

Unterprogramm zur Ermittlung charakteristischer Größen einer Flugbahn um einen
Himmelskörper unter Vorgabe des Zeitbedarfs, um z.B. einen Bahnwechsel
durchzuführen. Dabei wird auch hier die Dauer des notwendigen Antriebsmanövers
vernachlässigt.
Eingaben:
 Starthöhe (km)   h1  (der Ursprungskreisbahn)
 Zielhöhe (km)     hr
 Zeitbedarf           tu
Berechnungen:
 Anfangswerte:
  r1 = h1*1000+re
  r2 = hr*1000+re
  gh = (r1+r2)/2
  tm = pi*sqrt(gh*gh/mu*gh)
  v  = sqrt(mu*(2/r1-2/(r1+r2)))
 Schleife mit 30 Durchläufen (variiert wird v; t geht gegen 0):
  si = r1*v
  hh = v*v-2*mu/r1
  ex = sqrt(1+hh*si/mu*si/mu)
  pp = r1*r1*v*v/mu
  f  = (pp-r2)/ex/r2
  ph = pi/2-rasn(f)
  Für ex>1:
   gh = pp/(ex*ex-1)
   n  = hoch((pp/gh),1.5)/sqrt(mu/gh/gh/gh)
   t5 = (ln(abs(((ex-1)*rtan(ph/2)+sqrt(ex*ex-1))/((ex-1)*rtan(ph/2)
        -sqrt(ex*ex-1)))))/sqrt(ex*ex-1)
   t  = n*(ex*sin(ph)/(ex*ex-1)/(1+ex*cos(ph))-t5/(ex*ex-1))-tu
  Für ex<1:
   gh = pp/(1-ex*ex)
   n  = hoch((pp/gh),1.5)/sqrt(mu/gh/gh/gh)
   t5 = 2/sqrt(1-ex*ex)*arctan((1-ex)*rtan(ph/2)/sqrt(1-ex*ex))
   t  = n*(ex*sin(ph)/(ex*ex-1)/(1+ex*cos(ph))-t5/(ex*ex-1))-tu
  h2  = (gh*(1+ex)-re)/1000
  vt  = sqrt(hh+2*mu/r2)
  ps  = pi-arctan((1+ex*cos(ph))/ex/sin(ph))
Ausgaben:
 Startgeschwindigkeit           v
 Zeitdifferenz                          t (entspricht Rechenfehler)
 Exzentrizität                          ex
 Apsidenwinkel                      ph*180/pi
 Apo-höhe (km)                      h2
 Ortsgeschwindigkeit            vt
 Ortstangentenwinkel            ps*180/pi

7 Swingby-Manöver:

Unterprogramm zur Ermittlung charakteristischer Größen eines Swingby-Manövers
(gravity assist, fly-by).
Dabei wird ein ebener Vorbeiflug an einem Himmelskörper ohne Antriebsmanöver
und ohne Störung durch andere Gravitationsfelder vorausgesetzt. Als Nahbereich
wird eine Ortshöhe unterhalb der Apsidenhöhe festgelegt. Es wird außerdem ein
Koordinatensystem wie folgt festgelegt:
x-Richtung : Zielkörperflugrichtung
y-Richtung : in Zielkörperbahnebene senkrecht zu x nach außen
z-Richtung : auf Zielkörperbahnebene senkrecht nach oben');
Für Anflugbreiten über 90 Grad kommt es zu einer Beschleunigung, bei
Anflugbreiten unter 90 Grad zu einer Verzögerung der Satellitengeschwindig-
keit bzgl. dem gemeinsamen Gravitationszentrum (Hinterher- bzw. Vorausflug).
Eingaben:
 Anflugwinkel zu beider Gravitationszentrum (|fi|<=180)  f1
 Anflugortsgeschwindigkeit                                                v1
 Zielkörperflugwinkel zum Gravitationszentrum                fpl
 Zielkörpergeschwindigkeit                                                vpl
 Peri-höhe (km)                                                                   hp1
 Anflugbreitengrad (|in|<=180)                                           ink
Berechnungen:
 pr :=re+hp1*1000
 al :=f1-fpl
 vr :=sqrt(sqr(vpl)+sqr(v1)-2*v1*vpl*gcos(al))
 b  :=pr*sqrt(1+2*mu/pr/vr/vr)
 vp :=b/pr*vr
 vk :=sqrt(mu/pr)
 phi:=180-2*gatn(vr*vr*b/mu)
 bg :=gasn(gsin(phi)*gsin(-ink))
 Für ((gsin(-ink)>gsin(phi)) und (phi>90)): bg:=180-bg
 lg :=gacs(gcos(phi)/gcos(bg))
 x1 :=v1*gcos(al)-vpl
 y1 :=v1*gsin(al)
 x  :=(sqr(vpl)+sqr(vr)-sqr(v1))/(2*vpl*vr)
 Für ((abs(round(x*1e7))=1e7) und (x1>0)): bt:=180
 sonst für (abs(round(x*1e7))=1e7):        bt:=0
 sonst:                                    bt:=gacs(x)*sgn(y1)
 x2 :=-vr*gcos(bg)*gcos(lg+bt)
 y2 := vr*gcos(bg)*gsin(lg+bt)
 z2 := vr*gsin(bg)
 v2 :=sqrt(sqr(vpl+x2)+sqr(y2)+sqr(z2))
 psi:=gasn(z2/v2)
 xi :=gatn2(vpl+x2,y2)
 Für (vpl+x2)>0: f2 :=gacs(gcos(xi+fpl)*gcos(psi))
 sonst:          f2 :=-gacs(gcos(xi+fpl)*gcos(psi))
 dv :=sqrt(sqr(x2-x1)+sqr(y2-y1)+sqr(z2))
 si :=pr*vp
 hh :=vp*vp-2*mu/pr
 ex :=sqrt(1+hh/mu*si/mu*si)
 pp :=(1+ex)*pr
 gh :=pp/(ex*ex-1)
 t2 :=ln(abs(((ex-1)+sqrt(ex*ex-1))/((ex-1)-sqrt(ex*ex-1))))/sqrt(ex*ex-1)
 tu :=2*(ex/(ex*ex-1)-t2/(ex*ex-1))*sqrt((pp/mu*pp*pp))
 df :=f2-f1
 Für abs(df)>180: df:=360-abs(df)
Ausgaben:
 Relativgeschwindigkeit vorher (|v|,x,y)           vr  x1  y1
 Relativgeschwindigkeit nachher (|v|,x,y,z)      vr  x2  y2  z2
 Ortsgeschwindigkeit nachher                         v2
 Peri-geschwindigkeit                                      vp
 Kreisbahngeschwindigkeitsbedarf               vk-vp
 Exzentrizität                                                     ex
 Apsidenhöhe                                                  (pp-re)/1000 km
 Unbeeinflusste Vorbeiflughöhe                    (b-re)/1000 km
 Umlenkwinkel                                                 phi
 Flugwinkeländerung                                      df
 Abflugwinkel zum Zentrum                            f2
 Abflugneigungswinkel                                    psi
 Nahbereichaufenthaltsdauer                          tu
 Kinetischer Energiefaktor                               v2*v2/v1/v1 entspricht
  (sqr(v2)-sqr(v1))/2/sqrt(mu/re)/vpl*100 % des max. Energiegewinns
  oder v2-v1 m/s Geschwindigkeitsänderung
 Äquivalenter Antriebsbedarf                         dv
  entspricht dv/sqrt(mu/re)*100 % vom max. dv des Himmelskörpers
  oder dv/sqrt(mu/pr)*100 % vom max. dv bzgl. Peri-höhe

8 Bodenspur zeichnen:

Dieses Unterprogramm zeichnet für einen bestimmten Zeitraum die Bodenspur
einer geschlossenen Umlaufbahn um einen Himmelskörper. Für die Erde wird
zusätzlich eine geografische Karte gezeichnet. Die Satellitenspur kann
wahlweise punktiert oder als Polygonzug gezeichnet werden. Falls die
verwendete Grafik nicht lauffähig sein sollte, wird eine reduzierte
Grobgrafik aus ASCII-Zeichen initialisiert.
Eingaben:
 große Halbachse (km)                    grha
 Exzentrizität                                       ex
 Winkelabstand des Perizentrums   om
 Inklination                                           in1
 Rektaszension                                   rek
 Apsidenwinkel von (>=0)                 ff0
 Apsidenwinkel bis                            fm
 Apsidenschrittweite                          fi
Berechnungen:
 Die Berechnung läuft ohne Zutun von ff0 bis fm mit der Schrittweite fi
 Anfangswerte:
  lr   = rek/180*pi
  ink  = in1/180*pi
  ome  = om/180*pi
  beg  = round(ff0/fi)
  en   = round(fm/fi)
  tges = 2*pi*sqrt(gh*gh*gh/mu)
  f1   = 2*arctan(sqrt((1-ex)/(1+ex))*gtan(ff0/2))
  f2   = ex*sqrt(1-ex*ex)*gsin(ff0)/(1+ex*gcos(ff0))
  t2   = sqrt(gh*gh*gh/mu)*(f1-f2)
 Laufwerte:
  ip  = ip+1
  phi = fi*ip
  ph  = phi/180*pi
  ump = trunc(phi/360)
  bg  = rasn(sin(ink)*sin(ome+ph))
  l1  = racs(cos(ome+ph)/cos(bg))
  f1  = 2*arctan(sqrt((1-ex)/(1+ex))*rtan(ph/2))
  f2  = ex*sqrt(1-ex*ex)*sin(ph)/(1+ex*cos(ph))
  tu  = sqrt(gh*gh*gh/mu)*(f1-f2)
  tu  = tu+tges*ump
  lg  = l1-2*pi/lv[vc]*(tu-t2)+lr
  fp  = sqrt(mu/hoch((gh*(1-ex*ex)),3))*sqr(1+ex*cos(ph))
  vg  = sqrt(mu*(1+2*ex*cos(ph)+ex*ex)/(gh*(1-ex*ex)))
Ausgaben:
 Die Ausgaben werden für jeden Durchlauf aktualisiert und umfasst:
 Zeit                                                   tu
 Höhe (km)                                       (gh*(1-ex*ex)/(1+ex*cos(ph))-re)/1000
 Apsidenwinkel                                 ph*180/pi
 Breitengrad                                     bg*180/pi
 Längengrad                                     lg*180/pi
 Drehgeschwindigkeit (Grad/min)  fp*180/pi*60
 Geschwindigkeit                             vg
 Umlauf                                              ump

9 Seile:

Unterprogramm zur Ermittlung spezifischer Größen seilgefesselter Satelliten im
Einflussbereich eines Gravitationsfeldes mit Endmassekörpern.
Eingaben:
 Innenhöhe (km)                       h1 (des Seilendmassekörpers)
 Außenhöhe (km)                    h2
 Innenmasse                            m1
 Außenmasse                         m2
 Seiltragfähigkeit (N/mm^2)  sig (zulässiger Wert)
 Seildichte (g/cm^3)               dic
Berechnungen:
 Da die Masse des Seiles als Ergebnis der Berechnung einen z.T.
 beträchtlichen Einfluss auf die Berechnung selbst haben kann, wird die
 Lösung auf numerischem Wege in maximal 19 Schritten von oben  nach unten
 wie folgt ermittelt:
  r1  = re+1000*h1
  r2  = re+1000*h2
  rcm = (qi/2*(r2*r2-r1*r1)+m1*r1+m2*r2)/(qi*(r2-r1)+m1+m2)
  rcg = sqrt((qi*(r2-r1)+m1+m2)/(qi*(1/r1-1/r2)+m1/r1/r1+m2/r2/r2))
  rmz = hoch(((qi/2*(r2*r2-r1*r1)+m1*r1+m2*r2)/
             (qi*(1/r1-1/r2)+m1/r1/r1+m2/r2/r2)),(1/3))
  dga = mu*(1/r2/r2-r2/rmz/rmz/rmz)
  dgi = mu*(1/r1/r1-r1/rmz/rmz/rmz)
  kr  = max(dgi*m1,abs(dga*m2))
  kr  = kr+dic*1000*fl*mu*(1/r1-3/2/rmz+r1*r1/2/rmz/rmz/rmz)
  fl  = kr/1000000/sig
  ms  = fl*(r2-r1)*dic*1000
  om  = sqrt(mu/rmz/rmz/rmz)
  m   = min(m1,m2)
  dv  = sqrt(mu/r1)+om*(r2-r1)-sqrt(mu/r2)
  tu  = sqrt(4*pi*pi/mu*rmz*rmz*rmz)
  qi  = dic*1000*fl
Ausgaben:
 Nummer der Näherung
 Massenzentrumshöhe                                    rcm-re
 Gravitationszentrumshöhe                             rcg-re
 Metazentrumshöhe                                         rmz-re
 Abstand Massen- zu Gravitationszentrum   rcm-rcg
 Innenbeschleunigung                                      dgi
 Außenbeschleunigung                                  dga
 Seilzugkraft innen                                            dgi*m1
 Seilzugkraft außen                                         abs(dga)*m2
 Seilzugkraft im Metazentrum                           kr
 Seilmasse                                                         ms = ms/(m1+m2+ms)*100 %
 Seilmasse bzgl. kleinerer Endmasse            ms/m*100 %
 Konstanter Seildurchmesser                          2*sqrt(fl/pi)
 Transfergeschwindigkeitsbedarf                    dv
 Umlaufzeit                                                          tu
 Innengeschwindigkeit                                       om*r1
 Außengeschwindigkeit                                   om*r2
 Kreisbahngeschwindigkeitsbedarf innen       sqrt(mu/r1)-om*r1
 Kreisbahngeschwindigkeitsbedarf außen    sqrt(mu/r2)-om*r2

a Planetendaten:

Unterprogramm zur Ermittlung spezifischer Daten zu den verfügbaren
Himmelskörpern. Es kann eine Tabelle der wichtigsten Daten aller
Himmelskörper oder eines ganz bestimmten angefordert werden.
Eingaben:
 Himmelskörper       vc
Berechnungen:
 ip   = Zuweisung des Gravitationszentrums
 t0   = 2*rv[ip]/tv[vc]*180/pi*3600
 t1   = 2*rv[vc]/tv[vc]*180/pi*3600
 tt   = tb*lv[vc]/(max(lv[vc],tb)-min(lv[vc],tb))
 agz  = -2*gr*mv[vc]*rv[ip]/sqr(tv[vc])/tv[vc]
 ts1  = tb*31556926/(max(tb,31556926)-min(tb,31556926))
 ts2  = tb*lv[ip]/(max(tb,lv[ip])-min(tb,lv[ip]))
 me   = mv[ip]
 mm   = mv[vc]
 rges = tv[vc]
 r3   = rges
 r2   = rges
 r1   = rges
 Numerische Berechnung der Lagrange-Punkte L1,L2 und L3 mit 22 Durchläufen
  als Funktionen r1(l1=0), r2(l2=0) und r3(l3=0):
  l1 =me*(1-hoch((rges-r1)/rges,3))
      -mm*(sqr((rges-r1)/r1)+hoch((rges-r1)/rges,3))
  l2 =me*(1-hoch((rges+r2)/rges,3))
      +mm*(sqr((rges+r2)/r2)-hoch((rges+r2)/rges,3))
  l3 =me*(1-hoch(r3/rges,3))+mm*(sqr(r3/(rges+r3))-hoch(r3/rges,3))
 rm   = mv[vc]/(mv[ip]+mv[vc])*tv[vc]
 akts = tv[vc]*hoch((mv[vc]/mv[ip]/sqrt(2)),0.4)-rv[vc]
 neus = tv[vc]*hoch((mv[vc]/mv[ip]),0.5)-rv[vc]
 kos3 = sqrt(sqr(sqrt(2)-1)*gr*mv[ip]/tv[vc]+2*gr*mv[vc]/rv[vc])
 kos4 = sqrt(sqr(sqrt(sqr(sqrt(2)-1)*gr*mv[0]/tv[ip]
        +2*gr*mv[ip]/tv[vc])-sqrt(gr*mv[ip]/tv[vc]))+2*gr*mv[vc]/rv[vc])
Ausgaben:
 Gravitationszentrum                                    pls[ip]
 Masse                                                          mv[vc]
 Erdmassen                                                  mv[vc]/mv[3]
 Mittlerer Radius (km)                                  rv[vc]/1000
 Umfang (km)                                               2*pi*rv[vc]/1000
 Oberfläche (qkm)                                        4*pi*rv[vc]*rv[vc]/1e6
 Volumen (km^3)                                          4/3*pi*sqr(rv[vc])*rv[vc]/1e9
 Erdradien                                                     rv[vc]/rv[3]
 Erdoberflächen                                            sqr(rv[vc])/sqr(rv[3])
 Erdvolumina                                                 hoch(rv[vc],3)/hoch(rv[3],3)
 Scheinbarer mittl. Durchmesser des Gravitationszentrums                          t0
 Scheinbarer mittl. Himmelskörperdurchmesser auf Gravitationszentrum    t1
 Große Bahnhalbachse                            tv[vc]
 Gravitationszentrumsradien                     tv[vc]/rv[ip]
 Mittlere Dichte                                          mv[vc]/(4/3*pi*hoch(rv[vc],3))
 Bahnexzentrizität                                     ev[vc]
 Rotationsdauer                                         lv[vc]
 max. Rotationsgeschwindigkeit              2*pi*rv[vc]/lv[vc]
 max.Störpotential                                     -2/3/abv[vc]*100 %
 Abplattung                                                 1 : abv[vc]
 mittlere Oberflächenanziehung               gr*mv[vc]/sqr(rv[vc]
 Gezeitenbeschleunigung auf Gravitationszentrum     agz
 Oberflächenkreisbahngeschwindigkeit bzw. maximale
  Geschwindigkeitsänderung  beim Swingby               sqrt(gr*mv[vc]/rv[vc])
 Fluchtgeschwindigkeit                             sqrt(2*gr*mv[vc]/rv[vc])
 3. kosmische Geschwindigkeit               kos3
 4. kosmische Geschwindigkeit               kos4
 Gesamtenergie                                        -gr*mv[ip]/2/tv[vc]/3.6e6 kWh/kg
  bzgl. Gravitationszentrum
 Bahnneigung zur Ekliptik                          nwv[vc] Grad
 Bahnumlaufdauer                                      2*pi*sqrt(sqr(tv[vc])/(mv[ip]+mv[vc])*tv[vc]/gr)
 Mittlere Bahngeschwindigkeit                  sqrt(gr*mv[ip]/tv[vc])
 Himmelskörpertag                                     tt
 Synodische Periode zur Erde                  ts1
 Synodische Periode zur Gravitationszentrumsoberfläche      ts2
 Charakteristische Größen des Mehrkörpersystems Gravitationszentrum-
 Himmelskörper
 Schwerpunktsabstand                               rm/rv[ip]*100 %
  des Gravitationszentrumsradius
 Aktivsphäre (Höhe)                                   akts/1000 km
 Neutralsphäre (Höhe)                               neus/1000 km
 Maximale Störbeschleunigung                 hoch(4*mv[vc]/mv[ip],0.2)*100 %
  der Gravitationszentrumsanziehung
 Lagrangehöhe L1                                       (r1-rv[vc])/1000 km
 Lagrangehöhe L2                                       (r2-rv[vc])/1000 km
 Lagrangehöhe L3                                       (r3-rv[vc])/1000 km
 Himmelskörperdurchmesser in L1          2*rv[vc]/r2*180/pi*3600 ''
 Himmelskörperdurchmesser in L2          2*rv[vc]/r1*180/pi*3600 ''
 Gravitationszentrumsdurchmesser in L1
                  2*rv[ip]/(tv[vc]-r2)*180/pi*3600 Bogensekunden
 Gravitationszentrumsdurchmesser in L2
                  2*rv[ip]/(tv[vc]+r1)*180/pi*3600 Bogensekunden

b Erdabplattungseinflüsse:

Unterprogramm zur Ermittlung des Einflusses der Erdabplattung auf die
Bahnelemente Rektaszension und Argument des Perigäums von Ellipsenbahnen.
Eingaben:
 große Halbachse (km)  gh
 Inklination                         ink
 Exzentrizität                     ex
Berechnungen:
 dl  = hoch((6371.2/gh),3.5)*gcos(ink)/sqr(1-sqr(ex))*9.96
 dom = hoch((6371.2/gh),3.5)*(5*sqr(gcos(ink))-1)/sqr(1-sqr(ex))*4.98
Ausgaben:
 Änderung der Rektaszension                                     dl  in Grad pro Tag
 Änderung des Winkelabstandes des Perigäums    dom in Grad pro Tag

c Manövergeschwindigkeitsbedarf:

Unterprogramm zur Ermittlung des Antriebsbedarfs bei Änderung einer
Umlaufbahn um z.B. ein Rendevous durchzuführen.
Eingaben:
 Winkeländerung in der Flugebene      fi
 Neigungsänderung der Flugebene     inl
 Ortsgeschwindigkeit 1                         v1
 Ortsgeschwindigkeit 2                         v2
Berechnungen:
  ar = gcos(fi)*gcos(inl)
  w2 = pi/2-arctan(ar/sqrt(1-sqr(ar)))
  vk = sqrt(v1*v1+v2*v2-2*v1*v2*cos(w2))
Ausgaben:
 Rendevousbedarf                  vk
 Winkeldifferenz                      w2*180/pi

d Bahnelemente ermitteln:

Unterprogramm zur Ermittlung der fehlenden Bahnelemente Rektaszension und
Argument des Perizentrums, um z.B. Bodenspuren zu zeichnen.
Eingaben:
 Positionsbreitengrad       bg1
 Positionslängengrad        lg1
 Apsidenwinkel                  apswi
 Inklination                          in1
Berechnungen:
 om  = gasn(gsin(bg1)/gsin(in1)) - apswi
 rek = lg1 - gacs(gcos(om+apswi)/gcos(bg1))*sgn(gsin(2*in1))*sgn(bg1)
Ausgabe:
 Argument des Perizentrums  om
 Rektaszension                        rek

Numerische Unterprogramme:

Die nun folgenden numerischen Unterprogramme stellen Laufzeitprogramme dar.
Zur Programmsteuerung kann mit bestimmten Tasten Einfluss auf den Programmlauf
genommen werden.

e Raketengröße für Einstufenträger:

Unterprogramm zur Ermittlung der charakteristischen Größen einer einstufigen
Rakete unter Vorgabe des Antriebsbedarfs.
Eingaben:
 Geschwindigkeitsbedarf      ve
 Ausströmgeschwindigkeit   va
 Stufenmassenverhältnis      cc
 Nutzlast                                  nl
Berechnungen:
 sl = (exp(ve/va)*nl-nl)/(cc-exp(ve/va))
 s1 = sl+nl
 w1 = sqr(ve/va)/(exp(ve/va)-1)
 d  = cc*sl+nl
Ausgaben:
 Trockenmasse                                 s1
 Gesamtmasse                                 d
 Treibstoffmasse                              d-s1
 Gesamtmassenverhältnis              d/s1
 Nutzlastfaktor                                   nl/d*100  %
 Mechanischer Wirkungsgrad        sqr(ln(1+d/s1))/(d/s1)*100  %
 Mittlerer Vortriebswirkungsgrad    w1*100  %

f Raketengeschwindigkeitsvermögen:

Unterprogramm zur Ermittlung des Geschwindigkeitsvermögens von Trägersystemen
mit bis zu fünf Stufen.
Eingaben:
 Stufenanzahl (1-5)                                   sz
 Masse der vollen Stufen                         o1 bis o5
 Masse der leeren Stufen                         l1 bis l5
 Ausströmgeschwindigkeit der Stufen   v1 bis v5
 Masse der Nutzlast                                  nl
Berechnungen:
 mv1 = (o1+o2+o3+o4+o5+nl)/(l1+l2+l3+l4+l5+nl)
 g1  = v1*ln((o5+o4+o3+o2+o1+nl)/(l1+o5+o4+o3+o2+nl))
 g2  = v2*ln((o5+o4+o3+o2+nl)/(l2+o5+o4+o3+nl))
 g3  = v3*ln((o5+o4+o3+nl)/(l3+o5+o4+nl))
 g4  = v4*ln((o5+o4+nl)/(l4+o5+nl))
 g5  = v5*ln((o5+nl)/(l5+nl))
Ausgaben:
 Massenverhältnis                        mv1
 Nutzlastfaktor                              100*nl/(o1+o2+o3+o4+o5+nl) %
 Mechanischer Wirkungsgrad    sqr(ln(1+mv1))/mv1*100 %
 Geschwindigkeitsvermögen
  einstufig                                      g1
  zweistufig                                    g1+g2
  dreistufig                                     g1+g2+g3
  vierstufig                                      g1+g2+g3+g4
  fünffstufig                                     g1+g2+g3+g4+g5

g Raumschlepper:

Unterprogramm zur Ermittlung der Nutzlasttragfähigkeit eines einstufigen
Orbittransferfahrzeugs.
Eingaben:
 Ausströmgeschwindigkeit                     va (der Treibstoffe)
 einfacher Geschwindigkeitsbedarf       ve (für Orbitwechsel bzw. Rendevous)
 Schleppermasse voll                              vm (ohne Nutzlast)
 Schleppermasse leer                             tm (ohne Nutzlast)
Berechnungen:
 no  = (exp(ve/va)*tm-vm)/(1-exp(ve/va))
 nl  = (exp(2*ve/va)*tm-vm)/(1-exp(ve/va))
 ml  = exp(ve/va)*tm
 ql  = (vm+nl)/(tm+nl)
 znr = (exp(ve/va)*tm-vm/exp(ve/va))/(1-exp(ve/va))
 mr  = znr+vm/exp(ve/va)
 qvr = (vm+znr)/(znr+tm)
 nv  = (exp(2*ve/va)*tm-vm)/(1-exp(2*ve/va))
 mv1 = (nv+tm)*exp(ve/va)
 qv  = (vm+nv)/(nv+tm)
Ausgaben:
 Nutzlast ohne Rückflug                           no = no/(no+vm)*100 %
  Startmassenverhältnis                          (vm+no)/(tm+no)
  Geschwindigkeitsbedarf für Rückflug   va*ln(vm/tm)/2
 Nutzlast bei Leerrückflug                         nl = nl/(nl+vm)*100 %
  Masse vor Leerrückflug                          ml
  Treibstoffverbrauch hin,zurück               vm-ml  bzw.  ml-tm
  Startmassenverhältnis                            ql
  Mechanischer Wirkungsgrad                sqr(ln(1+ql))/ql*100 %
 Nutzlast bei Leerhinflug                           znr = znr/(znr+vm)*100 %
  Masse vor Vollrückflug                           mr
  Treibstoffverbrauch hin,zurück               vm-mr+znr  bzw.  mr-tm-znr
  Startmassenverhältnis                            qvr
  Mechanischer Wirkungsgrad                sqr(ln(1+qvr))/qvr*100 %
 Nutzlast ohne Leerflug                             nv = nv/(nv+vm)*100 %
  Masse vor Vollrückflug                            mv1
  Treibstoffverbrauch hin,zurück              vm-mv1+nv  bzw.  mv1-tm-nv
  Startmassenverhältnis                           qv
  Mechanischer Wirkungsgrad                sqr(ln(1+qv))/qv*100 %

h Raumfahrtantriebe:

Unterprogramm zur Ermittlung charakteristischer Größen von
Raumfahrt-Antrieben. Dabei wird differenziert nach Leistungsantrieben,
Schubantrieben und Sonnensegeln.
Für Leistungs- und Schubantriebe:
 Eingaben:
  Schub                                    kr 
        oder   Motorleistung       l
  Ausströmgeschwindigkeit  va
  Vollmasse                            m1
  Treibstoffmasse                  mt
 Berechnungen:
  ds = kr/va                   oder   ds = 2*l/sqr(va)
  l  = ds/2*sqr(va)            oder   kr = ds*va
  vv = va*ln(m1/m2)
  b1 = kr/m1
  b2 = kr/m2
 Ausgaben:
  Durchsatz        ds                              kg/s
                           ds*3600                    kg/h
                           ds*365*24*3600      kg/a
  Vorrat               mt/(ds*365*24*3600) Jahre
                            mt/(ds*3600)                   Stunden
  Leistung            l                            Watt
  Schub               kr                          Newton
  Anfangsbeschleunigung    b1                             m/sec^2
                                              365*24*3600*b1      m/sec pro a
                                              3600*b1                    m/sec pro h
  Endbeschleunigung           b2                             m/sec^2
                                              365*24*3600*b2      m/sec pro a
                                              3600*b2                    m/sec pro h
  Geschwindigkeitsvermögen         vv                  m/sec
Für Sonnensegel:
 Eingaben:
  Entfernung zur Sonne (km)                 hs
  Segelflächendichte (>=5 g/m^2)        sf
 Berechnungen:
  sl = 1370*sqr((696200+149.6e6)/(696200+hs))
  k  = 1.8*sl/299792458
  be = k/sf*1000
 Ausgaben:
  Strahlungsleistung                                       sl                            W/m^2
  Strahlungsdruck bei 90%-Reflexion          k*1e6                     N/qkm
  Beschleunigung                                           be                           m/sec^2
  Geschwindigkeitsvermögen pro Stunde   be*3600                m/sec pro h
             "             pro Tag                                be*3600*24          m/sec pro d
             "             pro Jahr                               be*3600*24*365  m/sec pro a

i Freiflug im Erde-Mond-System:

Dieses Unterprogramm simuliert einen Freiflug im Erde-Mond-System mit einem
Start aus einer kreisförmigen Umlaufbahn um die Erde. Mit einem beliebigen
Tastendruck kann die Ausführung unterbrochen werden. Es bestehen keine
Einflussmöglichkeiten (z.B. Antriebsmanöver). Bei extremen Vorgaben besteht
die Möglichkeit des Absturzes auf Erde oder Mond. Auch ein Entweichen aus
dem System ist berücksichtigt. Das Programm läuft mit Systemgeschwindigkeit.
Eingaben:
 Startgeschwindigkeit                                              v9
 Startkreisbahnhöhe                                                h9
 Abflugwinkel (0 <=> zwischen Erde und Mond)   ph9
  (in Gegen-Uhrzeigersinn)
 Anfangszeitschrittweite (10-1000 sec)                 d9
  (zur zeitlichen Steuerung)
Berechnungen:
 Anfangswerte:
  dd  = d9
  phi = ph9
  ho  = h9
  vo  = v9*100
  me  = 3.98546e20
  mo  = 4.9e18
  re  = 6.378e8
  rr  = re+ho*100000
  rm  = 3.844e10
  rx  = 7.688e10
  rl  = 3.476e8
  pp  = 2.66162e-6
  phi = phi/180*pi
  sinphi = sin(phi)
  cosphi = cos(phi)
  v2  = vo*cosphi
  vo  = -vo*sinphi
  x   = rr*cosphi
  y   = rr*sinphi
  t   = 0
  dt  = 100
  xm  = rm
  ym  = 0
  de  = rr
  dm  = sqrt(sqr(x-xm)+sqr(y-ym))
  r2  = dm
  de3 = sqr(de)*de
  dm3 = sqr(dm)*dm
  xb  = -me*x/de3-mo*(x-xm)/dm3
  yb  = -me*y/de3-mo*(y-ym)/dm3
 Schleife für aktuelle Zeit, Geschwindigkeit und Entfernung:
  xv = 0.5*(xb+xb)*dt
  yv = 0.5*(yb+yb)*dt
  v1 = vo+xv
  v3 = v2+yv
  dx = 0.5*(v1+vo)*dt
  dy = 0.5*(v3+v2)*dt
  x1 = dx+x
  y1 = dy+y
  r1 = sqrt(sqr(x1)+sqr(y1))
  Für r1<r2:
   dt = dd*r1/re
  sonst
   dt = dd*r2/rl
  t  = t+dt
  xm = rm*cos(pp*t)
  ym = rm*sin(pp*t)
  r2 = sqrt(sqr(x-xm)+sqr(y-ym))
  r4 = sqr(r1)*r1
  r5 = sqr(r2)*r2
  xc = -me*x1/r4-mo*(x1-xm)/r5
  yc = -me*y1/r4-mo*(y1-ym)/r5
  xv = 0.5*(xc+xb)*dt
  yv = 0.5*(yc+yb)*dt
  v1 = vo+xv
  v3 = v2+yv
  dx = 0.5*(v1+vo)*dt
  dy = 0.5*(v3+v2)*dt
  x1 = x+dx
  y1 = y+dy
  v  = sqrt(sqr(v1)+sqr(v3))
  ti = t/3600
  vi = v/100000
  ri = r1/100000
  rj = r2/100000
  xb = xc
  yb = yc
  vo = v1
  v2 = v3
  x  = x1
  y  = y1
Ausgaben:
 Zeit                                         t
 Geschwindigkeit                  v
 Entfernung zur Erde             r1-re
 Entfernung zum Mond          r2-rl

j Wiedereintrittsprogramm:

Dieses Unterprogramm simuliert einen Wiedereintrittskörper für Himmelskörper
mit und ohne Atmosphäre unter Annahme einfacher aerodynamischer und
kinematischer Zusammenhänge. Der freie Fall ohne Auftrieb aber mit
Luftwiderstand ist ebenfalls berechenbar.
Eingaben:
 Bodenatmosphärendruck                     p0
 Bodenatmosphärendichte                    rh
 Anfangsgeschwindigkeit                       v3
 Anfangsflugwinkel (zum Mittelpunkt)    w9
 Anfangshöhe                                           h0
 Masse                                                      mm
 Flügelfläche                                             ff
 Schrittweite pro Sekunde                      swp
Berechnungen:
 Die Berechnungen laufen in Real-Time oder mit Systemgeschwindigkeit.
 Mit den f/g/h/j-Tasten lässt sich der Anstellwinkel aw verändern.
 Anfangswerte:
  w1  = w9/180*pi
  wi3 = w1
  vh2 = sin(w1)*v3
  vv2 = cos(w1)*v3
  h4  = h0
  h3  = h0
  h2  = h0
 Schleife mit Laufvariable ep:
  g31 = gr*me/sqr(re+h3)
  g32 = vh2*vh2/(re+h3)
  g3  = g31-g32
  rh3 = rh*exp(-g31*h3*rh/p0)
  ca  = 0.05*aw+0.05           (Geradengleichung für Auftriebsbeiwert)
  cw  = 0.08/625*aw*aw+0.025   (Geradengleichung für Widerstandsbeiwert)
  bb  = cw/2*rh3*ff*v3*v3/mm
  cc  = ca/2*rh3*ff*v3*v3/mm
  vh3 = vh2-sin(wi3)*bb/swp+cos(wi3)*cc/swp
  vv3 = vv2+g31/swp-cos(wi3)*bb/swp-sin(wi3)*cc/swp
  vx  = vv2+g3/swp-cos(wi3)*bb/swp-sin(wi3)*cc/swp
  h3  = h2-vv3/swp
  h2  = sqrt(sqr(h3+re)+sqr(vh3/swp))-re
  v3  = sqrt(vv3*vv3+vh3*vh3)
  wi3 = arctan(vh3/vv3)+vh3/swp/(h2+re)
  vh2 = v3*sin(wi3)
  vv2 = v3*cos(wi3)
  v3  = sqrt(vv2*vv2+vh2*vh2)
  h4  = h4-vx/swp
  h2  = h4
  ghb = bb*mm*v3
  gb  = sqrt(bb*bb+cc*cc)
  ep  = ep+1
  e   = e+ep/swp
Ausgaben:
 Bremsbeschleunigung                  bb
 Auftriebsbeschleunigung              cc
 Anziehung                                       g31 - g32 = g3
 Lastvielfaches                                gb/9.81
 Maximum des Lastvielfachen und zugehörige Höhe
 Geschwindigkeit                            v3*3.6 km/h
 Horizontalgeschwindigkeit           vh2*3.6 km/h
 Sinkgeschwindigkeit                    vv2*3.6 km/h
 Höhe                                               h4
 Luftdruck                                         rh3
 Flugwinkel                                      wi3*180/pi
 Anstellwinkel                                  2.5*aw+2.5
 Bremsleistung                                ghb
 Maximum der Bremsleistung und zugehörige Höhe
 Zeit                                                  e

k Raketenstartprogramm:

Dieses Unterprogramm simuliert für beliebige Himmelskörper den Start eines
einstufigen Trägers. Dabei finden die Einflüsse der Eigendrehgeschwindigkeit
und des Atmosphärenwiderstands Berücksichtigung. Auch der Start aus einem
Orbit oder dem kräftefreien Weltraum ist möglich. Als Bezugspunkt für die
Endgeschwindigkeit ist der Mittelpunkt des Gravitationszentrums festgelegt.
Zur Erleichterung der Verwaltung der Eingabedaten besteht die Möglichkeit
diese Werte auf dem gültigen Laufwerk abzuspeichern bzw. einzulesen (in
Dateien mit der Extension ".RAK"). Dazu wird der COMMAND-Interpreter im
Root-Directory benötigt.
Die numerische Schrittweite beträgt eine Sekunde. Die Ausgabe erfolgt in
Real-Time oder in Systemzeit. Zur schnelleren Bearbeitung z.B. für
Parameterstudien kann die Ausgabe auf den Brennschluss beschränkt werden.
Eingaben:
 Steuervariable z3 für 1 Ausbrennen, sonst definierter Brennschluss
 Zielgeschwindigkeit: 1=Kreisbahngeschwindigkeit,
                      xxxxx=Endgeschwindigkeit         vj
 Bodendruck                                                       p0
 Bodendichte                                                      rh
 Startwinkel                                                         ws1
 Endwinkel                                                           wz
 Abflugrichtungswinkel (0=Osten)                      fr
 Starthöhe                                                             h0
 Breitengrad                                                           bg
 Raketenanfangsgeschwindigkeit bzgl.Startort   vv0
 Querschnittsfläche                                                   fl
 cw-Wert                                                                    cw
 Anfangsmasse                                                         ma
 Maximalbeschleunigung (*g)                                   gm
 Massenverhältnis (Bruch möglich)                          pp
 Richtungsänderungsfaktor (0-2)                              l
 Schubrichtungsunsicherheit (0.2-2 Grad)              su
 Düsenendluftwiderstandsfaktor (1-4)                     xx
 Vakuumausströmgeschwindigkeit                         va
 Brenndauer bei vollem Durchsatz                           bb
Berechnungen:
 Anfangswerte:
  wi   = ws1*pi/180
  we1  = wz*pi/180
  wb   = fr*pi/180
  ho   = h0
  w0   = wi-we1
  b3   = bg/180*pi
  abvc = abv[vc]
  aa3  = re*(1-1/abvc*bg/90)
  Für Erde:
   ge = 9.780501+(0.0178*bg/90)
  sonst
   ge = gr*me/aa3/aa3
  ve   = 2*pi*cos(b3)*aa3/te*cos(wb)
  vr   = vv0
  v1   = vr
  vw   = cos(wi)*vr
  m1   = ma
  m9   = ma
  gu   = gm*9.80665
  nl   = ma/pp
  ds   = (ma-nl)/bb
  ad   = ds
  v3   = va*(1-cos(su*pi/180))
 Schleife für die aktuelle Zeit:
  g1  = ge*sqr(aa3/(aa3+ho))
  g2  = g1-vw*vw/(aa3+ho)
  vg  = vg+g2
  v5  = v5+v3
  ds  = (gu+b1+b4)/va*m1
  wi  = wi-wa
  vo  = cos(wi)*ve
  v9  = vo-v8
  v7  = v7+v9
  v6  = v9
  v8  = vo
  m1  = m9-ds/2
  dc  = rh*exp(-ge*rh/p0*re*re*(1/re-1/(re+ho)))
  wa  = w0*l/(bb-ie)
  b1  = (cw/2*fl*dc*sqr((vr+v1)/2))/m1
  v1  = vr
  b4  = dc*xx
  be  = ds*va/m1-g2-b1-b4
  vr  = v1+be+v6-v3
  vw  = cos(wi)*(vr-v7)+ve-vn
  gz  = vw*vw/(aa3+ho)
  m9  = m1-ds/2
  vf  = (ds*va-b4*m9)/ds
  aaa = 2*vr/vf/(1+vr*vr/vf/vf)
  ho  = ho+be/2*sin(wi)+v1*sin(wi)
  bl  = bl+b1
  vv  = vv+b4
  ie  = ie+1
 Berechnungen nach Brennschluss:
  p   = ad/2*va*va
  bs1 = va*bb*(1-ln(pp)/(pp-1))
  wh  = gr*me*m9*(1/(re+h0)-1/(re+ho))
  wk  = m9/2*sqr(vr-v7-vv0)
  wt  = (ma-m9)/2*va*va
  nu  = (wh+wk)/wt*100
Ausgaben:
 Ausgaben in der Schleife:
  Zeit                                                    ie
  Anziehungsbeschleunigung           g1
  Flugwinkel                                        wi*180/pi
  Zentripedalbeschleunigung           gz
  Durchsatz                                        ds
  effekt.Ausströmgeschwindigkeit   vf
  Schubbeschleunigung                   be
  Schub                                              ds*va-b4*m9
  Bremsbeschleunigung                  b1
  Schubregelung                              ds/ad*100 %
  Geschwindigkeit                            vr
  Geschwindigkeitsverluste             vg+v5+vv+bl-v7
  Höhe                                                ho
  Luftdruck                                         dc/rh*p0
  Masse                                             m9
  Vortriebswirkungsgrad                 aaa*100 %
 Ausgaben nach Brennschluss:
  Anfangsdurchsatz                               ad
  Planetendrehgewinn                           v7
  Gravitationsverlust                              vg
  Luftwiderstandsverlust                        bl
  Lenkverlust                                          v5
  Triebwerksverlust                                vv
  Resttreibstoff                                       m9-ma/pp
  Kreisbahngeschwindigkeit                 sqrt(gr*me/(ho+re))
  Maximalleistung                                   p/1e6 MW
  ideale Brennschlusshöhe                    bs1/1000 km
  Maximalvakuumschub                         ad*va/1000 kN
  Geschw-vermögen (ohne Reserve)   vr+vg+v5+vv+bl-v7-vv0
  kinetischer Wirkungsgrad                   nu %
  mechan.Wirkungsgrad                        sqr(ln(1+pp))/pp*100 %

l Planetensystem:

Dieses Unterprogramm verdeutlicht die Größenverhältnisse des Planetensystems
bzgl. der Bahngröße und der Körpergröße und die relative Geschwindigkeit der
einzelnen Himmelskörper ohne Bezug auf die exakte Daten zu einem gegebenen
Zeitpunkt.

m Weltraumbahnhöfe:

Dieses Unterprogramm zeigt in einer Grafik die bedeutendsten Weltraumbahnhöfe
und deren zulässige Abschussrichtungen. Diese Darstellung veraltet im Gegensatz zu
den vorherigen Unterprogrammen mit fortschreitender Zeit. Der dargestellte Stand
kann dem Datum der Datei RFP.EXE grob entnommen werden.


Und hier liegt die Sammelprogramm RFP.EXE als ZIP-Version.



Zurück zum Anfang               Startseite