%% access_Earth2014_shcs10800.m
%
% This script shows how to access the surface spherical harmonic
% coefficients files to degree 10,800 of the Earth2014 models in Matlab. 
% 
% In detail, the script :
% 
%  a) reads the binary shc-files (.bshc)
%  b) plots the spectrum logarithmic in a triangle matrix
%  c) plots the power spectrum (degree variances) from SHCs

% dependencies: - read_SHCs_bshc2tri.m
%
% matlab version:  8.2.0.701 (R2013b)
% 
% authors : Christian Hirt, Moritz Rexer
% date : 2015 / 03 / 24
% contact: m.rexer@tum.de
% -------------------------------------------------------------------------

clear all; close all;
clc

% directory of the Earth2014 files 

input_dir='../data_1min/shcs_to10800/'

% Earth2014 models' file-names (.bshc)
mnames = [                                     % id
  'Earth2014.SUR2014.degree10800.bshc'; % 1  : SUR
  'Earth2014.BED2014.degree10800.bshc'; % 2  : BED
  'Earth2014.ICE2014.degree10800.bshc'; % 3  : ICE
  'Earth2014.RET2014.degree10800.bshc'; % 4  : RET
  'Earth2014.TBI2014.degree10800.bshc'; % 5  : TBI
];

% loop over all models
for id= 1 : 5 % choose desired models by id

% a) read SHC model from .bshc-file
filename_in = strcat(input_dir,strtrim(mnames(id,:)));
[SHC, nmax] =read_SHCs_bshc2tri(filename_in);

% b) plot triangle matrix
fig1=figure;
im=imagesc(log10(abs(SHC)));
cb=colormap('jet'); colormap([1 1 1; cb]); colorbar;
title({strcat('Earth2014 : ',strtrim(mnames(id,11:17)),' spectrum (log10)'); 'sin                   cos'});
ylabel('degree');xlabel('order');
set(gca,'XTick',[800 5800 10800 15800 21000]);
set(gca,'XTickLabel',[10000 5000 0 5000 10000]);
print(fig1,strcat(strrep(filename_in,'.bshc',''),'_spectrum.png'),'-r300','-dpng');

% c) plot the power spectrum (degree variances) from SHCs
cols=['r';'c';'m';'y';'g'];
fig2=figure;
semilogy(sum(SHC'.^2),cols(id)); hold on;
title({'Earth2014 degree variances (power spectrum)'}); grid on;
legend(strcat('Earth2014: ',' ',strtrim(mnames(id,11:17))));
xlabel('spherical harmonic degree');
ylabel('[m^2]');

print(fig2,strcat(strrep(filename_in,'.bshc',''),'_degvar.png'),'-r300','-dpng');
close all
end