%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute boundary layer top from vertical-beam profiler (60 m gates)
% from the height of maximum range-corrected signal-to-noise ratio (SNR) 
% above 800 m. Output 10-minutely 16-84th percentile-averaged daylong
% boundary layer height files.
%
% Simon de Szoeke  18 December 2007
% Adapted from profiler915_BLHT.m by Dan Wolfe (need to use Dan's version
% for different radar modes, i.e. multiple elevation angles)
% and from pick_mbl_snrcr.pro by Bill Otto.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
read_parameters_Stratus
%jd=input('Input Jd for raw file to plot (ex: rbc07291.raw, enter 291): ', 's');

for yday=340:355%:355; % -357? loop year day of files to read
%yday=340;
    fn=sprintf('rbc04%3d.r*w',yday); % read .raw or .rew
    adir=dir([way_raw_data_profiler fn]);

    cbfn = way_proc_data_ceilo;
    datev=datevec(datenum(2004,0,0,0,0,0)+yday);
    cbfn1=sprintf('%03d_%02d_%02d_%02d.txt',yday,datev(3:-1:2),datev(1)-2000);
    lf1 = length(way_proc_data_ceilo);
    lf2 = length(cbfn1);
    cbfn(lf1+1:lf1+lf2)=cbfn1;

    % rbclear                           0     0    13 91500
    %5  75  45  7 123  0  0 26
    %61  700 196 101  55 110 30  128
    %  -4.308  13.240  62.462   3.114  28.150 0
    %  -4.622  12.690  51.212   4.098  17.450 0

    fid=fopen([way_raw_data_profiler adir.name],'r');
    %frewind(fid);
    fgetl(fid);

    jd=[];hr=[];mn=[];sec=[];snr=[];width=[];
    k=0; % start with no scans, increment k for desired (vertical) scans
    while ~feof(fid) % loop over scans
        k=k+1;
        %for pp = 1:20  %short test loop DEW
        temp1=fgetl(fid);
        data=sscanf(temp1,'%d');
        %bb(k) = data(1); % beam index--arbitrary number
        el(k)=data(2); % elevation angle, degrees
        az(k)=data(3); % azimuth angle, degrees
        %yr(k)=data(4); % years since 1900
        jd(k)=data(5); % append each scan
        hr(k)=data(6);
        mn(k)=data(7);
        sec(k)=data(8);
        %clear temp1 data;

        temp2=fgetl(fid);
        data=sscanf(temp2,'%d');
        ipp(k)=data(1); % inter-pulse period (ms)
        pw(k)=data(2); % pulse width (ns)
        gate1(k)=data(3); % in m
        gate_space(k)=data(4);
        last_gate(k)=data(5);
        %clear temp2 data;

        for j=1:last_gate(k) % read a line for each gate
            temp=fgetl(fid);
            data=sscanf(temp,'%f');
            % DEW 5/17/2007 changed from beam test to EL(elevation angle) and PW
            %if el==90 && pw==400 % Simon eliminated conditions
                if j==1
                    %k=k+1;
                    ht(k,1:last_gate(k))=(gate1(k)+(0:1:last_gate(k)-1)*gate_space(k))/1000; % km % store gate heights for each scan
                end
                snr(k,j)=data(2);
                %pow(k,j)=data(3);
            %end
            %clear temp data;
        end
        %if el==90 && pw==400
        %    snr_cr=snr+20*log10(ht);
        %end
    end
    fclose(fid);
    
    % select only 400 ns pulses, not 700 ns
    ipulse=(pw(1:k)==400);
    iz=[true([1 52]) false([1 3])];
    snr_cr=NaN+zeros(sum(ipulse),sum(iz));
    snr_cr=snr(ipulse,iz)+20*log10(ht(ipulse,iz));
    time=hr(ipulse)'*3600+mn(ipulse)'*60+sec(ipulse)';
    
    % find BL top from max snr_cr
    M=spatialfilter(snr_cr,[1 2 3 2 1],[1 2 1]); % filtered SNR
    M(isnan(M))=snr_cr(isnan(M));

    cb=load(cbfn);  %Load Ceilometer cloud base txt file
    cb_time=cb(:,5)+((cb(:,6)+(cb(:,7)/60.))/60);
    [junk,ind]=unique(cb_time);
    iis=isfinite(cb(:,2)) & sparse(ind,1,ones(length(ind),1),length(cb),1);
    cbase=interp1(cb_time(iis)*3600,cb(iis,2)/1e3,time);
    %cbase(isnan(cbase)|cbase>1.500)=0.500; % default level to look above=500m
    
    % look above 800m for inversion height, except for when
    % cloud base is low, look above ceilometer cloud base
    h=ht(ipulse,iz);
    ii=(1:size(h,2))';
    ihmax=NaN+zeros(size(snr_cr,1),1); ihmaxtest=ihmax;
    bltop=ihmax; bltoptest=ihmax;
    for itime=1:size(snr_cr,1)
        % search for BL top echo above sea/ship clutter
        if cbase(itime)<1.000 % cloud base is defined and lower than 1 km.
            mask=h(itime,:)>cbase(itime) & h(itime,:)>0.500 & h(itime,:)<1.500;
            isrch=find(h(itime,:)>cbase(itime) & h(itime,:)>0.500,1,'first');
        else % Btw, reverts to this if cbase is NaN.
            mask=h(itime,:)>0.800 & h(itime,:)<2.500;
            isrch=find((diff(M(itime,:)).*mask(1:end-1))>0,1,'first'); % look above min (filtered) SNR
        end
        minin(itime)=h(itime,isrch);
        mask1=ii<=isrch; % exclude below SNR min. to exclude clutter
        [junk,ihmaxtest(itime)]=max(M(itime,:)-1000*(mask1'|~mask)); %filtered top
        [snrmax(itime),ihmax(itime)]=max(snr_cr(itime,:)-1000*(mask1'|~mask));
        bltoptest(itime)=h(itime,ihmaxtest(itime));
        bltop(itime)=h(itime,ihmax(itime));
    end
    bltop(abs(bltop-bltoptest)>.080)=NaN; % toss when filtered top disagrees by more than one gate

    % Form the mean every 10 minutes of the 16-84th percentile.
    minbase=(0:10:24*60)'; % marks beginning of time bins (last bin is empty of data)
    secbase=minbase*60;

    bltop10=NaN+zeros(length(minbase)-1,1);
    for tbin=1:length(minbase)-1
        iin=time>=secbase(tbin) & time<secbase(tbin+1) & isfinite(bltop);
        nobs(tbin)=sum(iin);
        if nobs(tbin)>=4; % 4 needed to toss away outliers
            ind=0.5:1:nobs(tbin);
            irange=ind>=0.16*nobs(tbin) & ind<=0.84*nobs(tbin);
            topsort=sort(bltop(iin));
            bltop10(tbin)=mean(topsort(irange));
        else
            bltop10(tbin)=NaN;
        end
    end

    % save 10-min BL top data daily
    ftop=fopen([way_proc_data_profiler cruise year '_top_' sprintf('%3d',yday) '.txt'],'w');
    fprintf(ftop,'%3d %5d %6.04f\n',[repmat(yday,[length(bltop10),1]) secbase(1:end-1,1) bltop10]');
    fclose(ftop);
    
    % plot to see how BL height looks
    %figure()
    clf
    set(gcf,'defaultaxesfonts',14)
    set(gca,'fontname','time','fontsize',14);
    [c,hcf]=contourf((repmat(time,[1 52]))/3600,h,snr_cr,16); set(hcf,'edgecolor','none');
    axis xy;
    hcb=colorbar; set(get(hcb,'ylabel'),'string','rcSNR','fontsize',14);
    xlim([0 24])
    ylim([0 2.5]) %ylim([0 4.29])
    set(gca,'xtick',0:4:24)
    ylabel('Height (km)')
    xlabel('Hour (UTC)')
    tit = {[cruise ' ' num2str(year) ' day ' num2str(jd(k))] 'Wind profiler rcSNR'};
    title(tit)
    hold on
    plot(secbase(1:end-1)/3600,bltop10,'m.')
    print_buffer=[way_images_profiler,  cruise,'_',year,'_',sprintf('%03d',yday)];
    print('-djpeg90', [print_buffer '.jpeg'])
    print('-depsc2',[print_buffer '.eps'])

    plot(cb_time,cb(:,2)/1000,'r.')
    axis([0,24,0,2.5])
    tit(2) = {'Wind profiler rcSNR, cloud top (magenta) and ceilometer cloud base (red)'};
    title(tit)
    print_buffer=[way_images_profiler,'wind_profiler_cloud_', cruise,year,'_',fn(7:9)];
    print('-djpeg', [print_buffer '.jpeg'])
    print('-depsc2', [print_buffer '.eps'])

end