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 solution (listed below in detail) is analytically !exact, but does not include Lorentz/Einstein "relativistic" !effects. All derived distances are measured relative to the !center of mass for this two-body program, and all units are !relative to solar numbers ("mass" is in solar units where one !solar mass equals 1.991E30 kg, and "distance" is in units of one !"Astronomical Unit" (AU) which is 1.495E11 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 !times, angles, and speeds (in units meters per second). ! !There are many references describing the classic two-body !gravitational problem and the analytical solutions between !times, angles, distances, and velocities. Most text books are !very thorough and lengthy, describing energetics, etc., and with !many re-defined variables and parameters. These are steps which !may be overtly confusing and sidetrack the basic issue of just !wanting to obtain the total four-parameter solution to this !problem. In the following I've written a description of this !solution in four key steps using only standard parameters from !physics and analytical geometry to keep things as simple as !possible. ! !The listing of speeds in this program enables an evaluation of !relativistic effects. The speed of light in vacuum ("c") is !close to 2.99792458E8 m/s. If orbital speed of either of the !two masses is comparable to "c", then the problem should be !re-addressed using General Relativity of Einstein. However, !even for extreme situations this program is still adequate for !fundamental analysis. Given a massive binary star system (try !this in this program) where both masses are 1000 solar masses, !a semi-major axis separation of only 0.01 AU, and eccentricity !of 0.5, the maximum speed for both masses is about 1.15E7 m/s !which is around 0.0384 c. It follows that relativistic mass and !time dilation effects then go as 1/sqrt(1 - 0.384^2) which is !about 1.000737 (i.e., close to 1, and may be characterized as !"non-relativistic"). ! !The full solution to the generalized two-body non-relativistic !problem is just fundamental physics regarding center of mass !along with multiple applications of the conservation of total !angular momentum of the two masses about their mutual center of !mass. The clue to analytically fully solving the generalized !two-body problem (i.e., times, angles, speeds, distances) about !center of mass is to solve for time as a function of angle. It !turns out that "distances" are not required to do this because !of central gravitational force between the two masses and !resulting conservation of total angular momentum. ! !Four steps are listed below to solve this classic problem. The !first (STEP 1) is to determine total angular momentum about the !center of mass and the central force equation for the two !masses, the second (STEP 2) is to derive the time-independent !equation of elliptical motion for the distance between the two !masses. The third (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 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 problem. Note that time is set to zero at' WRITE(*,*)'Perihelion (closest distance between the two masses).' WRITE(*,*)' ' WRITE(*,*)'Geometry of the two-body gravitational system:' WRITE(*,*)' ' WRITE(*,*)' (AT TIME=ZERO)' WRITE(*,*)' ' WRITE(*,*)' M2 cm M1' WRITE(*,*)' O-----+----------------------o' WRITE(*,*)' R2 R1' 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 between' WRITE(*,*)'the two masses.' WRITE(*,*)'"Time" in output printed table is in units of one' WRITE(*,*)'revolution period.' WRITE(*,*)'"Angle" in output tables (in degrees) is position' WRITE(*,*)'angle measured counterclockwise about cm of either M1' WRITE(*,*)'or M2 (relative to horizontal line shown above and is' WRITE(*,*)'set to zero at time = zero and corresponds to' WRITE(*,*)'"Perihelion" which is position of closet separation' WRITE(*,*)'between the masses M1 and M2).' WRITE(*,*)'V1 and V2 in output table are instantaneous speeds of' WRITE(*,*)'masses M1 and M2 relative to cm, respectively.' WRITE(*,*)' ' WRITE(*,*)'Input the following data in solar units:' WRITE(*,*)' ' WRITE(*,*)'Enter first mass (e.g., if the Sun, enter 1):' READ(*,*) M2 WRITE(*,*)'Enter second mass (e.g., if Earth, enter 3.0e-6):' READ(*,*) M1 WRITE(*,*)'Enter semi-major axis (i.e., time-averaged mean' WRITE(*,*)'distance between the two masses in AU - for the' WRITE(*,*)'Earth-Sun system enter 1):' READ(*,*) AX WRITE(*,*)'Enter eccentricity of orbit (e.g., for Earth-Sun' WRITE(*,*)'enter 0.0167):' READ(*,*) E WRITE(*,*)'Enter the increment for the angle sampling in the' WRITE(*,*)'output table in degrees (20 is a good number):' READ(*,*) INC N=IFIX(360.0/INC) C=PI/180.0 TSEC=SQRT(AX*AX*AX/(M1+M2))*3.155815E7 !<== Rev. Period (seconds) TYR=TSEC/3.155815E7 !<== 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.495E11 m, one solar mass = 1.991E30 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 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).' WRITE(*,*)' ' WRITE(*,*) & 'Time(period) Angle(deg) R1(AU) R2(AU) R(AU)' WRITE(*,*) & '--------------------------------------------------------' WRITE(12,*)' ' WRITE(12,*)'Fundamental parameters for your two-body system' WRITE(12,*)'(1 AU = 1.495E11 m, one solar mass = 1.991E30 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 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).' 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(*,*)' ' 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,*)' ' WRITE(12,*) & 'Time(period) Angle(deg) V1(m/s) V2(m/s) V(m/s)' WRITE(12,*) & '--------------------------------------------------------' AXM=AX*1.495E11 !<== 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