files=uipickfiles('FilterSpec','E:\My Data\Projects\DRI-BP\data\*.mat');
nfiles=length(files);
writedir=cellstr('E:\My Data\Projects\DRI-BP\Analysis');
for j=1:nfiles
    [datapath, filename, ext, versn] = fileparts(char(files(j)));
    load(char(files(j)));
    
    % remove 1st 16 lines of data form window contamination - adjust if
    % needed
    
    vad=lidar;
    vad.ranges=lidar.ranges(17:end);
    vad.number_of_gates=lidar.number_of_gates-18;
    vad.RWS=vad.RWS(1:16,17:end);
    vad.SNR=vad.SNR(1:16,17:end);
    vad.Beta=vad.Beta(1:16,17:end);
    vad.xd=vad.xd(1:16,17:end);
    vad.yd=vad.yd(1:16,17:end);
    vad.zd=vad.zd(1:16,17:end);
    vad.elevation=vad.elevation(1:16,17:end);
    vad.azimuth=vad.azimuth(1:16,17:end);
    [na,nr]=size(vad.RWS);
    for k=1:na
        locs=peakseek(vad.SNR(k,:),100,2);
        if ~isempty(locs)
           lay(k).locs=vad.zd(k,locs);
        else
            lay(k).locs=NaN;
        end
    end
    SNR_thresh=1.004;
    max_hvel=15;
%     [inda,indb]=size(vad.SNR);
%     numble=inda*indb;
%     snr=reshape(vad.SNR,numble,1);
%     rws=reshape(vad.RWS,numble,1);
%     zd=reshape(vad.zd,numble,1);
%     ele=reshape(vad.elevation,numble,1);
    ind=find(vad.SNR<SNR_thresh | abs(vad.RWS./cosd(vad.elevation))>max_hvel );
    vad.SNR(ind)=NaN;
    vad.RWS(ind)=NaN;
%     snr=reshape(snr,inda,indb);
%     rws=reshape(rws,inda,indb);
%     zd=reshape(zd,inda,indb);
%     vad.SNR=snr;
%     vad.RWS=rws;
    profile=vvp_regression(vad);
    profile.layer=lay;
    hfig=figure('Position',[100 100 900 500],'MenuBar','none','PaperOrientation','landscape')
  
    subplot(1,4,1);    
    plot(profile.Vh,profile.heights);
    maxx=(floor(max(profile.Vh)/5)+1)*5;
    set(gca,'XLim',[0 maxx])
    set(gca,'YLim',[0 1000])
    xlabel('Horizontal Velocity (m/s)')
    ylabel('Height AGL (m)')
    subplot(1,4,2);
    plot(profile.theta,profile.heights);
    set(gca,'XLim',[0 360])
    set(gca,'YLim',[0 1000])
    xlabel(gca,'Direction')
    set(gca,'XTick',[0 90 180 270 360])
    set(gca,'YTickLabel',[])
    subplot(1,4,3);
    plot(profile.w,profile.heights);
    maxx=(floor(max(abs(profile.w))/2)+1)*2;
    set(gca,'XLim',[-maxx maxx])
    set(gca,'YLim',[0 1000])
    xlabel(gca,'Vertical Velocity (m/s)')
    set(gca,'YTickLabel',[])
    subplot(1,4,4);
    plot(profile.turb,profile.heights);
    maxx=(floor(max(abs(profile.turb))/5)+1)*5;
    set(gca,'XLim',[0 maxx])
    set(gca,'YLim',[0 1000])
    xlabel(gca,'Residual (m/s)')
    set(gca,'YTickLabel',[])
    tt=datestr(lidar.start_time,'yyyy-mmm-dd HH:MM:SS');
    text(0,1025,tt,'HorizontalAlignment','center')
    pause(1)
    
     saveas(hfig,fullfile(char(writedir),[filename,'_profile.png']));
     close(hfig)
     textfile=fullfile(char(writedir),[filename,'_profile.dat']);
     
     prtid=print_file(textfile,profile,vad);
     textfile=fullfile(char(writedir),[filename,'_profile.mat']);
     save(textfile,'profile');
end
    
    