
%-------------------------------------------------------------------------
% MGM2011 read test 
%------------------------------------------------------------------------- 
%
% This Matlab-script reads and displays the content of the 12 MGM2011
% product files that can be downloaded from the MGM2011 website.
%
% 
%-------------------------------------------------------------------------
% Dr Christian Hirt
% Western Australian Centre for Geodesy
% Curtin University, Perth
%-------------------------------------------------------------------------
% Version of 29/09/2011
%-------------------------------------------------------------------------
clear all
close all

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

%-------------------------------------------------------------------------
% Paths
%-------------------------------------------------------------------------
%mgm2011data_dirname = 'MGM2011_data_files/'
%mgm2011img_dirname  = 'MGM2011_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 mgm2011 data files 
%-------------------------------------------------------------------------
% file name convention for publication:
% Productclass_functional_[remarks].[spectralband].unit.resolution[arcmin].bin

% filenames
vec_productfile_names = {
 'MGM2011_surface_gravity_minus370000mgal.unit_mgal.3m'
 'MGM2011_freeair_gravity.unit_mgal.3m'
 'MGM2011_NSDOV.unit_arcsec.3m'
 'MGM2011_EWDOV.unit_arcsec.3m'
 'MGM2011_normal_gravity_minus370000mgal.unit_mgal.3m'
 'MRTM85_gravity.unit_mGal.3m'
 'MRTM85_areoid.unit_dm.3m'
 'MRTM85_NSDOV.unit_arcsec.3m'
 'MRTM85_EWDOV.unit_arcsec.3m'
 'MRO110B2_gravity_atMOLA3D.band2_85.unit_mGal.3m'
 'MRO110B2_NSDOV_atMOLA3D.band2_85.unit_arcsec.3m'
 'MRO110B2_EWDOV_atMOLA3D.band2_85.unit_arcsec.3m'
}

  
% column 1: conversion factor to convert data from file unit to
%           basis units mgal, m, arc seconds
% column 2: conversion constant
% multiplication and constants
vec_conversion_constants = [
    1     370000  %to mgal  MGM2011 products
    1     0       %to mgal
    1     0       %to sec
    1     0       %to sec
    1     370000  %to mgal  MGM2011 normalgrav
    1     0       %to mgal  MRTM85 ----------
    0.1   0       %dm to m
    1     0       %to sec
    1     0       %to sec
    1     0       %to mgal  MRO110B2 --------
    1     0       %to sec 
    1     0       %to sec
    ]
   
   
% product class, functional name, unit
vec_class_funct_unit = {
    'MGM2011' 'surface gravity'   '[mgal]'
    'MGM2011' 'free-air gravity'  '[mgal]'      
    'MGM2011' 'NSDoV'             '[arc sec]'     
    'MGM2011' 'EWDoV'             '[arc sec]'   
    'Normal'  'gravity'           '[mgal]' 
    'MRTM85'  'gravity'           '[mgal]'
    'MRTM85'  'geoid'             '[m]'  
    'MRTM85'  'NSDoV'             '[arc sec]' 
    'MRTM85'  'EWDoV'             '[arc sec]'  
    'MRO110B2' 'gravity (band 2 to 85)'          '[mgal]' 
    'MRO110B2' 'NSDoV (band 2 to 85)'            '[arc sec]' 
    'MRO110B2' 'EWDoV (band 2 to 85)'            '[arc sec]'
    }


%-------------------------------------------------------------------------
% Load and display MGM2011 data files 
%-------------------------------------------------------------------------
%list_of_files = [1 2]  
list_of_files = 1:12
% This list can be customized 
% 1  = first file  (mgm2011 surface gravity)
% 2  = second file (mgm2011 free-air gravity)

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 = [classname ' ' functname ' ' unitname]
  title(titlestring)
  colorbar
  grid on
  ylim([-90 90])
  xlim([-180 180])
  print(strcat(inputfilename,'.png'),'-r200','-dpng')

  clear data
    
end

disp('MGM2011 read test finished.')












