%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
read_parameters_Stratus
%w=what('S:\GOMECC_2007\RHB\radar\profiler\Raw\raw\');
w = way_raw_data_profiler;
jd=input('Input Jd for raw file to plot (ex: rbc07291.raw, enter 291): ', 's');
%cbfn = 'S:\GOMECC_2007\RHB\ceilometer\Processed\';
cbfn = way_proc_data_ceilo;
cbfn1=input('Input ceilo CB txt filename (ex: 291_10_18_07.txt) ', 's');
lf1 = length(way_proc_data_ceilo);
lf2 = length(cbfn1);
cbfn(lf1+1:lf1+lf2)=cbfn1;
fn = 'rbc107120.r*w'; % read .raw or .rew
fn(7:9) = jd;
a=dir([w fn]);
i=1;
% 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([w a(i).name],'r');
fgetl(fid);

jd=[];hr=[];mn=[];sec=[];
k=1;

beam = 9;  % set beam to capture
while ~feof(fid)
%for pp = 1:20  %short test loop DEW
    temp1=fgetl(fid);
    data=sscanf(temp1,'%d');
   bb = data(1);
   el = data(2);
   az = data(3);
    jd=[jd data(5)];
    hr=[hr data(6)];
    mn=[mn data(7)];
    sec=[sec data(8)];
    data(6);

    clear temp1 data;

    temp2=fgetl(fid);
    data=sscanf(temp2,'%d');
    pw = data(2);

    gate1=data(3); % in m
    gate_space=data(4);
    last_gate=data(5);

    clear temp2 data;

    for j=1:last_gate

        temp=fgetl(fid);
        data=sscanf(temp,'%f');
%        if beam == bb
% hardwires SPdeS doesn't understand
%        if el == 90 & pw == 700 % DEW 5/17/2007 changed from beam test to EL and PW
        if el == 90 && pw == 400 % elevation angle, a
            if j == 1
                k = k + 1;
                ht=(gate1+(0:1:last_gate-1)*gate_space)/1000; % in km
            end
            snr(k,j)=data(2);
            %pow(k,j)=data(3);
        end  %  End test for beam 9 vertical high mode
        clear temp data;
    end
    %k=k+1;
end

fclose(fid);
clear fid;

time=hr*3600+mn*60+sec;
k;
clear hr mn sec k jd;
%Save correct ht scale for data
for j=1:last_gate
    snr_cr(:,j)=snr(:,j)+20*log10(ht(j));
    %pow_cr(:,j)=pow(:,j)+20*log10(height(j));
end

figure()
set(gcf,'defaultaxesfonts',14)
set(gca,'fontname','time','fontsize',14);
imagesc((time-time(1))/3600,ht,snr_cr',[-10 20])
axis xy;
hcb=colorbar; set(get(hcb,'ylabel'),'string','rcSNR','fontsize',14);
xlim([0 23])
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 ' fn(7:9)] 'Wind profiler rcSNR'};
title(tit)
print_buffer=[way_images_profiler,  cruise,'_',year,'_',fn(7:9)];
print('-djpeg90', print_buffer)
print('-depsc2',print_buffer)

cb=load(cbfn);  %Load Ceilometer cloud base txt file
% Ceilo cloud base file
%  

cb_time=cb(:,5)+((cb(:,6)+(cb(:,7)/60.))/60);
hold
plot(cb_time,cb(:,2)/1000,'k.')
axis([0,24,0,2.5])
tit(2) = {'Wind profiler rcSNR and ceilometer cloud base (.)'};
title(tit)

print_buffer=[way_images_profiler,'wind_profiler',  cruise,'_',year,'_',fn(7:9),'cb'];
print('-djpeg', print_buffer)
print('-depsc2', print_buffer)
