%-------------------------------------------------------------------------
% ReadTest_PotentialCoefficients_from_SpectralForwardModelling
%-------------------------------------------------------------------------
% Christian Hirt, TU Munich, & Michael Kuhn, Curtin University
%
% This archive contains 6 spectral models representing the external 
% gravitational potential implied by the Moon's topography. The spectral
% models were derived through spectral-domain gravity forward modelling
% as described in the Hirt and Kuhn (2017) paper. 
% The field coefficients are provided as solid spherical harmonics that
% refer to the reference sphere of radius 1,738,000 m and a GM value of
% 0.4902805823D+13 m^3 s^-2.

% The key difference between the six spectral models is the resolution
% of the topographic mass model used to define the field-generating mass
% distribution. The resolution levels are in terms of harmonic degrees:
% N = 90, 180, 360, 720, 1080 and 2160, with the lowest-resolution model
% (N = 90) reaching a spatial resolution of ~60 km, and the highest
% resolution model reaching a spatial resolution of ~2.5 km.

% The filenames further provide the following information, e.g.,
% dV_MTshape_in90_out900_R1738000_rho2500_pmax10
% |_ topgraphic potential
%    |_ of the Moon Topography represented as shape model
%            |_ input topography resolved to degree 90
%                  |_ implied gravity field modelled to degree 900
%                        |_relative to a sphere with radius 1738000 m
%                                   |_ and a mass density of 2500 kg m-^3
%                                      for the topography

% pmax = 10 denotes the order of the series expansion used to transform
% topography to gravity; here the contributions of the first 10 integer 
% powers of the topographic height function to the potential were
% considered.
%
% Each file contains (2pi) fully-normalized potential spherical harmonic 
% coefficients (SHCs) in our bshc format, a binary format to store the 
% SHCs. The format is explained in the readme-file to this archive and
% implemented in the read_SHCs_binary_bshc_format2014 function.
%
%
% Use of the SHCs in gravity syntheses:

% The six sets of SHCs were generated to study the effect of
% divergence and convergence in series expansions of gravity functionals
% computed from the SHCs. 
%
% We EXPLICITLY advise that invalid values for gravity functionals may be
% obtained in some cases when the SHCs are used to compute gravity values 
% inside the Brillouin-sphere (radius < 1749000 m). This is studied in 
% detail in the Hirt and Kuhn (2017) paper. Outside the Brillouin-sphere, 
% the series converge accurately.

% Full details are found 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.
% ------------------------------------------------------------------------
% Christian Hirt, 02 August 2017
% ------------------------------------------------------------------------

clear all
close all

% six spectral gravity forward models ------------------------------------
models = [ {'dV_MTshape_in90_out900_R1738000_rho2500_pmax10'}
           {'dV_MTshape_in180_out1800_R1738000_rho2500_pmax15'}
           {'dV_MTshape_in360_out3600_R1738000_rho2500_pmax50'}
           {'dV_MTshape_in720_out3600_R1738000_rho2500_pmax70'}
           {'dV_MTshape_in1080_out3600_R1738000_rho2500_pmax70'}
           {'dV_MTshape_in2160_out3600_R1738000_rho2500_pmax70'}]

% choose the desired file via the index ----------------------------------
fileName = strcat(models{1},'.bshc');       

% read the SHCs from the file --------------------------------------------
[n m C S]      = read_SHCs_binary_bshc_format2014(fileName);
% now the SHCs are loaded into Matlab and can be further processed.

% compute dimensionless potential degree variances -----------------------
[n_vec sigma2] = degree_variances2012(n,m,C,S,3600);

% plot the degree variances ----------------------------------------------
figure(1);
plot(n_vec,sigma2);
grid on
set(gca,'XMinorTick','off','YMinorTick','off');
set(gca,'XMinorGrid','off','YMinorGrid','off');
set(gca,'XScale','linear','YScale','log');
xlabel(' Harmonic degree');
ylabel('Potential degree variance');
ylim([10^(-35) 10^(-5)]); xlim([0 3600]);
title(['Degree variances of  SHCs from file ' strrep(fileName,'_','-')] )
infostring = 'For full details see Hirt & Kuhn (2017), JGR Planets 122, doi:10.1002/2017JE005298.'
text(10,10^-28,infostring);
print(strcat('Degree_variances_',fileName,'.png'),'-r300','-dpng')
 
    
  