// Program to find the time spent falling through a hole in the earth assuming
// a non-uniform earth. The density of the earth is modeled using several
// layers which model closely the known density of the earth.
// This program is written in Frink:
// https://frinklang.org/
//
// Alan Eliasen, eliasen@mindspring.com
use Grid.frink
// The mass of a whole shell with inner radius d0 and outer radius d1
// and inner density r0 and outer density r1
wholeShellMass[d0 is length, d1 is length, r0 is mass_density, r1 is mass_density] :=
{
((d1 - d0) pi (2 d0 d1 (r0 + r1) + d0^2 (3 r0 + r1) + d1^2 (r0 + 3 r1))) / 3
}
// The mass of the partial shell with inner radius d0 and outer radius d1
// and inner density r0 and outer density r1 which is below an object
// at radius d.
partShellMass[d is length, d0 is length, d1 is length, r0 is mass_density, r1 is mass_density] :=
{
(pi (3 d^4 (r0 - r1) + d^3 (-4 d1 r0 + 4 d0 r1) -
d0^3 (-4 d1 r0 + d0 (3 r0 + r1)))) / (3 (d0 - d1))
}
// The mass of a whole or partial shell with inner radius d0 and outer radius
// d1 and inner density r0 and outer density r1 which is below an object at
// radius d.
shellMass[d is length, d0 is length, d1 is length, r0 is mass_density, r1 is mass_density] :=
{
if (d < d0)
0 kg // inside inner surface of shell
else
if d > d1
wholeShellMass[d0, d1, r0, r1] // completely outside shell
else
partShellMass[d, d0, d1, r0, r1] // partway through shell
}
// The different layers of density in the Earth
shells = [[ 0 km, 1221 km, 13.1 g/cm^3, 12.8 g/cm^3],
[1221 km, 3480 km, 12.2 g/cm^3, 9.9 g/cm^3],
[3480 km, 5651 km, 5.6 g/cm^3, 4.4 g/cm^3],
[5651 km, 6341 km, 4.4 g/cm^3, 3.4 g/cm^3],
[6341 km, earthradius, 2.9 g/cm^3, 2.2 g/cm^3]]
// This finds the mass that's still below you at a given distance from the
// earth's center.
mass[d is length, shells] :=
{
var m is mass = 0 kg
for [d0, d1, r0, r1] shells
m = m + shellMass[d, d0, d1, r0, r1]
// Fudge factor to make integral come out
// with the known mass of the earth.
return 0.995775976555 m
}
// Find the acceleration at a given distance from the core.
a[dist is length, shells] := G mass[dist, shells]/dist^2
var v is velocity = 0 m/s // Velocity at end of timestep
var stepsize is time = 1/100 s
var d is length = earthradius
var t is time = 0 s
var a is acceleration = 0 gravity
ga = new graphics // Graphics for distance-vs-acceleration graph
pa = new polyline
while (d > 0 m)
{
t = t + stepsize
a = a[d, shells]
v = v + a stepsize
d = d - v stepsize
pa.addPoint[d/km, -a]
// Print results every second
if (t mod sec == 0 s)
println[(t -> sec) + "\t" + (d->km) + "\t" + (v->mph) + "\t" + (a->m/s^2)]
}
// Print final times
println["Time to core : " + (1. t -> ["min", "sec"]) + "\t, " + (v->"mph")]
println["Time through earth: " + (2. t -> ["min", "sec"])]
ga.add[pa]
grid = new Grid
grid.setUnits[1, -m/s^2]
grid.auto[ga, false, 10, 25]
ga.add[grid.getGrid[]]
ga.show[]