From 686c0f208d26553e73823aef7d44552e92deaef8 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Tue, 3 Aug 2021 14:13:23 -0700 Subject: [PATCH] feat: add harmonic resolution calculator (#44) increase version number --- doc/source/index.rst | 1 + .../user_guide/calc_harmonic_resolution.rst | 33 +++++++ scripts/calc_harmonic_resolution.py | 92 +++++++++++++++++++ version.txt | 2 +- 4 files changed, 127 insertions(+), 1 deletion(-) create mode 100644 doc/source/user_guide/calc_harmonic_resolution.rst create mode 100755 scripts/calc_harmonic_resolution.py diff --git a/doc/source/index.rst b/doc/source/index.rst index 284a3a20..a26de08f 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -99,6 +99,7 @@ missions user_guide/aod1b_geocenter.md user_guide/aod1b_oblateness.md user_guide/calc_mascon.rst + user_guide/calc_harmonic_resolution.rst user_guide/calc_sensitivity_kernel.rst user_guide/combine_harmonics.rst user_guide/convert_harmonics.rst diff --git a/doc/source/user_guide/calc_harmonic_resolution.rst b/doc/source/user_guide/calc_harmonic_resolution.rst new file mode 100644 index 00000000..57188717 --- /dev/null +++ b/doc/source/user_guide/calc_harmonic_resolution.rst @@ -0,0 +1,33 @@ +=========================== +calc_harmonic_resolution.py +=========================== + +- Calculates the spatial resolution that can be resolved by the spherical harmonics of a certain degree [Barthelmes2013]_ [HofmannWellenhof2006]_ +- Default method uses the smallest half-wavelength that can be resolved +- Secondary method calculates the smallest possible bump that can be resolved + +Calling Sequence +################ + +.. code-block:: bash + + python calc_harmonic_resolution --lmax 60 --cap + +`Source code`__ + +.. __: https://github.com/tsutterley/read-GRACE-harmonics/blob/main/scripts/calc_harmonic_resolution.py + +Command Line Options +#################### + +- ``--help``: list the command line options +- ``-l X``, ``--lmax X``: maximum spherical harmonic degree +- ``-R X``, ``--radius X``: Average radius of the Earth (km) +- ``-C``, ``--cap``: Calculate the smallest possible bump that can be resolved + +References +########## + +.. [Barthelmes2013] F. Barthelmes, "Definition of Functionals of the Geopotential and Their Calculation from Spherical Harmonic Models", GeoForschungsZentrum Scientific Technical Report, STR09/02, (2013). `doi: 10.2312/GFZ.b103-0902-26 `_ + +.. [HofmannWellenhof2006] B. Hofmann-Wellenhof and H. Moritz, *Physical Geodesy*, 2nd Edition, 403 pp., (2006). `doi: 10.1007/978-3-211-33545-1 `_ diff --git a/scripts/calc_harmonic_resolution.py b/scripts/calc_harmonic_resolution.py new file mode 100755 index 00000000..d8a72efb --- /dev/null +++ b/scripts/calc_harmonic_resolution.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python +u""" +calc_harmonic_resolution.py +Written by Tyler Sutterley (09/2020) + +Calculates the spatial resolution that can be resolved + by the spherical harmonics of a certain degree +Default method uses the smallest half-wavelength that can be resolved + (is equal to approximately 20000/lmax km) +Secondary method calculates the smallest possible bump that can be resolved + by dividing the area of a sphere by (lmax+1)^2 + +CALLING SEQUENCE: + python calc_harmonic_resolution.py --lmax 60 --cap + +COMMAND LINE OPTIONS: + -l X, --lmax X: maximum degree of spherical harmonics + -R X, --radius X: average radius of the Earth in kilometers + -C, --cap: calculate the smallest possible bump that can be resolved + +PYTHON DEPENDENCIES: + numpy: Scientific Computing Tools For Python (https://numpy.org) + +REFERENCES: + Hofmann-Wellenhof and Moritz, "Physical Geodesy" (2005) + http://www.springerlink.com/content/978-3-211-33544-4 + Barthelmes, "Definition of Functionals of the Geopotential and Their + Calculation from Spherical Harmonic Models", STR09/02 (2009) + http://icgem.gfz-potsdam.de/str-0902-revised.pdf + +UPDATE HISTORY: + Updated 09/2020: using argparse to set parameters + Updated 10/2019: changing Y/N flags to True/False + Updated 02/2014: minor update to if statement + Updated 08/2013: changed SPH_CAP option to (Y/N) + Written 01/2013 +""" +import sys +import argparse +import numpy as np + +#-- PURPOSE: Calculates minimum spatial resolution that can be resolved +#-- from spherical harmonics of a maximum degree +def calc_harmonic_resolution(LMAX, RADIUS=6371.0008, SPH_CAP=False): + """ + Calculates minimum spatial resolution that can be resolved from + spherical harmonics of a maximum degree + + Arguments + --------- + LMAX: maximum spherical harmonic degree + + Keyword arguments + ----------------- + RADIUS: average radius of the Earth in kilometers + SPH_CAP: calculate the smallest possible bump that can be resolved + """ + if SPH_CAP: + # Smallest diameter of a spherical cap that can be resolved by the + # harmonics. Size of the smallest bump, half-wavelength, which can + # be produced by the clm/slm + psi_min = 4.0*RADIUS*np.arcsin(1.0/(LMAX+1.0)) + else: + # Shortest half-wavelength that can be resolved by the clm/slm + # This estimation is based on the number of possible zeros along + # the equator + psi_min = np.pi*RADIUS/LMAX + return psi_min + +# Main program that calls calc_harmonic_resolution() +def main(): + # Read the system arguments listed after the program + parser = argparse.ArgumentParser() + parser.add_argument('--lmax','-l', metavar='LMAX', + type=int, nargs='+', + help='maximum degree of spherical harmonics') + parser.add_argument('--radius','-R', + type=float, default=6371.0008, + help='Average radius of the Earth in kilometers') + parser.add_argument('--cap','-C', + default=False, action='store_true', + help='Calculate smallest possible bump that can be resolved') + args,_ = parser.parse_known_args() + # for each entered spherical harmonic degree + for LMAX in args.lmax: + psi_min = calc_harmonic_resolution(LMAX, + RADIUS=args.radius, SPH_CAP=args.cap) + print('{0:5d}: {1:0.4f} km'.format(LMAX,psi_min)) + +# run main program +if __name__ == '__main__': + main() diff --git a/version.txt b/version.txt index 50610e00..67539e47 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -1.0.1.23 \ No newline at end of file +1.0.2.0 \ No newline at end of file