function [gh, t, params] = read_shc(filename)
% function [gh, t, params] = read_shc(filename);
% reads Gauss coefficients stored in Swarm SCARF SHC format

% February 2011, Nils Olsen, DTU Space


if exist(filename, 'file') == 0
    error(['coefficient file ', filename, ' does not exist'])
end


fid=fopen(filename, 'r');
line = fgetl(fid);

% skip comment lines at beginning of each block
while strcmp(line(1), '#')
    line = fgetl(fid);
end
% read header part (2 lines)
tmp = sscanf(line, '%d %d %d %d %d');
params.N_min = tmp(1);
params.N_max = tmp(2);
params.N_times = tmp(3);
params.sp_order = tmp(4);
params.N_steps = tmp(5);
t = sscanf(fgetl(fid), '%f', params.N_times);

gh = zeros(params.N_max*(params.N_max+2), params.N_times);
k = params.N_min^2;
for n = params.N_min:params.N_max
    tmp = sscanf(fgetl(fid), '%f');
    if n ~= tmp(1); error('error in n'); end
    if 0 ~= tmp(2); error('error in m'); end
    gh(k, :) = tmp(3:params.N_times+2);
    k = k + 1;
    for m = 1:n
        tmp = sscanf(fgetl(fid), '%f');
        if n ~= tmp(1); error('error in n'); end
        if m ~= tmp(2); error('error in m'); end
        gh(k, :) = tmp(3:params.N_times+2);
        k = k + 1;
        tmp = sscanf(fgetl(fid), '%f');
        if n ~= tmp(1); error('error in n'); end
        if m ~= -tmp(2); error('error in n'); end
        gh(k, :) = tmp(3:params.N_times+2);
        k = k + 1;
    end
end

fclose(fid);
