disp('rdson_07')   %('read_son_day_rhb_07')
fclose('all');
%Sonic model R3 used on AMMA06. Running at 10Hz
% S/N 000138

clear st1;clear st2;clear st3;clear jadson;clear sonraw;

jd=num2str(ddd);

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

plotit=1;
prtit=0;

hr='00';
jhold=0;
jax=1;
jaxx=1;
backchk=0;
st1=[];
st2=[];
st3=[];
jaxd=[];
%if 0 % DEW 5/15/2007
for jam=0:23,						%cycle thru 24 hourly stats files
    clear sonraw
    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;

    e=[way_raw_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
        % DEW commented out next line, because I son't see a blank space after SONIC title
        dmp2=fgets(flist);		%read header lines and skip blanks
        while feof(flist)==0,
                disp(['Reading Sonic file in rdson  ',e])
            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;
            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(2,:)/100;
v=st3(3,:)/100;
w=st3(4,:)/100;

if outputsonic==1
    Tsonic=st3(5,:)/100; %Temperature in degC
else
    speedson=st3(5,:)/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);
np=jhold;

if prtit==1;
    prt_jas_son_rhb_07;
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', year]);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', year]);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 (2006)'); ...
        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', year]);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