disp('rdson2_08')    %('eval_son_rhb_07')
fclose('all');
clear st1;clear st2;clear st3;clear jadson;

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

hr='00';
jhold=0;
jax=1;
jaxx=1;
end;

plotit=1;
backchk=0;
st1=[];
st2=[];
st3=[];
for jam=hhh:hhh,						%cycle thru 24 hourly stats files
    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;
            e=[way_raw_data, 'day', jd,'\son008' jd hr '_raw.txt'];
    flist=fopen(e,'r');
    if flist>0,	   				%if the file exists,
        disp(['Reading Sonic file in rdson_08 ',e])
        dmp=fgets(flist);			%read header lines and skip blanks
        dmp2=fgets(flist);		%read header lines and skip blanks
        while feof(flist)==0,
            st1=fscanf(flist,'%g,%g,%g,%g,%g',[5,inf]);
            % despike data  DEW 5/6/2007 fm CWK
            ss = size(st1)
            %put in because file P0_12812 was bombing after it look like it
            %read the file once after 64004 pts and then tried to read
            %again but was empty 71857 pts in file DEW 5/9/2007
            if ss(1) == 0
                break
            end
            sonraw = st1;   % load sonic array into temp array
            sonfix  % run despiking routine
            clear st1
            st1 = sonraw_dep; %copy despiked data into working array

            if backchk==0,
                tzilch=jam;
                backchk=1;
                tzero=st1(1);
            end;% end if
            [nr1,nl1]=size(st1);
            nll=nl1-1;
            if nr1==5
                st2(:,jax:jax+nl1-1)=st1(:,1:nl1);
                jax=jax+nl1;
            else
                x=fgets(flist);
            end;%end if
        end;%end while feof
        %      sonfac=1200 or 600; %sonic runs at 10Hz, 36000 values/hour...
        sonfac=3600*fsonic/60;
        if ~isempty(st2);
            jk=min(60,ceil(jax/sonfac));
            jade=str2num(jd)+((st2(1,:)-tzero)/1000/60/60+tzilch)/24;
            for j=1:jk
                jl=sonfac*(j-1)+1;jj=min(jax-1,sonfac*j);
                st3(:,j+jhold)=median(st2(:,jl:jj)')';
                jadson(j+jhold)=jade(jl);
            end;%end for
            stx=st2;
            jhold=jhold+jk;
            jax=1;
        end;  %if isempty
    end;    %end if flist
end;       %end for

if ~isempty(st3)
    clear u v w Tsonic U V dir
    u=st3(2,:)/100;
    v=st3(3,:)/100;
    w=st3(4,:)/100;

    if outputsonic==1
        Tsonic=st3(5,:)/100; %Temperature in degC
        Tslow=Tsonic;
    else
        speedson=st3(5,:)/50;
        Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15;
        Tslow=Tsonic;
    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);
    uu=stx(2,:)/100;
    vv=stx(3,:)/100;
    ww=stx(4,:)/100;
    tt=st3(5,:)/100;%-273.15;
    UU=uu;
    VV=vv;
    dir=atan2(-V,U)*180/pi+180;		%jeff coordinates, +boward, +portward, 180 means wind from
    dir=mod(dir+720,360);
    np=jhold;

    %clear w vul v uul u tzilch tzero sul st3 speedson prtit prt_it plotit nr1 npz np nl1 jl jk jj jhold jd jazx jaxx jax jam jadson jade j_st j_en j hul hr hll flist
    %clear f e dmp2 dmp dir backchk ans V U Tul Tsonic Tll S Rtul Rtll Rsul Rsll Rlul Rlll

    if 1   %hhh==0 | hhh==01 | hhh==12 || hhh==13 | hhh==15 | hhh==18 | hhh==22,
        disp(['Calculating spectra for ',num2str(hhh)])
        SS=sqrt(UU.^2+VV.^2);
        [Su,F]=psd2(detrend(SS),length(SS),20.);
        [Sus,Fsu]=specsmoo(Su,20.);
        [Sw,F]=psd2(detrend(ww),length(ww),20.);
        [Sws,Fsw]=specsmoo(Sw,20.);
        
        figure;
        plot(sonraw(2,:),'b')
        hold
        plot(sonraw_dep(2,:),'b.')
        plot(sonraw_dep(3,:),'g.')
        plot(sonraw_dep(4,:),'r.')
        plot(sonraw_dep(5,:),'k.')
        legend('raw','DS-U','DS-V','DS-W','DS-T')
        if 0
            if hhh == 1
                save F:\AMMA_2007\RHB\flux\Raw\spc01
            end
            if hhh == 13
                save F:\AMMA_2007\RHB\flux\Raw\spc13
            end
            if hhh == 15
                save F:\AMMA_2007\RHB\flux\Raw\spc15
            end
            if hhh == 22
                save F:\AMMA_2007\RHB\flux\Raw\spc22
            end
            if hhh == 0
                save F:\AMMA_2007\RHB\flux\Raw\spc00
                figure
                plot(sonraw(2,:),'b')
                hold
                plot(sonraw_dep(2,:),'b.')
                plot(sonraw_dep(3,:),'g.')
                plot(sonraw_dep(4,:),'r.')
                plot(sonraw_dep(5,:),'k.')
                legend('raw','DS-U','DS-V','DS-W','DS-T')
                title('Sonic Hr 00')
            end
            if hhh == 12
                save F:\AMMA_2007\RHB\flux\Raw\spc12
                figure
                plot(sonraw(2,:),'b')
                hold
                plot(sonraw_dep(2,:),'b.')
                plot(sonraw_dep(3,:),'g.')
                plot(sonraw_dep(4,:),'r.')
                plot(sonraw_dep(5,:),'k.')
                legend('raw','DS-U','DS-V','DS-W','DS-T')
                title('Sonic Hr 12')
            end
            if hhh == 18
                save F:\AMMA_2007\RHB\flux\Raw\spc18
                figure
                plot(sonraw(2,:),'b')
                hold
                plot(sonraw_dep(2,:),'b.')
                plot(sonraw_dep(3,:),'g.')
                plot(sonraw_dep(4,:),'r.')
                plot(sonraw_dep(5,:),'k.')
                legend('raw','DS-U','DS-V','DS-W','DS-T')
                title('Sonic Hr 18')
            end
        end
        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
        loglog(Fsu,Sus.*Fsu,'b');
        loglog(Fsw,Sws.*Fsw,'g');
        title(['Sonic Anemometer Spectra:   Day ',num2str(ddd),'   Hour ',num2str(hhh),]);
        legend('5/3','Spd','W' );
        xlabel('Frequency in Hz');
        ylabel('Spectrum * Frequency');
        print_buffer = ['F:\AMMA_2007\RHB\flux\Raw_Images\',jdr ,num2str(hhh), '_spc.jpg'];
        print('-djpeg90 ', print_buffer);
        if hhh==12,  % Print out spectra file for hr 12
            specfile=[way_spec_data, '\w_sp',jdr,num2str(hhh),'.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 UU VV Su F Sus Fsu Sw F Sws Fsw
    end;%end hhh=0 or 12 or 18
    clear ww Tsonic UU VV jade jadson tt uu vv ww SS
else
    U=[];V=[];w=[];dir=[];Tslow=[];S=[];
end;