%wsr=load([big_dir 'P3\WSRA\WSRAcat2.dat']);
WSRAheaders20150206;
mssw=mss;rrw=rr;
ssw=sod5;

hhw=ssw/3600;
msswp=interp1(hhw,mssw,hhm);
rrwp=interp1(hhw,rrw,hhm);
ii=find(abs(roll3m)>3);msswp(ii)=NaN;rrwp(ii)=NaN;

rainwsr=load('c:\data\cwf\DROPBOX\DOC\CLIMATE\CalWater2\aircraft\Wband_precip\WSRArr5feb2015x100.dat');
hodwsr=load('c:\data\cwf\DROPBOX\DOC\CLIMATE\CalWater2\aircraft\Wband_precip\WSRAhod5feb2015x1000.dat');
rainwsra=load('c:\data\cwf\DROPBOX\DOC\CLIMATE\CalWater2\aircraft\Wband_precip\WSRArr5feb2015x100a.dat');
k=length(rainwsra);rainwsr(1:k)=rainwsra;

im=1:length(rainwsr);dt=2/length(rainwsr);
hhws=hodwsr/1000;
rrwp2=interp1(hhws,rainwsr,hhm)/100;ii=find(rrwp2>20);rrwp2(ii)=NaN;


figure;
subplot(2,1,1);plot(hhw,mssw/1e4);ylabel('WSRA MSS');
subplot(2,1,2);plot(hhw,rrw/10);ylabel('WSRA Rainrate (mm/hr)');xlabel('Hour (GMT)');

kl=length(P3_rainm);
rng=nanmeanx(zzx);
rain_nc=-0.5*(nrcsm-nrcf)*1.3./sqrt(1.2)./P3_altrefm/1e-3;
%nrcf is wind speed model of nrcs; nrcsm is the measured value.  nrcsm=NaN
%if return from surface has SNR<-20
rain_nc=rain_nc-.0;ii=find(rain_nc<0);rain_nc(ii)=0;
rain_ncs=sgolayfilt(rain_nc,1,3);

