% Program = 'main_plot_P3_raw_data'
% updated = '2-Aug-2017'

% This routine reads P3 tail doppler radar

% **********************************************************************

%% Load the data
load C:\data\cwf\Dropbox\matlabstf\cwf\progs\wband_calwater\Feb5a
cd  c:\data\cwf\Dropbox\DOC\CLIMATE\CalWater2\aircraft\Wband_precip\TDR
load('p3_tail_aft_and_fore_below_aircraft_2015_02_05.mat');
% The variables are:
%   p3_tail_aft_altitude                   1295x1                  10360  double
%   p3_tail_aft_azimuth                    1295x1                  10360  double
%   p3_tail_aft_drift                      1295x1                  10360  double
%   p3_tail_aft_elevation                  1295x1                  10360  double
%   p3_tail_aft_fixed_angle                1295x1                  10360  double
%   p3_tail_aft_georefs_applied            1295x1                  10360  double
%   p3_tail_aft_heading                    1295x1                  10360  double
%   p3_tail_aft_latitude                   1295x1                  10360  double
%   p3_tail_aft_longitude                  1295x1                  10360  double
%   p3_tail_aft_norm_coh_pow               1295x1000            10360000  double
%   p3_tail_aft_nyquist_velocity           1295x1                  10360  double
%   p3_tail_aft_pitch                      1295x1                  10360  double
%   p3_tail_aft_prt                        1295x1                  10360  double
%   p3_tail_aft_prt_ratio                  1295x1                  10360  double
%   p3_tail_aft_pulse_width                1295x1                  10360  double
%   p3_tail_aft_radar_beam_width_h         1295x1                  10360  double
%   p3_tail_aft_radar_beam_width_v         1295x1                  10360  double
%   p3_tail_aft_range                      1295x1000            10360000  double
%   p3_tail_aft_reflectivity               1295x1000            10360000  double
%   p3_tail_aft_roll                       1295x1                  10360  double
%   p3_tail_aft_rotation                   1295x1                  10360  double
%   p3_tail_aft_spectrum_width             1295x1000            10360000  double
%   p3_tail_aft_tilt                       1295x1                  10360  double
%   p3_tail_aft_time_of_ray                1295x1                  10360  double
%   p3_tail_aft_time_start_scan            1295x6                  62160  double
%   p3_tail_aft_total_power                1295x1000            10360000  double
%   p3_tail_aft_unambiguous_range          1295x1                  10360  double
%   p3_tail_aft_velocity                   1295x1000            10360000  double
%   p3_tail_fore_aft_altitude              1295x1                  10360  double
%   p3_tail_fore_aft_azimuth               1295x1                  10360  double
%   p3_tail_fore_aft_drift                 1295x1                  10360  double
%   p3_tail_fore_aft_elevation             1295x1                  10360  double
%   p3_tail_fore_aft_georefs_applied       1295x1                  10360  double
%   p3_tail_fore_aft_heading               1295x1                  10360  double
%   p3_tail_fore_aft_latitude              1295x1                  10360  double
%   p3_tail_fore_aft_longitude             1295x1                  10360  double
%   p3_tail_fore_aft_nyquist_velocity      1295x1                  10360  double
%   p3_tail_fore_aft_pitch                 1295x1                  10360  double
%   p3_tail_fore_aft_prt                   1295x1                  10360  double
%   p3_tail_fore_aft_prt_ratio             1295x1                  10360  double
%   p3_tail_fore_altitude                     1x1295               10360  double
%   p3_tail_fore_azimuth                      1x1295               10360  double
%   p3_tail_fore_drift                        1x1295               10360  double
%   p3_tail_fore_elevation                    1x1295               10360  double
%   p3_tail_fore_fixed_angle               1295x1                  10360  double
%   p3_tail_fore_georefs_applied              1x1295               10360  double
%   p3_tail_fore_heading                      1x1295               10360  double
%   p3_tail_fore_latitude                     1x1295               10360  double
%   p3_tail_fore_longitude                    1x1295               10360  double
%   p3_tail_fore_norm_coh_pow              1295x1000            10360000  double
%   p3_tail_fore_nyquist_velocity             1x1295               10360  double
%   p3_tail_fore_pitch                        1x1295               10360  double
%   p3_tail_fore_prt                          1x1295               10360  double
%   p3_tail_fore_prt_ratio                    1x1295               10360  double
%   p3_tail_fore_pulse_width               1295x1                  10360  double
%   p3_tail_fore_radar_beam_width_h        1295x1                  10360  double
%   p3_tail_fore_radar_beam_width_v        1295x1                  10360  double
%   p3_tail_fore_range                     1295x1000            10360000  double
%   p3_tail_fore_reflectivity              1295x1000            10360000  double
%   p3_tail_fore_roll                      1295x1                  10360  double
%   p3_tail_fore_rotation                  1295x1                  10360  double
%   p3_tail_fore_spectrum_width            1295x1000            10360000  double
%   p3_tail_fore_tilt                      1295x1                  10360  double
%   p3_tail_fore_time_of_ray               1295x1                  10360  double
%   p3_tail_fore_time_start_scan           1295x6                  62160  double
%   p3_tail_fore_total_power               1295x1000            10360000  double
%   p3_tail_fore_unambiguous_range         1295x1                  10360  double
%   p3_tail_fore_velocity                  1295x1000            10360000  double

