Using Julia to calc distance to moon with UNCERTAINTIES

StevenSiew · August 25, 2019, 12:33pm

I am an old fart now. But once upon a time I was a high school kid.

And as a school kid I dream about measuring the distance to the moon. I fantasize about given a science project by my science teacher to measure the distance to the moon with a huge science budget of ONE DOLLAR.

Well, guess what it can be done. Julia is FREE.

The basic idea is to form a triangle between
A) Melbourne Australia
B) Johannesburg South Africa
C) Moon

Figure out the Angle between two Vectors using the DotProduct

Then use the “Law of Sin” to calc “distance from Johannesburg to the moon” given

Distance from Melbourne to Johannesburg
Angle_Melb_to_Moon_to_Johan
Angle_Johan_to_Melb_to_Moon

Here is the Source Code for the Julia progam

# Find the distance to the moon

using Measurements

# Radius of the Earth in metres
EarthRadius = 6371008.8 ± 1.0

# Spherical Coordinate Systems
#
#    SphericalGeoCoor  {Radius,Long, Lat}
#    Radius is EarthRadius
#    Long is  -180 to 180
#    Lat is -90 to 90

#    SphericalMathCoor {Radius,Long,CoLat}
#    Radius is EarthRadius
#    Long is -180 to 180
#    CoLat is 0 to 180  where 0 is North Pole

# Horizontal Coodinate system
#    HorizontalCoor   {Azimuth,Elevation}
#    Azimuth is 0 to 360  with 0 is North , 90 is East , 180 is South , 270 is West
#    Elevation is 0 to 90  with 0 being the horizon and 90 being Zenith (overhead)

# functions
GeoToMath(geo) = [geo[1], geo[2], 90.0 ± 0.0 - geo[3]]

function MathToCart(math)
  return [
    math[1] * cos(math[2] * π / 180.0) * sin(math[3] * π / 180.0),
    math[1] * sin(math[2] * π / 180.0) * sin(math[3] * π / 180.0),
    math[1] * cos(math[3] * π / 180.0)
  ]
end

CartDistance(a,b) = sqrt(  (a[1] - b[1])^2.0 + (a[2] - b[2])^2.0 + (a[3] - b[3])^2.0  )

CartLength(a) = sqrt(  a[1]^2.0 + a[2]^2.0 + a[3]^2.0  )

function CreateVectorFromMathCoorAndHorCoor(math, hor)
  radius = math[1]
  long = math[2]
  colat = math[3]
  azimuth = hor[1]
  elevation = hor[2]
  vertical = sin(elevation * π / 180.0)
  horizontal = cos(elevation * π / 180.0)
  "vertical is positive Z direction"
  horizontaleasting = horizontal * sin(azimuth * π / 180.0)
  horizontalnorthing = -1.0 * horizontal * cos(azimuth * π / 180.0)
  x = horizontalnorthing
  y = horizontaleasting
  z = vertical
  newx = x * cos(colat * π / 180.0) * cos(long * π / 180.0) -
         y * 1.0 * sin(long * π / 180.0) +
         z * sin(colat * π / 180.0) * cos(long * π / 180.0)
  newy = x * cos(colat * π / 180.0) * sin(long * π / 180.0) +
         y * 1.0 * cos(long * π / 180.0) +
         z * sin(colat * π / 180.0) * sin(long * π / 180.0)
  newz = -1.0 * x * sin(colat * π / 180.0) * 1.0 +
         z * cos(colat * π / 180.0) * 1.0
  return [newx, newy, newz]
end

NormaliseVector(v) = [ v[1], v[2], v[3] ] / sqrt( v[1]^2.0 + v[2]^2.0 + v[3]^2.0 )

DotProduct(a,b) = a[1]*b[1] + a[2]*b[2] + a[3]*b[3]

AngleBetweenTwoVectorInDegree(a,b) = acosd( DotProduct(NormaliseVector(a), NormaliseVector(b)) )

# Main program

# The coordinate for Melbourne Australia is
# Longtitude 144 degrees 57 minutes and 47 seconds East
# Latitude 37 degrees 48 minutes and 49 seconds South
MelbGeoCoor = [
  EarthRadius,
  (144.0 + 57.0 / 60.0 + 47.0 / (60.0^2.0)) ± 0.015,
  (-1.0 * (37.0 + 48.0 / 60.0 + 49.0 / (60.0^2.0))) ± 0.015
]

MelbMathCoor = GeoToMath(MelbGeoCoor)

MelbCartCoor = MathToCart(MelbMathCoor)

# The coordinate for Johannesburg South Africa is
# Longtitude 28 degrees 2 minutes and 44 seconds East
# Latitude 26 degrees 12 minutes and 16 seconds South
JohannesburgGeoCoor = [
  EarthRadius,
  (28.0 + 2.0 / 60.0 + 44.0 / (60.0^2.0)) ± 0.015,
  (-1.0 * (26.0 + 12.0 / 60.0 + 16.0 / (60.0^2.0))) ± 0.015
]

