function gps = read_gprm_CAPRICORN_2016(dfl,ddd,hhh)
%{

%%%% THIS VERSION FOR 1 HZ GPS DATA %%%%

reads 1Hz gps file specified by dfl and returns array with the
following columns:

1   jd_ref      10 Hz timestamp
2   cog         cog, deg
3   sog         sog, m/s
4   Ngps        filtered N component of ship speed, m/s
5   Egps        filtered E component of ship speed, m/s
6   lat         latitude (N positive)
7   lon         longitude (0-360 deg, east)

input parameters: dfl = string path to scs file
                  ddd = julian date
                  hhh = hour

gps files are reformatted to gprm using a python script.
GPRMC lines are stripped from the raw file and saved to a new gprm file.

Format for $GPRMC lines in the raw file:
0000078 $GPRMC,000000.00,A,5822.95924260,N,04515.92776422,W,1.58,223.96,151013,22.0,W,A*04
1 2 3                      4             5 6              7 8    9      10     11
PSD timecodes              lat             lon              spd  hdg    
%}
%% initialize
disp(['Reading gprm file for hour ',int2str(hhh)]);
flist = fopen(dfl);

temp = {};    % empty arrays for textscan
stx = [];
jd_ref = ddd+hhh/24+(0:1:3599)/60/60/24;   % exact 1Hz timestamp

%% read file
if feof(flist)==0,	% read entire file into cell array then scan for NMEA tags
    gpsdata = importdata(dfl);
    temp = gpsdata(cellfun(@(x) ~isempty(strfind(x,'$GPRMC')),gpsdata));
    
    for j=1:length(temp)
        test =textscan(temp{j,:},'%2f%2f%3f %*6c %*9c %*c %f %c %f %c %f %f %f %f %*c %*4c %*[^\n]','delimiter',', ','headerlines',0,'emptyvalue',NaN);
        try
            stx(:,j)=[test{1};test{2};test{3};test{4};double(test{5});test{6};double(test{7});test{8};test{9};test{10};test{11}];
        catch
            stx(:,j)=ones(11,1)*NaN;
        end;
    end;
    stx(:,isnan(stx(1,:)))=[];
    fclose(flist);
    jd_gps = ddd + (hhh + (stx(1,:)+(stx(2,:)+stx(3,:)/1000)/60)/60)/24;
else
    gps = NaN(3600,7);
    gps(:,1) = jd_ref';
    return
end;

%% format variables
glt = stx(4,:); 
gpslat_sign = stx(5,:);        % ascii code for "N" or "S"
gln = stx(6,:); 
gpslon_sign = stx(7,:);        % ascii code for "W" or "E"
gpsspeed1 = stx(8,:).*0.51444;  % convert to m/s
gpsspeed1(gpsspeed1>100) = NaN; % kill wierd values
gpshead1 = stx(9,:);          % in deg
gpshead1(gpshead1>361) = NaN; % kill wierd values
gpshead1 = unwrap(gpshead1*pi/180);  % unwrapped in radians for now...

bb = floor(glt/100);          % convert lat format to decimal degrees
glt = (glt-bb*100)/60;
glt = glt+bb;
ii = find(gpslat_sign==83); glt(ii) = -glt(ii);
gpslat1 = glt;

bbn = floor(gln/100);         % convert lon format to decimal degrees
gln = (gln-bbn*100)/60;
gln = gln+bbn;
ii = find(gpslon_sign==87);   % W longitudes
gln(ii) = -gln(ii) + 360;     % reformat to 0-360 deg
gpslon1 = gln;

%% despike variables - difference between successive points should be small
ii = find(abs(diff(gpslon1))>1); gpslon1(ii) = gpslon1(ii-1);
ii = find(abs(diff(gpslat1))>1); gpslat1(ii) = gpslat1(ii-1);

