Rain = 0;
kind = 37;
kind = 104;
freq = [14,35,94];
showonly_FLAG =0;
figure_FLAG = 0;


if showonly_FLAG == 0
   if Rain == 1
      Rlog = (-1.2):0.1:2;
      R = 10.^(Rlog);
      D = 0.1:0.1:2;
      S = 0.34:0.04:0.42;
      mode = 5;
   else
      Rlog = -4:0.1:1;
      R = 10.^(Rlog);
      D = [0.02:0.02:0.1,0.2:0.2:6];
      S = 0.1:0.2:0.5;
      mode = 6;
   end
   
   NR = length(R);
   ND = length(D);
   NS = length(S);
   NF = length(freq);
   
   ze1 = zeros(NR,ND,NS,NF);
   kdb1 = ze1;
   avtm1 = zeros(NR,ND,NS);
   avtz1 = ze1;
   
   

   [ze,kdb,PSDf,eta,zf,kkm,sb,avtm,avtz,Mx,Rx,Dx,Sx] = DFDPR(freq,[kind,mode,R(1),D(1),S(1)],280,1);
   
   for iR = 1:NR
      for iD = 1:ND
         for iS = 1:NS
            [ze,kdb,PSDf,eta,zf,kkm,sb,avtm,avtz,Mx,Rx,Dx,Sx] = DFDPR(freq,[kind,mode,R(iR),D(iD),S(iS)],280,0);
            ze1(iR,iD,iS,:) = ze(1,:);
            kdb1(iR,iD,iS,:) = kdb(1,:);
            avtm1(iR,iD,iS,:) = avtm(1);
            avtz1(iR,iD,iS,:) = avtz(1,:);
         end
      end
      disp(['Rain: ' num2str(R(iR)) ' done'])
   end
end


if 1 == 0
ze1 = 10*log10(ze1);
s = S;

% Calculation of the Jacobians
dRlog = 0.1;
dDlog = 0.1;
ds = 0.04;

[dzdD,dzdR,dzds] = gradient(ze1(:,:,:,1));
dzdD = dzdD / dDlog;
dzds = dzds / ds;
dzdR = dzdR / dRlog;

[dz2dD,dz2dR,dz2ds] = gradient(ze1(:,:,:,2));
dz2dD = dz2dD / dDlog;
dz2ds = dz2ds / ds;
dz2dR = dz2dR / dRlog;

[dvdD,dvdR,dvds] = gradient(avtz1(:,:,:,1));
dvdD = dvdD / dDlog;
dvds = dvds / ds;
dvdR = dvdR / dRlog;

JZV_RD = (dzdR .* dvdD - dzdD .* dvdR);
JZV_RS = (dzdR .* dvds - dzds .* dvdR);

JZZ_RD = (dzdR .* dz2dD - dzdD .* dz2dR);
JZZ_RS = (dzdR .* dz2ds - dzds .* dz2dR);

JZ2V_RD = (dz2dR .* dvdD - dz2dD .* dvdR);
end

if 1 == 0
figure
subplot(2,2,1)
x = dzdR;
pcolor(D,Rlog,x(:,:,2)); shading flat, caxis([-20 20]), colorbar
subplot(2,2,2)
x = dzdD;
pcolor(D,Rlog,x(:,:,2)); shading flat, caxis([-20 20]), colorbar
subplot(2,2,3)
x = dvdR;
pcolor(D,Rlog,x(:,:,2)); shading flat, caxis([-20 20]), colorbar
subplot(2,2,4)
x = dvdD;
pcolor(D,Rlog,x(:,:,2)); shading flat, caxis([-20 20]), colorbar
end

figure
subplot(2,2,1)
x = JZV_RD;
pcolor(D,Rlog,10*log10(abs(x(:,:,2)))); shading flat, caxis([0 20]), colorbar
subplot(2,2,2)
x = JZZ_RD;
pcolor(D,Rlog,10*log10(abs(x(:,:,2)))); shading flat, caxis([0 20]), colorbar
subplot(2,2,3)
x = JZV_RS;
pcolor(D,Rlog,10*log10(abs(x(:,:,2)))); shading flat, caxis([0 20]), colorbar
subplot(2,2,4)
x = JZZ_RS;
pcolor(D,Rlog,10*log10(abs(x(:,:,2)))); shading flat, caxis([0 20]), colorbar


