% read_gps_hourly
% Read the GPS ASCII data files.

disp(['read_gps: ',cruise]);

% gps data is 10 Hz
nsampl=3600*10;
ncols=11;
jd_gps=NaN+zeros(1,nsampl*1.01); % preallocate with tiny 1% buffer
st2=NaN+zeros(ncols,nsampl*1.01);

jax1=1;jd1=[];
jd=sprintf('%03i',ddd);


% for jam=00:23
    hr=sprintf('%02i',hhh);
    %dfl=[way_raw_data_flux, year(3:4),jd,'\gps0' year(3:4) jd hr '_raw.txt'];
    dfl=[way_raw_data_flux, year(3:4),jd,'/gps0' year(3:4) jd hr '_raw.txt']; %Mac file system
    disp(['Reading gps file  ',dfl])
    flist=fopen(dfl);
%     stx=[];temp=[];
   
    if flist>0
        headdata=importdata(dfl);
        temp=headdata(cellfun(@(x) ~isempty(strfind(x,'$GPRMC')),headdata));
    
        stx=ones(ncols,36000)*NaN;
        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', ', ',  '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,:)))=[];

        [nr1,nl1]=size(stx);
        st2(:,jax1:jax1+nl1-1)=stx;
        jd1=ddd + (hhh+(st2(1,jax1:jax1+nl1-1)+(st2(2,jax1:jax1+nl1-1)+st2(3,jax1:jax1+nl1-1)/1000)/60)/60)/24;
        jd_gps(:,jax1:jax1+nl1-1) = jd1;
        jax1=jax1+nl1;
    
    else
        disp('file not found')
    end;
% end;

% trim arrays to their filled size
st2=st2(:,1:jax1-1);
jd_gps=jd_gps(:,1:jax1-1);

[b,m,n]=unique(jd_gps);
jd_gps=jd_gps(m);st2=st2(:,m);

glt=st2(4,:); 
gpslat_sign=st2(5,:); 
gln=st2(6,:); 
gpslon_sign=st2(7,:); 
% gpsspeed1=st2(8,:).*.514;  %convert to m/s
gpsspeed1=st2(8,:);  %in knots
gpshead1=st2(9,:);  %deg

% jd_gps=ddd+((st2(1,:)-st2(1,1))/1000/60/60)/24;

bb=floor(glt/100);
glt=(glt-bb*100)/60;
glt=glt+bb;
ii=find(gpslat_sign==83);glt(ii)=-glt(ii);
% glt(char(gpslat_sign)=='S')=-glt(char(gpslat_sign)=='S');
gpslat1=glt;

bbn=floor(gln/100);
gln=(gln-bbn*100)/60;
gln=gln+bbn;
ii=find(gpslon_sign==87);gln(ii)=-gln(ii);
% gln(char(gpslon_sign)=='W')=-gln(char(gpslon_sign)=='W');
gpslon1=gln;

ii=find(abs(diff(gpslon1))>1);gpslon1(ii)=gpslon1(ii-1);
ii=find(abs(diff(gpslat1))>1);gpslat1(ii)=gpslat1(ii-1);
ii=find(gpsspeed1>100);gpsspeed1(ii)=NaN;
ii=find(gpshead1>400);gpshead1(ii)=NaN;

rdcon=pi/180;
gpsspeedx=gpsspeed1'.*cos(gpshead1'*rdcon);    %convert to N
gpsspeedy=gpsspeed1'.*sin(gpshead1'*rdcon);    %convert to E

if length(st2)< 36000
        disp('<<<<<<<<<<<<<<<<<<  achtung >>>>>>>>>>>>>>>>>>>>>>>>')
    if ~isempty(gpslon1) && ~isempty(gpsspeedx)
        jd_ref=ddd+hhh/24+(0:1/10:35999/10)/60/60/24;
    %     tref=0:1/10:35999/10;
        gpslon1=interp1(jd_gps,gpslon1,jd_ref); ii=find(isnan(gpslon1)); 
        try 
            gpslon1(ii)=gpslon1(ii(end)+1); 
        catch
            disp('')
        end
        gpslat1=interp1(jd_gps,gpslat1,jd_ref); ii=find(isnan(gpslat1));
        try 
            gpslat1(ii)=gpslat1(ii(end)+1); 
        catch
            disp ('');
        end
        gpsspeedxi=interp1(jd_gps,gpsspeedx,jd_ref); ii=find(~isnan(gpsspeedxi));msu=mean(gpsspeedxi(ii)); 
        try 
            ii=find(isnan(gpsspeedxi));gpsspeedxi(ii)=msu; 
        catch
            disp (''); 
        end
        gpsspeedyi=interp1(jd_gps,gpsspeedy,jd_ref); ii=find(~isnan(gpsspeedyi));msu=mean(gpsspeedyi(ii)); 
        try 
            ii=find(isnan(gpsspeedyi));gpsspeedyi(ii)=msu; 
        catch
            disp (''); 
        end
        
        
        gpsspeed1=sqrt(gpsspeedxi.*gpsspeedxi+gpsspeedyi.*gpsspeedyi);                         %1-min sog
        gpshead1=mod(atan2(gpsspeedyi,gpsspeedxi)/rdcon+360,360);

        gpslat1(1)=gpslat1(2);
        gpslon1(1)=gpslon1(2);
    end;   
% elseif length(gpsspeed1)< 36000
% % if length(gpsspeed1)< 36000
%     gpsspeed1=[gpsspeed1 (1:36000-length(st2)).*0+gpsspeed1(length(st2))];
%     gpshead1=[gpshead1 (1:36000-length(st2)).*0+gpshead1(length(st2))];
%     gpslat1=[gpslat1 (1:36000-length(st2)).*0+gpslat1(length(st2))];
%     gpslon1=[gpslon1 (1:36000-length(st2)).*0+gpslon1(length(st2))];
elseif length(gpsspeed1)> 36000
% if length(gpsspeed1)> 36000
    gpsspeed1=gpsspeed1(1:36000);
    gpshead1=gpshead1(1:36000);
    gpslat1=gpslat1(1:36000);
    gpslon1=gpslon1(1:36000);
end


gpslon=gpslon1; 
gpslat=gpslat1; 
gpsspeed=gpsspeed1;                         %1-min sog
gpshead=gpshead1;
    
sog=gpsspeed1;		%in knots as used in programs afterward...!!!
cog=gpshead1;

if clearit
    clear gpslon1 gpslat1 gpsspeed1 gpshead1 glt gln gpslat_sign gpslon_sign 
    clear gpsspeedx gpsspeedy jd_gps bb bbn headdata jd1 st2 stx temp
end