JohannesburgMathCoor = GeoToMath(JohannesburgGeoCoor)

JohannesburgCartCoor = MathToCart(JohannesburgMathCoor)

DistanceFromMelbToJohannesburg = CartDistance(MelbCartCoor, JohannesburgCartCoor)

# Melb to Johannesburg vector
VectorMelbToJohan = JohannesburgCartCoor - MelbCartCoor

CartLengthMelbToJohan = CartLength(VectorMelbToJohan)

# Johannesburg to Melb vector
VectorJohanToMelb = MelbCartCoor - JohannesburgCartCoor

CartLengthJohanToMelb = CartLength(VectorJohanToMelb)

# Melb to Moon normalized vector
# From Melb the moon is at 270 deg 37 min 05 sec azimuth and elevation of 26 deg 47 min 37 sec
HorizontalCoor = [ 270.6181±0.1, 26.7936±0.1 ]

VectorMelbToMoon = CreateVectorFromMathCoorAndHorCoor(MelbMathCoor, HorizontalCoor)

CartLengthMelbToMoon = CartLength(VectorMelbToMoon)

# Johan to Moon normalized vector
# From Johan the moon is at 88 deg 03 min 02 sec azimuth and elevation of 42 deg 02 min 04 sec
HorizontalCoor = [ 88.0506±0.1, 42.0344±0.1 ]

VectorJohanToMoon = CreateVectorFromMathCoorAndHorCoor(JohannesburgMathCoor, HorizontalCoor)

CartLengthJohanToMoon = CartLength(VectorJohanToMoon)

# Calc Angle for MoonToMelbToJohan

AngleMoonToMelbToJohanInDegree = AngleBetweenTwoVectorInDegree(VectorMelbToMoon, VectorMelbToJohan)

# Calc Angle for MoonToJohanToMelb

AngleMoonToJohanToMelbInDegree = AngleBetweenTwoVectorInDegree(VectorJohanToMoon, VectorJohanToMelb)

# Calc Angle for MelbToMoonToJohan

AngleMelbToMoonToJohanInDegree = 180.0 - AngleMoonToMelbToJohanInDegree - AngleMoonToJohanToMelbInDegree

# Finally calculate the distance

DistJohanToMoon = sin(AngleMoonToMelbToJohanInDegree  * π / 180.0) * CartLength(VectorMelbToJohan) / sin(AngleMelbToMoonToJohanInDegree * π / 180.0)

println("Distance from Johannesburg to moon is ",DistJohanToMoon," metres")
println("which is between ",Measurements.value(DistJohanToMoon)-Measurements.uncertainty(DistJohanToMoon)," and ",Measurements.value(DistJohanToMoon)+Measurements.uncertainty(DistJohanToMoon)," metres")

# At perigee — its closest approach — the moon comes as close as 363,104 kilometers.
# At apogee — the farthest away it gets — the moon is 405,696 km from Earth.
# On average, the distance from Earth to the moon is about 384,400 km.

with output

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 1.2.0 (2019-08-20)
 _/ |\__ _|_|_|\__ _|  |  Official https://julialang.org/ release
|__/                   |

[ Info: Precompiling Measurements [eff96d63-e80a-5855-80a2-b1b0885c5ab7]
Distance from Johannesburg to moon is 3.55e8 ± 3.3e7 metres
which is between 3.216740488208968e8 and 3.881716048152297e8 metres

longemen3000 · August 25, 2019, 9:02pm

Nice!

giordano · September 6, 2019, 10:41pm

One maybe little known feature of Measurements.jl is that it allows to have an idea of which measurement is contributing most to the total error. The analysis of the error budget can then tell you which measurement should be improved in order to reduce the error. This is described in the manual at Examples · Measurements.

In this case we have

julia> Measurements.uncertainty_components(DistJohanToMoon)
Dict{Tuple{Float64,Float64,UInt64},Float64} with 9 entries:
  (42.0344, 0.1, 0x000000000000000a)    => 2.07562e7
  (6.37101e6, 1.0, 0x0000000000000002)  => 55.709
  (270.618, 0.1, 0x0000000000000007)    => 9.17039e6
  (88.0506, 0.1, 0x0000000000000009)    => 9.13019e6
  (28.0456, 0.015, 0x0000000000000005)  => 3.47003e6
  (144.963, 0.015, 0x0000000000000003)  => 3.47003e6
  (-37.8136, 0.015, 0x0000000000000004) => 686652.0
  (-26.2044, 0.015, 0x0000000000000006) => 1.15949e6
  (26.7936, 0.1, 0x0000000000000008)    => 2.19387e7

The keys of this dictionary are the independent quantities on which DistJohanToMoon depends, the values of the dictionary indicate the corresponding contribution to the total uncertainty, which is the Euclidean norm of the values of the dictionary. Thus, in this case the elevations of Melbourne to the Moon (26.7936±0.1) and Johannesburg to the Moon (42.0344±0.1) give the largest contributions to the final error.