function [modelname, GM, a, Cnm, Snm] = readGRACE(filename);
%[GM, a, Cnm, Snm] = function readGRACE(filename);
%   filename is an optional argument
%
%Read the GRACE gravity field from one of the various files for the model.
%All of the files appear to have the same format.
%The file is assumed to be in the current directory unless the name
%of the file is passed as an argument.
%
%At present this function returns only the static components of the gravity
%field.
%
%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 'GGM02C.GEO' if it is not given.

if ( nargin > 1 )
    disp('Usage: readGRACE(filename), with filename optional');
elseif ( nargin == 0 )
    filename = 'GGM02C.GEO';
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(200,201);
Cnm = zeros(200,201);

% Open the file and read the coefficients

[fid, err] = fopen(filename, 'r');
if ( fid == -1 )
    disp(err);
end

% Read the three header lines

header_line = fgetl(fid);

tline = fgetl(fid);
[modelname, GM, a, epoch] = strread(tline, '%s %f %f %f')

header_line = fgetl(fid);

% Now read the coefficients from the file

A = fscanf(fid, '%6s %3d %3d %f %f %f %f %f');
fclose(fid);

% A now contains all the numerical values in the file, in the order
%   (string), n, m, Cnm, Snm, sig_Cnm, sig_Snm.
% A is stored as a vector. There are 13 elements for every line of the
%  file.

nlines = length(A)/13;

for i = 1:nlines,
    ioff = (i-1)*13;
    n = A(ioff+7);
    m = A(ioff+8);
    
    Cnm(n,m+1) = A(ioff+9);
    Snm(n,m+1) = A(ioff+10);
end
