% Plot UM has logger data
% compare 8 levels of UMASS data all N booms
% 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
% This is for the files before Dave renamed them.

clear
close all

read_parameters_UM; %Reading parameter files flux conf

pltit = 1;
%DEW 192008  switched 300 and 200 until we figure out if levels are
%switched
for x = 7:19 % loop through days
    close all
    xs = num2str(x);
    f300 = 'w:\data\umass\data\processed\SONIC_250m\sonic_6_20070713.dat';
    f250 = 'w:\data\umass\data\processed\SONIC_250m\sonic_4_20070813.dat';
    f200 = 'w:\data\umass\data\processed\SONIC_250m\sonic_2_20070813.dat';
    f150 = 'w:\data\umass\data\processed\SONIC_100m\sonic_2_20070813.dat';
    f100 = 'w:\data\umass\data\processed\SONIC_100m\sonic_4_20070813.dat';
    f50 =  'w:\data\umass\data\processed\SONIC_100m\sonic_6_20070813.dat';
    f22 =  'w:\data\umass\data\processed\SONIC_base\sonic_2_20070813.dat';
    f10 =  'w:\data\umass\data\processed\SONIC_base\sonic_4_20070813.dat';

    if x <= 9
        f300(55) = '0';
        f300(56) = xs;
        f250(55) = '0';
        f250(56) = xs;
        f200(55) = '0';
        f200(56) = xs
        f150(55) = '0';
        f150(56) = xs;
        f100(55) = '0';
        f100(56) = xs;
        f50(55) = '0';
        f50(56) = xs;
        f22(55) = '0';
        f22(56) = xs;
        f10(55) = '0';
        f10(56) = xs;
    end
    if x >= 10
        f300(55:56) = xs;
        f250(55:56) = xs;
        f200(55:56) = xs;
        f150(55:56) = xs;
        f100(55:56) = xs;
        f50(55:56) = xs;
        f22(55:56) = xs;
        f10(55:56) = xs;
    end
    f300
    f250;
    f200;
    f150;
    f100;
    f50;
    f22;
    f10;

    st300 = load(f300);
    st250 = load(f250);
    st200 = load(f200);
    st150 = load(f150);
    st100 = load(f100);
    st50 = load(f50);
    st22 = load(f22);
    st10 = load(f10);

    %load w:\data\umass\data\processed\SONIC_250m\sonic_2_20070827.dat  %1 day 1 sec averaged data
    %st200 = sonic_2_20070827;
    %load w:\data\umass\data\processed\SONIC_250m\sonic_4_20070827.dat  %1 day 1 sec averaged data
    %st250 = sonic_4_20070827;
    %load w:\data\umass\data\processed\SONIC_250m\sonic_6_20070827.dat  %1 day 1 sec averaged data
    %st300 = sonic_6_20070827;
    %load w:\data\umass\data\processed\SONIC_100m\sonic_2_20070827.dat  %1 day 1 sec averaged data
    %st150 = sonic_2_20070827;
    %load w:\data\umass\data\processed\SONIC_100m\sonic_4_20070827.dat  %1 day 1 sec averaged data
    %st100 = sonic_4_20070827;
    %load w:\data\umass\data\processed\SONIC_100m\sonic_6_20070827.dat  %1 day 1 sec averaged data
    %st50 = sonic_6_20070827;
    %load w:\data\umass\data\processed\SONIC_base\sonic_2_20070827.dat  %1 day 1 sec averaged data
    %st22 = sonic_2_20070827;
    %load w:\data\umass\data\processed\SONIC_base\sonic_4_20070827.dat  %1 day 1 sec averaged data
    %st10 = sonic_4_20070827;
    tit = 'BAO 2007';
    % 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 10m 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.808,-999.000,238.61,-999.0000,13.515,843.01,0.000,212.610

    st300(:,13) = datenum(st300(:,1),st300(:,2),st300(:,3),st300(:,4),st300(:,5),st300(:,6))-datenum(st300(:,1)-1,12,31);
    st250(:,13) = datenum(st250(:,1),st250(:,2),st250(:,3),st250(:,4),st250(:,5),st250(:,6))-datenum(st250(:,1)-1,12,31);
    st200(:,13) = datenum(st200(:,1),st200(:,2),st200(:,3),st200(:,4),st200(:,5),st200(:,6))-datenum(st200(:,1)-1,12,31);
    st150(:,13) = datenum(st150(:,1),st150(:,2),st150(:,3),st150(:,4),st150(:,5),st150(:,6))-datenum(st150(:,1)-1,12,31);
    st100(:,13) = datenum(st100(:,1),st100(:,2),st100(:,3),st100(:,4),st100(:,5),st100(:,6))-datenum(st100(:,1)-1,12,31);
    st50(:,13) = datenum(st50(:,1),st50(:,2),st50(:,3),st50(:,4),st50(:,5),st50(:,6))-datenum(st50(:,1)-1,12,31);
    st22(:,13) = datenum(st22(:,1),st22(:,2),st22(:,3),st22(:,4),st22(:,5),st22(:,6))-datenum(st22(:,1)-1,12,31);
    st10(:,13) = datenum(st10(:,1),st10(:,2),st10(:,3),st10(:,4),st10(:,5),st10(:,6))-datenum(st10(:,1)-1,12,31);

    u10 = st10(:,7);
    v10 = st10(:,8);
    w10 = st10(:,9);
    t10 = st10(:,10);
    yys10 = num2str(st10(1,1));
    mms10 = num2str(st10(1,2));
    dds10 = num2str(st10(1,3));
    hhs10 = num2str(st10(1,4));
    mns10 = num2str(st10(1,5));
    sss10 = num2str(st10(1,6));

    u22 = st22(:,7);
    v22 = st22(:,8);
    w22 = st22(:,9);
    t22 = st22(:,10);
    yys22 = num2str(st22(1,1));
    mms22 = num2str(st22(1,2));
    dds22 = num2str(st22(1,3));
    hhs22 = num2str(st22(1,4));
    mns22 = num2str(st22(1,5));
    sss22 = num2str(st22(1,6));

    u50 = st50(:,7);
    v50 = st50(:,8);
    w50 = st50(:,9);
    t50 = st50(:,10);
    yys50 = num2str(st50(1,1));
    mms50 = num2str(st50(1,2));
    dds50 = num2str(st50(1,3));
    hhs50 = num2str(st50(1,4));
    mns50 = num2str(st50(1,5));
    sss50 = num2str(st50(1,6));

    u100 = st100(:,7);
    v100 = st100(:,8);
    w100 = st100(:,9);
    t100 = st100(:,10);
    yys100 = num2str(st100(1,1));
    mms100 = num2str(st100(1,2));
    dds100 = num2str(st100(1,3));
    hhs100 = num2str(st100(1,4));
    mns100 = num2str(st100(1,5));
    sss100 = num2str(st100(1,6));

    u150 = st150(:,7);
    v150 = st150(:,8);
    w150 = st150(:,9);
    t150 = st150(:,10);
    yys150 = num2str(st150(1,1));
    mms150 = num2str(st150(1,2));
    dds150 = num2str(st150(1,3));
    hhs150 = num2str(st150(1,4));
    mns150 = num2str(st150(1,5));
    sss150 = num2str(st150(1,6));

    u200 = st200(:,7);
    v200 = st200(:,8);
    w200 = st200(:,9);
    t200 = st200(:,10);
    yys200 = num2str(st200(1,1));
    mms200 = num2str(st200(1,2));
    dds200 = num2str(st200(1,3));
    hhs200 = num2str(st200(1,4));
    mns200 = num2str(st200(1,5));
    sss200 = num2str(st200(1,6));

    u250 = st250(:,7);
    v250 = st250(:,8);
    w250 = st250(:,9);
    t250 = st250(:,10);
    yys250 = num2str(st250(1,1));
    mms250 = num2str(st250(1,2));
    dds250 = num2str(st250(1,3));
    hhs250 = num2str(st250(1,4));
    mns250 = num2str(st250(1,5));
    sss250 = num2str(st250(1,6));

    u300 = st300(:,7);
    v300 = st300(:,8);
    w300 = st300(:,9);
    t300 = st300(:,10);
    yys300 = num2str(st300(1,1));
    mms300 = num2str(st300(1,2));
    dds300 = num2str(st300(1,3));
    hhs300 = num2str(st300(1,4));
    mns300 = num2str(st300(1,5));
    sss300 = num2str(st300(1,6));

    S10=sqrt(u10.*u10+v10.*v10);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir10=(atan2(u10,v10)*180/pi);
    for j = 1:length(dir10)
        if dir10(j) < 0
            dir10(j) = 360 + dir10(j);
        end
        dir10(j) = dir10(j) - 26;
        if dir10(j) < 0
            dir10(j) = 360 + dir10(j);
        end
    end

    S22=sqrt(u22.*u22+v22.*v22);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir22=(atan2(u22,v22)*180/pi);
    for j = 1:length(dir22)
        if dir22(j) < 0
            dir22(j) = 360 + dir22(j);
        end
        dir22(j) = dir22(j) - 26;
        if dir22(j) < 0
            dir22(j) = 360 + dir22(j);
        end
    end

    S50=sqrt(u50.*u50+v50.*v50);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir50=(atan2(u50,v50)*180/pi);
    for j = 1:length(dir50)
        if dir50(j) < 0
            dir50(j) = 360 + dir50(j);
        end
        dir50(j) = dir50(j) - 26;
        if dir50(j) < 0
            dir50(j) = 360 + dir50(j);
        end
    end

    S100=sqrt(u100.*u100+v100.*v100);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir100=(atan2(u100,v100)*180/pi);
    for j = 1:length(dir100)
        if dir100(j) < 0
            dir100(j) = 360 + dir100(j);
        end
        dir100(j) = dir100(j) - 26;
        if dir100(j) < 0
            dir100(j) = 360 + dir100(j);
        end
    end

    S150=sqrt(u150.*u150+v150.*v150);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir150=(atan2(u150,v150)*180/pi);
    for j = 1:length(dir150)
        if dir150(j) < 0
            dir150(j) = 360 + dir150(j);
        end
        dir150(j) = dir150(j) - 26;
        if dir150(j) < 0
            dir150(j) = 360 + dir150(j);
        end
    end

    S200=sqrt(u200.*u200+v200.*v200);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir200=(atan2(u200,v200)*180/pi);
    for j = 1:length(dir200)
        if dir200(j) < 0
            dir200(j) = 360 + dir200(j);
        end
        dir200(j) = dir200(j) - 26;
        if dir200(j) < 0
            dir200(j) = 360 + dir200(j);
        end
    end

    S250=sqrt(u250.*u250+v250.*v250);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir250=(atan2(u250,v250)*180/pi);
    for j = 1:length(dir250)
        if dir250(j) < 0
            dir250(j) = 360 + dir250(j);
        end
        dir250(j) = dir250(j) - 26;
        if dir250(j) < 0
            dir250(j) = 360 + dir250(j);
        end
    end

    S300=sqrt(u300.*u300+v300.*v300);
    %dir=(atan2(-v,u)*180/pi) - 26;
    dir300=(atan2(u300,v300)*180/pi);
    for j = 1:length(dir300)
        if dir300(j) < 0
            dir300(j) = 360 + dir300(j);
        end
        dir300(j) = dir300(j) - 26;
        if dir300(j) < 0
            dir300(j) = 360 + dir300(j);
        end
    end

    %  Plot Temperatures
    figure(1)
    plot(st10(:,13),st10(:,10),'k')
    hold
    plot(st22(:,13),st22(:,10),'r')
    plot(st50(:,13),st50(:,10),'b')
    plot(st100(:,13),st100(:,10),'k+')
    plot(st150(:,13),st150(:,10),'c')
    plot(st200(:,13),st200(:,10),'g')
    plot(st250(:,13),st250(:,10),'y')
    plot(st300(:,13),st300(:,10),'r.')

    legend('10','22','50','100','150','200','250','300')
    ylabel('T (C)')
    xlabel('Year Day')
    xlim([st10(1,13) st10(1,13)+1])
    ylim([min(st10(:,10))-5 max(st10(:,10)+5)])
    title(tit);
    if pltit
        print_buffer = [way_images_sonic 'UMPSD_T_' yys10 mms10 dds10 '_8lvls.png'];
        print('-dpng ', print_buffer);
    end

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

    figure(2)
    plot(st10(:,13),S10,'k')
    hold
    plot(st22(:,13),S22,'r')
    plot(st50(:,13),S50,'b')
    plot(st100(:,13),S100,'k+')
    plot(st150(:,13),S150,'c')
    plot(st200(:,13),S200,'g')
    plot(st250(:,13),S250,'y')
    plot(st300(:,13),S300,'r.')

    legend('10','22','50','100','150','200','250','300')
    ylabel('Spd (m/s)')
    xlabel('Year Day')
    xlim([st10(1,13) st10(1,13)+1])
    ylim([0 20])
    title(tit);
    if pltit
        print_buffer = [way_images_sonic 'UM_SPD_' yys10 mms10 dds10 '_8lvls.png'];
        print('-dpng ', print_buffer);
    end

    figure(3)
    plot(st10(:,13),dir10,'k')
    hold
    plot(st22(:,13),dir22,'r')
    plot(st50(:,13),dir50,'b')
    plot(st100(:,13),dir100,'k+')
    plot(st150(:,13),dir150,'c')
    plot(st200(:,13),dir200,'g')
    plot(st250(:,13),dir250,'y')
    plot(st300(:,13),dir300,'r.')
    legend('10','22','50','100','150','200','250','300')

    ylabel('Direction')
    xlabel('Year Day')
    xlim([st10(1,13) st10(1,13)+1])
    ylim([0 360])
    if pltit
        print_buffer = [way_images_sonic 'UM_DIR_' yys10 mms10 dds10 '_8lvls.png'];
        print('-dpng ', print_buffer);
    end

    if 0
        figure(3)
        plot(st10(:,13),st10(:,7:9))
        legend('u','v','w')
        ylabel('Spd (m/s)')
        xlabel('Year Day')
        xlim([st10(1,13) st10(1,13)+1])
        ylim([-10 10])
        title('10');
        if pltit
            print_buffer = [way_images_sonic 'UMPSD_uvw_' yys10 mms10 dds10 '_10.png'];
            print('-dpng ', print_buffer);
        end
        figure(4)
        plot(st22(:,13),st22(:,7:9))
        legend('u','v','w')
        ylabel('Spd (m/s)')
        xlabel('Year Day')
        xlim([st22(1,13) st10(1,13)+1])
        ylim([-10 10])
        title('22');
        if pltit
            print_buffer = [way_images_sonic 'UMPSD_uvw_' yys22 mms22 dds22 '_22.png'];
            print('-dpng ', print_buffer);
        end

        figure(6)
        plot(st10(:,13),st10(:,11),'k.')
        hold
        plot(st22(:,13),st22(:,11),'r.')
        legend('10','22','Location','Southwest')
        title('BAO 2007 UMass sonic npts');
        ylim([0 34]);

    end  % skip uvw and pts plots

end % loop throug daily plots