Classical Free Fall Calculation

This Mathematica Notebook performs the calculation of the time of free fall from the orbit of the Earth to the Sun, using simple Newtonian functions.

Glen Langston,  NRAO GB, 2004 October 1

In[1]:=

Clear[ GM, Φ, r, t, m, v, A]

The Newtonian Acceleration of a particle in free fall is the derivative of the potential

In[2]:=

Φ[r_] := GM /r

During Free Fall, energy is conserved, so falling at rest from A to r yields

In[3]:=

SetAttributes[A, Constant]

In[4]:=

Energy = m v^2/2   - m Φ[r]  - m Φ[A]

Out[4]=

-(GM m)/r + (m v^2)/2 -(GM m)/A

In[5]:=

V = Solve[ Energy, v]

Out[5]=

{{v -2^(1/2) (-GM/A + GM/r)^(1/2)}, {v2^(1/2) (-GM/A + GM/r)^(1/2)}}

In[6]:=

f = 1/(v/.V[[1]])

Out[6]=

-1/(2^(1/2) (-GM/A + GM/r)^(1/2))

Now work with the special case of falling from rest at 1 AU towards the Sun

In[7]:=

A = 1

Out[7]=

1

Calculate the function, T[r], the time it takes a particle to fall from rest at 1 AU, to a ratius r. The Integral function of a freely falling particle is surprisingly complicated:

In[8]:=

T = Integrate[ f, r]

Out[8]=

-((-1 + r) r^(1/2) + (-1 + r)^(1/2) Log[(-1 + r)^(1/2) + r^(1/2)])/(2^(1/2) (GM (-1 + 1/r))^(1/2) r^(1/2))

Now calculate the duration of the fall from the limits at one AU and zero AU.

In[9]:=

AB = Simplify[ Limit[ T, r1]]

Out[9]=

0

In[10]:=

AA = Simplify[ Limit[ T, r0]]

Out[10]=

π/(2 2^(1/2) GM^(1/2))

In[11]:=

N[AA]

Out[11]=

1.11072/GM^(1/2)

Calculation Units: Conversion to Days

The Distance Units for this calculation are AU.  The time units for this calculation are days.  In these units, GM is a small number, indicating that the system is not relativistic.  GC = gravitational constant = 6.67300 × 10^-11 m3 kg^-1 s^-1; Astronomical Unit = 149 598 000 000;  One Day = 86400 s

In[12]:=

GC = 6.673 10^-11 m^3/kg/s^2

Out[12]=

(6.673*10^-11 m^3)/(kg s^2)

In[13]:=

MS = 1.9891  10^30 kg

Out[13]=

1.9891*10^30 kg

In[14]:=

AU = 1.49598 10^11 m

Out[14]=

1.49598*10^11 m

In[15]:=

day = 3600 24s

Out[15]=

86400 s

Calculate the Speed of Light in AU per day (Its a fairly small number)

In[40]:=

clight = (299782458 m/s )    (day / AU)

Out[40]=

173.139

The Distance Units for this calculation are AU.  The time units for this calculation are days.  In these units, GM (gm) is a small number, indicating that our solar system is not general-relativistic.

In[17]:=

gm = GC MS day^2/ AU^3

Out[17]=

0.000295956

Formally, the velocity of the particle goes to infinity as r goes to zero.  To avoid this, the calculation will stop at the Solar Radius (RS).   The Sun can not be treated as a point source inside a Solar Radius anyway.

In[18]:=

RS = 6.95987 10^8 m

Out[18]=

6.95987*10^8 m

In Astronomical Units (AU) the Solar Radius is 0.47 % of an AU.

In[19]:=

rs = RS / AU

Out[19]=

0.00465238

In[20]:=

MaxT = Limit[ AA, GM  -> gm]

General :: spell1 : Possible spelling error: new symbol name \"MaxT\" is similar to existing symbol \"Max\".  More…

Out[20]=

64.5641

Now numerically evaluate the time, T in days, it takes an object to fall from rest from 1 AU to the surface of the Sun .

In[21]:=

MaxT = Abs[Limit[ Limit[T, rrs], GM  -> gm]]

Out[21]=

64.5554

Calculation of the Speed: v[r]

In[22]:=

T

Out[22]=

-((-1 + r) r^(1/2) + (-1 + r)^(1/2) Log[(-1 + r)^(1/2) + r^(1/2)])/(2^(1/2) (GM (-1 + 1/r))^(1/2) r^(1/2))

t[r] is the time to fall a distance r(AU) in units of days

In[23]:=

t = Limit[ T, GM  gm]

Out[23]=

(-(-1 + r) r^(1/2) - (-1 + r)^(1/2) Log[(-1 + r)^(1/2) + r^(1/2)])/((-0.000591913 + 0.000591913/r)^(1/2) r^(1/2))

In[26]:=

Plot[ t, {r, rs, 1}, AxesLabel {"Distance/AU", "T/days"}]

[Graphics:HTMLFiles/index_47.gif]

Out[26]=

⁃Graphics⁃

Calculate the speed of the falling object in units of the speed of light.

In[28]:=

OneOverV = D[ t, r]

Out[28]=

(-((1/(2 (-1 + r)^(1/2)) + 1/(2 r^(1/2))) (-1 + r)^(1/2))/((-1 + r)^(1/2) + r^(1/2)) - (-1 + r ... ) - (-1 + r)^(1/2) Log[(-1 + r)^(1/2) + r^(1/2)])/(2 (-0.000591913 + 0.000591913/r)^(1/2) r^(3/2))

In[29]:=

v = 1/OneOverV

Out[29]=

1/((-((1/(2 (-1 + r)^(1/2)) + 1/(2 r^(1/2))) (-1 + r)^(1/2))/((-1 + r)^(1/2) + r^(1/2)) - (-1  ...  - (-1 + r)^(1/2) Log[(-1 + r)^(1/2) + r^(1/2)])/(2 (-0.000591913 + 0.000591913/r)^(1/2) r^(3/2)))

In[35]:=

Plot[ v, {r, rs, 1}, AxesLabel { "Distance/AU", "Velocity (AU/day)"}]

[Graphics:HTMLFiles/index_54.gif]

Out[35]=

⁃Graphics⁃

The velocity is negative because it is the infall velocity .   The magnitude of the velocity is :

In[41]:=

mag := Abs[v]/clight

In[42]:=

Plot[ mag, {r, rs, 1}, AxesLabel { "Distance/AU", "β / c"}] 

[Graphics:HTMLFiles/index_59.gif]

Out[42]=

⁃Graphics⁃

The maximum velocity is the velocity when reaching the surface of the Sun (in units of c).  

In[46]:=

MaxMag = Abs[Limit[ mag, rrs]]

Out[46]=

0.00205534

Converting  back to km/s

In[48]:=

MaxMag 299782. km/s

Out[48]=

(616.156 km)/s


Created by Mathematica  (October 1, 2004)