% read_son_day
% VOCALS 2008 :: 2008-10-16 :: Simon de Szoeke
% Read sonic anemometer data stream.

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

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

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

st1=[];st2=[];st3=[];
jaxd=[];jade=[];temp=[];

for jam=0:23,						%cycle thru 24 hourly stats files
    hr=sprintf('%02i',jam);
    e=[way_raw_data_flux, 'day',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(['Sonic Anemometer Spectra:   Day ',jd,'   Hour ',hr, ', ',cruise, '  ' yyyy]);
                legend('5/3','Spd','W' );       xlabel('Frequency in Hz');
                ylabel('Spectrum * Frequency');
                if prtit; print('-djpeg90 ',[way_images_flux 'son_spectrum' jd hr '.jpg']); end;
            end
            
            if jam==12,  % Print out spectra file for hr 12
                specfile=[way_spec_data_flux, 'w_sp',jd,num2str(hr),'.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;
        end;       %if st2 isempty
    fclose(flist); % SPdeS
    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 'VOCALS'
            U= u*cos(90+30/180*pi)-v*sin(90+30/180*pi);
            V= u*sin(90+30/180*pi)+v*cos(90+30/180*pi);
        otherwise
            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;
    prt_jas_son;
end;

if plotit;
    %plot T profiles
    j_st=jd_son(1); 
    j_en=jd_son(np);
    figure;
    plot(jd_son(1:np),Tsonic(1:np));xlabel(['JD', yyyy]);ylabel('Temperature (C)');title(['Sonic Temperature. ', cruise, yyyy]);
    if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_temp_' num2str(ddd) '.jpg']); 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(['JD ' yyyy]);ylabel('Sonic Winds (m/s)');title(['Relative Wind Components:  U=blue,  V=green,  W=red,   Speed=cyan. ', cruise, yyyy]);
    if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_winds_' num2str(ddd) '.jpg']); 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;
    plot(jd_son(1:np),dir(1:np),'b', ... 		%data points
        [ddd+4/24 ddd+4/24],[-180 180],'k', ...		%4 hour watch dividers
        [ddd+8/24 ddd+8/24],[-180 180],'k', ...
        [ddd+12/24 ddd+12/24],[-180 180],'k', ...
        [ddd+16/24 ddd+16/24],[-180 180],'k', ...
        [ddd+20/24 ddd+20/24],[-180 180],'k', ...
        [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(['JD ' yyyy]);
        ylabel('Sonic Wind Direction (deg)');
        title(['Relative Wind Direction. ', cruise, yyyy]);
        axis([j_st j_en -180 180]);
        if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_rel_wnd_dir_' num2str(ddd) '.jpg']); end;
        
        figure;
        plot(jd_son(1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));xlabel(['JD ' yyyy]);ylabel('Wind Speed m/s');title(['Sonic Relative Wind Speed. ', cruise, yyyy]);
        if prtit; print('-djpeg90 ',[way_images_flux '\Sonic_rel_wnd_speed_' num2str(ddd) '.jpg']); end;
end; % if plotit