C---+----+----+----+----+----+----+----+----+----+----+----+----+----+--
PROGRAM TWO_BODY
!Program: Calculates the times, angles, speeds, and distances of
!the classical two-body gravitational problem for two arbitrary
!celestial masses in bounded elliptical orbits about their center
!of mass. The system solution (listed below in detail) is exact
!analytically but does not include any relativistic effects. All
!derived distances are measured relative to the center of mass
!for this two-body system, and all units are relative to solar
!numbers (i.e., "mass" is in solar units where one solar mass
!equals about 1.989E30 kg, and "distance" is in units of one
!"Astronomical Unit" (AU) which is about 1.496E11 m). Time = zero
!corresponds to "Perihelion" position (closest position between
!the two masses). Times, angles, and distances are printed to
!a table (both to screen and to an ASCII file "two_body.out")
!following user input of the above two masses and their time
!averaged separation distance (i.e., "semi-major axis" of their
!relative elliptical motion) in AU. A similar second table is
!printed to screen and the ASCII file "two-body.out" listing
!the times, angles, and speeds (in units meters per second).
!
!Derivation of the two-body system is listed below. The first
!part (STEP 1) is to determine total angular momentum about the
!center of mass and the central force equation for the two
!masses, and STEP 2 derives the time-independent equation of
!elliptical motion for the distance between the two masses.
!The third part (STEP 3) is to combine STEP 1 and STEP 2
!and solve for time as a function of mutual angle of revolution,
!and finally STEP 4 derives the speeds for the two masses and
!their relative elliptical system speed.
!
!
!
!
!STEP 1: Conservation of total angular momentum L about center
!of mass yields a simple equation for time differential dt:
!dt = mR^2/L dPHI, where m is reduced mass m1*m2/(m1+m2), R is
!distance between masses m1 and m2, and PHI is rotation angle
!(counter-clockwise angle in standard polar coordinates) which is
!taken to be zero at Perihelion (closest position over time
!between the two masses). The derivation of this relation for dt
!(and this two-bocy problem in general) follows by placing the
!vector coordinate system origin at the center of mass (c.m.)
!between the two masses:
!
! m1 c.m. m2
! o------+------------------------o
! <----- ----------------------->
! R1 R2
! R (R = R2 - R1)
! ------------------------------>
!
!Vectors R1 and R2 (distance vectors from c.m. to masses m1 and
!m2, respectively): R = R2 - R1, where vector R1 points
!opposite vectors R and R2. From center of mass it follows that
!m1*R1 + m2*R2 = 0 for vectors R1 and R2 (symbol "*" denotes
!simple scalar multiplication). Since R = R2 - R1, it follows
!that vectors R1 and R2 are
!
! R1 = -R * m2/(m1 + m2) and R2 = R * m1/(m1 + m2) (1)
!
!It is now noted that the central force equation for the two
!combined masses can be readily determined from evaluating the
!central force equation for either mass m1 or m2. For mass m1
!(same final relation follows for mass m2, and symbol "^" above
!any following parameter denotes "unit vector"):
!
! 2
! d R1 -G m1 m2 ^ G m1 m2 ^
! m1 ---- = --------- R1 = --------- R
! 2 R^2 R^2
! dt
!
!From (1) it follows that
!
! 2 2
! d R G m1 m2 ^ d R -G m1 m2 ^
! -m1*m2/(m1+m2) --- = --------- R ---> m --- = --------- R (2)
! 2 R^2 2 R^2
! dt dt
!
!where m is reduced mass m1*m2/(m1+m2).
!
!Total angular momentum vector L = R1 X m1*V1 + R2 X m2*V2,
!where "X" denotes cross product, V1 = dR1/dt, and V2 = dR2/dt.
!From (1), L becomes
! dR
! L = R X m -- (3)
! dt
!
!There are two components to dR/dt, one that is parallel to R and
!one that is perpendicular. The one that is parallel does not
!contribute to L because of the cross-product. The perpendicular
!component yields (here only magnitudes are implied for vector
!parameters)
! L = R * m * (R * dPHI/dt) (4)
!
!where the value of the perpendicular component for dR/dt is
!R * dPHI/dt. The result is that dt = mR^2/L dPHI as stated
!in the beginning STEP 1 paragraph. It is noted that the same
!equivalent expression for dt may be obtained using (1) for the
!masses m1 and m2 and their relative distances of R1 and R2.
!
!
!
!
!STEP 2: The time-independent equation of ellipical motion is as
!follows (PHI = zero corresponds to Perihelion):
! ed
! Equation of orbital mean distance R: R = -------------- (5)
! 1 + e cos(PHI)
!
!where e = eccentricity and "d" is distance (a constant) between
!the "focus" and "directrix" of the ellipse (with directrix
!located to the right of the "focus" in polar coordinates). It
!is noted that the numerator "ed" can be replaced by a*(1 - e^2)
!where "a" is semi-major axis of the two-body ellipse which
!is equivalent to the time-averaged mean distance between m1 and
!m2. The derivation of equation (5) begins with (2) which is
!re-written as
!
! 2
! d R - G m1 m2
! --- = --------- (cos(PHI), sin(PHI), 0) (6)
! 2 m R^2
! dt
! ^
!where (cos(PHI), sin(PHI), 0) is the unit vector R in polar
!coordinates representing standard cartesian coordinates x, y, z.
!
!The relative velocity vector dR/dt is then determinined by
!integrating (6) over time differential dt:
!
! dR - G m1 m2
! -- = INTEGRAL --------- (cos(PHI), sin(PHI), 0) dt (7)
! dt m R^2
!
!From conservation of total angular momentum, dt = mR^2/L dPHI
!which changes (7) to a simple integration over the mutual angle
!differential dPHI:
!
! dR - G m1 m2
! -- = INTEGRAL --------- (cos(PHI), sin(PHI), 0) dPHI
! dt L
!
!The integration yields (C is a vector constant of integration)
!
! dR - G m1 m2
! -- = --------- (sin(PHI), -cos(PHI), 0) + C (8)
! dt L
!
!Next, introduce boundary condition where at PHI = 0 the velocity
!dR/dt equals (0, Vp, 0) where Vp is the mutual relative speed
!at Perihelion:
!
! - G m1 m2
! (0, Vp, 0) = --------- (0, -1, 0) + C
! L
!
!It follows then that C = (0, Vp - G*m1*m2/L, 0). Substituting
!this vector constant into (8) yields the following for dR/dt:
!
! dR - G m1 m2 L Vp
! -- = --------- (sin(PHI), -cos(PHI) + 1 - -------, 0) (9)
! dt L G m1 m2
!
!At this point, the solution for distance R is readily solved by
!going back again to conservation of total angular momentum L
!which can be written
! ^ dR
! L = m R R X --
! dt
! ^
!where R is again the unit vector (cos(PHI), sin(PHI), 0).
!
!Hence from (9) and the above expression for L (symbols "|"
!denote magnitude of cross product),
!
! L ^ dR
! --- = | R X -- |
! mR dt
!
! = | (cos(PHI), sin(PHI), 0)
!
! - G m1 m2 L Vp
! X --------- (sin(PHI), -cos(PHI) + 1 - -------, 0) | (10)
! L G m1 m2
!
!The magnitude of the cross product then yields
!
! L G m1 m2 L Vp
! --- = ------- ( 1 + (------- - 1) cos(PHI) )
! mR L G m1 m2
!
!Solving for distance R gives
!
! L^2 / (G*m1*m2*m)
! R = --------------------------------- (11)
! 1 + (L*Vp/(G*m1*m2) - 1) cos(PHI)
!
!This is the equation of a conic section (ellipse, parabola,
!hyperbola) in polar coordinates and is equivalent to (5) above
!with numerator "ed" and eccentricity "e" replaced by the
!respective quantities in (11). That is, ed = L^2 / (G*m1*m2*m)
!and eccentricity e = L*Vp/(G*m1*m2)-1.
!
!
!
!
!STEP 3: Step 3 involves solving for time t. Again, from the
!conservation of total angular momentum, dt = m R^2/L dPHI.
!Using the expression for R (either (5) or (11), but for
!simplicity (5) is invoked), time is then solved by integrating
!the time differential dt in angle PHI:
!
! m(ed)^2 dPHI
! dt = ------- * ------------------
! L [1 + e cos(PHI)]^2
!
!Integration in PHI yields time "t" for bounded elliptical motion
!(i.e., 0 <= e < 1, and t = zero for Perihelion):
!
! m(ed)^2 e sin(PHI)
! t = ------- * [ --------------------------
! L (e^2 - 1) (1 + e cos(PHI))
!
! 2
! + ------------- arctan (sqrt((1-e)/(1+e)) * tan (PHI/2)) ]
! (1 - e^2)^1.5
!
! m(ed)^2 2 PI
!Note that lim(t) = ------- * [ ------------- * --- ]
! PHI --> PI L (1 - e^2)^1.5 2
!
!which corresponds exactly to 1/2 the total revolution time
!period T of the two masses m1 and m2. It then follows that
!
! m(ed)^2 (1 - e^2)^1.5 T
! ------- = ------------- * --- (12)
! L PI 2
!
!Hence, "normalized" relative time (i.e., t/T) is then the sum
!of the following two relative time terms t1 and t2 (i.e.,
!t/T = t1 + t2, with 0 <= t/T <= 1):
!
! -e sqrt(1 - e^2) sin(PHI) (13a)
! t1 = -------------------------
! 2 PI (1 + e cos(PHI))
!
! 1
! t2 = --- * arctan(sqrt((1-e)/(1+e)) * tan (PHI/2)) (13b)
! PI
!
!This program lists normalized time from (13) in the tables to
!both screen and the ASCII text file "two_body.out". The output
!beginning header in the text file also lists time period of
!revolution in seconds and Earth years to easily convert the
!normalized time to these units.
!
!The final part of STEP 3 is to determine the time period "T".
!With central force, total angular momemtum is conserved and so
!also then is the incremental area swept out per unit time. This
!follows by writing the incremental area dA for the two-body
!system as
! 1
! dA = --- | R X dR/dt | dt.
! 2
!
!It follows that dA/dt = L/(2*m), where L is the magnitude of
!total angular momentum (a constant). Since the area swept out
!per unit time is constant, the value for dA/dt is simply the
!total area A of the ellipse divided by the total revolution
!time period T: dA/dt = A/T. For an ellipse, the area is given
!by A = PI*a*b where b = a*sqrt(1 - e^2) is the "semi-minor"
!axis. It follows then that
!
! L/(2*m) = PI a^2 sqrt(1 - e^2)/T. (14)
!
!From (5) and (11) above, ed = L^2/(G m1 m2 m) = a (1 - e^2),
!so that
! sqrt(1 - e^2) = L / sqrt(G m1 m2 m a) (15)
!
!Combining (14) and (15) we eliminate the terms L and (1 - e^2),
!and the revolution time period T for the generalized two-body
!system is then given by
!
! T^2 = 4 PI^2 a^3 / G(m1 + m2). (16)
!
!
!
!
!STEP 4: As noted in the first paragraph, this program also
!prints a second table to screen and the output text file which
!lists the magnitudes of the individual velocities as a function
!of angle and time (i.e., |dR1/dt|, |dR2/dt|, and |dR/dt|).
!These speeds are listed with units of meters per second.
!
!To derive speeds, first note that in (9) 1 - L*Vp/(G*m1*m2) is
!identically the negative value of eccentricity (i.e., "-e"), so
!that (9) yields
!
! |dR/dt| = (G*m1*m2/L)*sqrt(sin(PHI)^2 + (cos(PHI)+e)^2) (17)
!
!From (12) and the relation ed = a*(1 - e^2), it follows that
!
! L = m*a^2 sqrt(1 - e^2) 2*PI / T (18)
!
!Noting that reduced mass m = m1*m2/(m1+m2) and from (16) that
!G(m1 + m2) = 4 PI^2 a^3 / T^2, the relative speed |dR/dt| in
!(17) can then be written
!
! 2*PI*a
! |dR/dt| = ---------------- * sqrt(sin(PHI)^2 + (cos(PHI)+e)^2)
! T*sqrt(1 - e^2)
! (19)
!
!Note that when e = 0 (i.e., pure circular orbits for m1 and m2
!about center of mass) (19) becomes simply the circumference of
!the circle of semi-major axis "radius a" divided by the total
!revolution period - that is, |dR/dt| = 2*PI*a / T.
!
!By center of mass (i.e., from (1)) the speeds for masses m1 and
!m2 are then calculated from (19) as follows:
!
! |dR1/dt| = |dR/dt| * m2/(m1+m2) (20a)
!
! |dR2/dt| = |dR/dt| * m1/(m1+m2) (20b)
!
!Author: Dr. Jerry R. Ziemke
!
PARAMETER(PI=3.141592653589793)
REAL INC,M1,M2
WRITE(*,*)'TWO-BODY GRAVITATIONAL PROGRAM:'
WRITE(*,*)'This program calculates the times, angles, speeds,'
WRITE(*,*)'and distances of the generalized non-relativistic'
WRITE(*,*)'two-body gravitational problem.'
WRITE(*,*)' '
WRITE(*,*)'Geometry:'
WRITE(*,*)' '
WRITE(*,*)' (AT TIME=ZERO)'
WRITE(*,*)' '
WRITE(*,*)' M1 cm M2'
WRITE(*,*)' O-----+----------------------o'
WRITE(*,*)' R1 R2'
WRITE(*,*)' <------------ R ----------->'
WRITE(*,*)' '
WRITE(*,*)'NOTES:'
WRITE(*,*)' '
WRITE(*,*)'M1 = First mass (arbitrary).'
WRITE(*,*)'M2 = Second mass (arbitrary).'
WRITE(*,*)'"cm" = Center of mass of two-body system.'
WRITE(*,*)'R1 = Instantaneous distance from cm to M1.'
WRITE(*,*)'R2 = Instantaneous distance from cm to M2.'
WRITE(*,*)'R = R1+R2 is instantaneous separation distance'
WRITE(*,*)'between the two masses.'
WRITE(*,*)' '
WRITE(*,*)'"Time" in output printed tables is in units of one'
WRITE(*,*)'revolution period of the system and is set to zero at'
WRITE(*,*)'the "Perihelion" position (i.e., the closest position'
WRITE(*,*)'between the two masses M1 and M2).'
WRITE(*,*)' '
WRITE(*,*)'"Angle" in output tables is in degrees and is the'
WRITE(*,*)'angle measured counterclockwise about cm of either M1'
WRITE(*,*)'or M2 (i.e., relative to horizontal line shown above).'
WRITE(*,*)'"Angle" in output tables is set to zero at time = 0'
WRITE(*,*)'which again corresponds to "Perihelion" position.'
WRITE(*,*)' '
WRITE(*,*)'"V1" and "V2" in output tables are the instantaneous'
WRITE(*,*)'revolution speeds of masses M1 and M2 relative to cm.'
WRITE(*,*)' '
WRITE(*,*)'All tabulated calculations are made using input'
WRITE(*,*)'masses and distances in relative Earth-Sun units'
WRITE(*,*)'(e.g., masses are in units of one Sun mass = 1.0'
WRITE(*,*)'corresponding to 1.989e30 kg, and distances are in'
WRITE(*,*)'units of 1.0 AU corresponding to 1.496e11 meters).'
WRITE(*,*)' '
WRITE(*,*)' '
WRITE(*,*)' '
WRITE(*,*)'Enter mass M1 (e.g., if Sun enter 1, if Earth enter'
WRITE(*,*)'3.002e-6, if Moon enter 3.694e-8):'
READ(*,*) M1
WRITE(*,*)'Enter mass M2 (e.g., if Sun enter 1, if Earth enter'
WRITE(*,*)'3.002e-6, if Moon enter 3.694e-8):'
READ(*,*) M2
WRITE(*,*)'Enter semi-major axis "a" (i.e., time-averaged mean'
WRITE(*,*)'distance between the two masses in AU - e.g., for'
WRITE(*,*)'the Earth-Sun system enter 1.0, for the Earth-Moon'
WRITE(*,*)'system enter 0.002567):'
READ(*,*) AX
WRITE(*,*)'Enter orbital eccentricity "e" between M1 and M2'
WRITE(*,*)'(e.g., for Earth-Sun system enter 0.0167, for'
WRITE(*,*)'Earth-Moon enter 0.0549):'
READ(*,*) E
WRITE(*,*)'Enter the increment for the angle sampling in the two'
WRITE(*,*)'output tables in degrees (e.g., 10 is a good number):'
READ(*,*) INC
N=IFIX(360.0/INC)
C=PI/180.0
TSEC=SQRT(AX*AX*AX/(M1+M2))*3.15576E7 !<== Rev Period (seconds)
TDAYS=TSEC/86400.0 !<== Rev Period (days)
TYR=TSEC/3.15576E7 !<== Rev Period (Earth years)
OPEN(12,FORM='FORMATTED',STATUS='UNKNOWN',FILE='two_body.out')
WRITE(*,*)' '
WRITE(*,*)'Fundamental parameters for your two-body system'
WRITE(*,*)'(1 AU = 1.496E11 m, one solar mass = 1.989E30 kg):'
WRITE(*,*)' '
WRITE(*,*)'Mass M1 (in solar mass units): ',M1
WRITE(*,*)'Mass M2 (in solar mass units): ',M2
WRITE(*,*)'Semi-major axis (in AU): ',AX
WRITE(*,*)'Eccentricity: ',E
WRITE(*,*)'Revolution period (seconds):',TSEC
WRITE(*,*)'Revolution period (Earth days):',TDAYS
WRITE(*,*)'Revolution period (Earth years):',TYR
WRITE(*,*)' '
WRITE(*,*)'Below is a table listing times, angles, and'
WRITE(*,*)'distances. R1 and R2 are instantaneous distances'
WRITE(*,*)'between the mutual center of mass and masses M1 and'
WRITE(*,*)'M2, respectively. Distance R is the total distance'
WRITE(*,*)'between M1 and M2 (i.e., R = R1 + R2). "Time" and'
WRITE(*,*)'"Angle" are both = ZERO at Perihelion:'
WRITE(*,*)' '
WRITE(*,*)' Time(period) Angle(deg) R1(AU)'//
& ' R2(AU) R(AU)'
WRITE(*,*)'-------------------------------------------'//
& '--------------------------------------'
WRITE(12,*)'TWO-BODY GRAVITATIONAL PROGRAM:'
WRITE(12,*)'This program calculates the times, angles, speeds,'
WRITE(12,*)'and distances of the generalized non-relativistic'
WRITE(12,*)'two-body gravitational problem.'
WRITE(12,*)' '
WRITE(12,*)'Geometry of the two-body gravitational system:'
WRITE(12,*)' '
WRITE(12,*)' (AT TIME=ZERO)'
WRITE(12,*)' '
WRITE(12,*)' M1 cm M2'
WRITE(12,*)' O-----+----------------------o'
WRITE(12,*)' R1 R2'
WRITE(12,*)' <------------ R ----------->'
WRITE(12,*)' '
WRITE(12,*)' '
WRITE(12,*)'Fundamental parameters for your two-body system'
WRITE(12,*)'(1 AU = 1.496E11 m, one solar mass = 1.989E30 kg):'
WRITE(12,*)' '
WRITE(12,*)'Mass M1 (in solar mass units): ',M1
WRITE(12,*)'Mass M2 (in solar mass units): ',M2
WRITE(12,*)'Semi-major axis (in AU): ',AX
WRITE(12,*)'Eccentricity: ',E
WRITE(12,*)'Revolution period (seconds):',TSEC
WRITE(12,*)'Revolution period (Earth days):',TDAYS
WRITE(12,*)'Revolution period (Earth years):',TYR
WRITE(12,*)' '
WRITE(12,*)'Below is a table listing times, angles, and'
WRITE(12,*)'distances. R1 and R2 are instantaneous distances'
WRITE(12,*)'between the mutual center of mass and masses M1 and'
WRITE(12,*)'M2, respectively. Distance R is the total distance'
WRITE(12,*)'between M1 and M2 (i.e., R = R1 + R2). "Time" and'
WRITE(12,*)'"Angle" are both = ZERO at Perihelion:'
WRITE(12,*)' '
WRITE(12,*)' Time(period) Angle(deg) R1(AU)'//
& ' R2(AU) R(AU)'
WRITE(12,*)'-------------------------------------------'//
& '--------------------------------------'
!Evaluate times, distances, and angles, and print out:
DO I=0,N
ANG=I*INC
IF (ANG.GT.180.0) PHI=(360-ANG)*C
IF (ANG.LE.180.0) PHI=ANG*C
IF (ANG.NE.180.0) THEN
T1=-E*SQRT(1-E*E)*SIN(PHI)/(1+E*COS(PHI))/2/PI
T2=ATAN(SQRT((1-E)/(1+E))*TAN(PHI/2))/PI
T=T1+T2
ELSE
T=0.5
ENDIF
R=AX*(1-E*E)/(1+E*COS(PHI))
R1=R*M2/(M1+M2)
R2=R*M1/(M1+M2)
IF (ANG.GT.180.0) T=1-T
WRITE(*,'(F14.10,F17.8,E17.8,E17.8,E17.8)') T,ANG,R1,R2,R
WRITE(12,'(F14.10,F17.8,E17.8,E17.8,E17.8)') T,ANG,R1,R2,R
ENDDO
!Evaluate times, speeds, and angles, and print out:
WRITE(*,*)' '
WRITE(*,*)'Below is a table listing times, angles, and speeds.'
WRITE(*,*)'V1 and V2 are instantaneous revolution speeds of'
WRITE(*,*)'masses M1 and M2, respectively, about their mutual'
WRITE(*,*)'center of mass. Speed V is the relative speed of'
WRITE(*,*)'revolution between M1 and M2 (i.e., V = |dR/dt|).'
WRITE(*,*)'Again, "Time" and "Angle" are both = ZERO at'
WRITE(*,*)'Perihelion:'
WRITE(*,*)' '
WRITE(*,*)' Time(period) Angle(deg) V1(m/s)'//
& ' V2(m/s) V(m/s)'
WRITE(*,*)'--------------------------------------------'//
& '-------------------------------------'
WRITE(12,*)' '
WRITE(12,*)'Below is a table listing times, angles, and speeds.'
WRITE(12,*)'V1 and V2 are instantaneous revolution speeds of'
WRITE(12,*)'masses M1 and M2, respectively, about their mutual'
WRITE(12,*)'center of mass. Speed V is the relative speed of'
WRITE(12,*)'revolution between M1 and M2 (i.e., V = |dR/dt|).'
WRITE(12,*)'Again, "Time" and "Angle" are both = ZERO at'
WRITE(12,*)'Perihelion:'
WRITE(12,*)' '
WRITE(12,*)' Time(period) Angle(deg) V1(m/s)'//
& ' V2(m/s) V(m/s)'
WRITE(12,*)'--------------------------------------------'//
& '-------------------------------------'
AXM=AX*1.496E11 !<== Semi-major axis (units meters)
DO I=0,N
ANG=I*INC
IF (ANG.GT.180.0) PHI=(360-ANG)*C
IF (ANG.LE.180.0) PHI=ANG*C
IF (ANG.NE.180.0) THEN
T1=-E*SQRT(1-E*E)*SIN(PHI)/(1+E*COS(PHI))/2/PI
T2=ATAN(SQRT((1-E)/(1+E))*TAN(PHI/2))/PI
T=T1+T2
ELSE
T=0.5
ENDIF
DRDT=
& 2*PI*AXM/TSEC/SQRT(1-E**2)*SQRT(SIN(PHI)**2+(E+COS(PHI))**2)
DR1DT=DRDT*M2/(M1+M2)
DR2DT=DRDT*M1/(M1+M2)
IF (ANG.GT.180.0) T=1-T
WRITE(*,'(F14.10,F17.8,E17.8,E17.8,E17.8)')
& T,ANG,DR1DT,DR2DT,DRDT
WRITE(12,'(F14.10,F17.8,E17.8,E17.8,E17.8)')
& T,ANG,DR1DT,DR2DT,DRDT
ENDDO
CLOSE(12)
STOP
END