%% interp if necessary
if length(stx) < 3600 && ~isempty(gpslon1)
    gpslon1 = interp1(jd_gps,gpslon1,jd_ref);
    gpslat1 = interp1(jd_gps,gpslat1,jd_ref);
    gpshead1 = interp1(jd_gps,gpshead1,jd_ref);
    gpsspeed1 = interp1(jd_gps,gpsspeed1,jd_ref);   % 10 Hz sog
elseif length(gpsspeed1)> 3600     % else if long, truncate and screen for NaNs
    gpsspeed1 = gpsspeed1(1:3600);
    gpshead1 = gpshead1(1:3600);
    gpslat1 = gpslat1(1:3600);
    gpslon1 = gpslon1(1:3600);
end;

% Replace NaN with median of +/- 20 point range around NaN
ii=find(isnan(gpslon1));    % indices of NaNs
if ~isempty(ii)
    jj = (ii<21);       % NaNs in the first 20 points
    kk = (ii>3600-21); % NaNs in the last 20 points
    ll = ~xor(jj,kk);   % all others
    gpslon1(ii(jj)) = nanmedian1(gpslon1(1:ii(jj)+20)); % median of 1st 20 pnts only
    gpslon1(ii(ll)) = nanmedian1(gpslon1(ii(ll)-20:ii(ll)+20));
    gpslon1(ii(kk)) = nanmedian1(gpslon1(ii(kk)-20:end)); % median of last 20 pnts only
end;

ii=find(isnan(gpslat1));    % same for lat
if ~isempty(ii)
    jj = (ii<21);
    kk = (ii>3600-21);
    ll = ~xor(jj,kk);
    gpslat1(ii(jj)) = nanmedian1(gpslat1(1:ii(jj)+20));
    gpslat1(ii(ll)) = nanmedian1(gpslat1(ii(ll)-20:ii(ll)+20));
    gpslat1(ii(kk)) = nanmedian1(gpslat1(ii(kk)-20:end));
end;

ii=find(isnan(gpsspeed1));    % same for speed
if ~isempty(ii)
    jj = (ii<21);
    kk = (ii>3600-21);
    ll = ~xor(jj,kk);
    gpsspeed1(ii(jj)) = nanmedian1(gpsspeed1(1:ii(jj)+20));
    gpsspeed1(ii(ll)) = nanmedian1(gpsspeed1(ii(ll)-20:ii(ll)+20));
    gpsspeed1(ii(kk)) = nanmedian1(gpsspeed1(ii(kk)-20:end));
end;

ii=find(isnan(gpshead1));    % same for heading
if ~isempty(ii)
    jj = (ii<21);
    kk = (ii>3600-21);
    ll = ~xor(jj,kk);
    gpshead1(ii(jj)) = nanmedian1(gpshead1(1:ii(jj)+20));
    gpshead1(ii(ll)) = nanmedian1(gpshead1(ii(ll)-20:ii(ll)+20));
    gpshead1(ii(kk)) = nanmedian1(gpshead1(ii(kk)-20:end));
end;

gpslon = gpslon1; 
gpslat = gpslat1; 
gpsspeed = gpsspeed1;
gpshead = gpshead1.*180/pi;     % back to degrees
gpshead = mod(gpshead+360,360); % be sure range is 0-360

sog = gpsspeed;
cog = gpshead;

%% compute ship N & E velocity components from sog & cog
% in m/s
Ngps = sog.*cos(cog./180.*pi); % 1x36000 arrays
Egps = sog.*sin(cog./180.*pi);

%------filter ship velocities
fs = 1;          % GPS output @ 1 Hz
wp = 1/25/(fs/2);
ws = 0.8*wp;
% [n,wn]=buttord(wp,ws,10,25);
[n,wn] = buttord(wp,ws,3,7);
[bhi,ahi] = butter(n,wn,'low');
Ngps = filtfilt(bhi,ahi,Ngps);
Egps = filtfilt(bhi,ahi,Egps);

%% output
gps = [jd_ref' cog' sog' Ngps' Egps' gpslat' gpslon'];

