%%BAOPLT10100300
% Plot daily combined data for levels 10, 100 and 300 meters.

clear
close all

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

doplot = 1;  % do the plot 1=yes 0=no
pngflg = 0;  % save data to png's flag  1= yes 0 = no

dd = clock;

year = dd(1);
bd = datenum(dd(1),dd(2),dd(3),0,0,0)-datenum(dd(1)-1,12,31);
yms = num2str(dd(1));
pd = bd - 1;
pds = num2str(pd);
bds = sprintf('%03d',bd);

%% Get previous days data
name10p =   ['Z:\bao\Tower\Processed\daily\BAO_010_' yms pds '.dat'];
name100p =  ['Z:\bao\Tower\Processed\daily\BAO_100_' yms pds '.dat'];
name300p =  ['Z:\bao\Tower\Processed\daily\BAO_300_' yms pds '.dat'];

    fid10p = fopen( name10p );
    if( fid10p > 0 )
        fclose( fid10p );
    else
        fprintf(1,'No  10m: %s\n',name10p);
        continue;       
    end
    
    fid100p = fopen( name100p );
    if( fid100p > 0 )
        have100p = 1;
        fclose( fid100p );
    else
        have100p = 0;
        fprintf(1,'No 100m: %s\n',name100p);
    end

    fid300p = fopen( name300p );
    if( fid300p > 0 )
        fclose( fid300p );
    else
        fprintf(1,'No 300m: %s\n',name300p);
        continue;
    end
    
    data10p = load(name10p);

    if have100p
        data100p = load(name100p);
    end

    fprintf(1,' 10m: %s\n',name10p);
    if have100p
        fprintf(1,'100m: %s\n',name100p);
    end
    fprintf(1,'300m: %s\n',name300p);

%     nmax = 2880;
%     yearp = zeros(1,nmax);
%     monthp = zeros(1,nmax);
%     dayp = zeros(1,nmax);
%     hourp = zeros(1,nmax);
%     minutep = zeros(1,nmax);
%     secp = zeros(1,nmax);
%     ydayp = zeros(1,nmax);
%     data300p = zeros( nmax, 15);
% 
    f300p = fopen(name300p,'r');
    %load f:\BAO_CR10\data\bao10.dat
    %read in 4 header lines
    tmpp=fgetl(f300p);
    tmpp=fgetl(f300p);
    tmpp=fgetl(f300p);
    tmpp=fgetl(f300p);
    n = 0;
    while ~feof(f300p)
        %for x = 1:10000
        n = n + 1;
        tempp=fgetl(f300p);
        yearp(n)= str2double(tempp(2:5));
        monthp(n)= str2double(tempp(7:8));
        dayp(n)= str2double(tempp(10:11));
        hourp(n)= str2double(tempp(13:14));
        minutep(n)= str2double(tempp(16:17));
        secp(n)= str2double(tempp(19:20));
        ltp= length(tempp);
        % strip off date/time info from first part of each line
        s1p = tempp(23:ltp);
        % convert each , to a blank space
        s1p = strrep( s1p, ',', ' ' );

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

    end
    data300p(:,15) = datenum(yearp,monthp,dayp,hourp,minutep,secp)-datenum(yearp-1,12,31);

    hhmmp = data10p(:,4);
    hhp = fix(hhmmp/100);
    mnp = hhmmp - (hhp*100);
    data10p(:,3) = data10p(:,3) + (((mnp/60.) + hhp)/24); %Julian day
    clear hhmm hh mm

    if have100p
        hhmmp = data100p(:,4);
        hhp = fix(hhmmp/100);
        mnp = hhmmp - (hhp*100);
        data100p(:,3) = data100p(:,3) + (((mnp/60.) + hhp)/24); %Julian day
        clear hhmm hh mm
    end

    %10m prop vane was on the north boom until Aug 3 ~1600 UTC when it was moved to
    %the south boom
    % N is 154  to correct we add 206 and correct for > 360
    ld = length(data10p);
    jdb = fix(data10p(1,3));
    jde = fix(data10p(ld,3))+1;


