file=uipickfiles('FilterSpec','F:/C17FSR/0813-70/*.hpl');
nfiles=length(file);
for j=1:nfiles
    [datapath, filename, ext, versn] = fileparts(char(file(j)));
    if ~exist([datapath,'/analysis'],'dir')
        mkdir([datapath,'/analysis'])
    end
    writedir(j)=cellstr([datapath,'/analysis']);
    files(j)=cellstr(filename)
end
fid=fopen(char(file(1)));
    lidar=read_Proc_profiles(fid,75.0);
    fclose(fid);
    [naz,nr]=size(lidar.RWS);
    nn=naz-1;
    [na,nr]=size(lidar.azimuth);
    rws=lidar.RWS(1:nn,:);
    snr=lidar.SNR(1:nn,:);
    beta=lidar.Beta(1:nn,:);
    xd=lidar.xd(1:nn,:);
    yd=lidar.yd(1:nn,:);
    zd=lidar.zd(1:nn,:);
    azimuth=lidar.azimuth(1:nn,:);
    elevation=lidar.elevation(1:nn,:);
    lidar_decimal_time=lidar.lidar_decimal_time(1:nn);
    [sd sm sd sH sM sS]=datevec(lidar.start_time);
    for j=1:na-1;
        beamtime(j,1)=datenum([sd sm sd lidar.lidar_decimal_time(j) 0 0]);
    end
    ranges=lidar.ranges;
for jfiles=2:nfiles
        fid=fopen(char(file(jfiles)));
        lidar=read_Proc_profiles(fid,0.0);
        fclose(fid);
        [sd sm sd sH sM sS]=datevec(lidar.start_time);
        for j=1:na-1;
            beam(j,1)=datenum([sd sm sd lidar.lidar_decimal_time(j) 0 0]);
        end
        rws=cat(1,rws,lidar.RWS(1:nn,:));
        snr=cat(1,snr,lidar.SNR(1:nn,:));
        beta=cat(1,beta,lidar.Beta(1:nn,:));
        xd=cat(1,xd,lidar.xd(1:nn,:));
        yd=cat(1,yd,lidar.yd(1:nn,:));
        zd=cat(1,zd,lidar.zd(1:nn,:));
        elevation=cat(1,elevation,lidar.elevation(1:nn,:));
        azimuth=cat(1,azimuth,lidar.azimuth(1:nn,:));
        beamtime=cat(1,beamtime,beam);
        lidar_decimal_time=cat(1,lidar_decimal_time,lidar.lidar_decimal_time(1:nn));
end
[na,nr]=size(rws);
vad=lidar;
for j=1:nr-3
    rwsa(:,j)=mean(rws(:,j:j+3),2);
    snra(:,j)=mean(snr(:,j:j+3),2);
    betaa(:,j)=mean(beta(:,j:j+3),2);
    xda(:,j)=mean(xd(:,j:j+3),2);
    yda(:,j)=mean(yd(:,j:j+3),2);
    zda(:,j)=mean(zd(:,j:j+3),2);
    elevationa(:,j)=mean(elevation(:,j:j+3),2);
    azimutha(:,j)=mean(azimuth(:,j:j+3),2);
end
    
vad.RWS=rwsa;
vad.Beta=betaa;
vad.SNR=snra;
vad.beamtime=beamtime;
vad.lidar_decimal_time=lidar_decimal_time;
vad.elevation=elevationa;
vad.azimuth=azimutha;
vad.xd=xda;
vad.yd=yda;
vad.zd=zda;

    SNR_thresh=1.005;
    [na,nr]=size(vad.RWS);
    for k=1:na
        for j=1:nr-5
            p(k,j)=std(vad.RWS(k,j:j+5));
        end
    end
    [indax,inday]=find(vad.SNR<SNR_thresh);
    [indbx,indby]=find(p>0.5);

    numble=na*nr;
    snr=reshape(vad.SNR,na,nr);
    rws=reshape(vad.RWS,na,nr);
    for j=1:length(indax)
       snr(indax(j),inday(j))=NaN;
       rws(indax(j),inday(j))=NaN;
    end
    for j=1:length(indbx)
       snr(indbx(j),indby(j))=NaN;
       rws(indbx(j),indby(j))=NaN;
    end
    %snr=reshape(snr,na,nr);
    %rws=reshape(rws,na,nr);
    for i=1:10:na-48
        vad2.ranges=vad.ranges;
       vad2.SNR=snr(i:i+nn-1,:);
       vad2.RWS=rws(i:i+nn-1,:);
       vad2.elevation=vad.elevation(i:i+nn-1,:);
       vad2.azimuth=vad.azimuth(i:i+nn-1,:);
       vad2.Beta=vad.Beta(i:i+nn-1,:);
       vad2.xd=vad.xd(i:i+nn-1,:);
       vad2.yd=vad.yd(i:i+nn-1,:);
       vad2.zd=vad.zd(i:i+nn-1,:);
       %profile=vvpanalysis_3_olda(vad2);
       profile=vvp_regression(vad2,2);
       [year month day a b c]=datevec(lidar.start_time);
       profile.time=datenum([year month day mean(vad.lidar_decimal_time(i:i+11)) 0 0]);
       filename=['User1_70_',datestr(profile.time,'yyyymmdd_HHMMSS'),'_profile.mat'];
    
    %textfile=fullfile(char(writedir(jfiles)),[char(files(jfiles)),'_profile.dat']);
    %prtid=print_file(textfile,profile,vad);
       textfile=fullfile(char(writedir(jfiles)),filename);
       save(textfile,'profile');
       ind=find(profile.turb'./(profile.Vh+0.01)>0.7);
       Vhor=profile.Vh;
       %Vhor(ind)=NaN;
      
       plot(Vhor,profile.heights)
       set(gca,'XLim',[0 25])
       set(gca,'YLim',[0 700])
       pause(1)
    
    end
