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 jade;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:12:12,						%cycle thru 24 hourly stats files
    clear sonraw jade sonraw_dep st1 st2 st3 stx jdn
    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;

    e=[way_raw_data, 'day', jd,'\son008' jd hr '_raw.txt'];
    %03/20/2008 00:00 Gill-Sonic Anemometer  (SN: 070514)
    %0000039 B,+001.16,+000.02,-000.09,M,+014.66,00,2D
    flist=fopen(e,'r');
    if flist>0,	   				%if the file exists,
        disp(['Reading Sonic file in rdson  ',e])
        dmp=fgets(flist);			%read 1st header line
        YY = str2num(dmp(7:10));
        MM = str2num(dmp(1:2));
        DD = str2num(dmp(4:5));
        hh = str2num(dmp(12:13));
        ict = 0;
        while feof(flist)==0,
            % for p = 1:100
            dtmp=fgetl(flist);
            if length(dtmp) == 51
            ict = ict + 1;
            tim = dtmp(1:8);
            mm = str2num(tim(:,1:2));
            ss = str2num(tim(:,3:8))/1000;
            %calculate year day
            yd=datenum(YY,MM,DD,hh,0,0)-datenum(YY-1,12,31);
            jdn(ict,1) = yd + (((mm+(ss/60))/60))/24;
            st1(1,ict) = jdn(ict,1);
            % despike data  DEW 5/6/2007 fm CWK
            st1(2,ict) = str2num(dtmp(12:18));  %  U
            st1(3,ict) = str2num(dtmp(20:26));  %  V
            st1(4,ict) = str2num(dtmp(28:34));  %  W
            st1(5,ict) = str2num(dtmp(38:44));  %  T
            end % test for good data
        end  %WHILE LOOP
        ss = size(st1)
        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
        fclose(flist);
           jade=  st1(1,:);  %  Julian day already in col 1   str2num(jd)+((st2(1,:)-tzero)/1000/60/60+tzilch)/24;

    clear u v w Tsonic U V dir
    u=st1(2,:);%;
    v=st1(3,:);%;
    w=st1(4,:);%;
    %      sonfac=1200 or 600; %sonic runs at 10Hz, 36000 values/hour...
    sonfac=3600*fsonic/60;
    %st3 = st2;
 
    %% DEW 5/18/2007
    if 0
        plot(u,'k')
        hold
        plot(v,'b')
        plot(w,'r')
        title('Sonic raw components')
        legend('u','v','w')
    end  %skip
    if outputsonic==1
        Tsonic=st1(5,:);%     /100; %Temperature in degC
    else
        speedson=st1(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
        np = length(jade);
        j_st=jade(1);
        j_en=jade(np);

        figure;
        plot(jade(1:np),Tsonic(1:np));
        xlabel(['Year Day (UTC)', year]);ylabel('Temperature (C)');title('Sonic Temperature');
        xlim([j_st j_en]);
        print_buffer = ['c:\data\icealot\data\sonic\ICEALOT_Sonic_T_' jd '_' hr '.png'];
        print('-dpng ', print_buffer);

        figure;
        plot(jade(1:np),U(1:np),jade(1:np),V(1:np),jade(1:np),w(1:np),jade(1:np),S(1:np));
        xlabel(['Year Day (UTC)', year]);ylabel('Sonic Winds (m/s)');title('Relative Wind Components:  U=blue,  V=green,  W=red,   Speed=cyan');
        print_buffer = ['c:\data\icealot\data\sonic\ICEALOT_Sonic_UVW_' jd '_' hr '.png'];
        print('-dpng ', print_buffer);

        if 0  %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;
        subplot(2,1,1),plot(jade(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'); ...
        ylabel('Sonic Direction (deg)'); ...
            title('Relative Wind Direction'); ...
            axis([j_st j_en 0 360]);

        subplot(2,1,2),plot(jade(1:np),S);
        xlabel('Year Day (UTC)'); ...
            ylabel('Sonic Speed m/s');
        title('Sonic Relative Wind Speed')
        axis([j_st j_en 0 30]);
        print_buffer = ['c:\data\icealot\data\sonic\ICEALOT_Sonic_WIND_' jd '_' hr '.png'];
        print('-dpng ', print_buffer);

        figure
        plot(sonraw(2,:),'b')
        hold
        plot(sonraw_dep(2,:),'b.')
        plot(sonraw_dep(3,:),'g')
        plot(sonraw_dep(4,:),'r')
        legend('raw-U','DS-U','DS-V','DS-W','Location','Best')
        title('Despiking')
        print_buffer = ['c:\data\icealot\data\sonic\ICEALOT_Sonic_despk_' jd  '_' hr '.png'];
        print('-dpng ', print_buffer);
    end;
    sonraw = st1;   % load sonic array into temp array
    sonfix  % run despiking routine
    son_spc  % calculate sonic spectra was rdson2_08.m
        end;       %end flist

end;         %end of hour loop jam

%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 jade 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