parseGPS.m

function [lat, lon, elev, decl, decl_info, gain_char, error_flag, comment] = parseGPS(GPS)

% Usage:  [lat,lon,elev,decl, decl_info, gain_char,error_flag,comment] = parseGPS(GPS)

lat = NaN;
lon = NaN;
elev = NaN;
decl = NaN;
rdate = NaN;
decl_info.igrf = [];
decl_info.igrf_conf = [];
decl_info.file_decl = decl;
gain_char = 'U';
error_flag = 0;
check_elev = 0;
comment = '';    

if nargin == 0
    error_flag = 1;
    comment = 'No input GPS string. ';
    return
end

% Find gain flag H(igh) or L(ow) or E(X) high or E(Y) high and save it
iGain = [find(GPS==200); find(GPS==204); find(GPS==216); find(GPS==217)];
if ~isempty(iGain)    
    % this will ignore a short time with gain changed ('unique' will not!)
    gain_char = char( median(GPS(iGain)) - 128 );
    gain_char_start = char( GPS(min(iGain)) - 128 );
    if gain_char ~= gain_char_start
       comment = ['Gain problem (changed during run?). ' comment];
       gain_char = 'U';
    end
end

% Find and parse all active (fixed) $GPRMC and $GPGGA messages
[GPRMC, GPGGA] = gps_scan(char(GPS'));

% no GPS messages
if isempty(GPRMC)
    error_flag = 1;
    comment = 'No GPS. ';
    return
end

iFixed = find(strcmp(GPRMC(:,3),'A'));

% iFixed contains the indices of all valid GPS messages in GPRMC array
if isempty(iFixed)
    error_flag = 2;
    comment = ...
     ['GPS NEVER FIXED; Lat, Lon and Time may NOT be correct. ' comment];
    return
end

% Save the GPS information into arrays. Need a decimal year to use with the
% Aerospace Toolbox to compute the declination from the WMM. Need to
% replace with something that isn't toolbox dependent. Until this is done,
% keeping the old declination computation in gps_info.m for backwards
% compatibility with old GPS strings. Will produce NaN decl with new GPS.
% When updated, the decimal year decyear() could also be implemented.
j = 1;
for k = 1:length(iFixed)
    [info,err] = gps_info(GPRMC(iFixed(k),:),0);
    if err > 0
        continue
    end
    lat(j) = info.lat;
    lon(j) = info.lon;
    decl(j) = info.decl;
    rdate(j) = datenum(info.date);
    j = j + 1;
end

j = 1;
for k = 1:length(iFixed)
    [info,err] = gps_info(GPGGA(iFixed(k),:),0);
    if err > 0
        continue
    end
    elev(j) = info.elev;
    j = j+1;
end

% If we only have a few elevation values, have the user fix them!
if length(elev(~isnan(elev))) < 100
    check_elev = 1;
end

% Compute a robust average; use median but make it bulletproof
lat = median(lat(~isnan(lat)));
lon = median(lon(~isnan(lon)));
decl = median(decl(~isnan(decl)));
elev = median(elev(~isnan(elev)));

% At this stage we only obtain the declination from the GNSS sequence. Will
% later overwrite with IGRF if appropriate.
decl_info.igrf.version = 13;
decl_info.igrf.toolbox = 'GNSS sequence as implemented';
decl_info.file_decl = decl;

if isnan(elev)
    elev = 0;
    error_flag = 4;
    comment = ['No valid elevation recorded in the GPS. ' comment];
elseif check_elev
    error_flag = 4;
    comment = ['GPS elevation not accurate. Please check and fix location and elevation. ' comment];
end