
%------------------------------------------------------------------------
% LGM2011 read test 
%------------------------------------------------------------------------ 
%
% This Matlab-script reads and displays the content of the 14 LGM2011
% product files that can be downloaded from the LGM2011 website.
%
%------------------------------------------------------------------------
% Dr Christian Hirt
% Western Australian Centre for Geodesy
% Curtin University, Perth
%------------------------------------------------------------------------

clear all
close all

set(0,'DefaultAxesFontSize',12);
set(0,'DefaultAxesLineWidth',0.8);
set(0,'DefaultLineLineWidth',0.8); 

%-------------------------------------------------------------------------
% Paths
%-------------------------------------------------------------------------
%lgm2011data_dirname = 'LGM2011_data_files/'
%lgm2011img_dirname  = 'LGM2011_pictures/'
%-------------------------------------------------------------------------
% Grid definition 
%-------------------------------------------------------------------------

resolution = 3/60             % in degree

minlon = -180+resolution/2
maxlon = +180-resolution/2
minlat =  -90+resolution/2
maxlat =  +90-resolution/2

vec_lon = minlon:resolution:maxlon;
vec_lat = minlat:resolution:maxlat;
n = length(vec_lat);
m = length(vec_lon);

%-------------------------------------------------------------------------
% Description of the LGM2011 data files 
%-------------------------------------------------------------------------
% file name convention for publication:
% Productclass_functional_[remarks].[spectralband].unit.resolution[arcmin].bin

% filenames
vec_productfile_names = {
 'LGM2011_surface_gravity_minus162468.unit_mgal.3m'
 'LGM2011_freeair_gravity.unit_mgal.3m'
 'LGM2011_selenoid.unit_dm.3m'
 'LGM2011_NSDOV.unit_arcsec.3m'
 'LGM2011_EWDOV.unit_arcsec.3m'
 'LGM2011_normal_gravity_minus162468.unit_mgal.3m'
 'LRTM70_gravity.unit_mGal.3m'
 'LRTM70_selenoid.unit_dm.3m'
 'LRTM70_NSDOV.unit_arcsec.3m'
 'LRTM70_EWDOV.unit_arcsec.3m'
 'SGM100i_gravity_atLOLA3D.band2_70.unit_mGal.3m'
 'SGM100i_selenoid_atR1738km.band2_70.unit_dm.3m'
 'SGM100i_NSDOV_atLOLA3D.band2_70.unit_arcsec.3m'
 'SGM100i_EWDOV_atLOLA3D.band2_70.unit_arcsec.3m'
}
  
% column 1: conversion factor to convert data from file unit to
%           basis units mgal, m, arc seconds
% column 2: conversion constant
vec_conversion_constants = [
    1     162468  % mgal     LGM2011 surface gravity
    1     0       % mgal             freeair gravity
    0.1   0       % dm to m          selenoid
    1     0       % sec              NS DoV
    1     0       % sec              EW DoV
    1     162468  % mgal     LGM2011 normal gravity
    1     0       % mgal     LRTM70  gravity
    0.1   0       % dm to m          selenoid
    1     0       % sec              NS DoV
    1     0       % sec              EW DoV
    1     0       % mgal     SGM100i gravity
    0.1   0       % dm to m          selenoid
    1     0       % sec              NS DoV
    1     0       % sec              EW DoV
    ]
   
% product class, functional name, unit
vec_class_funct_unit = {
    'LGM2011' 'surface gravity'  '[mgal]'
    'LGM2011' 'free-air gravity'  '[mgal]'      
    'LGM2011' 'selenoid' '[m]' 
    'LGM2011' 'NSDoV'    '[arc sec]'     
    'LGM2011' 'EWDoV'    '[arc sec]'  
    'Normal'  'gravity'  '[mgal]'   
    'LRTM70'  'gravity'  '[mgal]'
    'LRTM70'  'selenoid' '[m]'  
    'LRTM70'  'NSDoV'    '[arc sec]' 
    'LRTM70'  'EWDoV'    '[arc sec]' 
    'SGM100i' 'gravity  (band 2 to 70)'  '[mgal]'
    'SGM100i' 'selenoid (band 2 to 70)' '[m]'  
    'SGM100i' 'NSDoV (band 2 to 70)'    '[arc sec]' 
    'SGM100i' 'EWDoV (band 2 to 70)'    '[arc sec]' 
    }


%-------------------------------------------------------------------------
% Load and display LGM2011 data files 
%-------------------------------------------------------------------------
list_of_files = [1 2 6 7 11]  
%list_of_files = 1:14
% This list can be customized 
% 1  = first file  (LGM2011 surface gravity)
% 2  = second file (LGM2011 free-air gravity)
% 14 = last file   (SGM100i EW DoV)

for i = 1:length(list_of_files)
    
  ix =  list_of_files(i)
    
  inputfilename  = char(vec_productfile_names(ix))
  fac            = vec_conversion_constants(ix,1)
  const          = vec_conversion_constants(ix,2)
  classname      = char(vec_class_funct_unit(ix,1))
  functname      = char(vec_class_funct_unit(ix,2))
  unitname       = char(vec_class_funct_unit(ix,3))

  
  % Read data and convert --------------------------
  disp('Read input file')
  fid = fopen(strcat(inputfilename,'.bin'),'rb')
  data = fread(fid,n*m,'int16','ieee-be');
  fclose(fid);
  data = data*fac+const;
  data = reshape(data,n,m);

  % Display ----------------------------------------
  figure(1)  
  imagesc(vec_lon,vec_lat,data);
  set(gca,'XDir','normal','YDir','normal');
  xlabel('Longitude [deg]')
  ylabel('Latitude [deg]')
  titlestring = strcat(classname,' ',functname,' ',unitname)
  title(titlestring)
  colorbar
  grid on
  ylim([-90 90])
  xlim([-180 180])
  print(strcat(inputfilename,'.png'),'-r200','-dpng')

  clear data
    
end

disp('LGM2011 read test finished.')












