% Plot UM and PSD 300m 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;

%jd = 212; %Jul 31 2007
% 1    2   3   4   5   6   7   8   9   10
% 213 214 215 216 217 218 219 220  221  222
% UMass data start Aug 7 JD 219
jd = 224;
for x = 13:14 % loop through days
    jd = jd + 1;
    close all
    xs = sprintf('%02d',x)
    jds = num2str(jd);
    psd = ['z:\BAO\Tower\processed\daily\BAO_300_2007' jds '.dat'];
    f3 = ['z:\BAO\Tower\UMASS\processed\sonic_300n_200708' xs '.dat'];
    
    st1 = load(f3);

%load w:\data\umass\data\processed\sonic_250m\sonic_6_20070819.dat  %1 day 1 sec averaged data
%st1 = sonic_6_20070819;
%fn = 'w:\data\bao\data\BAO_300_2007231.dat';
%fn = 'z:\bao\Tower\Processed\daily\BAO_300_2007240.dat';
%fn(36:38) = jds(1:3);  % daily directory

f300 = fopen(psd,'r');
%load f:\BAO_CR10\data\bao10.dat
%read in 4 header lines
tmp=fgetl(f300);
tmp=fgetl(f300);
tmp=fgetl(f300);
tmp=fgetl(f300);
x = 0;
while ~feof(f300)
    %for x = 1:1000
    x = x + 1;
    temp=fgetl(f300);
    yr(x)= str2num(temp(2:5));
    month(x)= str2num(temp(7:8));
    day(x)= str2num(temp(10:11));
    hour(x)= str2num(temp(13:14));
    minute(x)= str2num(temp(16:17));
    sec(x)= str2num(temp(19:20));
    lt= length(temp);
    % strip off date/time info from first part of each line
    s1 = temp(23:lt);
    % convert each , to a blank space
    % then read into an array
    for t=1:length(s1)
        if s1(t) == ','
            s1(t) = ' ';
        end
    end
    %"TIMESTAMP","RECORD","SE1_SONICWS_Avg","SE2_SONICWD_Avg","SE2_SONICWD_Std","SE1_SONICWS_Std","SE3_SONICTEMP_Avg","D3_PRTRES_Avg","D3_PRTTEMPAPPROX_Avg","SE4_HMPTEMP_Avg","SE9_HMPRH_Avg","P1_AIRFLOW_Avg","SE10_ACPOWER","BT_BATVOLT_Min","PRT_AVG","SONICWDT_AVG"
    %"TS","RN","ms","Deg","Deg","ms","mV","","Deg C","Deg C","%","Counts","ms","Volts","DegC","Deg True"

    [data]=sscanf(s1,'%f') ;
    psd300(x,1) = data(2,1); % WS ms Sonic
    psd300(x,2) = data(3,1); % WD deg Sonic
    psd300(x,3) = data(4,1); % WD  std deg Sonic
    psd300(x,4) = data(5,1); % WS  STD Sonic
    psd300(x,5) = data(6,1); % Temp deg C Sonic
    psd300(x,6) = data(7,1); % PRT temp resistance
    psd300(x,7) = data(8,1); % PRT temp deg C
    psd300(x,8) = data(9,1); % HMT Temp Deg C
    psd300(x,9) = data(10,1); % RH
    psd300(x,10) = data(11,1); % airflow
    psd300(x,11) = data(12,1); % VAC
    psd300(x,12) = data(13,1); % Bat VDC
    psd300(x,13) = data(14,1); % PRT corrected deg C
    psd300(x,14) = data(15,1); % Dir true


end

%Correct direction for pointing angle on tower
% Sonic is pointed 322 degres according to Jonathan Kofler
% This means N is reall 322 so we subtract 38
b300b = psd300;  % save data before any corrections to direction

juld = datenum(yr,month,day,hour,minute,sec)-datenum(yr-1,12,31);
psd300(:,15) = juld; % Bat VDC


tit = 'BAO 2007 300m';
% 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.810,-999.000,238.61,-999.0000,13.515,843.
%01,0.000,212.610

hhmm = psd300(:,4);
hh = fix(hhmm/100);
mn = hhmm - (hh*100);
psd300(:,3) = psd300(:,3) + (((mn/60.) + hh)/24); %Julian day
st1(:,13) = datenum(st1(:,1),st1(:,2),st1(:,3),st1(:,4),st1(:,5),st1(:,6))-datenum(st1(:,1)-1,12,31);

u = st1(:,7);
v = st1(:,8);
w = st1(:,9);
t = st1(:,10);
yys = num2str(st1(1,1));
mms = num2str(st1(1,2));
dds = num2str(st1(1,3));
hhs = num2str(st1(1,4));
mns = num2str(st1(1,5));
sss = num2str(st1(1,6));
lvl = '300m';

%  Plot Temperatures
figure(1)
clf
plot(st1(:,13),st1(:,10),'r')
hold
plot(psd300(:,15),psd300(:,13),'b')
legend('UMass','PSD')
ylabel('T (C)')
xlabel('Year Day')
mx = max(st1(:,13));
xlim([st1(1,13) mx])
ylim([min(psd300(:,13))-5 max(psd300(:,13)+5)])
title(tit);
if 1 %pltit
    print_buffer = [way_images_sonic 'UMPSD_T_' yys mms dds '_' lvl '.png'];
    print('-dpng ', print_buffer);
end

S=sqrt(u.*u+v.*v);
%dir=(atan2(-v,u)*180/pi) - 26;
dir=(atan2(u,v)*180/pi);
for j = 1:length(dir)
    if dir(j) < 0
        dir(j) = 360 + dir(j);
    end
                 dir(j) = dir(j) - 26;
        if dir(j) < 0
            dir(j) = 360 + dir(j);
        end
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)
subplot(2,1,1),plot(st1(:,13),S)
hold
subplot(2,1,1),plot(psd300(:,15),psd300(:,1),'r.')
legend('UMass','PSD')
ylabel('Spd (m/s)')
xlabel('Year Day')
xlim([st1(1,13) mx])
ylim([0 20])
title(tit);

subplot(2,1,2),plot(st1(:,13),dir)
hold
subplot(2,1,2),plot(psd300(:,15),psd300(:,14),'r.')
ylabel('Direction')
xlabel('Year Day')
xlim([st1(1,13) mx])
ylim([0 360])
if pltit
    print_buffer = [way_images_sonic 'UMPSD_SPDDIR_' yys mms dds '_' lvl '.png'];
    print('-dpng ', print_buffer);
end


figure(3)
plot(st1(:,13),st1(:,7:9))
legend('u','v','w')
ylabel('Spd (m/s)')
xlabel('Year Day')
xlim([st1(1,13) mx])
ylim([-10 10])
title(tit);
if pltit
    print_buffer = [way_images_sonic 'UM_uvw_' yys mms dds '_' lvl '.png'];
    print('-dpng ', print_buffer);
end

figure
plot(st1(:,13),st1(:,11),'.')
title('BAO 2007 UMass sonic npts');
ylim([0 34]);
if pltit
    print_buffer = [way_images_sonic 'UM_npts_' yys mms dds '_' lvl '.png'];
    print('-dpng ', print_buffer);
end

end %loop daily files