function [GM, Re, degree, C, S] = astReadSHAFile( filename )
%
%   Inputs are:
%   FILENAME :a string specifying the name of the ASCII SHA-file containing
%            the normalized coefficients of the spherical harmonic gravity
%            model.  
%
%   Output includes:
%   GM     :a scalar value specifying the planetary gravitational parameter
%          in meters cubed per second squared.  
%   RE     :a scalar value specifying the planetary equatorial radius in
%          meters. 
%   DEGREE :a scalar value specifying the maximum degree of the spherical
%          harmonic gravity model.
%   C      :a (DEGREE+1)-by-(DEGREE+1) matrix containing the normalized Cnm
%          coefficients of the spherical harmonic gravity model.
%   S      :a (DEGREE+1)-by-(DEGREE+1) matrix containing the normalized Snm
%          coefficients of the spherical harmonic gravity model.

%   Copyright 2009-2017 The MathWorks, Inc.

% open data file 
fid = fopen( filename, 'r' );
% read header line
header_line = fscanf(fid, '%e, %e, %e, %d, %d, %d, %e, %e', [8 1] );
% read data matrix
matrix = fscanf(fid, '%d, %d, %e, %e, %e, %e', [6 inf] );
% close data file
fclose(fid);

if header_line(6) ~=1
    error(message('aero:readSHAFileGravityModel:unnormalizedCoeff'));
end

% assign outputs and convert to expected units
Re     = header_line(1)*10^3; % meters
GM     = header_line(2)*10^9; % meters^3/s^2
degree = header_line(4);

% set degree index (zero vs one based indices)
n = matrix(1,:) + 1;
% set order index (zero vs one based indices)
m = matrix(2,:) + 1;

% put normalized coefficients into matrices
for k = length(n):-1:1
    C(n(k),m(k)) = matrix(3,k);
    S(n(k),m(k)) = matrix(4,k);
end
