disp(['read_son_day_', cruise, '_', year]);
fclose('all');

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

read_parameters_rhb

ddd=input('Input start yearday to evaluate.... ');
prtit=input('Print out text files? (Y=1 N=0) ');
plotit=input('Plot data? (Y=1 N=0) ');

jd=num2str(ddd);

if ddd<100
    jd=['0' num2str(ddd)];
end;
if ddd<10
    jd=['00' num2str(ddd)];
end;

jhold=0;
jax=1;

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

for jam=0:1,						%cycle thru 24 hourly stats files

    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;

    %e=[way_raw_data_flux, 'day',jd,'\P0_' jd hr];
    e=[way_raw_data_flux, 'DAY', jd,'\son0' year(3:4) jd hr '_raw.txt'];
    flist=fopen(e,'r');
    %10/06/07 Generic with ASCII Output (SN: 1234)
    %160000181 04,00,+01.51,-00.29,-00.80,+33.44,02
    %160000275 05,00,+01.09,-00.22,-00.60,+33.45,0A
    %160000384 06,01,+00.80,+00.03,-00.88,+33.44,0A

    if flist>0,	   				%if the file exists,
        disp(['Reading Sonic file  ',e])
        %             st1=fscanf(flist,'%g,%g,%g,%g,%g',[5,inf]);
        st1 =textscan(flist,'%2f%2f%2f%3f %*3c %f %f %f %f %f %*3c %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Sonic');
        st1=cell2mat(st1)';
        st1(5:9,any(isnan(st1),1))=NaN;      %put NaN into any rows containing NaNs (to remove bad lines, gap...), except JD

        [nr1,nl1]=size(st1);
        st2(:,jax:jax+nl1-1)=st1(:,1:nl1);
        uu=st1(6,:);  vv=st1(7,:);   ww=st1(8,:);
        jade=ddd + (st1(1,:)+(st1(2,:)+(st1(3,:)+st1(4,:)/1000)/60)/60)/24;
        jax=jax+nl1;

        % despiking used during GOMECC07...
        sonraw = st1;   % load sonic array into temp array
        sonfix  % run despiking routine
        clear st1
        st1 = sonraw_dep; %copy despiked data into working array

        [nr1,nl1]=size(st1);
        st2(:,jax:jax+nl1-1)=st1(:,1:nl1);
        uu=st1(2,:)/100;  vv=st1(3,:)/100;   ww=st1(4,:)/100;
        %                 jade=ddd + (jam+(st1(1,:)+(st1(2,:)+st1(3,:)/1000)/60)/60)/24;
        jax=jax+nl1;

        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);

            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, ' ', year]);
            legend('5/3','Spd','W' );       xlabel('Frequency in Hz');
            ylabel('Spectrum * Frequency');
            print('-djpeg90 ',[way_images_flux 'son_spectrum' jd hr '.jpg']);
            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 hhh=0 or 12 or 18

        sonfac=3600*fsonic/60;
        if ~isempty(st2);
            jk=min(60,ceil(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)=st2(1,jl);
                jd_son(j+jhold)=jade(jl);
            end;    %end for
            stx=st2;
            jhold=jhold+jk;
            st2=[]; st1=[];
            jaxd=[jaxd jax];
            jax=1;
        end;       %if isempty
    end;         %end if flist
end;            %end for
clear u v w Tsonic U V dir
u=st3(6,:);  %m/s
v=st3(7,:);  %m/s
w=st3(8,:);  %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(9,:)/100; %Temperature in degC
else
    speedson=st3(9,:)/50;
    Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15;
end;

if rotationsonic==1
    C=strcmp(sonicmodel,'R3');
    if C==1                                            %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
    else
        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);

if prtit==1;
    np=length(st3);
    prt_jas_son_Stratus_07;
end;

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