%-------------------------------------------------------------------------
% Testdriver for ERTM2160 access and use
%
% Requires ERTM2160_2013_v1.m for seamless access of ERTM2160 binary files
%
% Christian Hirt and Sten Claessens,
% Western Australian Centre for Geodesy,
% Curtin University, Perth, Australia
% 2013-06-05, last modified 2014-10-08
%-------------------------------------------------------------------------

clear all
close all

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

% Insert path where ERTM2160 files are located ----------------------------
pathERTM2160 = 'K:\ERTM2160\data';

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

% Define desired target area ---------------------------------------------
minlon1 =  10.5;  % 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 =  48.0;  % geodetic latitude of northern edge (negative for southern hemisphere)

% NB: In this example, the files N45E010.suffix are assumed to be available
%     at the following locations
%     'E:\ERTM2160\data\dg\N45E010.dg'     (gravity disturbances)
%     'E:\ERTM2160\data\ga\N45E010.ga'     (gravity accelerations)
%     'E:\ERTM2160\data\xi\N45E010.xi'     (NS DoVs)
%     'E:\ERTM2160\data\eta\N45E010.eta'   (EW DoVs)
%     'E:\ERTM2160\data\geoid\N45E010.ha'  (quasigeoid)
%     'E:\ERTM2160\data\dem\N45E010.dem'   (terrain height)

% Factor applied on ERTM2160 basis resolution of 0.002 deg ---------------- 
% Use facX=facY=1 to extract ERTM2160 data over target area *without*
% interpolation.
% Use facX=facY=2 to downsample ERTM2160 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
[x1,y1,z1 ] = ertm2160_2013_v2('gravity',      minlon1,maxlon1,minlat1,maxlat1,facX,facY,pathERTM2160,ndv,method); 
[x1,y1,z3 ] = ertm2160_2013_v2('xi',           minlon1,maxlon1,minlat1,maxlat1,facX,facY,pathERTM2160,ndv,method); 
[x1,y1,z4 ] = ertm2160_2013_v2('eta',          minlon1,maxlon1,minlat1,maxlat1,facX,facY,pathERTM2160,ndv,method); 
[x1,y1,z5 ] = ertm2160_2013_v2('geoid',        minlon1,maxlon1,minlat1,maxlat1,facX,facY,pathERTM2160,ndv,method); 
[x1,y1,z6 ] = ertm2160_2013_v2('dem',        minlon1,maxlon1,minlat1,maxlat1,facX,facY,pathERTM2160,ndv,method); 

% Select test point coordinates ------------------------------------------
% Select coordinates of a single point for which functionals will be
% calculated. The point must be inside the area selected above. The
% results will be displayed in Matlab's command window.
% Test point Zugspitze summit, German Alps.
test_phi = 47.421;  % geodetic latitude
test_lam = 10.984;  % longitude
methodtp = 'cubic'; % interpolation method for test point calculation

% ------------------------------------------------------------------------
%    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(' ERTM2160 Gravity disturbances (mGal)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
end

% Acceleration
if (exist('z2','var')==1) 
figure(2)
h=imagesc(x1(1,1:end),y1(1:end,1),z2); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' ERTM2160 Accelerations (mGal)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
end

% North-South DoV (xi)
if (exist('z3','var')==1) 
figure(3)
h=imagesc(x1(1,1:end),y1(1:end,1),z3); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' ERTM2160 North-South DoV (arc-sec)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
end

% East-West DoV (eta)
if (exist('z4','var')==1) 
figure(4)
h=imagesc(x1(1,1:end),y1(1:end,1),z4); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' ERTM2160 East-West DoV (arc-sec)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
end

% Geoid
if (exist('z5','var')==1) 
figure(5)
h=imagesc(x1(1,1:end),y1(1:end,1),z5); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' ERTM2160 Geoid (m)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
end

if (exist('z6','var')==1) 
figure(6)
h=imagesc(x1(1,1:end),y1(1:end,1),z6); set(gca,'XDir','normal','YDir','normal');
axis equal; axis ([minlon1 maxlon1 minlat1 maxlat1]); 
title(' Terrain height (m)'); colorbar(); hold on; grid on
xlabel(' Longitude [deg]'); ylabel('Latitude [deg]')
end

% ------------------------------------------------------------------------
%    CALCULATION AT A SINGLE POINT OF INTEREST
% ------------------------------------------------------------------------
disp('------------------------------------------------')
disp(' Interpolated ERTM2160 functionals at test point')
disp('------------------------------------------------')

% Interpolate ERTM2160 functionals for the test point------ ---------------
if (exist('z1','var')==1) 
    test_dg  = interp2(x1,y1,z1,test_lam,test_phi,methodtp); 
    fprintf('Gravity disturbance:  %10.4f mGal\n',test_dg)
end
if (exist('z2','var')==1)
    test_ga  = interp2(x1,y1,z2,test_lam,test_phi,methodtp); 
    fprintf('Acceleration:         %10.2f mGal\n',test_ga)
end
if (exist('z3','var')==1) 
    test_xi  = interp2(x1,y1,z3,test_lam,test_phi,methodtp);
    fprintf('North-South DoV:      %10.4f arcsec\n',test_xi)
end
if (exist('z4','var')==1) 
    test_eta = interp2(x1,y1,z4,test_lam,test_phi,methodtp);
    fprintf('East-West DoV:        %10.4f arcsec\n',test_eta)
end
if (exist('z5','var')==1) 
    test_ha  = interp2(x1,y1,z5,test_lam,test_phi,methodtp);
    fprintf('Geoid:                %10.4f m\n',test_ha)
end

% ------------------------------------------------------------------------
%    ASCII OUTPUT
% ------------------------------------------------------------------------
% A text file containing longitude, latitude and functional values in 
% ASCII format can be created for each functional using the code below.
% By default, this is commented out, as it may result in very large text
% files if a large area is selected with a low downsampling factor.
% This example (when uncommented) prints gravity disturbances (z1) in the
% selected area to the text file 'gravity.txt'.

%outputfile='K:\ERTM2160\gravity.txt'
%fid=fopen(outputfile,'wt');
%[m,n]=size(x1);
%for i=1:m
%     for j=1:n
%         fprintf(fid,'%8.3f %7.3f %6.1f \n',x1(i,j),y1(i,j),z1(i,j));
%     end
% end