%%  Process Aft beam

%% Construct a time array hr.fraction - aft beam

% Get the time of the start of the scan
time_aft_hrmin  = p3_tail_aft_time_start_scan(:,4) + p3_tail_aft_time_start_scan(:,5)./60 + ... 
   p3_tail_aft_time_start_scan(:,6)./(60*60);

% Add the seconds of the ray that is at rotation angle = 180
time_aft_hrmin  = time_aft_hrmin + p3_tail_aft_time_of_ray./(60*60);

%%  Fig 1. Plot lat, lon, & altitude
% 
% figure
% subplot(3,1,1)
% plot(time_aft_hrmin,p3_tail_aft_latitude,'k')
% hold on
% title('P3 TDR, Aft beam: Latitude [deg]')
% 
% subplot(3,1,2)
% plot(time_aft_hrmin,p3_tail_aft_longitude,'k')
% hold on
% title('P3 TDR, Aft beam: Longitude [deg]')
% 
% subplot(3,1,3)
% plot(time_aft_hrmin,p3_tail_aft_altitude,'k')
% hold on
% title('P3 TDR, Aft beam: Altitude [m]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig1_aft_aircraft_latlonalt.tif')
% close
% 
% %%    Fig 2. Plot pitch, roll, heading
% 
% figure
% subplot(3,1,1)
% plot(time_aft_hrmin,p3_tail_aft_pitch,'k')
% hold on
% title('P3 TDR, Aft beam: pitch [deg]')
% 
% subplot(3,1,2)
% plot(time_aft_hrmin,p3_tail_aft_roll,'k')
% hold on
% title('P3 TDR, Aft beam: roll [deg]')
% 
% subplot(3,1,3)
% plot(time_aft_hrmin,p3_tail_aft_heading,'k')
% hold on
% title('P3 TDR, Aft beam: heading [deg]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig2_aft_aircraft_pitch_roll_heading.tif')
% close
% 
% %%    Fig 3. Plot elevation, azimuth, and rotation
% 
% figure
% subplot(3,1,1)
% plot(time_aft_hrmin,p3_tail_aft_elevation,'k')
% hold on
% title('P3 TDR, Aft beam: Beam Pointing Elevation Angle [deg]')
% 
% subplot(3,1,2)
% plot(time_aft_hrmin,p3_tail_aft_azimuth,'k')
% hold on
% title('P3 TDR, Aft beam: Beam Pointing Azimuth Angle [deg]')
% 
% subplot(3,1,3)
% plot(time_aft_hrmin,p3_tail_aft_rotation,'k')
% hold on
% title('P3 TDR, Aft beam: Beam Pointing Rotation Angle [deg]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig3_aft_el_az_rotation.tif')
% close
% 
% %%    Fig 4. Plot Time-height Reflectivity
% 
% figure
% colormap('jet')
% 
% start_hour  = 19;
% end_hour    = 22;
% start_range =  0;
% end_range   =  5;
% 
% subplot(3,1,1)
% pcolor(time_aft_hrmin,p3_tail_aft_range(1,:)./1000,p3_tail_aft_reflectivity')
% shading flat
% colorbar
% hold on
% axis([start_hour end_hour start_range end_range])
% set(gca,'xtick',start_hour:0.5:end_hour);
% set(gca,'ytick',start_range:1:end_range);
% title('P3 TDR: Aft Beam, Reflectivity [dBZ]')
% ylabel('Range [km]')
% 
% subplot(3,1,2)
% pcolor(time_aft_hrmin,p3_tail_aft_range(1,:)./1000,p3_tail_aft_velocity')
% shading flat
% colorbar
% hold on
% axis([start_hour end_hour start_range end_range])
% set(gca,'xtick',start_hour:0.5:end_hour);
% set(gca,'ytick',start_range:1:end_range);
% title('P3 TDR: Aft Beam, Radial Velocity [m/s]')
% ylabel('Range [km]')
% 
% subplot(3,1,3)
% pcolor(time_aft_hrmin,p3_tail_aft_range(1,:)./1000,p3_tail_aft_spectrum_width')
% shading flat
% colorbar
% hold on
% axis([start_hour end_hour start_range end_range])
% set(gca,'xtick',start_hour:0.5:end_hour);
% set(gca,'ytick',start_range:1:end_range);
% title('P3 TDR: Aft Beam, Spectrum Width [m/s]')
% ylabel('Range [km]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig4_aft_time_range_ZVW.tif')
% close
% 
% %% Process fore beam

%% Construct a time array hr.fraction - fore beam

% Get the time of the start of the scan
time_fore_hrmin  = p3_tail_fore_time_start_scan(:,4) + p3_tail_fore_time_start_scan(:,5)./60 + ... 
   p3_tail_fore_time_start_scan(:,6)./(60*60);

% Add the seconds of the ray that is at rotation angle = 180
time_fore_hrmin  = time_fore_hrmin + p3_tail_fore_time_of_ray./(60*60);

%%  Fig 5. Plot lat, lon, & altitude

% figure
% subplot(3,1,1)
% plot(time_fore_hrmin,p3_tail_fore_latitude,'k')
% hold on
% title('P3 TDR, Fore beam: Latitude [deg]')
% 
% subplot(3,1,2)
% plot(time_fore_hrmin,p3_tail_fore_longitude,'k')
% hold on
% title('P3 TDR, Fore beam: Longitude [deg]')
% 
% subplot(3,1,3)
% plot(time_fore_hrmin,p3_tail_fore_altitude,'k')
% hold on
% title('P3 TDR, Fore beam: Altitude [deg]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig5_fore_aircraft_latlonalt.tif')
% close
% 
% %%    Fig 6. Plot pitch, roll, heading
% 
% figure
% subplot(3,1,1)
% plot(time_fore_hrmin,p3_tail_fore_pitch,'k')
% hold on
% title('P3 TDR, Fore beam: pitch [deg]')
% 
% subplot(3,1,2)
% plot(time_fore_hrmin,p3_tail_fore_roll,'k')
% hold on
% title('P3 TDR, Fore beam: roll [deg]')
% 
% subplot(3,1,3)
% plot(time_fore_hrmin,p3_tail_fore_heading,'k')
% hold on
% title('P3 TDR, Fore beam: heading [deg]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig6_fore_aircraft_pitch_roll_heading.tif')
% close
% 
% %%    Fig 7. Plot elevation, azimuth, and rotation
% 
% figure
% subplot(3,1,1)
% plot(time_fore_hrmin,p3_tail_fore_elevation,'k')
% hold on
% title('P3 TDR, Fore beam: Beam Pointing Elevation Angle [deg]')
% 
% subplot(3,1,2)
% plot(time_fore_hrmin,p3_tail_fore_azimuth,'k')
% hold on
% title('P3 TDR, Fore beam: Beam Pointing Azimuth Angle [deg]')
% 
% subplot(3,1,3)
% plot(time_fore_hrmin,p3_tail_fore_rotation,'k')
% hold on
% title('P3 TDR, Fore beam: Beam Pointing Rotation Angle [deg]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig7_fore_el_az_rotation.tif')
% close
% 
% %%    Fig 8. Plot Time-height Reflectivity
% 
% figure
% colormap('jet')
% 
% start_hour  = 19;
% end_hour    = 22;
% start_range =  0;
% end_range   =  2.8;
% 
% subplot(3,1,1)
% pcolor(time_fore_hrmin,p3_tail_fore_range(1,:)./1000,p3_tail_fore_reflectivity')
% shading flat
% colorbar
% hold on
% axis([start_hour end_hour start_range end_range])
% set(gca,'xtick',start_hour:0.5:end_hour);
% set(gca,'ytick',start_range:1:end_range);
% title('P3 TDR: Fore Beam, Reflectivity [dBZ]')
% ylabel('Range [km]')
% 
% subplot(3,1,2)
% pcolor(time_fore_hrmin,p3_tail_fore_range(1,:)./1000,p3_tail_fore_velocity')
% shading flat
% colorbar
% hold on
% axis([start_hour end_hour start_range end_range])
% set(gca,'xtick',start_hour:0.5:end_hour);
% set(gca,'ytick',start_range:1:end_range);
% title('P3 TDR: Fore Beam, Radial Velocity [m/s]')
% ylabel('Range [km]')
% 
% subplot(3,1,3)
% pcolor(time_fore_hrmin,p3_tail_fore_range(1,:)./1000,p3_tail_fore_spectrum_width')
% shading flat
% colorbar
% hold on
% axis([start_hour end_hour start_range end_range])
% set(gca,'xtick',start_hour:0.5:end_hour);
% set(gca,'ytick',start_range:1:end_range);
% title('P3 TDR: Fore Beam, Spectrum Width [m/s]')
% ylabel('Range [km]')
% xlabel('Time (UTC) of 5-Feb-2015')
% 
% print('-dtiff','fig8_fore_time_range_ZVW.tif')
% close
%%%%%%%%   add 3 dB to reflectivity
p3_tail_aft_reflectivity=p3_tail_aft_reflectivity+5;
p3_tail_fore_reflectivity=p3_tail_fore_reflectivity+5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p3_tail_fore_elevation=p3_tail_fore_elevation';p3_tail_fore_altitude=p3_tail_fore_altitude';

range_a=p3_tail_aft_altitude./sin((-p3_tail_aft_elevation)*pi/180);%range to the surface
range_f=p3_tail_fore_altitude./sin((-p3_tail_fore_elevation)*pi/180);
%%%%%%%%%%%%%%%%%   find reflectivity and Doppler at the surface
irange=[];ref1=[];ref2=[];dopsa=[];dospf=[];
for i=1:1295;
    ii=find(p3_tail_aft_range(i,:)>range_a(i));
    irange_a(i)=ii(1);
    ref1(i)=max(p3_tail_aft_reflectivity(i,ii(1)-1:ii(1)+1));
    kk=1;%kk=find(ref1(i)==p3_tail_aft_reflectivity(i,ii(1)-1:ii(1)+1));
    dopsa(i)=p3_tail_aft_velocity(i,ii(1)+kk-1);
   ii=find(p3_tail_fore_range(i,:)>range_f(i));
     irange_f(i)=ii(1);
    ref2(i)=max(p3_tail_fore_reflectivity(i,ii(1)-1:ii(1)+1));
   kk=1;% kk=find(ref2(i)==p3_tail_fore_reflectivity(i,ii(1)-1:ii(1)+1));
    dopsf(i)=p3_tail_fore_velocity(i,ii(1)+kk-1);
end;
%%%%%%%   rayleigh parameters at Xband
ksq=0.93;lam=0.0322;dr=75;
nrcsc=10*log10(pi^5*ksq/lam^4*dr);%NRCS constant at Xband
ii=find(abs(p3_tail_aft_roll)<1 & time_aft_hrmin>19 & time_aft_hrmin<22);figure;plot(time_aft_hrmin(ii),ref1(ii)-180+nrcsc,time_aft_hrmin(ii),ref2(ii)-180+nrcsc);xlim([19 22]);xlabel('Hour'); ylabel('Xband NRCS');
%%%%%%%%%%%%   plot xband nrcs at 20 zenith angle
%%%%%%%%%%%  sv2 = TRMM 9.4 GHz
%nrcf_wu
sv2=0.009+0.0276*log10(max(1,Us3m));
ii=find(Us3m>7);sv2(ii)=-0.084+0.138*log10(Us3m(ii));
svw=sv2;
nrcf_wu=10*log10(0.317./sv2);
%nrcf_cos
sv2=0.003+5.08e-3*max(1,Us3m);
svc=sv2;
nrcf_cox=10*log10(0.317./sv2);
%nrcf_TRMM
sv2=0.0036+0.028*log10(max(1,Us3m));
ii=find(Us3m>10);sv2(ii)=-0.0184+0.05*log10(Us3m(ii));
svt=sv2;
nrcf_trmm=10*log10(0.317./sv2);
ii=find(isnan(Us3m));nrcf_trmm(ii)=NaN;

ij=find(isnan(Us3m));sv2(ij)=NaN;
figure;plot(hhm,nrcf_trmm,'o',time_aft_hrmin,ref1-180+nrcsc,'.',hhm,10*log10(.8*.63*exp(-tan(20*pi/180)^2./sv2)./sv2./cos(20*pi/180)^4),'x');ylim([-20 10]);xlabel('Hour'); ylabel('NRCS');legend('NRCS W nadir TRMM','NRCS X 20 deg','NRCS X TRMM 20 deg');;

%%%%%%%%%%%%%%%   compute rainrate from Z-R
Za=10.^(p3_tail_aft_reflectivity/10);
Zf=10.^(p3_tail_fore_reflectivity/10);
a=210;b=1.47;%lee
a=300;b=1.40;%steiner
a=367;b=1.3;%tokay/short
a=275;b=1.49;%Iguchi 
a=297;b=1.42;%average
Ra=(nanmean(Za(:,6:20)')/a).^(1/b);
Rf=(nanmean(Zf(:,6:20)')/a).^(1/b);
Rr=(Ra+Rf)/2;%  combine aft and fore


%convert rangegates to altitude gates
htr_a=[];htr_f=[];
for i=1:length(dopsa)
htr_a(i,:)=p3_tail_aft_range(10,:).*sin(-p3_tail_aft_elevation(i,:)*pi/180);
htr_f(i,:)=p3_tail_fore_range(10,:).*sin(-p3_tail_fore_elevation(i,:)*pi/180);
end;


%%%%%%%  interpolate Xband range gates onto Wband altitude gates
ij=find(time_aft_hrmin>19 & time_aft_hrmin<22);%xband data in hrs 19 thru 21 correspoding to Wband interval
nj=length(ij);
ik=find(time_fore_hrmin>19 & time_aft_hrmin<22);
ik=[ik' 981]';%don't ask

%%%%%%%%%%%%%%%   interpolate Xband reflect and Doppl onto Wband height
%%%%%%%%%%%%%%%   gates
Ref_a=[];Ref_f=[];alt_a=[];hh_a=[];Dop_a=[];Dop_f=[];Dop_as=[];Dop_fs=[];Dop_as2=[];Dop_fs2=[];rollexa=[];rollexf=[];Dop_w=[];
for i=1:nj
    ka=find(~isnan(htr_a(ij(i),:)));kf=find(~isnan(htr_f(ik(i),:)));%time window
    Ref_a(i,:)=interp1(htr_a(ij(i),ka),p3_tail_aft_reflectivity(ij(i),ka),httot*1e3);
    Ref_f(i,:)=interp1(htr_f(ik(i),kf),p3_tail_fore_reflectivity(ik(i),kf),httot*1e3);
    Dop_a(i,:)=interp1(htr_a(ij(i),ka),p3_tail_aft_velocity(ij(i),ka),httot*1e3);
    Dop_f(i,:)=interp1(htr_f(ik(i),kf),p3_tail_fore_velocity(ik(i),kf),httot*1e3);
    alt_a(i)=p3_tail_aft_altitude(ij(i));
    hh_a(i)=time_aft_hrmin(ij(i));
    hh_f(i)=time_aft_hrmin(ik(i));
    Dop_as(i)=dopsa(ij(i))';
    Dop_fs(i)=dopsf(ij(i))';
    rollexa(i)=p3_tail_aft_roll(ij(i));
    rollexf(i)=p3_tail_fore_roll(ik(i));
    %%%%%%%%%%%%%%%  use ground speed from Wband files
    lk=find(hhm>time_aft_hrmin(ij(i)));lkx=min(length(sog3m),lk(1):lk+37);Dop_as2(i)=nanmean(sog3m(lkx)).*cos(pi/180*(p3_tail_aft_elevation(ij(i))));
    soga(i)=nanmean(sog3m(lkx));
    lk=find(hhm>time_fore_hrmin(ik(i)));lkx=min(length(sog3m),lk(1):lk+37);Dop_fs2(i)=nanmean(sog3m(lkx)).*cos(pi/180*(p3_tail_fore_elevation(ik(i))));
    sogf(i)=nanmean(sog3m(lkx));
    %%%%%%%%%%%%%%   combine aft and fore doppler to get estimate of
    %%%%%%%%%%%%%%   vertical component
    agh=(sin(pi/180*(p3_tail_aft_elevation(ij(i))))+sin(pi/180*(p3_tail_fore_elevation(ik(i)))));
    Dop_w(i,:)=-(Dop_a(i,:)+Dop_f(i,:))./agh+(soga(i).*cos(pi/180*(p3_tail_aft_elevation(ij(i))))-sogf(i).*cos(pi/180*(p3_tail_fore_elevation(ik(i)))))./agh;
 end;
 %***************   correct for velocity folding
Dop_as2=Dop_as2-2*p3_tail_aft_nyquist_velocity(1);
Dop_fs2=-(Dop_fs2-2*p3_tail_fore_nyquist_velocity(1));


%%%%%%%%%%%%%   average to 1 minute time resolution
[av bv]=size(Ref_a);Refa_min=[];Reff_min=[];Dopa_min=[];Dopf_min=[];Dopw_min=[];jdhh_a=36+hh_a/24;jdhh_f=36+hh_f/24;
ma=find(abs(rollexa)<2);  %don't use data where aircraft is rolled more than 2 deg
for i=1:bv
    [jdava,refmin1,sig]=timave_1_nonan(jdhh_a(ma),1,Ref_a(ma,i));
    mm=find(isnan(refmin1));refmin1(mm)=-25;
    Refa_min(:,i)=refmin1;
    [jdava,refmin1,sig]=timave_1_nonan(jdhh_f(ma),1,Ref_f(ma,i));
    mm=find(isnan(refmin1));refmin1(mm)=-30;
    Reff_min(:,i)=refmin1;    
   [jdava,refmin1,sig]=timave_1_nonan(jdhh_a(ma),1,Dop_a(ma,i));
    Dopa_min(:,i)=refmin1;
    [jdava,refmin1,sig]=timave_1_nonan(jdhh_f(ma),1,Dop_f(ma,i));
     Dopf_min(:,i)=refmin1;    
    [jdava,refmin1,sig]=timave_1_nonan(jdhh_f(ma),1,Dop_w(ma,i));
     Dopw_min(:,i)=refmin1;    
     
 end;
 %%%%%%%   surface Doppler
 [jdava,Dopas_min,sig]=timave_1_nonan(jdhh_a(ma),1,Dop_as(ma)');
 [jdavf,Dopfs_min,sig]=timave_1_nonan(jdhh_f(ma),1,Dop_fs(ma)');
 [jdava,Dopas2_min,sig]=timave_1_nonan(jdhh_a(ma),1,Dop_as2(ma)');
 [jdavf,Dopfs2_min,sig]=timave_1_nonan(jdhh_f(ma),1,Dop_fs2(ma)');
 [jdavf,Rollexa_min,sig]=timave_1_nonan(jdhh_f,1,rollexa');
 ii=find(Dopfs2_min<-16.9);Dopfs2_min(ii)=Dopfs2_min(ii)+2*16.9;%velocity folding
[jdavf,Rr_min,sig]=timave_1_nonan(jdhh_f,1,Rr(ij)');
 
hh_va=(jdava-36)*24;%convert jd to hrs
hh_vf=(jdavf-36)*24;


kk=25:60;  %select vertical range  25=2 km, 60=.93 km
vb=nanmean(Dopw_min(:,kk)');%% mean Xband vertical doppler
ii=find(vb<1);vb(ii)=NaN;
Zb=10.^(nanmean(Refa_min(:,kk)'/10));%mean aft reflect
figure;plot(10*log10(Zb),vb,'o');ylim([0 15]);title(['Altitude= 0.9 - 1.5 km']);xlabel('dBZe X band');ylabel('Vertical Doppler Xband (m/s)');mu=0;

%%  Rayleigh formulae
mu=0;
bu=1/gamma(7+mu)*(3.78*gamma(7.67+mu)/gamma(7+mu)).^((21+3*mu)/2);
au=6*pi/1e4*(3.78*gamma(7.67+mu)/gamma(7+mu)).^4.5*gamma(4.67+mu)/gamma(7.67+mu);
lb=(gamma(7.67+mu)*3.78/gamma(7+mu)./vb).^1.5;
D0=(3.67+mu)./lb;%rayleight D0
figure;plot(10*log10(Zb),D0,'o');ylim([0 3]);title(['Altitude= 0.9 - 1.5 km']);xlabel('dBZe X band');ylabel('Gamma D0 (mm)');
Rrx=au*Zb./vb.^3.5;%Rayleigh rainrate from Z and V
figure;loglog(Rrx,Rr_min,'.',[1e-3 1e2],[1e-3 1e2],'--');xlabel('Rain Gamma (mm/h)');ylabel('Rain Z-R (mm/h)');


kk=30;ik=find(hh_va>19 & hh_va<22);figure;plot(Refa_min(ik,kk:kk+30)',dBZe2t(kk:kk+30,ik),'.');xlabel('Xband dbZe');ylabel('Wband dBZe');

Zam=10.^((Refa_min)/10);
Zfm=10.^((Reff_min)/10);
Ram=(nanmean(Zam(:,40:60)')/297).^(1/1.42);
Rfm=(nanmean(Zfm(:,40:60)')/297).^(1/1.42);
Rrm=(Ram+Rfm)/2;
Refr_min=(Refa_min+Reff_min)/2;

Refa_min0=Refa_min(1201,:)-1;
Zmina0=10.^(Refa_min0/10);
Zefa_min=10.^(Refa_min/10);
Reff_min0=Reff_min(1201,:)-1;
Zminf0=10.^(Reff_min0/10);
Zeff_min=10.^(Reff_min/10);

kk=30;ik=find(hh_va>19 & hh_va<22);figure;plot(Refa_min(ik,kk:kk+30)'+5,dBZe2t(kk:kk+30,ik),'.',[-30 40],[-30 40]);xlabel('Xband dbZe');ylabel('Wband dBZe');grid;
ii=find(hh_va>20.9);figure;plot(Refmc_min(:,ii(1)),rng-ht,dBZe2t(:,ii(1)),rng-httot,10*log10(Zefa_min(ii(1),:))+7,rng-httot,'.',10*log10(Zeff_min(ii(1),:))+5,rng-httot,'o');ylim([0 2.5]);title(num2str(hh_va(ii(1))));xlim([-35 40]);
ik=find(hh_va>19 & hh_va<22);
kk=25:60;
Dlv=DVamc_min(kk-20,ik);
Snlv=SNRm_min(kk-20,ik);
Dlvx=Dopw_min(ik,kk);
for i=1:length(Dlv)
ii=find(Snlv(:,i)<-5);
Dlv(ii,i)=NaN;
end;

figure;
subplot(2,1,1);plot(-Dlv,Refa_min(ik,kk)','.');ylabel('Xband dBZe');xlim([-1 6]);ylim([-20 40]);grid;
subplot(2,1,2);plot(-Dlv,dBZe2t(kk,ik),'.');xlabel('Wband -Doppler (m/s)');ylabel('Wband dBZe');ylim([-20 40]);xlim([-1 6]);grid;
figure;
subplot(2,1,1);plot(Dlvx',Refa_min(ik,kk)','.');ylabel('Xband dBZe');xlim([-1 10]);ylim([-20 40]);grid;
subplot(2,1,2);plot(Dlvx',dBZe2t(kk,ik),'.');xlabel('Xband -Doppler (m/s)');ylabel('Wband dBZe');ylim([-20 40]);xlim([-1 10]);grid;

[.2 .6 1.5 3.1 5 15];
DDz=[-1 .25 .7 1.5 3 6 20]; 
DDz=[-1 0.1 .7 1.5 3 6 20]; 

%ii=find(rain_ncmin>DDz(i) & rain_ncmin<DDz(i+1));
i=1;ii=find(rain_za>DDz(i) & rain_za<DDz(i+1));
rmn(i)=nanmean(rain_za(ii));D0m(i)=nanmean(D0(ii))
% ii=find(Rr_min>DDz(i) & Rr_min<DDz(i+1));
% rmn=nanmean(Rr_min(ii));

figure;
subplot(1,2,1);plot(nanmean(dBZe2t(:,ii)'),rng-httot,'.',nanmean(Refa_min(ii,:)),rng-httot,nanmean(Reff_min(ii,:)),rng-httot);ylim([0 2]);xlim([-20 40]);title(['Rainrate= ' num2str(rmn(i))]);grid;
legend('dBZeW_c','dBZeX_a','dBZeX_f');
subplot(1,2,2);plot(nanmean(-DVamc_min(:,ii)'),rng-ht','.',nanmean(Dopw_min(ii,:)),rng-httot','.');ylim([0 2]);xlim([0 8]);grid;legend('Doppler_W','Doppler_X')

DVwx=DVamc_min'+Dopw_min(:,20:169);DVwx=max(DVwx,0);
dBZewx=(Refa_min(:,20:169)+Reff_min(:,20:169))/2-3.2*DVwx.^1.4;

i=6;ii=find(rain_za>DDz(i) & rain_za<DDz(i+1));rmn(i)=nanmean(rain_za(ii));D0m(i)=nanmean(D0(ii))
figure;plot(nanmean(Refmc_min(:,ii)'),rng-ht,'.',nanmean(dBZe2t(20:169,ii)'),rng-ht,'.',nanmean(dBZewx(ii,:)),rng-ht);ylim([0 2]);xlim([-20 40]);title(['Rainrate= ' num2str(rmn(i))]);grid;
xlabel('dBZe');ylabel('Altitude (km)');legend('dBZeW','dBZeW_c','dBZeX-3.2*DV^{1.4}');

figure;plot([-12 11 17 22 29 35],[-14 8 15 17 19 22],'x',[-12 11 17 22 29 35],[-14 8 15 17 19 22]+[.6 2 6 7 12 15],'o',[-12 11 17 22 29 35],[-12 11 17 22 29 35],'--');xlabel('Xband dBZe');ylabel('Wband dBZE');legend('dBZeW_c','dBZeW_c+\DeltadBZW');
figure;plot(10.^([11 14 17 19 22]/10),[.5 1 2 4 7.5],'o',10.^([11 14 17 19 22]/10),.048*10.^(1.0*[11 14 17 19 22]/10));xlabel('Ze mm^6/m^3');ylabel('Attenuation (dB/km)');title('Wband \gamma vs Ze using Xband reference');
kl=find(hh_va>19 & hh_va<22);
kk=25:60;slopx=[];intx=[];slopw=[];intw=[];
for i= 1:length(kl)
    [P,S]=polyfit(httot(kk)',Refr_min(kl(i),kk),1);
     slopx(i)=P(1);intx(i)=P(2);
    [P,S]=polyfit(httot(kk),Refmc_min(kk-19,kl(i)),1);
     slopw(i)=P(1);intw(i)=P(2);   
    
end;
       x=Dlvx'+Dlv;
       y=Refr_min(ik,kk)'-dBZe2t(kk,ik);
       binsize=0.5;xmin=0;xmax=4;
[X, Y, M, SD, SE, SDM, SEM] = binave_3(x,y,binsize,xmin,xmax);
figure;plot(X,Y+0,'-o',X,3.2*X.^1.4,'--');xlabel('DoppX-DoppW (m/s)');ylabel('dBZeX-dBZeWc');legend('Observations','3*DV^{1.4}');grid;