htl=6;htu=25;kkk=htl:htu;hhz=0.0131:.025:ht(htl);cc=length(hhz);ac=150;;Gh=zeros(1,ac+cc-htl);Iv=zeros(1,ac+cc-htl);
httot=[hhz(1:cc-htl) ht']';
rhoa=1.21./exp((nanmeanx(zzx)-httot)/9);rhov=13./exp((nanmeanx(zzx)-httot)/2.5);%rhov=10./exp((nanmeanx(zzx)-httot)/1.0);
gammav=-0.25+1.3*rhov/20+0.25./sqrt(1+0.5*rhov);%water vapor absorption coeff

for i=2:169;
    Gh(i)=Gh(i-1)+2*gammav(i)*.025;
    Iv(i)=Iv(i-1)+rhov(i).*25;
end;
Iv=Iv/1e4;%cm
Refmc=Refm;
for i=1:kl
    Refmc(:,i)=Refm(:,i)+Gh(20:169)'+1.6+2;
end

%Refm=Refmc;

ii=find(ht<rng);ji=ii(length(ii))-3;ji=25;jk=70;
kkk=5:ji;kkk2=5:jk;
slop=[];intr=[];slop2=[];intr2=[];sig_p=[];
for i=1:kl
    ij=find(SNRm(kkk,i)>-10 );
        if ~isempty(ij) & Refmc(5,i)>0
        %[PP,SS]=polyfit(ht(kkk(ij)),Refmc(kkk(ij),i),1);
        [PP,SS]=polyfit(ht(kkk(ij)),Refmc(kkk(ij),i),1);
        slop(i)=PP(1);
        intr(i)=PP(2);
        profp=polyval(PP,ht(kkk(ij)));
        sig_p(i)=std(profp-Refmc (kkk(ij),i));
       else
        slop(i)=NaN;
        intr(i)=NaN;
        sig_p(i)=NaN;
        end;
    if abs(roll3m)>3;slop(i)=NaN;end;
    ij=find(SNRm(kkk2,i)>-10 );
        if ~isempty(ij) & Refmc(5,i)>0
        [PP2,SS]=polyfit(ht(kkk2(ij)),Refmc(kkk2(ij),i),1);
        slop2(i)=PP2(1);
        intr2(i)=PP2(2);
        else
        slop2(i)=NaN;
        intr2(i)=NaN;
    end;
    if abs(roll3m)>3;slop2(i)=NaN;end;
end;
ij=find(slop>0 | slop<-50);slop(ij)=NaN;
ij=find(intr<-30);intr(ij)=NaN;
slops=sgolayfilt(slop,1,11);
intrs=sgolayfilt(slop,1,11);


ij=find(slop2>0 | slop2<-50);slop2(ij)=NaN;
ij=find(intr2<-30);intr2(ij)=NaN;
intr2s=sgolayfilt(intr2,1,11);
slop2s=sgolayfilt(slop2,1,11);

P3_rains=sgolayfilt(P3_rainm,1,3);rain_all=rain_nc*NaN;
rrwp2s=sgolayfilt(rrwp2,1,5);
rain_all=nanmedianx([rain_ncs' P3_rains' rrwp2s']');
rain_2=nanmeanx([ P3_rains' rrwp2s']');
%mm=find(isnan(rain_ncs));rain_ncs(mm)=rain_2(mm);
jdhhm=hhm/24;
[jdav,rain_nrcmin,sigr]=timave_1_nonan(jdhhm,1,rain_nc');


DDz=[-15 -5 0  7.5 10 15];
DDz=[-1 1 3 6 12 25 50]/2; 
DDz=[-1 .25 .7 1.5 3 6 13]; 

bq=length(DDz);
Yp=zeros(150,bq-1);Dyp=zeros(149,bq-1);np=[];zp=[];snp=Yp;Vp=Yp;ypz=[];nq=[];fq=[];rainp=[];
for j=1:bq-1
ij=find(rain_ncs>DDz(j) & rain_ncs<DDz(j+1));
ii=find(rain_nc>DDz(j) & rain_nc<DDz(j+1) & SNRm(5,:)>-20 & Refmc(5,:)-Refmc(75,:)>-5 );
ii=find(rain_ncs>DDz(j) & rain_ncs<DDz(j+1) & SNRm(5,:)>-20 & Refmc(5,:)-Refmc(75,:)>-5 );
%ii=find(rain_2>DDz(j) & rain_2<DDz(j+1) & SNRm(5,:)>-20 & Refmc(5,:)-Refmc(75,:)>-5 );


rainp(j)=nanmeanx(rain_ncs(ii));
nrcsp(j)=nanmeanx(nrcsm(ii));
nrczp(j)=nanmeanx(10.^(nrcsm(ii)/10));
nrcfp(j)=nanmeanx(nrcf(ii));
nrczf(j)=nanmeanx(10.^(nrcf(ii)/10));
zp(j)=nanmeanx(10.^(Refmc(5,ii)/10));
np(j)=length(ii);nq(j)=length(ij);
fq(j)=np(j)/length(ij);
yp=nanmeanx(Refmc(:,ii)');sp=nanmeanx(SNRm(:,ii)');vp=nanmedianx(DVamc(:,ii)');ypz=nanmeanx(10.^(Refmc(:,ii)'/10));
Yp(:,j)=yp;snp(:,j)=sp;Vp(:,j)=vp;Ypz(:,j)=ypz;
dyp=diff(yp)./diff(ht');
Dyp(:,j)=dyp';
end;
rainpp=rainp;
dzp=10*log10(zp);
alt=nanmeanx(zzx)-ht;
rainm1=num2str(rainpp(1));rainm1=rainm1(1:3);
rainm2=num2str(rainpp(2));rainm2=rainm2(1:3);
rainm3=num2str(rainpp(3));rainm3=rainm3(1:3);
rainm4=num2str(rainpp(4));rainm4=rainm4(1:3);
rainm5=num2str(rainpp(5));rainm5=rainm5(1:3);
rainm6=num2str(rainpp(6));rainm6=rainm6(1:3);

figure;fg=4:75;
subplot(2,1,1);plot(Yp(fg,:),alt(fg));ylim([0 2.35]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');
subplot(2,1,2);plot(Vp(fg,:),alt(fg));ylim([0 2.35]);xlim([-5 0]);xlabel('Doppler (m/s)');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);


htl=8;htu=67;Dyps=[];
for j=1:bq-1
    h=sgolayfilt(Dyp(htl:htu,j),1,11);
    Dyps(:,j)=h;
end;
figure;plot(Dyps(:,2)*fq(2),alt(htl:htu,:),Dyps(:,3)*fq(3),alt(htl:htu,:),Dyps(:,4),alt(htl:htu,:),'o',Dyps(:,5),alt(htl:htu,:),'--',Dyps(:,6),alt(htl:htu,:),'.');ylim([0 2.]);xlim([-20 5]);xlabel('\Delta dBZe /\DeltaRng');ylabel('Range from aircraft (km)');legend(rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');
figure;plot(-.5*1.3*Dyps(:,2)./sqrt(rhoa(htl:htu))*fq(2),alt(htl:htu,:),-.5*1.3*Dyps(:,3)*fq(3)./sqrt(rhoa(htl:htu)),alt(htl:htu,:),-.5*1.3*Dyps(:,4)./sqrt(rhoa(htl:htu)),alt(htl:htu,:),'o',-.5*1.3*Dyps(:,5)./sqrt(rhoa(htl:htu)),alt(htl:htu,:),'--',-.5*1.3*Dyps(:,6)./sqrt(rhoa(htl:htu)),alt(htl:htu,:),'.');ylim([0 2.]);xlim([-2 12]);xlabel('Rain Rate (mm/hr)');ylabel('Altitude (km)');legend(rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');
for i=1:6
    Dypsf(:,i)=Dyps(:,i)*fq(i);
end;

[av bv]=size(Refmc);Refmc_min=[];
for i=1:av
    [jdav6,refmin1,sig]=timave_1_nonan(jdhhm,1,Refmc(i,:)');
    Refmc_min(i,:)=refmin1;
end;
rg=-0.5*(Refmc_min(5,:)-Refmc_min(70,:))/(alt(70)-alt(5))*1.3/sqrt(1.25);rain_dzmin=max(0,rg);
ii=find(24*jdav6>19 & 24*jdav6<22); m=length(ii);
slopmin=NaN*ones(1,1440);intrmin=NaN*ones(1,1440);
for i=1:m;
    [PP,SS]=polyfit(ht(kkk2),Refmc_min(kkk2,ii(i)),1);
        slopmin(ii(i))=PP(1);
        intrmin(ii(i))=PP(2);
end
rain_slopmin=max(0,-.5*slopmin*1.3./sqrt(1.25));
rain_dz=rain_nc*NaN;rain_slop=rain_nc*NaN;rain_slop=-.5*slops*1.3./sqrt(1.25);
ii=find(Refmc(5,:)>0 & Refmc(5,:)-Refmc(70,:)>0);rain_dz(ii)=-0.5*(Refmc(5,ii)-Refmc(70,ii))/(alt(70)-alt(5))*1.3/sqrt(1.25);
rain_dzs=sgolayfilt(rain_dz,1,11);rain_slops=sgolayfilt(rain_slop,1,11);
%figure;plot(hhw,10*log10(.13./(mssw)),hhm,nrcf,'.');xlabel('Hr (JD 036)');ylabel('NRCS');ylim([0 10]);legend('WSRA 0.13/MSS','SFMR Wind Formulla');
ij=find(intrmin<0);rain_slopmin(ij)=0;rain_dzmin(ij)=0;

Refmc_minx=max(Refmc_min(5:70,:));
%ij=find(Refmc_minx<0);rain_slopmin(ij)=0;

rain_mat=[rain_ncs' P3_rains' rrwp2s' rain_dzs' rain_slops'];
[jdav,rain_ncmin,sigr]=timave_1_nonan(jdhhm,1,rain_ncs');
%[jdav1,rain_dzmin,sigr]=timave_1_nonan(jdhhm,1,rain_dzs');
%[jdav2,rain_slopmin,sigr]=timave_1_nonan(jdhhm,1,rain_slops');
[jdav3,rain_P3min,sigr]=timave_1_nonan(jdhhm,1,P3_rains');
[jdav4,rain_WSRmin,sigr]=timave_1_nonan(jdhhm,1,rrwp2s');
figure;plot(jdav*24,rain_ncmin,jdav*24,rain_P3min,jdav*24,rain_WSRmin);xlim([19 22]);xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');legend('Z_{integral}','NRCS');
legend('NRCS','SFMR','WSRA');

figure;plot(jdav*24,rain_ncmin,jdav*24,rain_slopmin,jdav*24,rain_dzmin);xlim([19 22]);xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');ylim([0 14]);
legend('NRCS','Slope Method1','\DeltadBZ/ \Deltaz Method2');

figure;plot(jdav*24,rain_ncmin,jdav*24,rain_dzmin);xlim([19 22]);xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');ylim([0 14]);
legend('NRCS','Slope Method1','\DeltadBZ/ \Deltaz Method2');


rain2_min=nanmeanx([rain_P3min rain_WSRmin]');

ii=find(~isnan(rain_ncmin));corrcoef(rain_ncmin(ii),rain2_min(ii))
ii=find(~isnan(rain_WSRmin));corrcoef(rain_WSRmin(ii),rain2_min(ii))
ii=find(~isnan(rain_P3min));corrcoef(rain_P3min(ii),rain2_min(ii))
ii=find(~isnan(rain_dzmin) & ~isnan(rain2_min));corrcoef(rain_dzmin(ii),rain2_min(ii))
ii=find(~isnan(rain_slopmin) & ~isnan(rain2_min));corrcoef(rain_slopmin(ii),rain2_min(ii))



Ypm=[];
htl=6;htu=25;kkk=htl:htu;hhz=0.0131:.025:ht(htl);cc=length(hhz);[ac,bc]=size(Yp);Ypm=zeros(ac+cc-htl,bc);
for i=1:bq-1
     [PP,SS]=polyfit(ht(kkk),Yp(kkk,i),1);
        slopm(i)=PP(1);
        intrm(i)=min(30,PP(2));
        Ypm(1:cc,i)=polyval(PP,hhz);
        Ypm(cc+1:ac+cc-htl,i)=Yp(htl+1:ac,i);
end;
httot=[hhz(1:cc-htl) ht']';

Zintrm=10.^(intrm/10);
Zpm=10.^(Ypm/10);
for i=1:6
    Zpmq(:,i)=Zpm(:,i).*fq(i);
end;
Ypmq=10*log10(Zpmq);
figure;
subplot(2,1,1);plot(Ypm(1:20,:),rng-httot(1:20),'.',Ypm(21:95,:),rng-httot(21:95));ylim([0 2.5]);xlim([-25 20]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');
subplot(2,1,2);plot(Ypmq(1:20,:),rng-httot(1:20),'.',Ypmq(21:95,:),rng-httot(21:95));ylim([0 2.5]);xlim([-25 20]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');

alp=.033;bet=0.97;%matrosov
alp=.058;bet=.85;%PIP
alp=.05;bet=1.;

Gampm=alp*Zpm.^bet;Spm=Ypm*0;
alp2=.12;bet2=1;
Gampm2=alp2*Zpm.^(bet2);
qq=0.2*log(10)*bet;%
%ii=find(Zpm<60);Gampm(ii)=Gampm2(ii);
ii=find(httot>rng);htlm=ii(1)-1;
[aq bq]=size(Gampm);
for i=1:bq
    for j=2:htlm;
        Spm(j,i)=Spm(j-1,i)+Gampm(j,i)*0.025;
    end;
end;
Zpe=Zpm./(1-qq*Spm).^(1/bet);
dBZe=10*log10(Zpm./(1-qq*Spm).^(1/bet));

As=nrczp./nrczf;
%As=10.^-((2*rainpp*rng*sqrt(1.2)/1.3)/10);
As=10.^-((2*(rainpp)*rng*sqrt(1.2)/1.3)/10);
Spmm=0*Spm;
Spmm(1:htlm,:)=Spm(1:htlm,:);
for i=1:6
    Spmm(1:htlm,i)=Spmm(htlm,i)-Spmm(1:htlm,i);
end;
for i=1:6
    Zpe2(:,i)=Zpm(:,i)./(As(i).^bet+qq*Spmm(:,i)).^(1/bet);
    Zpe2q(:,i)=Zpe2(:,i)*fq(i);
end;
dBZe2=10*log10(Zpe2);
dBZe2q=10*log10(Zpe2q);

ii=find(httot>rng);htlm=ii(1)-1;ib=length(dBZe);

figure;
subplot(2,1,1);plot(Ypm(1:20,:),rng-httot(1:20),'.',Ypm(21:95,:),rng-httot(21:95));ylim([0 2.5]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');
subplot(2,1,2);plot(Ypmq(1:20,:),rng-httot(1:20),'.',Ypmq(21:95,:),rng-httot(21:95));ylim([0 2.5]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('CALWATER2');

figure;
%subplot(2,1,1);plot(dBZe(25:95,:),alt(5:75));ylim([0 2.35]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('Top Down');
subplot(2,1,1);plot(Ypmq(1:20,:),rng-httot(1:20),'.',Ypmq(21:94,:),rng-httot(21:94));ylim([0 2.5]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('Mean Zem');
subplot(2,1,2);plot(dBZe2q(1:20,:),rng-httot(1:20),'.',dBZe2q(21:94,:),rng-httot(21:94));ylim([0 2.5]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('Retrieved Ze');
%plot(10*log10(Zpm./(1-0.21*Spm).^(1/0.46)),rng- httot);ylim([0 2.5]);
figure;
%subplot(2,1,1);plot(dBZe(25:95,:),alt(5:75));ylim([0 2.35]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('Top Down');
plot(Ypmq(1:20,:),rng-httot(1:20),'.',Ypmq(21:94,:),rng-httot(21:94));ylim([0 2.5]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('Mean Zem');
hold on;plot(dBZe2q(1:20,:),rng-httot(1:20),'--',dBZe2q(21:94,:),rng-httot(21:94),'--');ylim([0 2.5]);xlim([-20 25]);xlabel('dBZe');ylabel('Altitude (km)');legend(rainm1,rainm2,rainm3,rainm4,rainm5,rainm6);title('Retrieved Ze');


figure;plot((Zpe2q(25:95,2:6)/15).^(1/1.1),alt(5:75));ylim([0 2.5]);xlim([-2 12]);xlabel('Rainrate (mm/h)');ylabel('Altitude (km)');legend(rainm2,rainm3,rainm4,rainm5,rainm6);title('Rainrate from retrieved Ze');
altm1=num2str(alt(18));altm1=altm1(1:4);
altm2=num2str(alt(38));altm2=altm2(1:4);
altm3=num2str(alt(58));altm3=altm3(1:4);
figure;loglog(rainpp,Zpmq(1,:),'d',rainpp,Zpe2(18,:),'x',rainpp,Zpe2(38,:),'o',rainpp,Zpe2(58,:),'d',rainpp,15*rainpp.^1.1);xlabel('Rainrate (mm/hr)');ylabel('Zei, Ze');
legend('Zei',[altm1 ' km'],[altm2 ' km'],[altm3 ' km'],'Z=15*Rr^{1.1}');

figure;subplot(2,1,1);semilogy(rainpp,Zpmq(1,:),'square',rainpp,Zpe2(28,:),'x',rainpp,Zpe2(38,:),'o',rainpp,Zpe2(58,:),'d',rainpp,15*rainpp.^1.1);xlabel('Rainrate (mm/hr)');ylabel('Ze, Zei');
legend('2.50 km (Zei)',[altm1 ' km'],[altm2 ' km'],[altm3 ' km'],'Z=15*Rr^{1.1}');
subplot(2,1,2);plot(rainpp,-.5*(Dyps(18,:)+.1).*fq,'x',rainpp,-.5*(Dyps(38,:)+.1).*fq,'o',rainpp,-.5*(Dyps(55,:)+.1).*fq,'d',rainpp,1*rainpp.^1.0);xlabel('Rainrate (mm/hr)');ylabel('Attenuation (dB/km)');
legend([altm1 ' km'],[altm2 ' km'],[altm3 ' km'],'\gamma=1.0*Rr^{1.0}')


figure;loglog(Zpe2(18,:),-.5*(Dyps(18,:)+.1).*fq,'x',Zpe2(38,:),-.5*(Dyps(38,:)+.1).*fq,'o',Zpe2(58,:),-.5*(Dyps(58,:)+.1).*fq,'d',Zpe2(58,:),.05*Zpe2(58,:).^1);xlabel('Zei');ylabel('Atten (dB/km)');
legend([altm1 ' km'],[altm2 ' km'],[altm3 ' km'],'\gamma=0.05*Ze^{1.0}')
figure;plot(Zpmq(1,:),-.5*Dyps(5,:).*fq,'square',Zpe2(18,:),-.5*(Dyps(18,:)+.1).*fq,'x',Zpe2(38,:),-.5*(Dyps(38,:)+.1).*fq,'o',Zpe2(58,:),-.5*(Dyps(55,:)+.1).*fq,'d',[0 160],.05*[0 160].^1);xlabel('Ze, Zei');ylabel('Atten (dB/km)');
legend('2.50 km (Zei)',[altm1 ' km'],[altm2 ' km'],[altm3 ' km'],'\gamma=0.05*Ze^{1.0}')


figure;
subplot(2,1,1);plot(hhws,rainwsr/100,'or',hhm,rain_ncs,'ok',[18 22],[0 0 ],'--');xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');;ylim([-1 20]);xlim([19 22]);legend('WSRA','NRCS');
subplot(2,1,2);plot(hhm,P3_rains,'or',hhm,rain_ncs,'ok',[18 22],[0 0 ],'--');xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');;ylim([-1 20]);xlim([19 22]);legend('SFMR','NRCS');

ii=find(hhm<23 & hhm>17);nanmeanx([rain_nc(ii)' P3_rains(ii)' rrwp(ii)'/10])
kk=find(SNRm(5,ii)>-0 & rain_nc(ii)>1 & ~isnan(slops(ii)) );Zrint=10.^(intr(ii(kk))/10);Zrfm=10.^((Refmc(5,ii(kk))-slops(ii(kk))*ht(5))/10);figure;loglog(Zrfm,-.5*slops(ii(kk)),'o',Zrfm,.065*(Zrfm).^1.1,'.',Zrfm,.8*(Zrfm).^.47,'.');xlim([1 300]);ylim([1 30]);xlabel('Zbei');ylabel('Attenuation (dBZe/km)');legend('\gamma Measaured','\gamma=0.065*Zfrm^{1.1}','\gamma=0.8*Zfrm^{0.47}');
kk=find(SNRm(5,ii)>-0 & rain_nc(ii)>1 & ~isnan(slops(ii)) );Zrint=10.^(intr(ii(kk))/10);Zrfm=10.^((Refmc(5,ii(kk))-slops(ii(kk))*ht(5))/10);figure;loglog(Zrfm,-.5*slops(ii(kk)),'o',Zrfm,.033*(Zrfm).^0.97,'.',Zrfm,.24*(Zrfm).^.57,'.',Zrfm,.035*(Zrfm).^1,'.',Zrfm,.012*(Zrfm).^1.12,'.');xlim([10 1000]);ylim([1 30]);xlabel('Zbei');ylabel('Attenuation (dBZe/km)');legend('\gamma Measaured','\gamma=0.033*Zfrm^{0.97} HMT','\gamma=0.24*Zfrm^{0.57} log-log','\gamma=0.035*Zfrm^{1.0} Bootstrap','\gamma=0.012*Zfrm^{1.12} Llermitte');
kk=find(sig_p(ii)<1 & ~isnan(slops(ii)) );Zrint=10.^(intr(ii(kk))/10);Zrfm=10.^((Refmc(5,ii(kk))-slops(ii(kk))*ht(5))/10);figure;loglog(Zrfm,-.5*slops(ii(kk)),'o',Zrfm,.033*(Zrfm).^0.97,'.',Zrfm,.24*(Zrfm).^.57,'.',Zrfm,.035*(Zrfm).^1,'.',Zrfm,.012*(Zrfm).^1.12,'.');xlim([10 1000]);ylim([1 30]);xlabel('Zbei');ylabel('Attenuation (dBZe/km)');legend('\gamma Measaured','\gamma=0.033*Zfrm^{0.97} HMT','\gamma=0.24*Zfrm^{0.57} log-log','\gamma=0.035*Zfrm^{1.0} Bootstrap','\gamma=0.012*Zfrm^{1.12} Llermitte');
%kk=find(sig_p(ii)<1 & ~isnan(slops(ii)) );Zrint=10.^(intr(ii(kk))/10);Zrfm=10.^((Refmc(5,ii(kk))-slops(ii(kk))*ht(5))/10);figure;loglog(Zrfm,-.5*slops(ii(kk)),'o',Zrfm,.033*(Zrfm).^0.97,'.',Zrfm,.24*(Zrfm).^.57,'.',Zrfm,.035*(Zrfm).^1,'.',Zrfm,.012*(Zrfm).^1.12,'.');xlim([10 1000]);ylim([1 30]);xlabel('Zbei');ylabel('Attenuation (dBZe/km)');legend('\gamma Measaured','\gamma=0.033*Zfrm^{0.97} HMT','\gamma=0.24*Zfrm^{0.57} log-log','\gamma=0.035*Zfrm^{1.0} Bootstrap','\gamma=0.012*Zfrm^{1.12} Llermitte');

figure;
subplot(2,1,1);plot(rain_nc(ii(kk)),10*log10(Zrfm),'.',rain_nc(ii(kk)),10*log10(25*rain_nc(ii(kk)).^0.91),'.',[1.07 2.15 4.18 7.79],[11 15 18 22.5],'or');ylim([5 35]);xlim([1 11]);ylabel('dBZei');legend('dBZei measured','Z=25*Rr^{0.91}','Mean Profile')
subplot(2,1,2);plot(rain_nc(ii(kk)),-.5*slops(ii(kk)),'.',rain_nc(ii(kk)),0.93*rain_nc(ii(kk)).^1,'.',[1.07 2.15 4.18 7.79],[1.8 2.25 3.9 7.2],'ro');ylim([0 15]);xlim([1 11]);xlabel('Rainrate (hh/hr)');ylabel('Attenuation (dBZ/km)');legend('\gamma Measured','\gamma=0.93*Rr^{1.0}','Mean Profiles')


figure;ij=find(P3_altrefm<3e3);
subplot(2,1,1);plot(hhm,nrcf,hhm(ij),nrcsm(ij));xlim([19 20]);xlabel('Hrs (UTC)');ylabel('NRCS (dB)');legend('NRCS','NRCS_m');
subplot(2,1,2);plot(hhm(ij),rain_ncs(ij),'.',[19 20],[0 0],'--');xlim([19 20]);ylim([-2 12]);xlabel('Hrs (UTC)');ylabel('Rainrate NRCS (mm/hr)');


Ym=[];bq=length(Refmc);
htl=6;htu=25;kkk=htl:htu;hhz=0.0131:.025:ht(htl);cc=length(hhz);[ac,bc]=size(Yp);Ym=-26*ones(ac+cc-htl,bq);slopm=NaN*hhm;intrm=NaN*hhm;
for i=1:bq
      ij=find(SNRm(kkk,i)>-10); %
      if length(ij)>5
        [PP,SS]=polyfit(ht(kkk(ij)),Refmc(kkk(ij),i),1);
        if PP(1)>0;
            PP(1)=0;
            PP(2)=Refmc(kkk(1),i);
        end;
        slopm(i)=PP(1);
        intrm(i)=min(30,PP(2));
        Ym(1:cc,i)=min(30,polyval(PP,hhz));
        Ym(cc+1:ac+cc-htl,i)=Refmc(htl+1:ac,i);
      end;
end;
%alp=.8;bet=0.47;
% alp=.12;bet=1.0;
%  alp=0.033;bet=0.97;
%alp=0.042;bet=.94;
%alp=0.06;bet=1.0;
Zm=10.^(Ym/10);
Gamm=alp*Zm.^bet;Sm=Zm*0;
%Gampm=0.2*Zpm.^(bet+.22);Spm=Ypm*0;
qm=0.2*log(10)*bet;%
alp2=.12;bet2=1;
Gamm2=alp2*Zm.^(bet2);
qm2=0.2*log(10)*bet2;%
%ii=find(Zm<60);Gamm(ii)=Gamm2(ii);

ii=find(httot>rng);htlm=ii(1)-1;
bq=length(Zm);atten2=[];
for i=1:bq
    for j=2:htlm;
        Sm(j,i)=Sm(j-1,i)+Gamm(j,i)*0.025;
        %Sm(j,i)=min(Sm(j,i),1/qm-.01);
        atten2(j,i)=-1/qm/httot(j)*log(1-qm*Sm(j,i));
    end;
end;
%Sm=Sm;
Ast=10.^((nrcsm-nrcf)/10);Ast=min(Ast,1);
Smm=0*Sm;htlm=94;
Smm(1:htlm,:)=Sm(1:htlm,:);atten=[];
for i=1:length(Sm);
    Smm(1:htlm,i)=Sm(htlm,i)-Sm(1:htlm,i);
    for j=2:htlm;
        nmx=Ast(i).^bet+qm*Smm(j,i);nmx=min(nmx,1);
        atten(j,i)=-1/qm/httot(j)*log(nmx);
    end;
end;
for i=1:length(Sm)
    nmx=Ast(i).^bet+qm*Smm(j,i);nmx=min(nmx,1);
    Zpe2t(:,i)=Zm(:,i)./(nmx).^(1/bet);
end;
dBZe2t=10*log10(Zpe2t);
rain_z=[];
[jdav4,Zpe2t_min,sigr]=timave_1_nonan(jdhhm,1,Zpe2t(60,:)');
ii=find(SNRm(60,:)>-10);
[jdav4,Ze_min,sigr]=timave_1_nonan(jdhhm(ii),1,Refm(60,ii)');

ii=find(SNRm(5,:)>-10);
[jdav4,DV_min,sigr]=timave_1_nonan(jdhhm(ii),1,DVamc(40,ii)');
rain_zc=(Zpe2t_min/15).^(0.91);
rain_zac=rain_zc;
ik=find(rain_zc>2);rain_zac(ik)=rain_dzmin(ik);
Zmc_min=10.^(Refmc_min/10);
rain_z=(Zmc_min(4,:)/15).^(0.91);
rain_za=rain_z;
ik=find(DV_min<-3.0);rain_za(ik)=rain_dzmin(ik);
figure;plot(jdav*24,rain_ncmin,jdav*24,rain_za);xlim([19 22]);xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');ylim([0 14]);
legend('NRCS','Hybrid Method3');


rain_in=[];
%rain_in=-2/qm*log(1-qm*Sm(70,:))*.5*1.3./sqrt(1.2)./P3_altrefm/1e-3;
pkh=30;ddh=rng-httot(pkh);ddh=double(ddh);
rain_in=-1/qm*log(Ast.^bet+qm*Smm(pkh,:))*.5*1.3./sqrt(1.2)/ddh;
rain_ins=sgolayfilt(rain_in,1,11);
[jdav4,rain_zimin,sigr]=timave_1_nonan(jdhhm,1,rain_in');

pkh2=70;
rain_in2=-1/qm*log(max(1-qm*Sm(70,:),.03))*1.3./sqrt(1.2)/httot(pkh2);;%./P3_altrefm/1e-3;
%rain_in2=atten(pkh2,:)*1.3./sqrt(1.2)/ddh;
rain_in2s=sgolayfilt(double(rain_in2),1,11);
[jdav4,rain_zimin2,sigr]=timave_1_nonan(jdhhm,1,rain_in2');

figure;plot(jdav*24,rain_ncmin,jdav*24,rain_slopmin,jdav*24,rain_dzmin,jdav*24,rain_zimin2);xlim([19 22]);xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');legend('Z_{integral}','NRCS');ylim([0 14]);
legend('NRCS','Slope Method1','\DeltadB Z\Deltaz Method2','Z Integral Method');

ii=find(~isnan(rain_zimin)' & ~isnan(rain2_min));corrcoef(rain_zimin(ii)',rain2_min(ii))
ii=find(~isnan(rain_zimin) & ~isnan(rain_ncmin));corrcoef(rain_zimin(ii)',rain_ncmin(ii))


figure;plot(hhm,rain_ins,'or',hhm,rain_ncs,'ok',[19 22],[0 0],'--');ylim([-2 12]);xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');legend('Z_{integral}','NRCS');
ii=find(rain_ncs>.1 & SNRm(5,:)>-10 & sig_p<1);figure;loglog(Zpe2t(40,ii),.5*atten(40,ii),'.',Zpe2t(30,ii),.5*atten(30,ii),'.',Zpe2t(40,ii),.03*Zpe2t(40,ii).^.94,'.');ylim([.1 10]);xlim([1 1000]);xlabel('Ze');ylabel('\gamma');

ii=find(hhm<23 & hhm>17);
kk=find(SNRm(5,ii)>-0 & rain_ncs(ii)>1 & ~isnan(slopm(ii)) );Zrintm=Zpe2t(50,ii(kk));%Zrintm=10.^(intrm(ii(kk))/10);
figure;
subplot(2,1,1);plot(rain_ncs(ii(kk)),10*log10(Zrintm),'.',rain_ncs(ii(kk)),10*log10(18*rain_ncs(ii(kk)).^1.1),'.');ylim([5 25]);xlim([1 11]);ylabel('dBZei');legend('dBZei','Z=15*Rr^{0.91}');
subplot(2,1,2);plot(rain_ncs(ii(kk)),-.5*slopm(ii(kk)),'.',rain_ncs(ii(kk)),1.2*rain_ncs(ii(kk)).^.91,'.');ylim([0 15]);xlim([1 11]);xlabel('Rainrate (hh/hr)');ylabel('Attenuation (dBZ/km)');legend('\gamma Measured','\gamma=1.2*Rr^{0.91}')

k=9500;kk=k:k+20;figure;plot(atten(:,kk)*1.3/sqrt(1.2),rng-httot(1:94),'.');ylim([0 2.]);xlabel('Rainrate (mm/hr)');ylabel('Altitude (km)');
rain_ab=atten*1.3/sqrt(1.2);
nq=floor(bq/33);di=33;rainp=zeros(94,nq);rain_ncp=[];
for i=1:nq
    kk=(i-1)*33+1:i*33;
    rainp(:,i)=nanmeanx(atten(:,kk)'*1.3/sqrt(1.2));
    hhmp(i)=nanmeanx(hhm(kk));
    rain_ncp(i)=nanmeanx(rain_ncs(kk));
end;
raini=atten(60,:)*1.3/sqrt(1.2);rainis=sgolayfilt(raini,1,11);
rain_all2=nanmedianx([rain_ncs' P3_rains' rrwp2' raini']');
rain_alldz=nanmedianx([rain_ins' rain_dzs' rain_slops' ]');



ii=find(~isnan(rain_ncs) & ~isnan(rain_alldz));[length(ii)/1000 mean([rain_ncs(ii)' rain_alldz(ii)'])]
corrcoef(rain_ncs(ii),rain_alldz(ii))

ii=find(~isnan(rain_ncs) & ~isnan(rain_2));[length(ii)/1000 mean([rain_nc(ii)' rain_2(ii)'])]
corrcoef(rain_ncs(ii),rain_2(ii))

ii=find(~isnan(rain_ncs) & ~isnan(P3_rains));[length(ii)/1000 mean([rain_ncs(ii)' P3_rains(ii)'])]
corrcoef(rain_ncs(ii),P3_rains(ii))

ii=find(~isnan(rain_ncs) & ~isnan(rrwp2s));[length(ii)/1000 mean([rain_ncs(ii)' rrwp2s(ii)'])]
corrcoef(rain_ncs(ii),rrwp2s(ii))

ii=find(~isnan(rain_all) & ~isnan(rain_alldz));[length(ii)/1000 mean([rain_all(ii)' rain_alldz(ii)'])]
corrcoef(rain_all(ii),rain_alldz(ii))

ii=find(~isnan(rain_ncs) & ~isnan(rain_ins));[length(ii)/1000 mean([rain_ncs(ii)' rain_ins(ii)'])]
corrcoef(rain_ncs(ii),rain_ins(ii))


rain_dzs0=rain_dzs;ii=find(isnan(rain_dzs));rain_dzs0(ii)=0;
ii=find(~isnan(rain_ncs) & ~isnan(rain_dzs));[length(ii)/1000 mean([rain_ncs(ii)' rain_dzs(ii)'])]
corrcoef(rain_ncs(ii),rain_dzs(ii))

ii=find(~isnan(rain_ncs) & ~isnan(rain_dzs0));[length(ii)/1000 mean([rain_ncs(ii)' rain_dzs0(ii)'])]
corrcoef(rain_ncs(ii),rain_dzs0(ii))


rain_slops0=rain_slops;ii=find(isnan(rain_slops));rain_slops0(ii)=0;
ii=find(~isnan(rain_ncs) & ~isnan(rain_slops));[length(ii)/1000 mean([rain_ncs(ii)' rain_slops(ii)'])]
corrcoef(rain_ncs(ii),rain_slops(ii))

ii=find(~isnan(rain_ncs) & ~isnan(rain_slops0));[length(ii)/1000 mean([rain_ncs(ii)' rain_slops0(ii)'])]
corrcoef(rain_ncs(ii),rain_slops0(ii))

ii=find(atten(70,:)>.5);nanmeanx([rain_ncs(ii)' P3_rains(ii)' rrwp2s(ii)' raini(ii)' rain_slops0(ii)' rain_dzs0(ii)'])

ii=find(~isnan(rain_ncs) & ~isnan(rain_2));[length(ii)/1000 mean([rain_ncs(ii)' rain_2(ii)'])]
corrcoef(rain_ncs(ii),rain_2(ii))

ii=find(~isnan(rain_slops) & ~isnan(rain_2));[length(ii)/1000 mean([rain_ncs(ii)' rain_2(ii)'])]
corrcoef(rain_slops(ii),rain_2(ii))

ii=find(~isnan(rain_dzs) & ~isnan(rain_2));[length(ii)/1000 mean([rain_ncs(ii)' rain_2(ii)'])]
corrcoef(rain_dzs(ii),rain_2(ii))

ii=find(~isnan(rain_ins) & ~isnan(rain_2));[length(ii)/1000 mean([rain_ncs(ii)' rain_2(ii)'])]
corrcoef(rain_ins(ii),rain_2(ii))

ii=find(~isnan(rrwp2s) & ~isnan(rain_2));[length(ii)/1000 mean([rain_ncs(ii)' rain_2(ii)'])]
corrcoef(rrwp2(ii),rain_2(ii))


figure;
subplot(2,1,1);plot(hhm,rain_slops0,'or',hhm,rain_nc,'xk',[19 22],[0 0],'--k');xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');title('CALWATER2');legend('Slope Method1','NRCS');ylim([-1 15]);
subplot(2,1,2);plot(hhm,rain_dzs0,'or',hhm,rain_nc,'xk',[19 22],[0 0],'--k');xlabel('Hrs (UTC)');ylabel('Rainrate (mm/hr)');legend('\DeltadBZ\Deltaz Method4','NRCS');ylim([-1 15]);


%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);


 
            lk=400;lm=2070;ln=2994;
            Refn=-32+20*log10(ht);rng=nanmedian(zzx);
            rng=2.43;figure;plot(Refmc(:,lk),rng-ht,'r',Refmc(:,lm),rng-ht,'b',Refmc(:,ln),rng-ht,'k',Refn,rng-ht,'g');ylim([0 2]);xlim([-40 20])
            xlabel('Rectivity (dBZe)');ylabel('Altitude (km)');
            legend('19.04','19.20','19.49','noise');
          


r=[.4 1 2.1 4.1 7.7];
vw=[3.4 3.7 3.9 4.1 4.3];dzw=[6 11 16 19 22];
vb=(vw/1.65).^2;%vb=(vw/1.48).^(1/0.58);
ii=find(vw<3);vb(ii)=vw(ii)+1/14*vw(ii);vb=min(vb,9);
dvb=vb-vw;
dzb=dzw+3*(dvb).^1.4;dzb=dzw+5.35*dvb.^1.17;
zb=10.^(dzb/10);
mu=0;
lb=(gamma(7.67+mu)*3.78/gamma(7+mu)./vb).^1.5;
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);
n0=bu*zb./vb.^((21+3*mu)/2);
rb=au*zb./vb.^3.5;

DV_mn=-DV_min;
vb_mn=(DV_mn/1.65).^2;
ii=find(DV_mn<3);vb_mn(ii)=DV_mn(ii)+1/14*DV_mn(ii);
ii=find(vb_mn<1);vb_mn(ii)=NaN;
dvb=vb_mn-DV_mn;
zb_mn=Zpe2t_min.*(3*(dvb).^1.4);

Rzmn=1.39*zb_mn./vb_mn.^3.5;




















