disp('read_son_day_sjohnson_05')
fclose('all');

clear st1;clear st2;clear st3;clear jadson;
jd=num2str(ddd);

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

plotit=1;
prtit=1;

hr='00';
jhold=0; 
jax=1;
jaxx=1;
backchk=0;
st1=[];
st2=[];
st3=[];
for jam=0:23,						%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;
        st2=[]; st1=[];
        jax=1;
     end;       %if isempty
     fclose(flist);
   end;         %end if flist   
end;            %end for
clear u v w Tsonic U V dir speedson
u=-st3(2,:)/100;
v=st3(3,:)/100;
w=st3(4,:)/100;
%speedson=st3(5,:);

%for tst=1:length(speedson)
%    if speedson(1,tst)<-30000    %Fixes SOS to proper value
%       speedson(1,tst)=(speedson(1,tst)+32768*2)/100;
%    else
%        %speedson(1,tst)=sqrt((speedson(1,tst)/100)*403); %Fixes temp K
%        speedson(1,tst)=sqrt(((speedson(1,tst)/100)+273.15)*403);
%        %FixesTemp C
%    end
%end

Tsonic=st3(5,:)/100;    
%Tsonic=(speedson.*speedson)/403-273.15;
%Tsonic=st3(5,:)/100-273.15; %From read_son_day_epi_etl_03.m
%Tsonic=(speedson.*speedson)/403+(u.*u/403+ w.*w/403+v.*v/403)/3+273.15; %From read_son_day_dana.m
%Tsonic=(speedson.*speedson+(u.*u + w.*w)/2. + v.*v)/403-273.15;  %Original from NAME
%U=u;    
U=-u*cos(30/180*pi)-v*sin(30/180*pi);		   %to North
%V=v;    
V=-u*sin(30/180*pi)+v*cos(30/180*pi);  		   %to West
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=jhold;

if prtit==1;
  prt_jas_son_sjohnson_05;
end;

if plotit==1;
   %set some plot limits
   Tll=-25;
   hll=70;
   uul=10;
   vul=10;
   Tul=-0;
   hul=100;
   sul=12;
   Rlul=30;
   Rlll=-100;
   Rtul=300;
   Rtll=125;
   Rsul=1000;
   Rsll=0;
   
   %plot T profiles
   j_st=jadson(1);
   j_en=jadson(np);
   figure;plot(jadson(1:np),Tsonic(1:np));xlabel('Jday (2005)');ylabel('Temperature (C)');title('Sonic Temperature'); 
   figure;plot(jadson(1:np),U(1:np),jadson(1:np),V(1:np),jadson(1:np),w(1:np),jadson(1:np),S(1:np));xlabel('Jday (2005)');ylabel('Sonic Winds (m/s)');title('Relative Wind Components:  U=blue,  V=green,  W=red,   Speed=cyan'); 

   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(jadson(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('Jday (2005)'); ...
	ylabel('Sonic Wind Direction (deg)'); ...
	title('Relative Wind Direction'); ...
	axis([j_st j_en -180 180]);       
    figure;plot(jadson(1:np),sqrt(U(1:np).^2+V(1:np).^2+w(1:np).^2));xlabel('Jday (2005)');ylabel('Wind Speed m/s');title('Sonic Relative Wind Speed');%axis([j_st j_en 0 360]);
end;

%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