% NTAS 2009 :: 2009-06-12 :: LB
% Read Crescent GPS data (differential GPS)
%  Note:  hhh=hour

disp(['read_head_pitch_day_', cruise]);

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

jax1=1;jd1=[];

for jam=0:23
    
    hr=sprintf('%02i',jam);   
    jd=sprintf('%03i',ddd);  
    
    nmang='pitch';
    dfl=[way_raw_data_flux, yyyy(3:4),jd,'\hed0' yyyy(3:4) jd hr '_raw.txt'];
    disp(['Reading head file  ',dfl])
    flist=fopen(dfl);
    if flist>0
        
    %     temp =textscan(flist,'%f %6c %6c %c %f %c %f %c %f %f %f %f %4c %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN);
        temp =textscan(flist,'%2f%2f%3f %5c %3c %f %f %f %f %4c %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN);
         %%%%  check for pitch or roll
         xy=temp{1,8};
        if ~isnan(xy(1));
            try
                stx=[temp{1,1}';temp{1,2}';temp{1,3}';(temp{1,6}');(temp{1,7}');(temp{1,8}')];
            catch
                [l10,c10]=size(temp{1,10});
                if l10==0 || c10==0; continue; end % read next batch of lines
                for ll=1:10
                    if length(temp{1,ll})~=l10
                        temp{1,ll}(l10+1)=[];
                    end;
                end;
                if isempty(temp{1,1});
                    temp{1,10}=temp{1,10}';
                end;
                stx=[temp{1,1}';temp{1,2}';temp{1,3}';(temp{1,6}');(temp{1,7}');(temp{1,8}')];
            end;    
        else
            try
                stx=[temp{1,1}';temp{1,2}';temp{1,3}';(temp{1,6}');(temp{1,7}');(temp{1,9}')];
            catch
                [l10,c10]=size(temp{1,10});
                if l10==0 || c10==0; continue; end % read next batch of lines
                for ll=1:10
                    if length(temp{1,ll})~=l10
                        temp{1,ll}(l10+1)=[];
                    end;
                end;
                if isempty(temp{1,1});
                    temp{1,10}=temp{1,10}';
                end;
                stx=[temp{1,1}';temp{1,2}';temp{1,3}';(temp{1,6}');(temp{1,7}');(temp{1,9}')];
            end;    
            nmang='Pitch';
        end;
        stx(4:end,any(isnan(stx),1))=NaN;      %put NaN into any rows containing NaNs (to remove bad lines, gap...), except JD
% stx(:,isnan(stx(1,:)))=[];

        [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_hed(:,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_hed_pitch=jd_hed(:,1:jax1-1);

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

time_hed_pitch=single(st2(1,:));
headx_pitch=single(st2(5,:)); 
pitchx_pitch=single(st2(6,:)); 

tPSD=(jd_hed_pitch-316)*24*3600;
GPSh=floor(st2(4,:)/1e4);
GPSm=floor((st2(4,:)/1e4 - GPSh)*100);
GPSs=([(st2(4,:)/1e4 - GPSh)*100]-GPSm)*100;
tPSDGPS=GPSh*3600+GPSm*60+GPSs;

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

if plotit &&  ~isempty(st2(~isnan(st2)))
figure;
    if exist('gpstim','var') && exist('lrg','var')
        plot(jd_hed_pitch,headx_pitch,'.',gpstim./24+ddd,lrg);legend('PSD','SCS')
    else
        plot(jd_hed_pitch,headx_pitch,'.');legend('PSD')
    end    
    
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD GPS heading',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Heading (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
    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_heading_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;

    ii=find(isnan(pitchx_pitch));pitchx_pitch(ii)=nanmean(pitchx_pitch);
    %lowpass filter time series at 2 Hz for signal sampled at 10 Hz

    wp=[.01 3]/5;ws= [.003 4]/5;rp=2;rs=8;
    wp=[.02/5 2/5];ws= [.01/5 3/5];rp=3;rs=15;
    [n,wmn]=cheb1ord(wp,ws,rp,rs);
    [b,a]=cheby1(n,rp,wmn);
    yxp_pitch=filter(b,a,pitchx_pitch-mean(pitchx_pitch));%deg
    figure;plot(jd_hed_pitch,pitchx_pitch,'b-',jd_hed_pitch,yxp_pitch,'g');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD GPS %s angle ',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,nmang),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel([nmang ' (deg)']);  xlim([ddd ddd+1])
    legend('Raw signal','Pass-band filtered (0.02 to 2Hz)','Location','Best')
    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
    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_' sprintf('%s_%04i_%02i_%02i_%03i',nmang,Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
end;

if ~isempty(headx_pitch) && ~isempty(pitchx_pitch)
    rdcon=pi/180;
    headx_pitchc=cos(headx_pitch'*rdcon);        %convert to N
    headx_pitchs=sin(headx_pitch'*rdcon);        %convert to E
%     jd_ref=ddd+(0:1439)/60/24;
    jd_ref=jd_pc;
    %headxi_pitchc=interp1(jd_hed_pitch,headx_pitchc,jd_ref); %deg
    headxi_pitchc=intavg(jd_hed_pitch,headx_pitchc,jd_ref); %deg
    headxi_pitchs=intavg(jd_hed_pitch,headx_pitchs,jd_ref); %deg
    pitchxi_pitch=intavg(jd_hed_pitch,pitchx_pitch,jd_ref); %deg
    
    headxi_pitch=mod(atan2(headxi_pitchs,headxi_pitchc)/rdcon+360,360);   

    np=length(headxi_pitch); 
    
    if saveit
        eval([ 'prt_jas_gps_', currentMname])
    end;
end;

% figure;plot(jd_ref,headxi,'.');title(['PSD GPS Heading. ', cruise, yyyy],'FontWeight','Bold');xlabel(['Yearday ', yyyy]);ylabel('Heading (deg)');  
% 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_heading_' num2str(ddd) graphformat]); end;


% prt_jas_gps_Xcruise_2009
% if saveit && exist('gpslat1i','var') && exist('headxi','var') && ~isempty(gpslat1i) && ~isempty(headxi)
%     eval([ 'prt_jas_gps_', currentMname])
% end;
% if saveit && (~isempty(gpslon1) && ~isempty(gpsspeedx)) && (~isempty(headx) && ~isempty(pitchx))
%     eval([ 'prt_jas_gps_', currentMname])
% end;

if clearit
    clear headx pitchx jd_hed jdx time_hed st2 yxp jaz jazz jd1 stx xy yxp_pitch time_hed_pitch pitchx_pitch jd_hed_pitch
    clear headx_pitchc headx_pitchs headx_pitch
%clear also big size GPS variables
%     clear gpslon1 gpslat1 gpsspeed1 gpshead1 glt gln gpslat_sign gpslon_sign 
%     clear gpsspeedx gpsspeedy jd_gps
end
