% Plot UM and PSD 50m data
% UMass sonic 2D Model 81000
% u is E/W + fm E
% v is N/S + fm N
% w -s up/dn + is up
% on N boom +v = 334
% on S boom +v = 154

clear
close all

read_parameters_UM; %Reading parameter files flux conf

pltit = 1;

% UMASS data start Jul 7 2007 JD 188
%PSD first 100m data Jul 17 JD 198
%jd = 212; %Jul 31 2007
%jd = 197; %Jul 6 2007
jd = 224; %must be day before loop day
for x =13:28    %31 % loop through days
    jd = jd + 1;
    close all
    xs = sprintf('%02d',x);
    jds = num2str(jd);
    %psd = ['z:\BAO\Tower\processed\daily\BAO_50_2007' jds '.dat'];
    f50s = ['z:\BAO\Tower\UMASS\processed\sonic_050s_200708' xs '.dat']
    f50n = ['z:\BAO\Tower\UMASS\processed\sonic_050n_200708' xs '.dat']
    st1 = load(f50n);
    st2 = load(f50s);

    %load w:\data\umass\data\processed\sonic_base\sonic_4_20070811.dat  %1 day 1 sec averaged data
    %st1 = sonic_4_20070811;
    %psd10 = load(psd);
    %psd10 = BAO_010_2007223;
    %clear BAO_010_2007223
    tit = 'BAO 2007 50m';
    % yy  mm dd hh mn ss    u       v     w       t     npt bat
    %2007 08 20 00 00 00    3.89   0.50  -0.40   27.10  32 11.1
    %2007 08 20 00 00 01    3.51   0.85  -0.17   27.09  32 11.1

    %%%%%%% PSD 50m data
    %1    2   3   4     5     6      7      8       9      10        11     12     13    14
    %ID  YR  YD  hhmm  T      RH     WSv   WSs     Dir    WDt       Bat    press  precip WDT
    %10,2007,165,1405,15.095,50.858,1.810,-999.000,238.61,-999.0000,13.515,843.
    %01,0.000,212.610
if 0  % skip until 50m gets sensors
    hhmm = psd10(:,4);
    hh = fix(hhmm/100);
    mn = hhmm - (hh*100);
    psd10(:,3) = psd10(:,3) + (((mn/60.) + hh)/24); %Julian day
end
    st1(:,13) = datenum(st1(:,1),st1(:,2),st1(:,3),st1(:,4),st1(:,5),st1(:,6))-datenum(st1(:,1)-1,12,31)-212;
    st2(:,13) = datenum(st2(:,1),st2(:,2),st2(:,3),st2(:,4),st2(:,5),st2(:,6))-datenum(st2(:,1)-1,12,31)-212;

    un = st1(:,7);
    vn = st1(:,8);
    wn = st1(:,9);
    tn = st1(:,10);
    yysn = num2str(st1(1,1));
    mmsn = num2str(st1(1,2));
    ddsn = num2str(st1(1,3));
    hhsn = num2str(st1(1,4));
    mnsn = num2str(st1(1,5));
    sssn = num2str(st1(1,6));

    us = st2(:,7);
    vs = st2(:,8);
    ws = st2(:,9);
    ts = st2(:,10);
    yyss = num2str(st2(1,1));
    mmss = num2str(st2(1,2));
    ddss = num2str(st2(1,3));
    hhss = num2str(st2(1,4));
    mnss = num2str(st2(1,5));
    ssss = num2str(st2(1,6));
    lvl = '50m';

    %  Plot Temperatures
    figure(1)
    plot(st1(:,13),st1(:,10),'k')
    hold
    plot(st2(:,13),st2(:,10),'b')
    %plot(psd10(:,3),psd10(:,5),'r')
    legend('UMass N','UMass S')      %,'PSD')
    ylabel('T (C)')
    xlabel('August')
    mx = max(st1(:,13));
    xlim([st1(1,13) mx])
    ylim([min(st1(:,10))-5 max(st1(:,10)+5)])
    title(tit);
    if pltit
        print_buffer = [way_images_sonic 'UM_T_' yyss mmss ddss '_' lvl '.png'];
        print('-dpng ', print_buffer);
    end

    %Correction for Norht boom is 26 deg and fgor South boom is 206 deg
    Sn=sqrt(un.*un+vn.*vn);
    dirn=(atan2(un,vn)*180/pi);
    for j = 1:length(dirn)
        if dirn(j) < 0
            dirn(j) = 360 + dirn(j);
        end
        dirn(j) = dirn(j) - 26;
        if dirn(j) < 0
            dirn(j) = 360 + dirn(j);
        end
    end

    Ss=sqrt(us.*us+vs.*vs);
    dirs=(atan2(us,vs)*180/pi);
    for j = 1:length(dirs)
        if dirs(j) < 0
            dirs(j) = 360 + dirs(j);
        end
        dirs(j) = dirs(j) - 206;
        if dirs(j) < 0
            dirs(j) = 360 + dirs(j);
        end
    end

    % on N boom +v = 334
    % on S boom +v = 154
    % 50m is on N boom so we need to add 26 degrees to directions

    if 1  % skip speed and direction plots
        figure(2)
        subplot(2,1,1),plot(st1(:,13),Sn,'k')
        hold
        subplot(2,1,1),plot(st2(:,13),Ss,'b')
        %subplot(2,1,1),plot(psd10(:,3),psd10(:,7)*.51444,'r.')
        legend('UMass N','UMass S')     %,'PSD')
        ylabel('Spd (m/s)')
        xlabel('August')
        xlim([st1(1,13) mx])
        ylim([0 20])
        title(tit);

        subplot(2,1,2),plot(st1(:,13),dirn,'k')
        hold
        subplot(2,1,2),plot(st2(:,13),dirs,'b')
        %subplot(2,1,2),plot(psd10(:,3),psd10(:,14),'r.')
        ylabel('Direction')
        xlabel('August')
        xlim([st1(1,13) mx])
        ylim([0 360])
        if pltit
            print_buffer = [way_images_sonic 'UM_SPDDIR_' yyss mmss ddss '_' lvl '.png'];
            print('-dpng ', print_buffer);
        end

        %plot UVW
        figure(3)
        plot(st1(:,13),st1(:,7:9))
        legend('u','v','w')
        ylabel('Spd (m/s)')
        xlabel('August')
        xlim([st1(1,13) mx])
        ylim([-10 10])
        tit = 'BAO 2007 N 50m';
        title(tit);
        if pltit
            print_buffer = [way_images_sonic 'UM_uvwS_' yyss mmss ddss '_' lvl '.png'];
            print('-dpng ', print_buffer);
        end

         figure(4)
        plot(st2(:,13),st2(:,7:9))
        legend('u','v','w')
        ylabel('Spd (m/s)')
        xlabel('August')
        xlim([st2(1,13) mx])
        ylim([-10 10])
        tit = 'BAO 2007 S 50m';
        title(tit);
        if pltit
            print_buffer = [way_images_sonic 'UM_uvwN_' yyss mmss ddss '_' lvl '.png'];
            print('-dpng ', print_buffer);
        end
        
    end % skip direction plots until 50m gets a propvane

    figure
    plot(st1(:,13),st1(:,11),'o')
    hold
    plot(st2(:,13),st2(:,11),'r.')    
    title('BAO 2007 UMass  50m sonic npts');
    legend('N','S')
    xlabel('August');
    ylabel('Number of Pts');
    ylim([0 34]);
    if pltit
        print_buffer = [way_images_sonic 'UM_npts_' yyss mmss ddss '_' lvl '.png'];
        print('-dpng ', print_buffer);
    end

end  %loop daily files
