-- -- package Four_Body -- -- The package defines a data structure for the N-Body -- gravitational problem in two dimensions. -- -- N = 4, but can be set to anything. -- -- N-Body equation: -- Acceleration(Body_j) -- = SumOver(Body_k) { Mass(Body_k) * G * Delta_R / NORM(Delta_R)**3 } -- -- where Delta_R = (X(Body_k) - X(Body_j)) is a vector, -- and where G is the gravitaional constant. -- -- Natural units: -- Unit of time is earth_orbital_period : t_0. -- Unit of distance is earth's orbital semi-major axis : a_0. -- Unit of mass is sum of Sun and Earth mass : m_0. -- The 3 quantities are related to each other by Kepler's law: -- -- t_0 = 2 * Pi * Sqrt (a_0**3 / (G * m_0)) -- -- where G is the gravitaional constant. If you divide variable of -- time in the N_Body equation by t_0 to get natural units, the equation -- becomes: -- -- Acceleration(Body_j) = 2*2*Pi*Pi* -- SumOver(Body_k) { Mass(Body_k) * Delta_R / NORM(Delta_R)**3 } -- -- where Delta_R = (X(Body_k) - X(Body_j)). -- -- Mass, time, and distance are in the natural units given above, -- which gets rid of the gravitational constant G. In the Equation -- as implemented below, the masses are multiplied by 2*2*pi*pi -- when they are defined in order to avoid excessive multiplication -- during evaluation of F(t, Y). -- -- Body1 is the larger, mass = 1.0. -- Body2 is the smaller, mass = 0.6. -- Assume that the orbital period of the 2 stars is 80.0 years. -- Say the planet is body 3: one earth distance -- from the larger star, with 1 earth mass. -- From these we get the semimajor axis "a" in -- Kepler's law given above ... then put the three bodies in -- near circular orbits and observe stability, and -- remember to use the reduced-mass formulas to get distance -- from center of mass: -- -- First_Body_Radius r1 = a*m2 / (m1+m2) -- Second_Body_Radius r2 = a*m1 / (m1+m2) -- -- Planet_Orbital_Radius : constant Real := 1.0;-- earth's orbital period -- Planet_Period : constant Real := 1.0; -- -- Planet_Orbital_Radius : constant Real := 2.0;-- Twice earth's orbit -- Planet_Period : constant Real := 2.82842712474619; -- -- Planet_Orbital_Radius : constant Real := 3.0;-- Thrice earth's orbit -- Planet_Period : constant Real := 5.196152422706632; -- -- Planet_Orbital_Radius : constant Real := 4.0;-- 4 times earth's orbit -- Planet_Period : constant Real := 8.0;-- 4**1.5 from Kepler's law -- -- a : constant Real := 21.71534093275925; -- OrbitalPeriod : constant Real := 80.0; -- generic type Real is digits <>; package Four_Body is No_Of_Bodies : constant := 4; subtype Bodies is Integer range 0 .. No_Of_Bodies-1; subtype XYUV_Index is Integer range 0 .. 3; type CartesianCoordinate is array(XYUV_Index) of Real; subtype State_Index is Integer range 0 .. 4*No_Of_Bodies-1; type Dynamical_Variable is array(State_Index) of Real; function F (Time : Real; Y : Dynamical_Variable) return Dynamical_Variable; -- This defines the equation to be integrated as -- dY/dt = F (t, Y). Even if the equation is t or Y -- independent, it must be entered in this form. procedure Update_State (Y : in out Dynamical_Variable; Body_id : in Bodies; X, Z, U, W : in Real); function State_Val (Y : Dynamical_Variable; Body_id : Bodies; XYUV_id : XYUV_Index) return Real; function "*" (Left : Real; Right : Dynamical_Variable) return Dynamical_Variable; function "+" (Left : Dynamical_Variable; Right : Dynamical_Variable) return Dynamical_Variable; function "-" (Left : Dynamical_Variable; Right : Dynamical_Variable) return Dynamical_Variable; function Norm (Y : Dynamical_Variable) return Real; pragma Inline (F, "*", "+", "-", Norm); R_Pi : constant := 3.14159_26535_89793_23846_26433_83279_50288; Pii : constant Real := R_Pi; TwoPii : constant Real := 2.0 * R_Pi; TwoPii2 : constant Real := 4.0 * R_Pi * R_Pi; Mass : constant array(Bodies) of Real := (1.0*TwoPii2, 0.6*TwoPii2, 3.0E-6 * TwoPii2, 1.0E-22 * TwoPii2); -- body1 has sun mass, body2 has 0.6 sun mass, and -- body3 has earth mass, body4 has negligible mass. end Four_Body;