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

function [retval] = UMPSD10_loop_func( by, bm, bd, ey, em, ed )

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

pinpath = 'z:\bao\Tower\Processed\daily\';
uinpath = 'z:\bao\Tower\UMASS\Processed\';
outpath  = 'z:\bao\Tower\UMASS\Processed_Images\';

if( nargin == 0 )
    bys = input('Input year to process (2007) ','s');
    bms = input('Input the month to process (1-12) ','s');
    bds = input('Input the beginning day to process (1-31) ','s');
    eys = input('Input year to process (2007) ','s');
    ems = input('Input the month to process (1-12) ','s');
    eds = input('Input end day to process (1-31) ','s');
    by = str2double( bys );
    bm = str2double( bms );
    bd = str2double( bds );
    ey = str2double( eys );
    em = str2double( ems );
    ed = str2double( eds );
end

bdn = datenum([by bm bd 0 0 0]);
edn = datenum([ey em ed 0 0 0]);

% 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

for i = bdn:edn
    dv=datevec(i);
    yday = md2yd( dv(1), dv(2), dv(3) );
    fnpsd10 = sprintf('%sBAO_010_%4d%03d.dat', pinpath, dv(1), yday );
    fnumass10 = sprintf('%ssonic_010n_%4d%02d%02d.dat',uinpath,dv(1),dv(2),dv(3));

    fid10 = fopen( fnpsd10 );
    if( fid10 > 0 )
        fclose( fid10 );
    else
        fprintf(1,'No PSD 10m: %s\n',fnpsd10);
        continue;
    end

    fid10 = fopen( fnumass10 );
    if( fid10 > 0 )
        fclose( fid10 );
    else
        fprintf(1,'No UMASS 10m: %s\n',fnumass10);
        continue;
    end

    [pathp, namep, extp, versnp] = fileparts(fnpsd10);
    [pathu, nameu, extu, versnu] = fileparts(fnumass10);

    fprintf(1,'*** Processing %s%s and %s%s ***\n', namep,extp, nameu,extu );

    p10 = load(fnpsd10);
    hhmm = p10(:,4);
    hh = fix(hhmm/100);
    mm = hhmm - (hh*100);
    p10(:,3) = p10(:,3) + (((mm/60.) + hh)/24); %decimal year day
    clear hhmm hh mm

    u10 = load(fnumass10);
    uyd10 = datenum(u10(:,1),u10(:,2),u10(:,3),u10(:,4),u10(:,5),u10(:,6))-datenum(u10(:,1)-1,12,31);
    us10=sqrt(u10(:,7).*u10(:,7) + u10(:,8).*u10(:,8));
    ud10=(atan2(u10(:,7),u10(:,8))*180./pi) - 26.;
    indx = find(ud10 < 0. & ud10 >= -180.);
    if length(indx)
        ud10(indx) = 360. + ud10(indx);
    end

    %  Do plots if requested
    
    if pltflg
        %  Plot Temperatures
        figure('Position',[10 750 650 650]);
        plot(uyd10,u10(:,10))
        hold
        plot(p10(:,3),p10(:,5),'r')
        legend('UMass','PSD','Location','Best','Orientation','horizontal');
        legend('boxoff');
        ylabel('T (C)')
        xlabel('Year Day')
        mx = max(uyd10);
        xlim([uyd10(1) mx])
        blims(1) = min(u10(:,10))-5;
        elims(1) = max(u10(:,10))+5;
        blims(2) = min(p10(:,5))-5;
        elims(2) = max(p10(:,5))+5;
        ylim([min(blims) max(elims)])
        ylim([min(p10(:,5))-5 max(p10(:,5)+5)])
        stitle = sprintf('BAO UMASS vs PSD North 10m %4d-%02d-%02d', dv(1),dv(2),dv(3));
        title(stitle);
        if pngflg
            pngpath = sprintf('%sUMPSD_%4d%02d%02d_10_T.png',outpath,dv(1),dv(2),dv(3));
            f=getframe(gcf);
            imwrite(f.cdata,pngpath,'png');
        end

        figure('Position',[10 750 650 650])
        subplot(2,1,1),plot(uyd10,us10)
        hold
        subplot(2,1,1),plot(p10(:,3),p10(:,7),'r.')
        legend('UMass','PSD','Location','Best','Orientation','horizontal');
        legend('boxoff');
        ylabel('Spd (m/s)')
        xlabel('Year Day')
        xlim([uyd10(1) mx])
        blims(1) = 0; %min(us10)-5;
        elims(1) = max(us10)+5;
        blims(2) = 0; %min(p10(:,7))-5;
        elims(2) = max(p10(:,7))+5;
        ylim([min(blims) max(elims)])
        title(stitle);

        subplot(2,1,2),plot(uyd10,ud10)
        hold
        subplot(2,1,2),plot(p10(:,3),p10(:,14),'r.')
        ylabel('Direction')
        xlabel('Year Day')
        xlim([uyd10(1) mx])
        ylim([0 360])
        if pngflg
            pngpath = sprintf('%sUMPSD_%4d%02d%02d_10_WS.png',outpath,dv(1),dv(2),dv(3));
            f=getframe(gcf);
            imwrite(f.cdata,pngpath,'png');
        end

        figure('Position',[10 750 650 650])
        plot(uyd10,u10(:,7:9))
        legend('u','v','w','Location','Best','Orientation','horizontal');
        legend('boxoff');
        ylabel('Spd (m/s)')
        xlabel('Year Day')
        xlim([uyd10(1) mx])
        blims = min(u10(:,7:9))-5;
        elims = max(u10(:,7:9))+5;
        ylim([min(blims) max(elims)])
        title(stitle);
        if pngflg
            pngpath = sprintf('%sUMPSD_%4d%02d%02d_10_UVW.png',outpath,dv(1),dv(2),dv(3));
            f=getframe(gcf);
            imwrite(f.cdata,pngpath,'png');
        end

        figure('Position',[10 750 650 650])
        plot(uyd10,u10(:,11),'.')
        stitle = sprintf('BAO UMASS North 10m sonic npts %4d-%02d-%02d', dv(1),dv(2),dv(3));
        title(stitle);
        ylim([0 34]);
        if pngflg
            pngpath = sprintf('%sUMPSD_%4d%02d%02d_10_NPTS.png',outpath,dv(1),dv(2),dv(3));
            f=getframe(gcf);
            imwrite(f.cdata,pngpath,'png');
        end

    end %if plots are wanted

    clear p10 u10 uyd10 us10 ud10 indx
    close all
    
end  %loop daily files

retval = 0;
