%BAOPLT07
% Plot data from GMD 300m BAO data loggers
% Dave has split these data into daily files and corrected the PRT and Wind Dir
%You are correct on all counts.  yes, m/s on the wind speed.
%The PRT temperature is (accurate to within 0.1 deg C) is:  a0 +
%a1*(D3_PRTRES_Avg) + a2*(D3_PRTRES_Avg)^2 + a3*(D3_PRTRES_Avg)^3,
%a0 =    -246.058
%a1 =      2.36626
%a2 =  0.000923807
%a3 =  2.43906e-007

%The sonic anemometer south to north pole point to 322 degrees.
%The the true wind direction should then be    (SE2_SONICWD_Avg) - 38 degrees
%
% 300m
%"TOA5","CR1000","CR1000","2846","CR1000.Std.11","CPU:060516b_Met.CR1","10044","wkt30m"
%"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"
%"2007-06-14 14:22:00",0,3.37,230.5,0.502,0.012,13.42,105.8471,14.24,15.09,42.15,0,110.6,13.31,15.0,192.5
%"2007-06-14 14:22:30",1,3.418,231.1,0.269,0.026,13.38,105.8718,14.3,15.11,42.26,0,110.6,13.31,15.1,193.1
clear
close all

jds = input('Input Julian day to plot    ','s');
jd = str2num(jds);
jdb = jd;
jde = jd+1;

fn = 's:\bao\Tower\Processed\BAO_300_2007226.dat';
fn(36:38) = jds(1:3);  % daily directory

f300 = fopen(fn,'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:10000
    x = x + 1;
    temp=fgetl(f300);
    year(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') ;
    bao300(x,1) = data(2,1); % WS ms Sonic
    bao300(x,2) = data(3,1); % WD deg Sonic
    bao300(x,3) = data(4,1); % WD  std deg Sonic
    bao300(x,4) = data(5,1); % WS  STD Sonic
    bao300(x,5) = data(6,1); % Temp deg C Sonic
    bao300(x,6) = data(7,1); % PRT temp resistance
    bao300(x,7) = data(8,1); % PRT temp deg C
    bao300(x,8) = data(9,1); % HMT Temp Deg C
    bao300(x,9) = data(10,1); % RH
    bao300(x,10) = data(11,1); % airflow
    bao300(x,11) = data(12,1); % VAC
    bao300(x,12) = data(13,1); % Bat VDC
    bao300(x,13) = data(14,1); % PRT corrected deg C
    bao300(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 = bao300;  % save data before any corrections to direction

juld = datenum(year,month,day,hour,minute,sec)-datenum(2006,12,31);
bao300(:,15) = juld; % Bat VDC

if 1

    figure(1)
    clf
    subplot(2,1,1),plot(bao300(:,15),bao300(:,5),'k')
    hold
    subplot(2,1,1),plot(bao300(:,15),bao300(:,7),'r')
    subplot(2,1,1),plot(bao300(:,15),bao300(:,8),'g')
    subplot(2,1,1),plot(bao300(:,15),bao300(:,13),'b')
    
    %xlabel('Year Day')
    ylabel('T (C)')
    title('BAO 2007 300m')
    ylim([0,40])
    %axis([110,116,-10,30])
    legend('Sonic T','PRT','HMT','PTRc','Location','Best')
    %legend('10m','50m')
    xlim([jdb,jde])

    subplot(2,1,2),plot(bao300(:,15),bao300(:,9))
    %hold
    %subplot(2,1,2),plot(data50(:,3),data50(:,6),'r')
    xlabel('Year Day')
    ylabel('Rh (%)')
    %axis([110,116,0,100])
    %legend('300m')
    %legend('10m','50m')
    ylim([0,109])
    xlim([jdb,jde])

    figure(2)
    clf
    subplot(2,1,1),plot(bao300(:,15),bao300(:,1))
    %hold
    %subplot(2,1,1),plot(data50(:,3),data50(:,7),'r')
    %xlabel('Year Day')
    ylabel('Spd (m/s)')
    title('BAO 2007 300m')
    %axis([110,116,0,30])
    ylim([0,max(bao300(:,1))+5])
    xlim([jdb,jde])
    %legend('300m')
    %legend('10m','50m')
    subplot(2,1,2),plot(bao300(:,15),bao300(:,2))
    %hold
    %subplot(2,1,2),plot(data50(:,3),data50(:,8),'r')
    xlabel('Year Day')
    ylabel('Dir (true)')
    xlim([jdb,jde])
    %axis([110,116,0,360])
    %legend('300m')
    %legend('10m','50m')
    xlim([jdb,jde])
    
figure(4)
clf
plot(bao300(:,15),bao300(:,5),'k')
hold
plot(bao300(:,15),bao300(:,7),'r')
plot(bao300(:,15),bao300(:,8),'b')
plot(bao300(:,15),bao300(:,13),'g')
xlabel('Year Day (UTC)')
ylabel('T (C)')
title('BAO 2007 300m')
ylim([min(bao300(:,13))-3,max(bao300(:,13))+3])
%axis([110,116,-10,30])
legend('Sonic T','PRT est','HMT T','PRT corrected','Location','Best')
xlim([jdb,jde])

end
save bao300d.dat bao300 -ascii
