%BAOPLT300_month_loop
% Plot data from GMD 300m BAO data loggers
% Reads in 10 min file which has a different format than the 30 sec data
% 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","Public"
%"TIMESTAMP","RECORD","BT_BATVOLT","SE1_SONICWS","SE2_SONICWD","SE3_SONICTEMP","D3_PRTTEMPAPPROX","SE4_HMPTEMP","SE9_HMPRH","SE10_ACPOWER","P1_AIRFLOW","SW12_12VFLAG","D3_PRTRES"
%"TS","RN","Volts","mV","mV","mV","Deg C","Deg C","%","mV","Counts","",""
%"2007-06-14 14:18:58.79",4479243,13.32056,3.517115,229.1796,13.54004,14.05173,14.93604,42.78085,110.567,0,1,105.7731
%"2007-06-14 14:28:56.39",4479841,13.2959,3.516633,231.2331,13.42691,14.39247,15.26545,40.26095,110.634,0,1,105.9062

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

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

cd(inpath);

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

nmax = 90000;
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);

for ii = 1:lm
    disp(['Reading ',m(ii).name]);
    f300 = fopen(m(ii).name,'rt');

    %read in 4 header lines
    tmp=fgetl(f300);
    tmp=fgetl(f300);
    tmp=fgetl(f300);
    tmp=fgetl(f300);
    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

    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,6),'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', data300(1,15));
        title(stitle)
        ylim([min(bao300(1:n,13))-3,max(bao300(1:n,13))+3])
        %axis([110,116,-10,30])
        legend('Sonic T','HMP','PRT corrected','Location','Best')
        %legend('10m','50m')
        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])
        ylim([0,109])
        xlim([jdb,jde])

        if pngflg
            print_buffer = [outpath m(ii).name(1:14) '_TRH.png'];
            print('-dpng ', print_buffer);
        end

        figure(2)
        clf
        subplot(2,1,1),plot(bao300(1:n,15),bao300(1:n,3))
        %hold
        %subplot(2,1,1),plot(data300(1:n,3),data300(1:n,7),'r')
        %xlabel('Year Day')
        ylabel('Spd (m/s)')
        stitle=sprintf('BAO %d 300m', data300(1,15));
        title(stitle)
        %axis([110,116,0,30])
        ylim([0,max(bao300(1:n,3))+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])
        ylim([0,360])

        if pngflg
            print_buffer = [outpath m(ii).name(1:14) '_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)')
        stitle=sprintf('BAO %d 300m', data300(1,15));
        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:14) '_T.png'];
            print('-dpng ', print_buffer);
        end

    end % if doing plots

end % loop through all files

