
%-------------------------------------------------------------------------
% ReadTest_GravityGrids_from_NewtonianIntegration
%-------------------------------------------------------------------------
% Christian Hirt, TU Munich, & Michael Kuhn, Curtin University
%
% This archive contains 24 gravity grids computed from various topographic
% mass models of the Moon. The gravity grids were derived through numerical
% integration (Newtonian integration) in the spatial domain. The grids
% contain gravity disturbances (first radial derivatives of the topographic
% potential) implied by the topographic mass model. The gravity values
% are sometimes denoted "gravity-from-topography" or topographic gravity.
%
% The grid resolution is 5 arc-min; the grid points are equally spaced 
% in terms of planetocentric latitudes and longitudes (4320 x 2160 nodes).
% The unit of the gravity values is mGal.
% The exact grid layout is defined below in Sect. A: of this Matlab script.

% The gravity values are stored in double precision ieee-big endian binary 
% format; each grid file is of 72.9 MB size.

% The gravity values differ by a) four different evaluation heights (e.g., 
% at the surface of the topography or at the Brillioun-sphere), and b) six 
% spectral resolutions of the input topography model. See Sect. B of this
% Matlab script. 
%
% Zipped files: Please note that the gravity grids of the same
% spectral resolution of the input topography have been compressed into
% zip files. As such, each of these zip-files contains four different
% gravity grids that differ in the evaluation heights (see Sect. B).

% This archive is to document the results published in the publication
%
% Hirt, C. and M. Kuhn (2017), Convergence and divergence in 
% spherical harmonic series of the gravitational field generated by 
% high-resolution planetary topography – a case study for the Moon. 
% Journal of Geophysical Research – Planets, 122, doi:10.1002/2017JE005298.
%
% For full details on the computations (e.g., exact definition of 
% computation heights, constants and procedures used) the reader is 
% referred to the  above reference and references cited therein.
% ------------------------------------------------------------------------
% Christian Hirt, 02 August 2017
% ------------------------------------------------------------------------



clear all
close all

%-------------------------------------------------------------------------
% A: Grid information (resolution and grid latitude/longitude base vectors)
%-------------------------------------------------------------------------
res = 5; %in arcmin
res = res/60; r2 = res/2;
lam_vec = -180+r2:res:180-r2;
phi_vec =  -90+r2:res:90-r2;
n = length(lam_vec); m = length(phi_vec);

%-------------------------------------------------------------------------
% B: Gravity grids from numerical (Newtonian) integration
%-------------------------------------------------------------------------
% A total of 24  5arc-min resolution grids of topographic gravity grids
% are available, that differ by 
% a) the evaluation height (at the topography Hsur, at the positive 
% topography Hpos, and at heights of 11 km (H11km) and 23 km (H23km)
% b) the resolution of the input topography model that generates the 
% gravity field; resolution levels are degrees 90,180,360,720,1080, 2160

% Construct filename -----------------------------------------------------
fnA = 'Vz_MoonTop'; fnB = '_180W180E90N90S_30s05m.bin'
fnResolution = {'N90' 'N180' 'N360' 'N720' 'N1080' 'N2160'}
fnHeight     = {'Hsur' 'Hpos' 'H11km' 'H23km'}
% Use indices to specify the desired resolution level and evaluation
% height, as per the example below.
fileName = strcat(fnA,fnHeight{1},'MERG_',fnResolution{3},fnB);

% Read the gravity grid --------------------------------------------------
fid = fopen(fileName,'r');
[data, cnt] = fread(fid, inf,'double','ieee-be');
fclose(fid);
data = reshape(data, m, n);

% Display the gravity grid -----------------------------------------------
figure(1)
imagesc(lam_vec,phi_vec,data)
set(gca,'XDir','normal','YDir','normal');grid on
colorbar; xlabel('Longitude [deg]'); ylabel('Latitude [deg]')
title(['Gravity grid from file' strrep(fileName,'_','-') ' in mGal'] )
infostring = 'For full details see Hirt & Kuhn (2017), JGR Planets 122, doi:10.1002/2017JE005298.'
text(-170,-80,infostring);
print(strcat('Gravity_grid_',fileName,'.png'),'-r300','-dpng')


