function [GM, a, Cnm, Snm] = readEGM96(filename);
%[GM, a, Cnm, Snm] = function readEGM96(filename);
%   filename is an optional argument
%
%Read the EGM96 gravity field from the file egm96_to360.ascii
%The file is assumed to be in the current directory unless the name
%of the file is passed as an argument
%
%values returned:
% GM = 3986004.415*10^8 m^3/s^2
% a  = 6378136.3        m
% Cnm(N,N) = coefficients for cosine longitude term
% Snm(N,N) = coefficients for sine longitude term
% maximum degree is given by length(Cnm(:,1))

% Check the number of arguments to see if filename is given. Set filename
% to 'PGM2000A.TO360.ASCII' if it is not given.

if ( nargin > 1 )
    disp('Usage: readEGM96(filename), with filename optional');
elseif ( nargin == 0 )
    filename = 'egm96_to360.ascii';
end

% Define the constants used in this gravity field:
GM = 3.986004415*10^14;    % m^2/sec^3
a  = 6378136.3;            % m

%     Start by setting Cnm and Snm to zeros (degree and order 360)
Snm = zeros(360,360);
Cnm = zeros(360,360);

% Open the file and read the coefficients

[fid, err] = fopen(filename, 'r');
if ( fid == -1 )
    disp(err);
end
A = fscanf(fid, '%d %d %g %g %g %g');
fclose(fid);

% A now contains all the numerical values in the file, in the order
%   n, m, Cnm, Snm, sig_Cnm, sig_Snm.
% A is stored as a vector. So we find the number of lines in the file (6
% values per line), and extract n, m, Cnm, Snm from each line.

nlines = length(A)/6;

for i = 1:nlines,
    ioff = (i-1)*6;
    n = A(ioff+1);
    m = A(ioff+2);
    
    Cnm(n,m+1) = A(ioff+3);
    Snm(n,m+1) = A(ioff+4);
end
