@MarkWieczorek I have a very peculiar motivation. I need to properly calculate hyperspherical harmonics to high values of k. Hyperspherical Harmonics have three quantum numbers k,l,m. In principle this would be simple since the hyperspherical harmonics is just a gebenbauer polynomial evaluated at cos(ksi) times the Associate Legendre Polynomial evaluates at cos(theta) times cos(m*phi). The issue, I think, is the normalization factors. They involve large number factorials. I used mpmath. I think I did it rigth...yet... the results tells me that I did it wrong.
So, I suspect there is a lot of room for mistakes and I am sure I am making them. So, in that vein, I would welcome if shtools were to implement Hyperspherical Harmonics. I need the normalization factors properly done to be able to calculate spectra. Considering how much difficult I am having to get this done properly, I wouldn't trust myself to add my work to shtools.
mark-wieczorek I probably won't have time to do this myself (I've never had to use them before...), but perhaps someone else would be interested in working on this.
I am using pyshtools to find the Earth's moon gravity vector to maintain an orbiter vehicle. I have a python script that uses MakeGravGridDH:
degree_order_max = 660 # TEC
degree_order = 330 # TEC
rad, theta, phi, total, pot = gravmag.MakeGravGridDH(
clm, gm, r0, lmax=degree_order_max, a=constant.a_orbit_moon.value,
#f=constant.f_moon.value,
lmax_calc=degree_order,
omega=constant.omega_moon.value, normal_gravity=0, extend=1) My orbiter has a semi major axis around the moon at 88000 km. Can I simply use 88000 km in place of r0? I want to get the magnitude of the gravity vector for each lat/long as the vehicle orbits the moon. I am using GRAIL660 data set.
@timcurry1 Hi there. You can do this, but there is one thing to keep in mind: The routine MakeGravGrid was written to compute the gravity on the surface of a planet in a coordinate system rotating with the planet. By specifying omega, the routine will include the centrifugal potential. If you are a spacecraft that is not attached to the planet, then you would want to set omega to zero. This way, the computed gravity will only be from the mass within the Moon, and won't have any rotational components.
If you are really doing orbit computations, then you will probably want to add a gravitational contribution arising from Earth. Also, as orbits are usually eccentric, I'm not sure how you would use this routine in that situation. If the altitude of the spacecraft is really ~80000 km above the surface, I also really doubt that you need to compute the gravity up to degree 600 (maybe ~60 is good enough?), but this is something that you will need to test.
Lastly, you should check out the python classes SHGravCoeffs and SHGravGrid. This might simplify things for you, and will allow you to easily plot maps.
And to answer your question, you should set
And to answer your question, you should set
a to the semimajor axis of the spacecraft, and r0 should be the reference radius of the spherical harmonic coefficients.
Thanks, Mark. I have been comparing our outputs over the last few weeks. My simulation's gravity magnitude is off by the 5th decimal places from your SHTools at a semi major axis of 3999.9 km. Our values are about 1.0054132 ft/sec^2. I expect that my sim has a different frame. You mentioned that MakeGravGrid uses a fixed frame rotating with the moon. Is there a way to get the Cartesian coordinates and the moon's gravity using SHTools? I am not understanding the flatness attribute while using gravmag.MakeGravGridDH. My thoughts are that my frame does not match yours because of the ecliptic plane, elliptic plane, inertial frame or fixed frame. I am not modeling the earth's influence on the moon's gravity.
I also read in your documentation about errors and anomalies. Can I toggle errors while using the gravmag.MakeGravGridDH?
Thanks for your suggestions.
mark-wieczorek If the result is off in the 5th decimal place, this is possibly a single vs double precision problem. All computations in SHtools are in double.
mark-wieczorek As for your units of ft/sec^2: What's up with that? 😲
mark-wieczorek And for the reference frame: shtools uses the reference frame of the gravity model. As long as you don't include the normal gravity of centrifugal potential, MakeGravGrid gives you the gravitational potential from the mass of the Moon in the coordinate system used for the gravity model.
mark-wieczorek * And for the reference frame: shtools uses the reference frame of the gravity model. As long as you don't include the normal gravity or the centrifugal potential, MakeGravGrid gives you the gravitational potential from the mass of the Moon in the coordinate system used for the gravity model.
mark-wieczorek * And for the reference frame: shtools uses the reference frame of the gravity model. As long as you don't include the normal gravity or the centrifugal potential, MakeGravGrid gives you the gravitational acceleration from the mass of the Moon in the coordinate system used for the gravity model.
Mark: Our Sim also uses double floating points. I agree with your disappointment in using feet. Our newer sims use metric.
I was looking at your ellipsoid model but I realized that it is not affecting gravity because normal gravity is not used. I assume the reference frame is Z-up through the moon's North pole. Is there a way to extract the reference frame's Cartesian coordinates X, Y and Z? I want to understand how the reference frame is transformed back into earth's ECEF and into the Moon's center inertial (or fixed). I currently have a quaternion between the moon's inertial and fixed and I want to validate it, too.
Thanks for helping me with my validation.
mark-wieczorek You will have to get an ephemeris for the moon for what you are asking. The standard one is run by JPL. Also be aware that there are different ways to define lunar coordinates. The gravity models use a principal axis frame, but image data sets use a mean sub earth frame. You should look at some of the papers by konopliv et al to see how they define the coordinate system.
Hi Everyone. I'm working on adding datasets for the gravity field of Earth to pyshtools, but I'm not sure which ones to include. If you use these, could you let me know what you would like to see? I was thinking of about 4 in total, ranging from models based on spacecraft only data, up to something like egm2008. I looked on the ICGEM website, but there are too many models available...
Hello, Mark. I think, the EGM2008 should definitely be included because it is the de facto standard for now and recommended by IERS Conventions 2010. Although, EGM2020 is coming soon. But we need to be careful with the tidal system. Either it must be specified in the documentation, or in the future we need to add a function to switch from one to another. It is not difficult
Probably, another good choice is GOCO06S - it is a modern satellite-only model. The previous one GOCO05s is recommended for experiments for the International Height Reference Frame (IHRF). As well as XGM2016
Someone also suggested GGMplus and WGM2012, but I'm not familiar with either of those...
Also a good choice. All of them are alternatives for the EGM2008. But there is no very big difference between them. In my opinion, as an example, one high-degree model is enough for now. Probably, in the future, we need another repo for working with different data sources, including all ICGEM models
Well, or inside SHTOOLS itself
For Topo, I've added Earth2012, and am working on Earth2014 now [the later doesn't have spherical harmonic coefficients for shape models..., so I need both]
Hi @ioshchepkov, I have a quick question about EGM2008: On the ICGEM web site, there is only a single for that is "tide free". But if you download the original archive from NGA, there are two files for "tide free" and "zero tide". Would you be ok with me only including the single file on ICGEM as a pyshtools dataset, or should I include the two models? ICGEM is probably better for me to use because it takes forever to download anything from the NGA website.
Hi @ioshchepkov. I was just looking at the read_icgem_gfc() function that you wrote for pyshtools, and after looking at the gfc files, it looks like there are a few header items that are read, but that aren't returned. This really isn't important, but there is the model name, the normalization, tide_system, and product type. I am just wondering if there is an easy way to return these that won't break your code. For example, I could perhaps add an optional argument like
return_header=True that returns would append a dictionary of the full header to the output.
...and you might want to take a look at the last couple of commits here:
SHTOOLS/SHTOOLS#231
I added support for reading zip, gzip, and urls in the
SHTOOLS/SHTOOLS#231
I added support for reading zip, gzip, and urls in the
read_icgem_gfc() function. Let me know if you have any comments.
For info, I just merged a big PR onto the develop branch for "datasets". With a single line of code, you can download a spherical harmonic dataset from the original source, and return it to the user as a
I'm happy to add more datasets if I've forgotten any, so just let me know. As part of this, I've also updated some of the spherical harmonic io routines to deal with gz files, zip files, and urls. I also added two new shio routines for binary data and [degree, order, value] formats.
SHCoeffs class instance. This makes extensive use of pooch (Thanks @leouieda for all the help!), and I have included almost all the research grade topography, gravity, and magnetic field models that you would probably ever want.I'm happy to add more datasets if I've forgotten any, so just let me know. As part of this, I've also updated some of the spherical harmonic io routines to deal with gz files, zip files, and urls. I also added two new shio routines for binary data and [degree, order, value] formats.
And here is the long list of datasets
Mercury
Topography
- GTMES150: Topography constructed using laser altimeter and occultation data.
Gravity
- JGMESS160A : Konolpiv et al. (2020)
- JGMESS160A_ACCEL : Konolpiv et al. (2020)
- JGMESS160A_TOPOSIG : Konolpiv et al. (2020)
- GGMES100V08 : Genova et al. (2019)
Venus
Topography
- VenusTopo719 : Wieczorek (2015)
Gravity
- MGNP180U : Konopliv et al. (1999)
The Moon
Topography
- MoonTopo2600p : Wieczorek (2015)
Gravity
- GRGM900C : Lemoine et al. (2014)
- GRGM1200B : Goossens et al. (2020)
- GRGM1200B_RM1_1E0 : Goossens et al. (2020)
- GL0900D : Konopliv et al. (2014)
- GL1500E : Konopliv et al. (2014)
Magnetic Field
- T2015_449 : Wieczorek (2018), using data from Tsunakawa et al. (2015)
Mars
Topography
- MarsTopo2600 : Wieczorek (2015)
Gravity
- GMM3 : Genova et al. (2016)
- GMM3_RM1_1E0 : Goossens et al. (2017)
- MRO120D : Konopliv et al. (2016)
Magnetic Field
- Langlais2019 : Langlais et al. (2019)
- Morschhauser2014 : Morschhauser et al. (2014)
Vesta
Gravity
- VESTA20H : Konopliv et al. (2014)
Ceres
Gravity
- CERES18D : Konopliv et al. (2018)
Earth
Topography
- Earth2012 : Hirt et al. (2012)
- Earth2014 : Hirt and Rexer (2014)
Gravity
- EGM2008 : Pavlis et al. (2012); degree 2190
- EIGEN_6C4 : Förste et al. (2014); degree 2190
- GGM05C : Ries et al. (2016); degree 360
- GOCO06S : Kvas et al. (2019); degree 300, satellite only
- EIGEN_GRGS_RL04_MEAN_FIELD : Lemoine et al. (2019); degree 300, satellite only
- XGM2019E : Zingerle et al. (2019); degree 2190
Magnetic field
- IGRF_13 : Thébault et al. (2015); degree 13
- SWARM_MLI_2D_0501 : Thébault et al. (2013); degree 133
- NGDC_720_V3 : Maus (2010); degree 740
- WDMAM2_800 : Lesur et al. (2016); degree 800
Hi @MarkWieczorek According to the regulation of the SHTOOLS, the maximum spherical harmonic degree lmax=nlat/2-1, where nlat is the number of points in latitude. However, as far as I know, the maximum spherical harmonic degree of gravity feld lmax should be equal to nlat. When I do gravity inversion in spherical harmonic domain, I first create spherical harmonics of gravity kernel with degrees 180, and my observation data is 180*360. So in the inversion, the maximum spherical harmonic degree is 89 in this test, and degrees larger than 89 are not involved in the inversion. How can I solve this problem?
We will be putting out a huge new release in about 2 weeks time. If you want to test this version, please checkout out this PR: SHTOOLS/SHTOOLS#239
Everything is done, and I am just waiting for GMT and pygmt to fix a couple of bugs on their end, and as soon as that is done, I'll merge this into master.
There are lots of new features, of which the addition of datasets is probably the most significant!