"""
According to this website
https://en.wikipedia.org/wiki/Wikipedia:Reference_desk/Archives/Science/2006_September_27
I only need to use 64 decimal digits of Pi.
However I need an additional 4 decimal digits for guard digits
So I need a total of 68 decimal digits of Pi
I need to use 3.32 binary digits of accuracy for each 1 decimal digit of accuracy
Thus I need 226 binary digits of Pi.
So I need to use BigFloat with 226 binary digits of accuracy
""";
println("Start of program")
using Printf
setprecision(226) # Set BigFloat to 226 binary digits of accuracy
# Check the accurary of pi using mathematica
# according to Mathematica N[Pi, 68]
# the 68 decimal digits of pi is
# 3.1415926535897932384626433832795028841971693993751058209749445923078
pi68 = BigFloat("3.1415926535897932384626433832795028841971693993751058209749445923078")
println("accurary of pi is ",BigFloat(pi) - pi68)
"""
Now consider that earth is a ellipsoid with
A = 6378137.000000 # Semi-major Axis for Earth in metres
B = 6356752.314245 # Semi-minor Axis for Earth in metres
Then the length of equator of earth is
circumference = 2 * pi * r # where r = A
""";
A = BigFloat("6378137.000000") # Semi-major Axis for Earth in metres
B = BigFloat("6356752.314245") # Semi-minor Axis for Earth in metres
r = A
println("typeof(r) is ",typeof(r))
circumference = BigFloat(2) * BigFloat(pi) * r
@printf("equatorial circumference is %76.68f",circumference)
"""
Calculating the polar circumference of an ellipsoid is not EASY
https://www.mathsisfun.com/geometry/ellipse-perimeter.html
A much easier problem is to calculate the quarter circumference
which turns the whole problem into calculating the quarter meridian
https://en.wikipedia.org/wiki/Meridian_arc#Quarter_meridian
The polar circumference is just 4 times the quarter meridian
""";
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _ | |
| | |_| | | | (_| | | Version 1.9.2 (2023-07-05)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
accurary of pi is 3.709206150687421385731735261547639513367564778757791002453039058917581e-68
typeof(r) is BigFloat
equatorial circumference is 40075016.68557848615317681776140035737460936288284427863135134035429705953315
Now the next step with both of these results from the point of view of a physicist (or actually any scientifically thinking person bar pure mathematicians) would be estimation of the deviation of the real world from the mathematical model. E.g. due to tidal motion. How many digits of precision are “physical”?
using GMT
G = earthtide(T="2022-07-07T12:00:00");
C = grd2cpt(G, cmap=:geo);
imshow(G, proj=:guess, cmap=C, shade=true, coast=true, colorbar=true, title="Earth tide at noon 7-July-2022")
using QuadGK
using Printf
println("Start of program")
setprecision(226) # Set BigFloat to 226 binary digits of accuracy
A = BigFloat("6378137.000000")
B = BigFloat("6356752.314245")
e_square = ( (A+B)*(A-B) )/(A^2)
integral, error = quadgk(x -> sqrt(1 - e_square * sin(x)), BigFloat(0), BigFloat(pi)/2,rtol=1e-68)
println("integral is ",integral)
quartermeridian = A * integral
circumference = 4 * quartermeridian
@printf("equatorial circumference is %76.68f",circumference)
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _ | |
| | |_| | | | (_| | | Version 1.9.2 (2023-07-05)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
Start of program
integral is 1.567444724577919670791384090705783384345415625384518766551994250585625
equatorial circumference is 39989508.77314095534120938460056765247071486472257255348855854781378953208270
It’s hard to even define what is meant here. I mean if you rolled a surveyor wheel along the surface of the earth frozen in time it probably couldn’t be more than 3 digits, there are even smallish mountains that would add kilometers of distance just in the up and down undulations not to mention ocean waves.
It’d be an interesting calculation to just add a sinusoidal perturbation of height 3m and period 300m (emulating ocean waves) and calculate the difference via integration.
Considering a right triangle setup with 300m base, height 3m, the hypotenuse is x1.00005 longer than 300m. So overall, 5 significant digits are to be expected. The land portion of the circumference may be more crucial as it can support wilder curves, so we can assume only 4 significant digits are “physical”.
Another physical issue is the length of the measuring rod (like in the famous coastline length paradox). A km long rod can overlook the ocean waves, while a micrometer long rod will measure molecular crevices.
Well, 3 digits is presumably a reasonable lower limit of the precision estimation. I’d also maintain it makes little sense to give the result in micrometer: that could be an upper limit estimation. Thus 3 <= p <= 11 , where p is the number of decimal digits.
Sure, “estimation by feeling” is not a mathematically rigorous thing.
Since the French originally defined the meter as ten-million’th of the distance between equator and pole. My feeling is the only accurate circumference of the Earth is the mathematically exact 40,000km.
Surprisingly, despite the Earth’s unspherical shape, this is a really good value (40,075.017 km i.e. 3 significant digits).
The only greater surprise might be the accurate estimation by Eratosthenes in ancient times.
The nonspherical shape of the Earth was known by that time, and to my knowledge was taken into account. By the way, the history of triangulation is quite fascinating, as it combined the most modern mathematical and technical achievements of the time with quite adventurous expeditions.
I find it quite enjoyable though. For other ways of thinking about it consider an equipotential surface. The gravitational acceleration at the surface of the earth varies already in the 2nd decimal place
Then there is another name for it, which my employer normally uses when selling our services: Expert opinion
As to the Earth’s shape: A friend of mine, who used to work for a navigation software company once told me, on some map the sea level had a step like a meter or two (can’t remember that exactly). No, if you travel by sea, you wouldn’t feel it - it was just different countries using different models of the Earth, so that non-existing step “marked” the sea border.
Many countries (like mine) have the Hydrographic zero as > 2 meters bellow the mean sea level for navigation reasons in harbors (so that even in low tide there is water to navigate). Not all countries follow this so to say convention. But this makes a difference only in shallow waters. On the other hand it is a permanent source of troubles when depths and height data are join.
And if one looks at the Geoid anomaly with respect to the ellipsoid (the geoid is the equipotential corresponding to the mean sea level) we can see that the south of India is > 100 meters close to center of the Earth, whilst other places may be as ~80 further when comparing to the ellipsoid.
using GMT
G = gmtread("@earth_geoid_06m_g");
viz(G, region=:global, proj=:guess, shade=true, coast=true, colorbar=true, title="Geoid (meters)")
It’d be an interesting calculation to just add a sinusoidal perturbation of height 3m and period 300m (emulating ocean waves) and calculate the difference via integration.
#=
(360.0 deg)/(39989508.77314 m) = 9.002361145326282e-6 deg/m
300 m * 9.002361145326282e-6 deg/m = 0.0027007083435978845 deg
convert to rad = 4.713625273186541e-5 rad
=#
integral2, error2 = quadgk(x -> sqrt(1 - e_square * sin(x)) + 1.5 * sin(2*pi*x/4.713625273186541e-5), BigFloat(0), BigFloat(pi)/2,rtol=1e-68)
quartermeridian2 = A * integral2
circumference2 = 4 * quartermeridian2
@printf("equatorial circumference2 is %76.68f",circumference2)
println()
equatorial circumference2 is 39990037.63922746734945825554723989362469830257394170332800238910744266711623
compared to equatorial circumference is 39989508.77314095534120
if I replace a “straight line” with a curvy line, I just need to multiply the length of the straight line with a constant. And use the formula above to find the constant.
I think my suggestion corresponds to making the “radius” undulate. So in polar coordinates, instead of r = constant, you do r = constant * (1+epsilon*sin(k*x)). then derive the arc length from that and integrate that arc length.
if I replace a “straight line” with a curvy line, I just need to multiply the length of the straight line with a constant. And use the formula above to find the constant.
