From aa2eb01dbfd776b398ef72ad7e74a61bb9e5c65c Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 19 Jul 2021 22:09:49 -0400 Subject: [PATCH] Python ObserverGravity function. --- demo/python/astronomy.py | 32 ++++++++++ demo/python/gravity.py | 37 +++++++++++ demo/python/test/.gitignore | 1 + demo/python/test/gravity_correct.txt | 91 ++++++++++++++++++++++++++++ demo/python/test/test | 8 ++- generate/template/astronomy.py | 32 ++++++++++ source/python/README.md | 24 ++++++++ source/python/astronomy.py | 32 ++++++++++ 8 files changed, 256 insertions(+), 1 deletion(-) create mode 100755 demo/python/gravity.py create mode 100644 demo/python/test/gravity_correct.txt diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index 428ef6eb..1ae73465 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -4404,6 +4404,38 @@ def VectorObserver(vector, ofdate): ovec = _nutation(ovec, vector.t, _PrecessDir.From2000) return _inverse_terra(ovec, gast) +def ObserverGravity(latitude, height): + """Calculates the gravitational acceleration experienced by an observer on the Earth. + + This function implements the WGS 84 Ellipsoidal Gravity Formula. + The result is a combination of inward gravitational acceleration + with outward centrifugal acceleration, as experienced by an observer + in the Earth's rotating frame of reference. + The resulting value increases toward the Earth's poles and decreases + toward the equator, consistent with changes of the weight measured + by a spring scale of a fixed mass moved to different latitudes and heights + on the Earth. + + Parameters + ---------- + latitude : float + The latitude of the observer in degrees north or south of the equator. + By formula symmetry, positive latitudes give the same answer as negative + latitudes, so the sign does not matter. + height : float + The height above the sea level geoid in meters. + No range checking is done; however, accuracy is only valid in the + range 0 to 100000 meters. + + Returns + ------- + float + The effective gravitational acceleration expressed in meters per second squared [m/s^2]. + """ + s2 = math.sin(math.radians(latitude)) ** 2 + g0 = 9.7803253359 * (1.0 + 0.00193185265241*s2) / math.sqrt(1.0 - 0.00669437999013*s2) + return g0 * (1.0 - (3.15704e-07 - 2.10269e-09*s2)*height + 7.37452e-14*height*height) + @enum.unique class Refraction(enum.Enum): """Selects if/how to correct for atmospheric refraction. diff --git a/demo/python/gravity.py b/demo/python/gravity.py new file mode 100755 index 00000000..90e7a63f --- /dev/null +++ b/demo/python/gravity.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python3 +import sys +from astronomy import ObserverGravity + +UsageText = r''' + + USAGE: + + gravity.py latitude height + + Calculates the gravitational acceleration experienced + by an observer on the surface of the Earth at the specified + latitude (degrees north of the equator) and height + (meters above sea level). + The output is the gravitational acceleration in m/s^2. + +''' + +if __name__ == '__main__': + if len(sys.argv) != 3: + print(UsageText) + sys.exit(1) + + latitude = float(sys.argv[1]) + if latitude < -90.0 or latitude > +90.0: + print("ERROR: Invalid latitude '{}'. Must be a number between -90 and +90.".format(sys.argv[1])) + sys.exit(1) + + height = float(sys.argv[2]) + MAX_HEIGHT_METERS = 100000.0 + if height < 0.0 or height > MAX_HEIGHT_METERS: + print("ERROR: Invalid height '{}'. Must be a number between 0 and {}.".format(sys.argv[1], MAX_HEIGHT_METERS)) + sys.exit(1) + + gravity = ObserverGravity(latitude, height) + print("latitude = {:8.4f}, height = {:6.0f}, gravity = {:8.6f}".format(latitude, height, gravity)) + sys.exit(0) diff --git a/demo/python/test/.gitignore b/demo/python/test/.gitignore index 200dc141..10f65388 100644 --- a/demo/python/test/.gitignore +++ b/demo/python/test/.gitignore @@ -11,3 +11,4 @@ lunar_angles.txt jupiter_moons.txt galactic.txt triangulate.txt +gravity.txt diff --git a/demo/python/test/gravity_correct.txt b/demo/python/test/gravity_correct.txt new file mode 100644 index 00000000..2a65cbff --- /dev/null +++ b/demo/python/test/gravity_correct.txt @@ -0,0 +1,91 @@ +latitude = 0.0000, height = 0, gravity = 9.780325 +latitude = 1.0000, height = 0, gravity = 9.780341 +latitude = 2.0000, height = 0, gravity = 9.780388 +latitude = 3.0000, height = 0, gravity = 9.780467 +latitude = 4.0000, height = 0, gravity = 9.780577 +latitude = 5.0000, height = 0, gravity = 9.780718 +latitude = 6.0000, height = 0, gravity = 9.780889 +latitude = 7.0000, height = 0, gravity = 9.781092 +latitude = 8.0000, height = 0, gravity = 9.781325 +latitude = 9.0000, height = 0, gravity = 9.781589 +latitude = 10.0000, height = 0, gravity = 9.781882 +latitude = 11.0000, height = 0, gravity = 9.782205 +latitude = 12.0000, height = 0, gravity = 9.782558 +latitude = 13.0000, height = 0, gravity = 9.782939 +latitude = 14.0000, height = 0, gravity = 9.783348 +latitude = 15.0000, height = 0, gravity = 9.783785 +latitude = 16.0000, height = 0, gravity = 9.784249 +latitude = 17.0000, height = 0, gravity = 9.784740 +latitude = 18.0000, height = 0, gravity = 9.785258 +latitude = 19.0000, height = 0, gravity = 9.785800 +latitude = 20.0000, height = 0, gravity = 9.786368 +latitude = 21.0000, height = 0, gravity = 9.786960 +latitude = 22.0000, height = 0, gravity = 9.787575 +latitude = 23.0000, height = 0, gravity = 9.788213 +latitude = 24.0000, height = 0, gravity = 9.788873 +latitude = 25.0000, height = 0, gravity = 9.789554 +latitude = 26.0000, height = 0, gravity = 9.790256 +latitude = 27.0000, height = 0, gravity = 9.790976 +latitude = 28.0000, height = 0, gravity = 9.791716 +latitude = 29.0000, height = 0, gravity = 9.792473 +latitude = 30.0000, height = 0, gravity = 9.793247 +latitude = 31.0000, height = 0, gravity = 9.794037 +latitude = 32.0000, height = 0, gravity = 9.794842 +latitude = 33.0000, height = 0, gravity = 9.795661 +latitude = 34.0000, height = 0, gravity = 9.796492 +latitude = 35.0000, height = 0, gravity = 9.797336 +latitude = 36.0000, height = 0, gravity = 9.798191 +latitude = 37.0000, height = 0, gravity = 9.799055 +latitude = 38.0000, height = 0, gravity = 9.799928 +latitude = 39.0000, height = 0, gravity = 9.800809 +latitude = 40.0000, height = 0, gravity = 9.801697 +latitude = 41.0000, height = 0, gravity = 9.802590 +latitude = 42.0000, height = 0, gravity = 9.803488 +latitude = 43.0000, height = 0, gravity = 9.804389 +latitude = 44.0000, height = 0, gravity = 9.805293 +latitude = 45.0000, height = 0, gravity = 9.806198 +latitude = 46.0000, height = 0, gravity = 9.807103 +latitude = 47.0000, height = 0, gravity = 9.808007 +latitude = 48.0000, height = 0, gravity = 9.808909 +latitude = 49.0000, height = 0, gravity = 9.809808 +latitude = 50.0000, height = 0, gravity = 9.810702 +latitude = 51.0000, height = 0, gravity = 9.811591 +latitude = 52.0000, height = 0, gravity = 9.812474 +latitude = 53.0000, height = 0, gravity = 9.813349 +latitude = 54.0000, height = 0, gravity = 9.814216 +latitude = 55.0000, height = 0, gravity = 9.815073 +latitude = 56.0000, height = 0, gravity = 9.815919 +latitude = 57.0000, height = 0, gravity = 9.816754 +latitude = 58.0000, height = 0, gravity = 9.817576 +latitude = 59.0000, height = 0, gravity = 9.818384 +latitude = 60.0000, height = 0, gravity = 9.819177 +latitude = 61.0000, height = 0, gravity = 9.819955 +latitude = 62.0000, height = 0, gravity = 9.820715 +latitude = 63.0000, height = 0, gravity = 9.821459 +latitude = 64.0000, height = 0, gravity = 9.822183 +latitude = 65.0000, height = 0, gravity = 9.822889 +latitude = 66.0000, height = 0, gravity = 9.823574 +latitude = 67.0000, height = 0, gravity = 9.824238 +latitude = 68.0000, height = 0, gravity = 9.824880 +latitude = 69.0000, height = 0, gravity = 9.825499 +latitude = 70.0000, height = 0, gravity = 9.826095 +latitude = 71.0000, height = 0, gravity = 9.826666 +latitude = 72.0000, height = 0, gravity = 9.827213 +latitude = 73.0000, height = 0, gravity = 9.827734 +latitude = 74.0000, height = 0, gravity = 9.828229 +latitude = 75.0000, height = 0, gravity = 9.828697 +latitude = 76.0000, height = 0, gravity = 9.829137 +latitude = 77.0000, height = 0, gravity = 9.829550 +latitude = 78.0000, height = 0, gravity = 9.829934 +latitude = 79.0000, height = 0, gravity = 9.830289 +latitude = 80.0000, height = 0, gravity = 9.830614 +latitude = 81.0000, height = 0, gravity = 9.830910 +latitude = 82.0000, height = 0, gravity = 9.831176 +latitude = 83.0000, height = 0, gravity = 9.831411 +latitude = 84.0000, height = 0, gravity = 9.831616 +latitude = 85.0000, height = 0, gravity = 9.831789 +latitude = 86.0000, height = 0, gravity = 9.831931 +latitude = 87.0000, height = 0, gravity = 9.832042 +latitude = 88.0000, height = 0, gravity = 9.832121 +latitude = 89.0000, height = 0, gravity = 9.832169 +latitude = 90.0000, height = 0, gravity = 9.832185 diff --git a/demo/python/test/test b/demo/python/test/test index 5964ed6c..a824e981 100755 --- a/demo/python/test/test +++ b/demo/python/test/test @@ -5,7 +5,7 @@ Fail() exit 1 } -rm -f test/{jupiter_moons,camera,moonphase,positions,riseset,seasons,culminate,horizon,lunar_eclipse,lunar_angles}.txt +rm -f test/{jupiter_moons,camera,constellation,moonphase,positions,riseset,seasons,culminate,horizon,lunar_eclipse,lunar_angles,galactic,triangulate,gravity}.txt echo "Testing example: jupiter_moons.py" ./jupiter_moons.py 2021-04-16T00:26:18Z > test/jupiter_moons.txt || Fail "Error testing jupiter_moons.py." @@ -59,5 +59,11 @@ echo "Testing example: triangulate.py" ./triangulate.py 48.16042 24.49986 2019 18 7 48.27305 24.36401 662 83 12 > test/triangulate.txt || Fail "Error running triangulate.py." diff test/triangulate.txt test/triangulate_correct.txt || Fail "Error comparing triangulate.py output." +echo "Testing example: gravity.py" +for latitude in {0..90}; do + ./gravity.py ${latitude} 0 >> test/gravity.txt || Fail "Error running gravity.py." +done +diff test/gravity.txt test/gravity_correct.txt || Fail "Error comparing gravity.py output." + echo "PASS: Python examples" exit 0 diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index fbb4a34d..e2956393 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -2374,6 +2374,38 @@ def VectorObserver(vector, ofdate): ovec = _nutation(ovec, vector.t, _PrecessDir.From2000) return _inverse_terra(ovec, gast) +def ObserverGravity(latitude, height): + """Calculates the gravitational acceleration experienced by an observer on the Earth. + + This function implements the WGS 84 Ellipsoidal Gravity Formula. + The result is a combination of inward gravitational acceleration + with outward centrifugal acceleration, as experienced by an observer + in the Earth's rotating frame of reference. + The resulting value increases toward the Earth's poles and decreases + toward the equator, consistent with changes of the weight measured + by a spring scale of a fixed mass moved to different latitudes and heights + on the Earth. + + Parameters + ---------- + latitude : float + The latitude of the observer in degrees north or south of the equator. + By formula symmetry, positive latitudes give the same answer as negative + latitudes, so the sign does not matter. + height : float + The height above the sea level geoid in meters. + No range checking is done; however, accuracy is only valid in the + range 0 to 100000 meters. + + Returns + ------- + float + The effective gravitational acceleration expressed in meters per second squared [m/s^2]. + """ + s2 = math.sin(math.radians(latitude)) ** 2 + g0 = 9.7803253359 * (1.0 + 0.00193185265241*s2) / math.sqrt(1.0 - 0.00669437999013*s2) + return g0 * (1.0 - (3.15704e-07 - 2.10269e-09*s2)*height + 7.37452e-14*height*height) + @enum.unique class Refraction(enum.Enum): """Selects if/how to correct for atmospheric refraction. diff --git a/source/python/README.md b/source/python/README.md index 42e75d27..c93c6292 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -1637,6 +1637,30 @@ Keep calling this function as many times as you want to keep finding more transi --- + +### ObserverGravity(latitude, height) + +**Calculates the gravitational acceleration experienced by an observer on the Earth.** + +This function implements the WGS 84 Ellipsoidal Gravity Formula. +The result is a combination of inward gravitational acceleration +with outward centrifugal acceleration, as experienced by an observer +in the Earth's rotating frame of reference. +The resulting value increases toward the Earth's poles and decreases +toward the equator, consistent with changes of the weight measured +by a spring scale of a fixed mass moved to different latitudes and heights +on the Earth. + +| Type | Parameter | Description | +| --- | --- | --- | +| `float` | `latitude` | The latitude of the observer in degrees north or south of the equator. By formula symmetry, positive latitudes give the same answer as negative latitudes, so the sign does not matter. | +| `float` | `height` | The height above the sea level geoid in meters. No range checking is done; however, accuracy is only valid in the range 0 to 100000 meters. | + +### Returns: `float` +The effective gravitational acceleration expressed in meters per second squared [m/s^2]. + +--- + ### ObserverVector(time, observer, ofdate) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 428ef6eb..1ae73465 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -4404,6 +4404,38 @@ def VectorObserver(vector, ofdate): ovec = _nutation(ovec, vector.t, _PrecessDir.From2000) return _inverse_terra(ovec, gast) +def ObserverGravity(latitude, height): + """Calculates the gravitational acceleration experienced by an observer on the Earth. + + This function implements the WGS 84 Ellipsoidal Gravity Formula. + The result is a combination of inward gravitational acceleration + with outward centrifugal acceleration, as experienced by an observer + in the Earth's rotating frame of reference. + The resulting value increases toward the Earth's poles and decreases + toward the equator, consistent with changes of the weight measured + by a spring scale of a fixed mass moved to different latitudes and heights + on the Earth. + + Parameters + ---------- + latitude : float + The latitude of the observer in degrees north or south of the equator. + By formula symmetry, positive latitudes give the same answer as negative + latitudes, so the sign does not matter. + height : float + The height above the sea level geoid in meters. + No range checking is done; however, accuracy is only valid in the + range 0 to 100000 meters. + + Returns + ------- + float + The effective gravitational acceleration expressed in meters per second squared [m/s^2]. + """ + s2 = math.sin(math.radians(latitude)) ** 2 + g0 = 9.7803253359 * (1.0 + 0.00193185265241*s2) / math.sqrt(1.0 - 0.00669437999013*s2) + return g0 * (1.0 - (3.15704e-07 - 2.10269e-09*s2)*height + 7.37452e-14*height*height) + @enum.unique class Refraction(enum.Enum): """Selects if/how to correct for atmospheric refraction.