disp('eval_son_sjohnson_05')
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)];
end;

plotit=1;
prtit=1;

hr='00';
jhold=0; 
jax=1;
jaxx=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=['d:\data\day',jd,'\P0_' jd hr];
   flist=fopen(e,'r');
   if flist>0,	   				%if the file exists, 
      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]);
         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
     disp(['Reading Sonic file from hour ',num2str(jam)])   
     sonfac=1200;
     if ~isempty(st2);
        jk=min(60,ceil(jax/sonfac));
        jade=str2num(num2str(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;
    Tsonic=st3(5,:)/100;
    %speedson=st3(5,:)/50;
    %Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15;
    Tslow=Tsonic;
    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
    S=sqrt(u.*u+v.*v+w.*w);
    uu=st1(2,:)/100;
    vv=st1(3,:)/100;
    ww=st1(4,:)/100;
    tt=st3(5,:)/100;
    UU=uu;
    VV=vv;
    dir=atan2(-V,U)*180/pi;		%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 hhh==0 | hhh==12,
        disp('Computing Spectra')
        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;loglog(Fsu,Sus.*Fsu,Fsw,Sws.*Fsw,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));
        title(['Sonic Anemometer Spectra:   Day ',num2str(ddd),'   Hour ',num2str(hhh),'    U=blue;   W=green' ]);
        xlabel('Frequency in Hz');
        ylabel('Spectrum * Frequency');
        if hhh==12,
            specfile=['d:\data\spec_files\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
    clear ww Tsonic UU VV jade jadson tt uu vv ww SS
else
    U=[];V=[];w=[];dir=[];Tslow=[];S=[];
end;