% Xcruise 2009 :: 2009-09-07 :: LB
% Read sonic anemometer data stream.

disp(['read_son_day: ', cruise]);
fclose('all');

clear st1;clear st2;clear st3; clear jd_son;

jd=sprintf('%03i',ddd);
jhold=0;
jax=1;

%declare arrays
st1=[];st2=[];
st3=ones(8,1440)*NaN;
jd_son=ones(1,1440)*NaN;
jaxd=[];jade=[];temp=[];

for jam=0:23,						%cycle thru 24 hourly stats files
    hr=sprintf('%02i',jam);
    e=[way_raw_data_flux, yyyy(3:4),jd,'\son0' yyyy(3:4) jd hr '_raw.txt'];
    flist=fopen(e,'r');
    st1=[];temp=[];
    if flist>0,	   				%if the file exists,
        disp(['Reading Sonic file  ',e])
        k=0; %SPdeS keep count
        while ~feof(flist) % end-of-file not reached
            k=k+1;
            try
                temp=textscan(flist,'%2f%2f%3f %*4c %f %f %f %f %f %*3c %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Sonic');
                st1=[st1 cell2mat(temp)'];
            catch
                % if one line is cut short, blank it out totally
                [l8,c8]=size(temp{1,8});
                if l8==0 || c8==0; continue; end % read next batch of lines
                for ll=1:8
                    if length(temp{1,ll})~=l8
                        temp{1,ll}(l8+1)=[];
                    end;
                end;
                if isempty(temp{1,1});
                    temp{1,9}=temp{1,9}';
                end;
                st1=[st1 cell2mat(temp)'];
            end; %end try/catch
        end;
        
        [nr1,nl1]=size(st1);     
        st2(:,jax:jax+nl1-1)=st1(:,1:nl1);

        ii=find(~isnan(st1(1,:)));
        uu=st1(5,ii);
        vv=st1(6,ii);
        ww=st1(7,ii);
        jade=ddd+(jam+(st1(1,:)+(st1(2,:)+st1(3,:)/1000)/60)/60)/24; % time in decimal yearday

        if jam==0 || jam==12 || jam==18,
            disp(['Calculating sonic spectra for ',hr])
            SS=sqrt(uu.^2+vv.^2); %SS=sqrt(UU.^2+VV.^2);
            [Su,F]=psd2(detrend(SS),length(SS),fsonic);[Sus,Fsu]=specsmoo(Su,fsonic);
            [Sw,F]=psd2(detrend(ww),length(ww),fsonic);[Sws,Fsw]=specsmoo(Sw,fsonic);

            if plotit;
                figure;
                loglog(Fsu(length(Fsu)-31:length(Fsu)),((Fsu(length(Fsu)-31:length(Fsu))).^(-2/3))./(Fsu(length(Fsu)-18).^(-2/3)).*Sus(length(Fsu)-18).*Fsu(length(Fsu)-18),'r');
                hold on;
                loglog(Fsu,Sus.*Fsu,'b',Fsw,Sws.*Fsw,'g');
                title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Sonic anemometer spectra: hour %02i',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,jam),'FontWeight','Bold','Interpreter','none')
                legend('5/3','Spd','W','Location','Best');       xlabel('Frequency (Hz)');
                ylabel('Spectrum * Frequency (m^2 s^-^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(graphdevice,[way_images_flux cruise '_PSD_Sonic_UWspectra_hour' hr '_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
            end
            
            if jam==12,  % Print out spectra file for hr 12
                specfile=[way_spec_data_flux, cruise '_PSD_Sonic_w_sp_hour' hr '_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt'];
                fspec=fopen(specfile,'w');
                specs=[Fsu;Sus;Fsw;Sws];
                fprintf(fspec,'%11.5f %11.5f %11.5f %11.5f \r\n',specs);
                fclose(fspec);
                clear fspec specs specfile
            end;%end hhh=12
            clear SS ww Su F Sus Fsu Sw Sws Fsw
        end;%end if hour==0 or 12 or 18

        jax=jax+nl1;
                
        sonfac=3600*fsonic/60;
        if ~isempty(st2);
            jk=min(60,floor(jax/sonfac));
            for j=1:jk,
                jl=sonfac*(j-1)+1;
                jj=min(jax-1,sonfac*j);
                st3(:,j+jhold)=median(st2(:,jl:jj)')';
                jd_son(j+jhold)=jade(jl);
            end;
            stx=st2;
            jhold=jhold+jk;
            st2=[]; st1=[];
            jaxd=[jaxd jax];
            jax=1;
        else
            jhold=jhold+60;
        end;       %if st2 isempty
    fclose(flist); % 
    else
        jhold=jhold+60;
    end;         %end if flist
end;            %end for jam=0:23 % cycle through hour

clear u v w Tsonic U V dir
u=st3(5,:);  %m/s
v=st3(6,:);  %m/s
w=st3(7,:);  %m/s

% % % Time in flux data files is the number of 100-nano second intervals since
% % % Jan 1, 1601 divided by 10000. It was chosen to maintain compatibility
% % % with an HP-UNIX version of the code form long ago
% % %  0             'dd-mmm-yyyy HH:MM:SS'  
% % %01-Mar-2000 15:45:17
% % dates = datestr(datenum(1601,1,1) + datenum(jd_son(1,:)*10000*100e-9/86400));
% % hh = str2num(dates(:,13:14));
% % mm = str2num(dates(:,16:17));
% % ss = str2num(dates(:,19:20));
% % clear jd_son
% % jd_son= ddd + (hh+((mm+(ss/60))/60))/24;
% % % test to see if first point is from the previous day
% % if str2num(dates(1,1:2)) ~= str2num(dates(2,1:2))
% %     jd_son(1) = jd_son(1)-1;
% % end    

if outputsonic==1
    Tsonic=st3(8,:); %Temperature in degC
else
    speedson=st3(8,:)/50;
    Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15;
end;

if rotationsonic
    switch sonicmodel
        case 'R3'                                       %means sonicmodel==R3
            U=u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
            V=u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
        case 'R2'    
            U=-u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
            V=-u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
    end  
else
    U=u;
    V=v;
end;
S=sqrt(u.*u+v.*v+w.*w);
dir=atan2(-V,U)*180/pi+180;		%jeff coordinates, +boward, +portward, 180 means wind from
dir=mod(dir+720,360);

np=length(st3);
if saveit;
    if length(isnan(jd_son)==1)==1440;   %means no data files detected during the 24H period
        jd_son=ddd+(0:1439)/60/24;       %creates a fake decimal JD base so that timatch programs work properly
    end;
    eval([ 'prt_jas_son_', currentMname])
end;

if plotit &&  ~isempty(st3(~isnan(st3)))
    %plot T profiles
    j_st=jd_son(1); 
    j_en=jd_son(np);
    figure;
    plot(jd_son(1:np),Tsonic(1:np));xlabel('Hour (UTC)');ylabel('Temperature (degC)');xlim([ddd ddd+1])
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Sonic air temperature',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    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(graphdevice,[way_images_flux cruise '_PSD_Sonic_temperature_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;
    plot(jd_son(1:np),U(1:np),jd_son(1:np),V(1:np),jd_son(1:np),w(1:np),jd_son(1:np),S(1:np));
    xlabel('Hour (UTC)');ylabel('Wind speed (m/s)');xlim([ddd ddd+1])
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Sonic relative wind components',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    legend('U','V','W','speed','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(graphdevice,[way_images_flux  cruise '_PSD_Sonic_winds_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
%     for i=1:np % Jeff Otten's mod to display direction -180 to +180.
%         if dir(i) > 180
%             dir(i) = dir(i) - 360;
%         end
%     end
    dir(dir>180)=dir(dir>180)-360;

    figure;subplot(2,1,1)
    plot(jd_son(1:np),dir(1:np),'b.', ... 		%data points
        [ddd ddd+1],[45 45],'g', ...			%45 degree good zone
        [ddd ddd+1],[-45 -45],'g', ...
        [ddd ddd+1],[60 60],'r', ...			%60 degree questionable zone
        [ddd ddd+1],[-60 -60],'r');
        xlabel('Hour (UTC)');
        ylabel('Wind direction (deg)');xlim([ddd ddd+1])
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Sonic relative wind direction',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        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 ~isnan(j_st) && ~isnan(j_en)
            axis([j_st j_en -180 180]);
        end;
       
%         if prtit; print(graphdevice,[way_images_flux  cruise '_' yyyy '_PSD_Sonic_rel_wnd_dir_' num2str(ddd) graphformat]); end;
        
        subplot(2,1,2)
        plot(jd_son(1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));xlabel('Hour (UTC)');ylabel('Wind speed (m/s)');xlim([ddd ddd+1])
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Sonic relative wind speed',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
        orient tall
        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(graphdevice,[way_images_flux  cruise '_PSD_Sonic_rel_wnd_speed_dir_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
end; % if plotit


if clearit 
    clear stx uu vv ww
end