disp('rd_sonic_plt')
close all;
clear
%Licro 7500
% S/N 000138

bd = 1;
ed = 0;
ddd=input('Input yearday to evaluate.... ');

read_parameters_rhb; %Reading parameter files flux conf from D:/data/lb/flux_eval_rhb_06

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 files
    clear h2o co2 co ho lt p code jdn jad jade sonraw sonraw_dep st1
    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
    %mmss      U         V      W     T      QC
    %  1       2         3      4     5       6
    flist=fopen(e,'r');
   
   if flist>0,	   				%if the file exists,
        disp(['Reading Sonic file ',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
                if dtmp(12) ~= 'd' & dtmp(12) ~= 'C'
            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
            end % test for good data
        end  %WHILE LOOP
        ss = size(st1)
        if ss(1) == 0
            break
        end
        fclose(flist);
        %      sonfac=1200 or 600; %sonic runs at 10Hz, 36000 values/hour...
        sonfac=3600*fsonic/60;
        jade=  st1(1,:);  %  Julian day already in col 1   str2num(jd)+((st2(1,:)-tzero)/1000/60/60+tzilch)/24;
        % save only every second
        r = 0;
        for k = 1:100:length(st1);
            r = r + 1;
            son(bd+r-1,1) = jade(1,k);
            son(bd+r-1,2) = st1(2,k);
            son(bd+r-1,3) = st1(3,k);
            son(bd+r-1,4) = st1(4,k);
            son(bd+r-1,5) = st1(5,k);
        end
        ed = bd + r - 1;
        bd = ed + 1;
    end;  %while loop
end;       %end jam
if 1
    jds = num2str(fix(jade(1,1)));
    f=['c:\data\icealot\data\sonic\ICEALOT_sonic_' jds '.txt'];
    flist=fopen(f,'w');
    save( f, 'son', '-ascii','-tabs')
    fclose('all');
end;

if 1
   
    u = son(:,2); %m/s
    v = son(:,3);  %m/s
    w = son(:,4); %m/s
    Tsonic = son(:,5);  %C
    U=u;
    V=v;
    jad = son(:,1); % year day
    jds = num2str(fix(jad(1)));
    s = 'c:\data\icealot\data\sonic\Raw_Images\';
    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);
        j_st = str2num(jd);
        j_en= j_st+1;

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

        figure;
        plot(son(:,1),U,son(:,1),V);
        xlim([j_st j_en]);
        xlabel(['Year Day (UTC)', year]);ylabel('Speed (m/s)');title('Sonic Relative Wind Components');
        legend('U','V');
        print_buffer = ['c:\data\icealot\data\sonic\Raw_Images\ICEALOT_Sonic_UV_' jd '.png'];
        print('-dpng ', print_buffer);
        
        figure;
        plot(son(:,1),w);
        xlim([j_st j_en]);
        xlabel(['Year Day (UTC)', year]);ylabel('Sonic Winds (m/s)');title('Sonic W');
        print_buffer = ['c:\data\icealot\data\sonic\Raw_Images\ICEALOT_Sonic_W_' jd '.png'];
        print('-dpng ', print_buffer);
        
        figure;
        subplot(2,1,1),plot(son(:,1),dir)
        ylabel('Sonic Direction (deg)'); ...
            title('Sonic Relative Wind Direction'); ...
            axis([j_st j_en 0 360]);

        subplot(2,1,2),plot(son(:,1),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\Raw_Images\ICEALOT_Sonic_WIND_' jd '.png'];
        print('-dpng ', print_buffer);
end %skip any plots


