Frederik's page of Matlab routines is located here. Pages providing more detailed information about these routines are provided below.
This program reads in the Level-2 GRACE geoid products from either the CSR or GFZ data centers, does some processing, and saves them as a plmt matrix in a .mat file. This program is the best avenue if you desire geopotential as a timeseries of spherical harmonics and/or if you want the calibrated errors from the GRACE data centers. They remain potential coefficients, not geoid coefficients. In particular, the coefficients have several common processings applied to them:
Reordered to our prefered lmcosi format, for easy use with our existing codes.
The geopotential has been referenced to the WGS84 ellipsoid. Currently this is only done up to degree l=4.
The C2,0 coefficients of GRACE are more variable due to a large semiannual signal over Antarctica in the C2,0 coefficients of all processing centers. Don Chambers and colleagues at the U. of Texas pinpointed it as a 161 day fluctuation, the GRACE alias for the S2 tide. Because of this known error in the C2,0 coefficient, all C2,0 coefficients have been replaced with more accurate measurements from satellite laser ranging (SLR). The SLR data are available from the GRACE Tellus website: http://grace.jpl.nasa.gov/data/J2/
) The data come from Cheng, M.K., Tapley, B.D., Variations in the Earth's oblateness during the past 28 years, J. Geophys. Res., 109(B09402), 2004.
KERNELC calculates the localization matrix for some domain on the sphere. It creates the matrix from Simons et al. (2006)(eqns. 4.2–4.6) where
and is a real surface spherical harmonic in the form
KERNELC goes about this calculation in several steps, the main ones being:
Setup the domain and decide how to integrate
The function gausslegendrecof determines the absissas (array x) and weights (array w) of the integration.
Calculate the Legendre polynomials (in array Xlm)
Calculate the products of these polynomials for combinations of l and m (array XlmXlmp).
This is part of what gets integrated. As it turns out, there is much redundancy here. Obviously the product lml'm' will be the same as l'm'lm, so the array XlmXlmp will be of size (101,()) where is of length (L+1)*(L+2)/2. This length occurs because the functions are as m>=0. That is, there is redundancy for the negative m values.
101 is the number of integration points from north to south in the domain. Longitudinal integrals across the domain (E-W) are calculated at these points and then the weighted sum finishes the integration.
Assemble the matrix (array Klmlmp)
The longitudinal integrals are calculated by the routines sinsin, coscos, and sincos. These are then multiplied by the polynomial products XlmXlmp and the weights w. The array bigo indexes XlmXlmp to take into account the reuse of the products. Thus the final square size of Klmlmp is .