%%


name10 =   ['Z:\bao\Tower\Processed\daily\BAO_010_' yms bds '.dat'];
name100 =  ['Z:\bao\Tower\Processed\daily\BAO_100_' yms bds '.dat'];
name300 =  ['Z:\bao\Tower\Processed\daily\BAO_300_' yms bds '.dat'];


    %%%%%%%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
    %10,2007,165,1406,15.115,50.448,1.692,-999.000,245.47,-999.0000,13.514,842.98,0.000,219.470

    %%%%%%%100m 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
    %100,2007,244,5,29.848,19.414,0,0,0,0,13.132
    %100,2007,244,6,29.849,19.333,0,0,0,0,13.132


    %%%%300m data
    %"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"
    %"TS","RN","mV","Deg","Deg","mV","mV","","Deg C","Deg C","%","Counts","mV","Volts"
    %"","","WVc","WVc","WVc","Std","Avg","Avg","Avg","Avg","Avg","Avg","Smp","Min"
    %"2007-10-01 00:00:00",312159,5.06,201.1,8.43,0.523,15.68,106.923,17,17.91,23.68,0,110.7,13.26,17.8,163.1
    %"2007-10-01 00:00:30",312160,5.008,183.5,3.734,0.294,15.49,106.9227,17,17.9,25.13,0,110.7,13.26,17.8,145.5

    fid10 = fopen( name10 );
    if( fid10 > 0 )
        fclose( fid10 );
    else
        fprintf(1,'No  10m: %s\n',name10);
        continue;       
    end
    
    fid100 = fopen( name100 );
    if( fid100 > 0 )
        have100 = 1;
        fclose( fid100 );
    else
        have100 = 0;
        fprintf(1,'No 100m: %s\n',name100);
    end

    fid300 = fopen( name300 );
    if( fid300 > 0 )
        fclose( fid300 );
    else
        fprintf(1,'No 300m: %s\n',name300);
        continue;
    end
    
    data10 = load(name10);

    if have100
        data100 = load(name100);
    end

    fprintf(1,' 10m: %s\n',name10);
    if have100
        fprintf(1,'100m: %s\n',name100);
    end
    fprintf(1,'300m: %s\n',name300);