s1 = 2;
s0 = 1;
s2 = 3;


iRex = 7;
iDex = 10;

R = R * 0.65;

for ifr = 1:3
   if figure_FLAG == 1
      figure(ifr)
   else
      figure(1)
      subplot(3,1,ifr)
   end
   Zlevs = -20:10:40;
   [cs1,h1] = contour(D',R',10*log10(ze1(:,:,s1,ifr)),Zlevs,'k'); 
   set(h1,'linewid',2);
   hold on
   set(gca,'yscale','log')
   [cs,h] = contour(D',R',10*log10(ze1(:,:,s0,ifr)),Zlevs,'k:');                                                   
   [cs,h] = contour(D',R',10*log10(ze1(:,:,s2,ifr)),Zlevs,'k--');
   
   kdBlevs = [0.1,1,10];
   [cs2,h2] = contour(D',R',(kdb1(:,:,s1,ifr)),kdBlevs,'b'); 
   set(h2,'linewid',2);
   [cs,h] = contour(D',R',(kdb1(:,:,s0,ifr)),kdBlevs,'b:');      
   [cs,h] = contour(D',R',(kdb1(:,:,s2,ifr)),kdBlevs,'b--');      
   
   [cs3,h3] = contour(D',R',(avtz1(:,:,s1,ifr)),'r'); 
   set(h3,'linewid',2);
   % set(h,'color','r','fontsize',10,'fontw','bold')
   [cs,h] = contour(D',R',(avtz1(:,:,s0,ifr)),'r:');                                                  
   [cs,h] = contour(D',R',(avtz1(:,:,s2,ifr)),'r--');
   
   h = clabel(cs1,h1,'color','k','fontw','bold');  
   h = clabel(cs2,h2,'color','b','fontw','bold'); 
   if ~isempty(h3)
   h = clabel(cs3,h3,'manual','color','r','fontw','bold','fontsize',10,'hor','left'); 
   end
   
   if ifr == NF | figure_FLAG == 1
      if mode == 5
         xlabel('D" [mm] - Mean Drop Size (uncorrelated from R)');
      else
         xlabel('D_m [mm] - Mean Drop Size');
      end
   end
   
   if mode == 5
      ylabel('R [mm/hr] - Rain Rate')
   else
      ylabel('M [g/m^3] - Mass content')
   end
   
   if ifr == 0 | figure_FLAG == 2
      ht = title(['Black=Z [dBZ], Blue=k [dB/km], Red=v_p [m/s]';...
            'solid:s_m=0.34 , dot :s_m=0.34,dash:s_m=0.34']);
      set(ht,'fontsize',8)
      ht = text(1.5,80,'dotted  s* = 0.34');
      ht = text(1.5,60,'solid   s* = 0.38');
      ht = text(1.5,50,'dashed  s* = 0.42');
   end
   
   % grid on
   ht = text(D(ND-5),10*log10(R(NR-1)),[num2str(freq(ifr)) ' GHz']);
   set(ht,'fontsize',10,'fontw','bold')
   
   
   
   if 1 == 1
      figure(7)
      subplot(3,1,ifr)
      Zenow = 10*log10(ze1(:,:,:,ifr));
      Zex1 = (Zenow(iRex,iDex,s1));
      Zex0 = (Zenow(iRex,iDex,s0));
      Zex2 = (Zenow(iRex,iDex,s2));
      Zexmax = max([Zex0,Zex1,Zex2]);
      Zexmin = min([Zex0,Zex1,Zex2]);
      
      [cs,h] = contour(D',R',10*log10(ze1(:,:,s1,ifr)),[Zexmin-1,Zexmax+1],'k');
      hold on
      set(h,'linew',2)  
      hp = plot(D(iDex),R(iRex),'ro');
      set(hp,'linewid',3)
      
      know = (kdb1(:,:,:,ifr));
      kex1 = (know(iRex,iDex,s1));
      kex0 = (know(iRex,iDex,s0));
      kex2 = (know(iRex,iDex,s2));
      kexmax = max([kex0,kex1,kex2]);
      kexmin = min([kex0,kex1,kex2]);
      
      [cs,h] = contour(D',R',(kdb1(:,:,s1,ifr)),[kexmin-0.1,kexmax+0.1],'b'); 
      set(h,'linew',2)  
      
      vnow = (avtz1(:,:,:,ifr));
      vex1 = (vnow(iRex,iDex,s1));
      vex0 = (vnow(iRex,iDex,s0));
      vex2 = (vnow(iRex,iDex,s2));
      vexmax = max([vex0,vex1,vex2]);
      vexmin = min([vex0,vex1,vex2]);
      
      [cs,h] = contour(D',R',(avtz1(:,:,s1,ifr)),[vexmin-0.5,vexmax+0.5],'r'); 
      set(h,'linew',2)  
      set(gca,'yscale','log')
      grid on
      if ifr == NF, xlabel('D* [mm] - Average Drop Size'); end
      ylabel('R [mm/hr] - Rain Rate')
      if ifr == 1, title(['Black=Z [dBZ], Blue=k [dB/km], Red=v_p [m/s]']); end
      ht = text(D(ND-5),10*log10(R(NR-1)),[num2str(freq(ifr)) ' GHz']);
      set(ht,'fontsize',10,'fontw','bold')
      
      
   end
   
