function  s = SwarmL2_MIO_SHA_Read(filename)
% s = SwarmL2_MIO_SHA_Read(filename);
%
% Purpose:
% Read Swarm L2 MIO_SHA model from "*.DBL" file.
%
% Input:
%    filename           name of file to read
%
% Output: s is structure with elements
%             m_e_b(:,:)            primary Sq model for r < rho_Sq (REAL coefficients)
%             m_i(:,:)              secondary Sq model (REAL coefficients)
%             nmax, mmax            max degree and order in dipole coord.
%             p_vec(:)              diurnal wavenumbers
%             s_vec(:)              seasonal wavenumbers
%             theta_NGP             co-latitude of Northern dipole (in deg)
%             phi_NGP               longitude of Northern dipole (in deg)
%             h                     altitude of Sq current layer (in km)
%             N                     coefficient of the F10.7 dependence
%
% A. Chulliat, 20/05/2011

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

fid = fopen(filename);

% skip comment lines

tline = fgetl(fid);
while strcmp(tline(1), '#')
    tline = fgetl(fid);
end

% read header part (1 line)

f = sscanf(tline, '%i %i %i %i %i %i %f %f %f %f');
s.nmax = f(1);
s.mmax = f(2);
pmin = f(3);
pmax = f(4);
smin = f(5);
smax = f(6);
s.theta_NGP = f(7);
s.phi_NGP = f(8);
s.h = f(9);
s.N = f(10);


s.p_vec = pmin:pmax;
s.s_vec = smin:smax;

% read model coefficients

N_nm = s.mmax*(s.mmax+2) + (s.nmax-s.mmax)*(2*s.mmax+1);
N_sp = length(s.p_vec)*length(s.s_vec);

for i_model = 1:2

    m_MIO = zeros(N_nm,2*N_sp);
        
    for i_nm = 1:N_nm
        tline = fgetl(fid);
        if ~ischar(tline)
            ErrorMessage = 'Wrong number of model coefficients!';
            break
        end
        f = sscanf(tline,'%f',2*N_sp+2);
        m_MIO(i_nm,:) = f(3:2*N_sp+2);
        ErrorMessage = '';
    end
    
    if (i_model == 1)
        s.m_e_b = m_MIO;
    elseif (i_model == 2)
        s.m_i = m_MIO;
    end
    
end

fclose(fid);


