%-------------------------------------------------------------------------
% Testdriver for SRTM2gravity access and use
%-------------------------------------------------------------------------
%
% Requires 
% 1) s2g2018_v2.m for seamless access of SRTM2gravity binary files
% 2) data tiles available from the data server
%
% Please see the SRTM2gravity readme file for directory structures
% and all further details.
%
% Christian Hirt 
% TU Munich
% 2018-11-01, last edited 2010-04-24
%-------------------------------------------------------------------------


clear all
close all

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

% Insert path where GGMplus files are located ----------------------------
pathS2g = 'D:\SRTM2gravHarvest\SRTM2gravity_Release\data';
%
% Choose no-data value for points where GGMplus is unavailable -----------
ndv     = NaN;  %NaN produces best results for plotting

% Define desired target area ---------------------------------------------
minlon1 =  7;  % longitude of western edge (negative for western hemisphere)
maxlon1 =  9;  % longitude of eastern edge (negative for western hemisphere)
minlat1 =  45;  % geodetic latitude of southern edge (negative for southern hemisphere) 
maxlat1 =  47;  % geodetic latitude of northern edge (negative for southern hemisphere)


% Factor applied on GGMplus basis resolution of 0.002 deg ---------------- 
% Use facX=facY=1 to extract GGMplus data over target area *without*
% interpolation.
% Use facX=facY=2 to downsample GGMplus 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 -----------------------------------------
% Functionals not of interest can be commented out

disp(' Zuvor')
%[x1,y1,z1] = s2g2018_v1('gravity_residual', minlon1,maxlon1,minlat1,maxlat1,...
%    facX,facY,pathS2g,NaN,'linear');
[x1,y1,z1] = s2g2018_v1('gravity_residual', minlon1,maxlon1,minlat1,maxlat1,...
    facX,facY,pathS2g,NaN,'linear');

disp(' Danach')


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

figure(2)
h=imagesc(x1(1,1:end),(y1(1:end-1,1)+y1(2:end,1))/2,(z1(1:end-1,:)-z1(2:end,:))/90*100); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' S2g full-scale topographc gravity Hz gradient (mGal/100m)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
set(gca,'clim',[-2 2]);
end