%     nmax = 2880;
%     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);
%     data300 = zeros( nmax, 15);

    f300 = fopen(name300,'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);
    n = 0;
    while ~feof(f300)
        %for x = 1:10000
        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') ;
        data300(n,1) = data(2,1); % WS ms Sonic
        data300(n,2) = data(3,1); % WD deg Sonic
        data300(n,3) = data(4,1); % WD  std deg Sonic
        data300(n,4) = data(5,1); % WS  STD Sonic
        data300(n,5) = data(6,1); % Temp deg C Sonic
        data300(n,6) = data(7,1); % PRT temp resistance
        data300(n,7) = data(8,1); % PRT temp deg C
        data300(n,8) = data(9,1); % HMT Temp Deg C
        data300(n,9) = data(10,1); % RH
        data300(n,10) = data(11,1); % airflow
        data300(n,11) = data(12,1); % VAC
        data300(n,12) = data(13,1); % Bat VDC
        data300(n,13) = data(14,1); % PRT corrected deg C
        data300(n,14) = data(15,1); % Dir true

    end
    data300(:,15) = datenum(year,month,day,hour,minute,sec)-datenum(year-1,12,31);

    hhmm = data10(:,4);
    hh = fix(hhmm/100);
    mn = hhmm - (hh*100);
    data10(:,3) = data10(:,3) + (((mn/60.) + hh)/24); %Julian day
    clear hhmm hh mm

    if have100
        hhmm = data100(:,4);
        hh = fix(hhmm/100);
        mn = hhmm - (hh*100);
        data100(:,3) = data100(:,3) + (((mn/60.) + hh)/24); %Julian day
        clear hhmm hh mm
    end

    %10m prop vane was on the north boom until Aug 3 ~1600 UTC when it was moved to
    %the south boom
    % N is 154  to correct we add 206 and correct for > 360
    ld = length(data10);
    jdb = fix(data10(1,3));
    jde = fix(data10(ld,3))+1;

%    if have100
%        ld2 = length(data100);
%        jdb2 = fix(data100(1,3));
%        jde2 = fix(data100(ld2,3))+1;
%    end
    %CSD experiment
    %jdb = 208;
    %jde = 222;
    
    % combine 2 days and plot out only last 24 hours
    ln10 = length(data10);
    ln10p = length(data10p);
    ln100 = length(data100);
    ln100p = length(data100p);
    ln300 = length(data300);
    ln300p = length(data300p);
    ep = data300(ln300,15); %find last year day in file
    bp = ep - 1;
    bao10 = data10p;
    bao100 = data100p;
    bao300 = data300p;
    
    l1 = ln10p+1;
    l2 = (ln10p+ln10);
    bao10(l1:l2,:)  = data10(1:ln10,:);
    
    l1 = ln100p+1;
    l2 = (ln100p+ln100);
    bao100(l1:l2,:)  = data100(1:ln100,:);
    
    l1 = ln300p+1;
    l2 = (ln300p+ln300);
    bao300(l1:l2,:)  = data300(1:ln300,:); 
    
    yrs = num2str(year(1));
    ms = sprintf('%02d',(month(1)));
    dsb = num2str(dayp(1));
    dse = num2str(day(ln300));
    tit = ['BAO ' yrs '-' ms ' ' dsb '-' dse];

    if doplot %skip plotting

        figure(1)
        clf
        subplot(2,1,1),plot(bao10(:,3),bao10(:,7))
        hold
        %subplot(2,1,1),plot(data100(:,3),data100(:,7),'k')
        subplot(2,1,1),plot(bao300(:,15),bao300(:,1),'r')
        title(tit)
        ylabel('Speed (m/s)')
        xlim([bp,ep])        
        %xlim([jdb,jde])
        elim = max(data300(:,1))+5;
        ylim([0 elim])
        legend('10m','300m','Location','Best','Orientation','horizontal');
        %legend('boxoff');

        %axis([110,116,820,850])
        subplot(2,1,2),plot(bao10(:,3),bao10(:,14))
        hold
        %subplot(2,1,1),plot(data100(:,3),data100(:,14),'k')
        subplot(2,1,2),plot(bao300(:,15),bao300(:,14),'r')
        xlabel('Year Day (UTC)')
        ylabel('Dir (True)')
        xlim([bp,ep])        
        ylim([0 360])
        %axis([110,116,0,.5])
        if pngflg
            print_buffer = [outpath 'BAO_10300_' ym is '_WS.png'];
            print('-dpng ', print_buffer);
        end

        if have100
            figure(2)
            clf
            subplot(2,1,1),plot(bao10(:,3),bao10(:,5))
            hold
            subplot(2,1,1),plot(bao100(:,3),bao100(:,5),'k')
            subplot(2,1,1),plot(bao300(:,15),bao300(:,13),'r')
            title(tit)
            ylabel('Temp (C)')
            xlim([bp,ep])        
            blim = min(data10(:,5))-10;
            elim = max(data10(:,5))+10;
            ylim([blim elim])
            legend('10m','100m','300m','Location','Best','Orientation','horizontal');
            %legend('boxoff');
            %legend('10m','100m','300m','Location','SouthEast');

            subplot(2,1,2),plot(bao10(:,3),bao10(:,6))
            hold
            subplot(2,1,2),plot(bao100(:,3),bao100(:,6),'k')
            subplot(2,1,2),plot(bao300(:,15),bao300(:,9),'r')
            xlabel('Year Day (UTC)')
            ylabel('RH (%)')
            xlim([bp,ep])        
            ylim([0 105])
            if pngflg
                print_buffer = [outpath 'BAO_10100300_' ym is '_TRH.png'];
                print('-dpng ', print_buffer);
            end
        end % if have 100 data
    end % skip plotting