%BAOPLT300_loop
% 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 is on the south boom
% The north pole points 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

pltflg = 1; %flag to plot the data: 1=yes, 0=no
pngflg = 1; %flag to make a png of the plot: 1=yes, 0=no

%jds = input('Input Julian day to plot    ','s');
%jd = str2num(jds);

inpath = 'z:\bao\Tower\Processed\daily\';
outpath  = 'z:\bao\Tower\Processed_Images\daily\';

cd(inpath);

m=dir('BAO_300_200802*');
lm = length(m);
fprintf(1,'number of files to process %d\n', lm);

nmax = 3000;
year = zeros(1,nmax);
month = zeros(1,nmax);
day = zeros(1,nmax);
hour = zeros(1,nmax);
minute = zeros(1,nmax);
sec = zeros(1,nmax);
yday = zeros(1,nmax);
bao300 = zeros( nmax, 15);

% 76 days 69 = 233
for ii = 1:lm
    disp(['Reading ',m(ii).name]);
    f300 = fopen(m(ii).name,'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);
    bao300 = zeros( 3000, 15);
    n = 0;
    while ~feof(f300)
        n = n + 1;
        temp=fgetl(f300);
        year(n)= str2double(temp(2:5));
        month(n)= str2double(temp(7:8));
        day(n)= str2double(temp(10:11));
        hour(n)= str2double(temp(13:14));
        minute(n)= str2double(temp(16:17));
        sec(n)= str2double(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
        s1 = strrep( s1, ',', ' ' );

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

    yday = datenum(year,month,day,hour,minute,sec)-datenum(year-1,12,31);
    bao300(1:n,15) = yday(1:n); % Year Day
    jdb = fix(bao300(1,15));
    jde = fix(bao300(n,15)) + 1;

    if pltflg % do plots if pltflg=1

        figure(1)
        clf
        subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,5),'k')
        hold
        %subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,7),'r')
        subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,8),'g')
        subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,13),'b')

        %xlabel('Year Day')
        ylabel('T (C)')
        stitle=sprintf('BAO %d 300m', year(1));
        title(stitle)
        ylim([0,40])
        %axis([110,116,-10,30])
        legend('Sonic T','HMP','PRTc','Location','Best')
        ylim([min(bao300(1:n,13))-3,max(bao300(1:n,13))+3])
        xlim([jdb,jde])

        subplot(2,1,2),plot(bao300(1:n,15),bao300(1:n,9))
        %hold
        %subplot(2,1,2),plot(data300(1:n,3),data300(1:n,6),'r')
        xlabel('Year Day')
        ylabel('Rh (%)')
        %axis([110,116,0,100])
        %legend('300m')
        %legend('10m','50m')
        ylim([0,109])
        xlim([jdb,jde])
        if pngflg
            print_buffer = [outpath, m(ii).name(1:15),'_TRH.png'];
            print('-dpng ', print_buffer);
        end

        figure(2)
        clf
        subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,1))
        %hold
        %subplot(2,1,1),plot(data300(1:n,3),data300(1:n,7),'r')
        %xlabel('Year Day')
        ylabel('Spd (m/s)')
        title(stitle)
        %axis([110,116,0,30])
        ylim([0,max(bao300(1:n,1))+5])
        xlim([jdb,jde])
        %legend('300m')
        %legend('10m','50m')
        subplot(2,1,2),plot(bao300(1:n,15),bao300(1:n,14))
        %hold
        %subplot(2,1,2),plot(data300(1:n,3),data300(1:n,8),'r')
        xlabel('Year Day')
        ylabel('Dir (true)')
        xlim([jdb,jde])
        %axis([110,116,0,360])
        %legend('300m')
        %legend('10m','50m')
        xlim([jdb,jde])
        if pngflg
            print_buffer = [outpath, m(ii).name(1:15),'_WS.png'];
            print('-dpng ', print_buffer);
        end

        figure(3)
        clf
        plot(bao300(1:n,15),bao300(1:n,5),'k')
        hold
        plot(bao300(1:n,15),bao300(1:n,7),'r')
        plot(bao300(1:n,15),bao300(1:n,8),'b')
        plot(bao300(1:n,15),bao300(1:n,13),'g')
        xlabel('Year Day (UTC)')
        ylabel('T (C)')
        
        title(stitle)
        ylim([min(bao300(1:n,13))-3,max(bao300(1:n,13))+3])
        %axis([110,116,-10,30])
        legend('Sonic T','PRT est','HMP T','PRT corrected','Location','Best')
        xlim([jdb,jde])
        if pngflg
            print_buffer = [outpath, m(ii).name(1:15),'_Tall.png'];
            print('-dpng ', print_buffer);
        end
    % Scatter plot of HMP vs corrected PRT at 300m
        figure(4)
        clf
        plot(bao300(1:n,8),bao300(1:n,13),'k.')
        xlabel('HMP (C)')
        ylabel('PRT (C)')
        title(stitle)
        axis('square')
        %ylim([min(bao300(1:n,15))-3,max(bao300(1:n,15))+3])
        %axis([110,116,-10,30])
        %xlim([jdb,jde])
        if pngflg
            print_buffer = [outpath, m(ii).name(1:15),'_HMPPRT.png'];
            print('-dpng ', print_buffer);
        end
    end % if doing plots

fclose(f300);
end % loop though each day