end

ifr = 1;
vd = zeros(NR,ND,3);
vd(:,:,1) = (avtz1(:,:,s1,ifr)-avtm1(:,:,s1));
vd(:,:,2) = (avtz1(:,:,s0,ifr)-avtm1(:,:,s0));
vd(:,:,3) = (avtz1(:,:,s2,ifr)-avtm1(:,:,s2));
stvd = std(vd,0,3);


figure
[cs,h] = contour(D',R',vd(:,:,1),'r'); set(h,'linewid',2), h = clabel(cs,h); 
set(h,'color','r','fontsize',10,'fontw','bold')
hold on
set(gca,'yscale','log')
xlabel('D* [mm] - Average Drop Size')
ylabel('R [mm/hr] - Rain Rate')
title(['Backscattering wgt. avg. Velocity - Mass wgt. avg. Velocity [m/s]'])
grid on
ht = text(2.5,80,[num2str(freq(ifr)) ' GHz']);
set(ht,'fontsize',10,'fontw','bold')


figure
[cs,h] = contour(D',R',stvd,'r'); set(h,'linewid',2), h = clabel(cs,h); 
set(h,'color','r','fontsize',10,'fontw','bold')
set(gca,'yscale','log')
xlabel('D* [mm] - Average Drop Size')
ylabel('R [mm/hr] - Rain Rate')
title(['Backscattering wgt. avg. Velocity - Mass wgt. avg. Velocity (std) [m/s]'])
grid on
ht = text(2.5,80,[num2str(freq(ifr)) ' GHz']);
set(ht,'fontsize',10,'fontw','bold')



%[cs,h] = contour(D',R',(avtz1(:,:,s0,ifr)-avtm1(:,:,s0)),'b'); h = clabel(cs,h); 
%set(h,'color','b','fontsize',10)                                                 
%[cs,h] = contour(D',R',(avtz1(:,:,s2,ifr)-avtm1(:,:,s2)),'k'); h = clabel(cs,h); 
%set(h,'color','k','fontsize',10)


figure
[cs,h] = contour(D',R',10*log10(ze1(:,:,s1,1))-10*log10(ze1(:,:,s1,2)),'k'); set(h,'linewid',2), h = clabel(cs,h); 
set(gca,'yscale','log')

figure
[cs,h] = contour(D',R',(kdb1(:,:,s1,1))-(kdb1(:,:,s1,2)),[-10,-5,-2,-1,-0.2,0],'k'); set(h,'linewid',2), h = clabel(cs,h); 
set(gca,'yscale','log')