%-------------------------------------------------------------------------
% Testdriver to access the SRTM2gravity model of terrain corrections
% and short-scale gravity effects
%-------------------------------------------------------------------------
%
% The use of this testdriver requires
%
% 1) th s2g2018_v2.m access function enabling seamless access to the
%    SRTM2gravity binary files
% 2) the s2g2018.mat index file
% 3) all 1° x 1° binary files over your area of interest,      
%    stored in terms of 30 deg x 30 deg regions
%
% Important: Please consult the readme-file for further details on the model,
% the specifities of the two components "full-scale" and "residual" gravity,
% possible applications, data handling and storage formats.
%
%-------------------------------------------------------------------------
% Christian Hirt 
% Institute for Astronomical and Physical Geodesy (IAPG) &
% Institute for Advanced Study (IAS)
% Technical University Munich, Munich, Germany
% 2019-01-18, last edited 2019-04-25
%-------------------------------------------------------------------------

clear all
close all

% ------------------------------------------------------------------------
%    INPUT SELECTION
% ------------------------------------------------------------------------

% Insert path where the SRTM2gravity files are located -------------------
% Please use preferably the absolute path as in the example
pathS2g = 'E:\SRTM2gravity_Release\data\';

% the expected sub-directory structure is 
%      \SRTM2gravity_Release\data\FullScaleGravity
%      \SRTM2gravity_Release\data\ResidualGravity
%      \SRTM2gravity_Release\data\FullScaleGravity\N30E000\N47E015_full.bin
%                                                         \N47E016_full.bin
%                                                          ...
%      \SRTM2gravity_Release\data\FullScaleGravity\N30E000\N47E015_res.bin
%                                                         \N47E016_res.bin
%                                                          ...

% Choose no-data value for points where SRTM2gravity is unavailable ------
ndv     = NaN;  %NaN produces best results for plotting

% Define desired target area ---------------------------------------------
minlon1 =  10.0;  % longitude of western edge (negative for western hemisphere)
maxlon1 =  12.0;  % longitude of eastern edge (negative for western hemisphere)
minlat1 =  47.0;  % geodetic latitude of southern edge (negative for southern hemisphere) 
maxlat1 =  49.0;  % geodetic latitude of northern edge (negative for southern hemisphere)

% Factor applied on SRTM2gravity basis resolution of 3/3600 degrees ------ 
% Use facX=facY=1 to extract SRTM2gravity data over target area *without*
% interpolation (recommanded setting)
% Use facX=facY=2 to downsample data by factor 2 using the
% interpolation method specified below. Other factors can also be used. 
facX =1.0;
facY =1.0;
% Interpolation methods available are 'nearest', 'linear', 'spline' and
% 'cubic'. 
method = 'cubic';

% Select functionals of interest -----------------------------------------
functional1 = 'gravity_full'  % full-scale gravity signal of the topography
                              % = complete spherical terrain correction
functional2 = 'gravity_residual' % residual gravity signal of the topography
                                % at spatial scales less than 10 km
                                % also known as RTM gravity in geodesy

% ------------------------------------------------------------------------
% SRTM2gravity DATA ACCESS via function s2g2018_v2
% ------------------------------------------------------------------------              
% please comment out if the functional is not required
[x1,y1,z1] = s2g2018_v2(functional1, minlon1,maxlon1,minlat1,maxlat1,...
    facX,facY,pathS2g,ndv,method);

[x2,y2,z2] = s2g2018_v2(functional2, minlon1,maxlon1,minlat1,maxlat1,...
    facX,facY,pathS2g,ndv,method);

lamvec = x1(1,1:end); clear x1
phivec = y1(1:end,1); clear y1;

% ------------------------------------------------------------------------
%    MAP PRODUCTION
% ------------------------------------------------------------------------
% Gravity disturbance
figure(1)
h=imagesc(lamvec,phivec,z1); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' SRTM2gravity full-scale topographic gravity (mGal)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')

figure(2)
h=imagesc(lamvec,phivec,z2); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' SRTM2gravity residual topographic gravity (mGal)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')


