
%-------------------------------------------------------------------------
% read SHCs in binary format (bshc)
%-------------------------------------------------------------------------
% Reading large SHC files in binary format. The data is stored
% as double and big endian. The file structure is:
%
% n_min n_max   (minimum and maximum harmonic degree)
% C-coefficients ascending in degree, followed by the order,
% S-coefficients ascending in degree, followed by the order,
%
% Example for a degree-10800 file:
% 0 10800
% C(0,0), C(1,0), C(1,1), C(2,0), C(2,1), ... C(10800,10799), C(10800,10800),
% S(0,0), S(1,0), S(1,1), S(2,0), S(2,1), ... S(10800,10799), S(10800,10800).
%
% Input: filename
% Output: triangle of SHCs and maximum degree nmax
%
%-------------------------------------------------------------------------
% Christian Hirt, Curtin University
% Moritz Rexer, TU Munich
% ------------------------------------------------------------------------

function [SHCs,nmax] = read_SHCs_bshc2tri(filename_in)

 fid = fopen(filename_in,'rb')
 [data, cnt] = fread(fid, inf,'double');
 fclose(fid); 
 
 % extract min and max SHC
 input_nmin = data(1);
 input_nmax = data(2);
 
 %work out number of SHC records for C and S
 
 n_down = ((input_nmin-1)+1)*((input_nmin-1)+2)/2;
 n_up = ( input_nmax+1)*( input_nmax+2)/2;
 n_rows = n_up-n_down;
 offset = 2; % for correct computation of indices to access C, S blocks
 
 % access the C and S data blocks
 C = data(offset+1:offset+n_rows);
 S = data(offset+n_rows+1:offset+2*n_rows);
 
 % build new n, m vectors
 
 % build matrices with degrees and orders
 a = 0:input_nmax;          % vector from min to to max_degree
 b = a'*ones(1,input_nmax+1); % extend this into matrix
 degree_matrix = tril(b);            % set upper triangle to zero
 order_matrix = tril(b');
 clear b n m

% disp('Build n m vectors')
 ix1 = 1; ix2 = 1+input_nmin;
 
 for i = input_nmin+1: input_nmax+1
   % i
    n(ix1:ix2) = degree_matrix(i,1:i);
    m(ix1:ix2) = order_matrix(i,1:i);
    ix1 = ix2 +1; ix2 = ix1 + i;
 end

 % bring coefficients into triangle matrix
 [SHCs]=vec2tri(n,m,C,S);
 
nmax=max(n);
 
end


function [tri]=vec2tri(l,m,Cnm,Snm)

% Function creates a matrix with the coefficients Cnm and 
% Snm ordered in triangle shape. Sine-associated coefficients are assigned
% to the left, cosine-associated coefficients are assigned to the right
% side of the matrix.
%
% Input: l - degree vector 
%        m - order vector
%        Cnm,Snm - coefficient's vector
%
% Out  : tri - triangle matrix
%
% Moritz Rexer, IAPG/IAS - 23-03-2015

% maximum degree and order 
L=max(l);M=max(m);
% initialize triangle matrix
tri=zeros(L+1,(2*M)+1);

% set values to triangle
for i=1:length(l)
    % SNM
    tri(l(i)+1,M+1-m(i))=Snm(i);
    % CNM
    tri(l(i)+1,M+1+m(i))=Cnm(i);
end
end