% read_gps_daily
% NTAS 2009 :: 2009-06-12 :: LB
% Read the GPS ASCII data files.

%  Note:  hhh=hour

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

% gps data is 10 Hz
nsampl=86400*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',jam);
    dfl=[way_raw_data_flux, yyyy(3:4),jd,'\gps0' yyyy(3:4) jd hr '_raw.txt'];
    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,:)))=[];
        stx(4:end,any(isnan(stx),1))=NaN;      %put NaN into any rows containing NaNs (to remove bad lines, gap...), except JD


        [nr1,nl1]=size(stx);
        st2(:,jax1:jax1+nl1-1)=stx;
        jd1=ddd + (jam+(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;
        
    fclose(flist); % 
    else
        disp(sprintf('\b<--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);

% time_gpspsd=st2(1,:);
% jd_gps=jd_gps(ii);
glt=single(st2(4,:)); 
gpslat_sign=single(st2(5,:)); 
gln=single(st2(6,:)); 
gpslon_sign=single(st2(7,:)); 
gpsspeed1=single(st2(8,:).*.514);  %convert to m/s
gpshead1=single(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 plotit &&  ~isempty(st2(~isnan(st2)))
    figure;
    plot(gpslon1,gpslat1);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD GPS Lat/Lon',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none');
    xlabel('Longitude (deg)');ylabel('Latitude (deg)');grid
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
    if prtit; print('-djpeg90 ',[way_images_flux cruise '_PSD_GPS_track_LAT_LON_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
    
    figure;
    subplot(2,1,1)
    if exist('gpstim','var') && exist('sog','var')
        plot(jd_gps,gpsspeed1,gpstim./24+ddd,sog);legend('PSD','SCS')
    else
        plot(jd_gps,gpsspeed1);legend('PSD')
    end    
    xlabel('Hour (UTC)');ylabel('SOG (m/s)');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD GPS SOG/COG',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlim([ddd ddd+1]); ylim([0 8]);set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))
    grid
    subplot(2,1,2)
    if exist('gpstim','var') && exist('cog','var')
         plot(jd_gps,gpshead1,'.',gpstim./24+ddd,cog);legend('PSD','SCS')
    else
         plot(jd_gps,gpshead1,'.');legend('PSD')
    end    
    xlabel('Hour (UTC)');ylabel('COG (deg)'); axis([ddd ddd+1 0 360]); 
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    orient tall
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
    if exist('gpstim','var') && exist('cog','var')
        if prtit; print('-djpeg90 ',[way_images_flux cruise '_PSD-SCS_GPS_track_COG_SOG_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
    else
        if prtit; print('-djpeg90 ',[way_images_flux cruise '_PSD_GPS_track_COG_SOG_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
    end 
    
    figure;
    eval(['map_', currentMname,'(gpslon1,gpslat1,cruise,ship,[Lonmin Lonmax Latmin Latmax],Fluxroot)'])
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  R/V %s track',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,ship),'FontWeight','Bold','Interpreter','none')
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
    if prtit; print(graphdevice,[way_images_flux cruise '_shiptrack_map_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]);  end

end


if ~isempty(gpslon1) && ~isempty(gpsspeedx)
    jd_ref=ddd+(0:1439)/60/24;
%     jd_ref=jd_pc;
    gpslon1i=single(intavg(jd_gps,gpslon1,jd_ref)); 
    gpslat1i=single(intavg(jd_gps,gpslat1,jd_ref)); 
    gpsspeedxi=single(intavg(jd_gps,gpsspeedx,jd_ref)); 
    gpsspeedyi=single(intavg(jd_gps,gpsspeedy,jd_ref)); 

    gpsspeedi=single(sqrt(gpsspeedxi.*gpsspeedxi+gpsspeedyi.*gpsspeedyi));                         %1-min sog
    gpsheadi=single(mod(atan2(gpsspeedyi,gpsspeedxi)/rdcon+360,360));
else
    gpslat1i=ones(length(qvais),1)*NaN;
    gpslon1i=ones(length(qvais),1)*NaN;
end;

% if plotit
%     figure;
%     subplot(2,1,1)
%     plot(gpslon1i,gpslat1i);title(['DOY',jd,' UTC  Longitude/Latitude from PSD GPS. ', cruise, yyyy],'FontWeight','Bold');xlabel('Longitude (deg)');ylabel('Latitude (deg)');
%     subplot(4,1,3)
%     plot(jd_pc,gpsspeedi);xlabel(['Yearday ', yyyy]);title(['DOY',jd,' UTC  Speed-Over-Ground and Course-Over-Ground from PSD GPS. ' , cruise, yyyy],'FontWeight','Bold');ylabel('SOG (m/s)');
%     subplot(4,1,4)
%     plot(jd_pc,gpsheadi);xlabel(['Yearday ', yyyy]);ylabel('COG (deg)');
%     orient tall
%     annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
%     if prtit; print('-djpeg90 ',[way_images_flux cruise '_' yyyy '_PSD_GPS_track_COG_SOG_' num2str(ddd) graphformat]); end
% end
% 
% if saveit;
%   np=length(jd_pc);  
%   eval([ 'prt_jas_gps_', currentMname])
% %   [r,c] = size(st2');
% end;

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