%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Statistics of Convective, Stratiform, Other, and No Rain 
%                     regions of a tropical cyclone.
%
% PURPOSE:  1) Creates statistical graphs and output of reflectivity for
%          convective, stratiform, other, and no rain regions of a cyclone.
% 
% DATA:  TRMM Precipitation Radar Version 6 Reflectivity/Geolocation
%       
% INPUT: TClocation txt file: Data, Time, Center Location, Orbit Number,
%        Storm Number, Storm Name.
%
% OUTPUT: Percent of areal coverage of each region as a function of raining
% areas. Contour Frequency by Altitude Diagrams (CFADs).  Vertical profiles 
% of mean reflectivity. Cumulative Frequency Diagram (CDFs). 
% Probablity Density Function (PDFs)
%         
%
% Created October 24, 2006
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                        % Initializes Variables %

% count variables for percentage                        
convect_cnt=0;strat_cnt=0;other_cnt=0;total_cnt=0;

% count variables for mean reflectivity
meanindex=ones(1,80);Z_binsize=5;
convect_Z=meanindex;strat_Z=meanindex;other_Z=meanindex;tot_Z=meanindex;
convect_Z(:)=0;strat_Z(:)=0;other_Z(:)=0;tot_Z(:)=0;
convect_Zcnt=meanindex;strat_Zcnt=meanindex;other_Zcnt=meanindex;
tot_Zcnt=meanindex;
convect_Zcnt(:)=0;strat_Zcnt(:)=0;other_Zcnt(:)=0;tot_Zcnt(:)=0;

% count variables for PDF
convect_ZPDF=(0:Z_binsize:70);strat_ZPDF=(0:Z_binsize:70);
other_ZPDF=(0:Z_binsize:70);tot_ZPDF=(0:Z_binsize:70);
convect_ZPDF(:)=0;strat_ZPDF(:)=0;other_ZPDF(:)=0;tot_ZPDF(:)=0;
convect_PDFcnt=0;strat_PDFcnt=0;other_PDFcnt=0;tot_PDFcnt=0;

% count variables for CFAD
CFADi=ones(1,80); CFADj=ones(1,15);
convect_CFADZ=CFADi'*CFADj;strat_CFADZ=CFADi'*CFADj;other_CFADZ=CFADi'*CFADj;
tot_CFADZ=CFADi'*CFADj;
convect_CFADcnt=CFADi;strat_CFADcnt=CFADi;other_CFADcnt=CFADi;
tot_CFADcnt=CFADi;
convect_CFADZ(:,:)=0;strat_CFADZ(:,:)=0;other_CFADZ(:,:)=0;tot_CFADZ(:,:)=0;
convect_CFADcnt(:)=0;strat_CFADcnt(:)=0;other_CFADcnt(:)=0;tot_CFADcnt(:)=0;

% count variables for CDF
CDFi=ones(1,15);
convect_CDF=CDFi;strat_CDF=CDFi;other_CDF=CDFi;tot_CDF=CDFi;
convect_CDF(:)=0;strat_CDF(:)=0;other_CDF(:)=0;tot_CDF(:)=0;

% count variables for CCFAD
CFADi=ones(1,80); CFADj=ones(1,15);
convect_CCFAD=CFADi'*CFADj;strat_CCFAD=CFADi'*CFADj;other_CCFAD=CFADi'*CFADj;
tot_CCFAD=CFADi'*CFADj;
convect_CCFAD(:,:)=0;strat_CCFAD(:,:)=0;other_CCFAD(:,:)=0;tot_CCFAD(:,:)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                           % Loading Data %

% Loads Tropical Cyclone Location text file
[YR,MN,DY,TM,LAT,LON,DIST,ORB,NM]=textread(...
      '/home/meteo/abh5000/non_RI_Atlantic.txt'...
      ,'%q %q %q %q %q %q %q %q %q');

  andy='non-RI';
  
% sets size of TC location text file to a variable
[TCi,TCj]=size(YR);

% Beginning of for loop to loop through all storms in TC location txt file
for snum=15
  
    % Creates filename to open text file with reflectivity/geolocation data
          if (str2num((char(YR(snum))))<10)
            if (str2num((char(MN(snum))))<10)
                if (str2num((char(DY(snum))))<10)
                    file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.0',...
                        char(YR(snum)),'0',char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                else
                    file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.0',...
                        char(YR(snum)),'0',char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                end
            else if (str2num((char(DY(snum))))<10)
                     file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.0',...
                        char(YR(snum)),char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                else
                     file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.0',...
                        char(YR(snum)),char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                end
            end
        else if (str2num((char(MN(snum))))<10)
                    if (str2num((char(DY(snum))))<10) 
                        file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.',...
                        char(YR(snum)),'0',char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                    else
                        file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.',...
                        char(YR(snum)),'0',char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                    end

             else  if (str2num((char(DY(snum))))<10) 
                        file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.',...
                        char(YR(snum)),'0',char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                   else
                        file=strcat(...
                        '/home/meteo/abh5000/HDFfiles/2A25.',...
                        char(YR(snum)),'0',char(MN(snum)),...
                        char(DY(snum)),'.',char(ORB(snum)),...
                        '.6Z80levels.txt');
                   end 
            end
         end
    
    % Loads above file
    data=load(file);
    
    % sets size of data to a variable
    [di,dj]=size(data);
    di=(di/82);
    
    % Sets geolocation data to variables
    lon=data(1:49,2:402);
    lat=data(50:98,2:402);
    
    % Sets dBZ data to variables
    for v=1:80
        for y=2:402
            for x=1:49
                dBZ(x,y-1,v)=data(98+(x)+(v-1)*49,y);
            end
        end
    end
    
    % Takes only reflectivity from within 200 km of center
rLAT = (str2num((char(LAT(snum)))));
rLON = (str2num((char(LON(snum)))));
for iloc = 1:49
    for jloc = 1:401
        dist=distance(rLAT,rLON,lat(iloc,jloc),lon(iloc,jloc))*111.1218;
        if dist >= 300
            for kloc = 1:80
                dBZ(iloc,jloc,kloc)=NaN;
%           rLAT = (str2num((char(LAT(snum))))*(pi/180));
 %          rLON = (str2num((char(LON(snum))))*(pi/180));
 %          for iloc = 1:49
  %             for jloc = 1:401
   %                rLAT_1 = (lat(iloc,jloc)*(pi/180));
    %               rLON_1 = (lon(iloc,jloc)*(pi/180));
     %              dlat = rLAT - rLAT_1;
      %             dlon = rLON - rLON_1*(cos(rLAT));
       %            dlat_km = rad2km(dlat);
        %           dlon_km = rad2km(dlon);
         %          distance = sqrt(((dlat_km)^2)+((dlon_km)^2));
          %         if distance >= 300
           %            for kloc = 1:80
            %               dBZ(iloc,jloc,kloc)=NaN; 
            end
        end
    end
end
    
    % cleans dBZ variable (removes any negative numbers)
    aek=dBZ;
    k1=find(aek<0);
    dBZ(k1)=NaN;
   
    % Converts dBZ to Z
    for x=1:49    
        for y=1:401
            for v=1:80
                Z(x,y,v)=10^((dBZ(x,y,v))/10);
            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   % Defining/pre-setting variable %
       
    qq=(1:di); % Creates 1-D array w/ size of 'zonal' index
    rr=(2:dj); % Creates 1-D array w/ size of 'meridonal' index
    qq(:)=1;   
    rr(:)=1;
    Ltot=qq'*rr; % Creates a 2-D array-Lower Layer Total Reflect. variab. 
    Lcnt=qq'*rr; % Creates a 2-D array-Lower Layer Total count variab. 
    Mtot=qq'*rr; % Creates a 2-D array-Middle Layer Total Reflect. variab. 
    Mcnt=qq'*rr; % Creates a 2-D array-Middle Layer Total count variab. 
    Utot=qq'*rr; % Creates a 2-D array-Upper Layer Total Reflect. variab. 
    Ucnt=qq'*rr; % Creates a 2-D array-Upper Layer Total count variab. 
    Dtot=qq'*rr; % Creates a 2-D array-Deep Layer Total Reflect. variab. 
    Dcnt=qq'*rr; % Creates a 2-D array-Deep Layer Total count variab. 
    Ttot=qq'*rr; % Creates a 2-D array-Total Reflect. variab. 
    Tcnt=qq'*rr; % Creates a 2-D array-Total count variab. 
% Initially set counts to zero
    Ltot(:,:)=0;Lcnt(:,:)=0;Mtot(:,:)=0;Mcnt(:,:)=0;Utot(:,:)=0;Ucnt(:,:)=0;
    Dtot(:,:)=0; Dcnt(:,:)=0;Ttot(:,:)=0; Tcnt(:,:)=0;
    regime=qq'*rr; % Creates a 2-D array-variable indicating regions of 
                   %   convective(1)/stratiform(2)/other(3)/no rain(4)
                   %   based on index. 
    regime(:,:)=0; % Initially set to 0
                    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                        % Computing Layer Averages %

    % Lower Layer Total Reflectivity/Count (Height: 0.5-4km/Index: 2-16)
    for x=1:49
        for y=1:401
            for v=2:16
                if ~isnan(Z(x,y,v))
                    Ltot(x,y)=Ltot(x,y)+Z(x,y,v);
                    Lcnt(x,y)=Lcnt(x,y)+1;
                end
            end
        end
    end
    
    % Lower Layer Reflectivity Average- Converts back to dBZ
    for x=1:49
        for y=1:401
            Lavg(x,y)=10*log10((Ltot(x,y)/Lcnt(x,y)));
        end
    end
    
    % Middle Layer Total Reflectivity/Count (Height: 5.5-8km/Index: 22-32)
    for x=1:49
        for y=1:401
            for v=22:32
                if ~isnan(Z(x,y,v))
                    Mtot(x,y)=Mtot(x,y)+Z(x,y,v);
                    Mcnt(x,y)=Mcnt(x,y)+1;
                end
            end
        end
    end
    
    % Middle Layer Reflectivity Average- Converts back to dBZ
    for x=1:49
        for y=1:401
            Mavg(x,y)=10*log10((Mtot(x,y)/Mcnt(x,y)));
        end
    end
    
    % Upper Layer Total Reflectivity/Count (Height: 10-13km/Index: 40-52)
    for x=1:49
        for y=1:401
            for v=40:52
                if ~isnan(Z(x,y,v))
                    Utot(x,y)=Utot(x,y)+Z(x,y,v);
                    Ucnt(x,y)=Ucnt(x,y)+1;
                end
            end
        end
    end
    
    % Upper Layer Reflectivity Average- Converts back to dBZ
    for x=1:49
        for y=1:401
            Uavg(x,y)=10*log10((Utot(x,y)/Ucnt(x,y)));
        end
    end
    
    % Deep Layer Total Reflectivity/Count (Height: 0.5-8km/Index: 2-32)
    for x=1:49
        for y=1:401
            for v=2:32
                if ~isnan(Z(x,y,v))
                    Dtot(x,y)=Dtot(x,y)+Z(x,y,v);
                    Dcnt(x,y)=Dcnt(x,y)+1;
                end
            end
        end
    end
    
    % Deep Layer Reflectivity Average- Converts back to dBZ
    for x=1:49
        for y=1:401
            Davg(x,y)=10*log10((Dtot(x,y)/Dcnt(x,y)));
        end
    end
    
    % Total Reflectivity/Count (Height: 0.5-8km/Index: 2-32)
    for x=1:49
        for y=1:401
            for v=1:80
                if ~isnan(Z(x,y,v))
                    Ttot(x,y)=Ttot(x,y)+Z(x,y,v);
                    Tcnt(x,y)=Tcnt(x,y)+1;
                end
            end
        end
    end
    
    % Total Reflectivity Average- Converts back to dBZ
    for x=1:49
        for y=1:401
            Tavg(x,y)=10*log10((Ttot(x,y)/Tcnt(x,y)));
        end
    end
    
    % Cleaning up Layer Averages. Average < 1 = NaN
%     aek=Lavg;
%     LL=find(aek<1);
%     Lavg(LL)=0;
%     aek=Mavg;
%     MM=find(aek<1);
%     Mavg(MM)=0;
%     aek=Uavg;
%     UU=find(aek<1);
%     Uavg(UU)=0;
%     aek=Davg;
%     DD=find(aek<1);
%     Davg(DD)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Beginning of Sorting Algorithm %

     for x=1:49
        for y=1:401
            if Lavg(x,y)>=17
                if Lavg(x,y)>=30
                    if (Mavg(x,y)>=27) || (Uavg(x,y)>=20)
                        regime(x,y)=1;
                    else
                        if Davg(x,y)>=40
                            regime(x,y)=3;
                        else
                            regime(x,y)=2;
                        end
                    end
                else
                    if (Mavg(x,y)>=25) || (Uavg(x,y)>=20) 
                        regime(x,y)=3;                  
                    else
                        if Davg(x,y)>=30
                            regime(x,y)=3;  
                        else
                            regime(x,y)=2;
                        end
                    end
                end
            else
                if Mavg(x,y)>=20
                    regime(x,y)=3;
                else
                    regime(x,y)=4;
                end
            end
        end
     end
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Counts number of points that are convective, stratiform %
               % other and total raining points % 
  
    for x=1:49
        for y=1:401
            for v=1:80
                if ~isnan(dBZ(x,y,v))
                  if dBZ(x,y,v)>1
                      if regime(x,y)==1;
                          convect_cnt=convect_cnt+1;
                      end
                      if regime(x,y)==2;
                          strat_cnt=strat_cnt+1;
                      end
                      if regime(x,y)==3;
                          other_cnt=other_cnt+1;
                      end
                      total_cnt=total_cnt+1;
                  end
                end
            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Calculates Mean reflectivity components for all area types %

    % Adds Z values in convective regions and counts number of convective
    % points
    for x=1:49
        for y=1:401
            if regime(x,y)==1;
                for v=1:80
                    if ~isnan(dBZ(x,y,v));
                        if dBZ(x,y,v)>1
                            convect_Z(v)=convect_Z(v)+Z(x,y,v);
                            convect_Zcnt(v)=convect_Zcnt(v)+1;
                        end
                    end
                end
            end
        end
    end
    
    % Adds Z values in stratiform regions and counts number of stratiform
    % points
    for x=1:49
        for y=1:401
            if regime(x,y)==2;
                for v=1:80
                    if ~isnan(dBZ(x,y,v));
                        if dBZ(x,y,v)>1
                            strat_Z(v)=strat_Z(v)+Z(x,y,v);
                            strat_Zcnt(v)=strat_Zcnt(v)+1;
                        end
                    end
                end
            end
        end
    end
    
    % Adds Z values in other regions and counts number of other
    % points
    for x=1:49
        for y=1:401
            if regime(x,y)==3;
                for v=1:80
                    if ~isnan(dBZ(x,y,v));
                        if dBZ(x,y,v)>1
                            other_Z(v)=other_Z(v)+Z(x,y,v);
                            other_Zcnt(v)=other_Zcnt(v)+1;
                        end
                    end
                end
            end
        end
    end
    
     % Adds Z values in total and counts number of other
    % points
    for x=1:49
        for y=1:401
            for v=1:80
                if ~isnan(dBZ(x,y,v));
                    if dBZ(x,y,v)>1
                        tot_Z(v)=tot_Z(v)+Z(x,y,v);
                        tot_Zcnt(v)=tot_Zcnt(v)+1;
                    end
                end
            end
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Calculates PDF components %
      
    % Calculates Convective PDF components
    for x=1:49
        for y=1:401
            if regime(x,y)==1;
                for v=1:80
                    if ~isnan(dBZ(x,y,v))
                        if dBZ(x,y,v)>1
                            k=floor(dBZ(x,y,v)/Z_binsize)+1;
                            if k<16
                                convect_ZPDF(k)=convect_ZPDF(k)+1;
                            else
                                convect_ZPDF(15)=convect_ZPDF(15)+1;
                            end
                            convect_PDFcnt=convect_PDFcnt+1;
                        end
                    end
                end
            end
        end
    end
 
    % Calculates Stratiform PDF components
     for x=1:49
        for y=1:401
            if regime(x,y)==2;
                for v=1:80
                    if ~isnan(dBZ(x,y,v))
                        if dBZ(x,y,v)>1
                            k=floor(dBZ(x,y,v)/Z_binsize)+1;
                            if k<16
                                strat_ZPDF(k)=strat_ZPDF(k)+1;
                            else
                                strat_ZPDF(15)=strat_ZPDF(15)+1;
                            end
                            strat_PDFcnt=strat_PDFcnt+1;
                        end
                    end 
                end
            end
        end
     end

     % Calculates Other PDF components
     for x=1:49
        for y=1:401
            if regime(x,y)==3;
                for v=1:80
                    if ~isnan(dBZ(x,y,v))
                        if dBZ(x,y,v)>1
                            k=floor(dBZ(x,y,v)/Z_binsize)+1;
                            if k<16
                                other_ZPDF(k)=other_ZPDF(k)+1;
                            else
                                other_ZPDF(15)=other_ZPDF(15)+1;
                            end
                            other_PDFcnt=other_PDFcnt+1;
                        end
                    end 
                end
            end
        end
     end
     
     % Calculates Total PDF components
     for x=1:49
        for y=1:401
            for v=1:80
                if ~isnan(dBZ(x,y,v))
                    if dBZ(x,y,v)>1
                        k=floor(dBZ(x,y,v)/Z_binsize)+1;
                        if k<16
                            tot_ZPDF(k)=tot_ZPDF(k)+1;
                        else
                            tot_ZPDF(15)=tot_ZPDF(15)+1;
                        end
                        tot_PDFcnt=tot_PDFcnt+1;
                    end
                end 
            end
        end
     end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Calculates CFAD components %
                   
     % Calculates CFAD convective components
     for x=1:49
         for y=1:401
             if regime(x,y)==1;
                for v=1:80
                    if ~isnan(dBZ(x,y,v))
                        if dBZ(x,y,v)>1
                            k=floor(dBZ(x,y,v)/Z_binsize)+1;
                            if k<16
                                convect_CFADZ(v,k)=convect_CFADZ(v,k)+1;
                            else
                                convect_CFADZ(v,15)=convect_CFADZ(v,15)+1;
                            end
                            convect_CFADcnt(v)=convect_CFADcnt(v)+1;
                        end
                    end 
                end
            end
        end
     end
     
      % Calculates CFAD stratiform components
     for x=1:49
         for y=1:401
             if regime(x,y)==2;
                for v=1:80
                    if ~isnan(dBZ(x,y,v))
                        if dBZ(x,y,v)>1
                            k=floor(dBZ(x,y,v)/Z_binsize)+1;
                            if k<16
                                strat_CFADZ(v,k)=strat_CFADZ(v,k)+1;
                            else
                                strat_CFADZ(v,15)=strat_CFADZ(v,15)+1;
                            end
                            strat_CFADcnt(v)=strat_CFADcnt(v)+1;
                        end
                    end 
                end
            end
        end
     end
     
      % Calculates CFAD other components
     for x=1:49
         for y=1:401
             if regime(x,y)==3;
                for v=1:80
                    if ~isnan(dBZ(x,y,v))
                        if dBZ(x,y,v)>1
                            k=floor(dBZ(x,y,v)/Z_binsize)+1;
                            if k<16
                                other_CFADZ(v,k)=other_CFADZ(v,k)+1;
                            else
                                other_CFADZ(v,15)=other_CFADZ(v,15)+1;
                            end
                            other_CFADcnt(v)=other_CFADcnt(v)+1;
                        end
                    end 
                end
            end
         end
     end
     
      % Calculates CFAD Total components
     for x=1:49
         for y=1:401
            for v=1:80
                if ~isnan(dBZ(x,y,v))
                    if dBZ(x,y,v)>1
                        k=floor(dBZ(x,y,v)/Z_binsize)+1;
                        if k<16
                            tot_CFADZ(v,k)=tot_CFADZ(v,k)+1;
                        else
                            tot_CFADZ(v,15)=tot_CFADZ(v,15)+1;
                        end
                        tot_CFADcnt(v)=tot_CFADcnt(v)+1;
                    end
                end 
            end
         end
     end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% End for loop- includes all 34 storms
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     % Calculates Areal Percentage %

percent_conv=(convect_cnt/total_cnt)*100;
percent_strat=(strat_cnt/total_cnt)*100;
percent_other=(other_cnt/total_cnt)*100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Calculates Mean reflectivity components for all area types %

% Calculates Convective Mean Reflectivity
for v=1:80;
    mean_convect(v)=10*log10((convect_Z(v))/(convect_Zcnt(v)));
end

% Calculates Stratiform Mean Reflectivity
for v=1:80
   mean_strat(v)=10*log10((strat_Z(v))/(strat_Zcnt(v)));
end

% Calculates Other Mean Reflectivity
for v=1:80
  mean_other(v)=10*log10((other_Z(v))/(other_Zcnt(v)));
end

% Calculates Total Mean Reflectivity
for v=1:80
  mean_tot(v)=10*log10((tot_Z(v))/(tot_Zcnt(v)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            % Calculates PDF %
      
% Calculates Convective PDF 
for pdfcnt=1:15
    PDF_convect(pdfcnt)=(convect_ZPDF(pdfcnt)/convect_PDFcnt)*100;            
end

% Calculates Stratiform PDF 
for pdfcnt=1:15
    PDF_strat(pdfcnt)=(strat_ZPDF(pdfcnt)/strat_PDFcnt)*100;            
end

% Calculates Other PDF 
for pdfcnt=1:15
    PDF_other(pdfcnt)=(other_ZPDF(pdfcnt)/other_PDFcnt)*100;            
end

% Calculates Total PDF 
for pdfcnt=1:15
    PDF_tot(pdfcnt)=(tot_ZPDF(pdfcnt)/tot_PDFcnt)*100;            
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            % Calculates CDF %

% Calculates Convective CDF
convect_CDF(1)=PDF_convect(1);
for cdfcnt=2:15
    convect_CDF(cdfcnt)=convect_CDF(cdfcnt-1)+PDF_convect(cdfcnt);          
end

% Calculates Stratiform CDF 
strat_CDF(1)=PDF_strat(1);
for cdfcnt=2:15
   strat_CDF(cdfcnt)=strat_CDF(cdfcnt-1)+PDF_strat(cdfcnt);       
end

% Calculates Other CDF
other_CDF(1)=PDF_other(1);
for cdfcnt=2:15
    other_CDF(cdfcnt)=other_CDF(cdfcnt-1)+PDF_other(cdfcnt);            
end

% Calculates Total CDF
tot_CDF(1)=PDF_tot(1);
for cdfcnt=2:15
    tot_CDF(cdfcnt)=tot_CDF(cdfcnt-1)+PDF_tot(cdfcnt);            
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                % Calculates CFAD %

 % Calculates CFAD convective 
for v=1:80
    for cfadcnt=1:15
        CFAD_convect(v,cfadcnt)=(convect_CFADZ(v,cfadcnt)/convect_CFADcnt(v))*100;
    end
end

  % Calculates CFAD stratiform 
for v=1:80
    for cfadcnt=1:15
        CFAD_strat(v,cfadcnt)=(strat_CFADZ(v,cfadcnt)/strat_CFADcnt(v))*100;  
    end
end

  % Calculates CFAD other 
for v=1:80
    for cfadcnt=1:15
        CFAD_other(v,cfadcnt)=(other_CFADZ(v,cfadcnt)/other_CFADcnt(v))*100; 
    end
end

  % Calculates CFAD Total
for v=1:80
    for cfadcnt=1:15
        CFAD_tot(v,cfadcnt)=(tot_CFADZ(v,cfadcnt)/tot_CFADcnt(v))*100; 
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                        % Calculates Cumulative CFAD %
                    
% Sets first bin to first percentage for all levels
for v=1:80
    convect_CCFAD(v,1)=CFAD_convect(v,1);
    strat_CCFAD(v,1)=CFAD_strat(v,1);
    other_CCFAD(v,1)=CFAD_other(v,1);
    tot_CCFAD(v,1)=CFAD_tot(v,1);
end
    
% Calculates Convective Cumulative CFAD
for v=1:80
    for cdfcnt=2:15
        if ~isnan(CFAD_convect(v,cdfcnt));
            convect_CCFAD(v,cdfcnt)=convect_CCFAD(v,cdfcnt-1)+CFAD_convect(v,cdfcnt);   
        else
            convect_CCFAD(v,cdfcnt)=convect_CCFAD(v,cdfcnt-1)+0;
        end
    end
end

% Calculates Stratiform Cumulative CFAD 
for v=1:80
    for cdfcnt=2:15
        if ~isnan(CFAD_strat(v,cdfcnt));
            strat_CCFAD(v,cdfcnt)=strat_CCFAD(v,cdfcnt-1)+CFAD_strat(v,cdfcnt);
        else
            strat_CCFAD(v,cdfcnt)=strat_CCFAD(v,cdfcnt-1)+0;
        end
    end
end

% Calculates Other Cumulative CFAD
for v=1:80
    for cdfcnt=2:15
        if ~isnan(CFAD_other(v,cdfcnt));
            other_CCFAD(v,cdfcnt)=other_CCFAD(v,cdfcnt-1)+CFAD_other(v,cdfcnt);
        else
            other_CCFAD(v,cdfcnt)=other_CCFAD(v,cdfcnt-1)+0;
        end
    end
end

% Calculates Total Cumulative CFAD
for v=1:80
    for cdfcnt=2:15
        if ~isnan(CFAD_tot(v,cdfcnt));
            tot_CCFAD(v,cdfcnt)=tot_CCFAD(v,cdfcnt-1)+CFAD_tot(v,cdfcnt);
        else
            tot_CCFAD(v,cdfcnt)=tot_CCFAD(v,cdfcnt-1)+0;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                 % Graphics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Initializing graphic Variables %
              
dBZbin=([0:14])*5;
hgt=([1:80])*0.25;
ii=dBZbin;
jj=hgt;
ii(:)=1;
jj(:)=1;
CFADx=jj'*dBZbin;
CFADy=hgt'*ii;
meanx=[1:49];
meany=[1:401];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            % Defining colorbar for CFAD %
                            
vv=[0,0.1,0.5,1,2,5,15,30,40,45];
minv=min(diff(vv));
new_v=min(vv):minv:max(vv);
map1=([0.644 0.118 0.800;0.074 0.261 0.800;0.439 0.673 1.00;0.371 0.728 0.579;...
        0.020 0.800 0.000;0.665 0.825 0.318;...
        0.939 0.664 0.280; 1.0 0.578 0.00;0.860 0.016 0.394]);
d_colorbar=zeros(length(new_v),3);

for ax=1:2
    d_colorbar(ax,:)=map1(1,:);
end
for ax=2:3
    d_colorbar(ax,:)=map1(2,:);
end
for ax=3:10
    d_colorbar(ax,:)=map1(3,:);
end
for ax=11:20
    d_colorbar(ax,:)=map1(4,:);
end
for ax=21:50
    d_colorbar(ax,:)=map1(5,:);
end
for ax=51:150
    d_colorbar(ax,:)=map1(6,:);
end
for ax=151:300
    d_colorbar(ax,:)=map1(7,:);
end
for ax=301:400
    d_colorbar(ax,:)=map1(8,:);
end
for ax=401:451
    d_colorbar(ax,:)=map1(9,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 % Defining colorbar for Cumulative CFAD %
                            
ww=[0,0.01,0.1,1,5,25,50,75,95,99,99.9,99.99,105];
minw=min(diff(ww));
new_w=min(ww):0.001:max(ww);
map2=([0.800 0.255 0.769;0.644 0.118 0.800;0.074 0.261 0.800;...
        0.439 0.673 1.00;0.371 0.728 0.579;...
        0.020 0.800 0.000;0.665 0.825 0.318;0.996 0.996 0.000;...
        0.939 0.664 0.280; 1.0 0.578 0.00;0.996 0.255 0.351;...
        0.860 0.016 0.394]);
dw_colorbar=zeros(length(new_w),3);

[wx,wy]=size(ww);
wcnt=1;
while wcnt~=wy+1;
    for i=1:length(new_w)
        if new_w(i)==ww(wcnt);
            index(wcnt)=i;
            wcnt=wcnt+1;
            break;
        end
    end
end

for i=index(1):index(2)
    dw_colorbar(i,:)=map2(1,:);
end

colorcnt=2;
indcnt=2;
while indcnt~=length(index)
    for jj=(index(indcnt)+1):index(indcnt+1)
        dw_colorbar(jj,:)=map2(colorcnt,:);
    end
    indcnt=indcnt+1;
    colorcnt=colorcnt+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                        % Graphing Mean Reflectivity %
                        
hf=figure;
ha0=axes;
set(ha0,'Xlim',[0 50],'xtick',[0:5:50],'Ylim',[1 20],'Ytick',[1:20]);
hold on;
grid on;
plot(mean_convect,hgt,'k-');
plot(mean_strat,hgt,'c-');
plot(mean_other,hgt,'r-');
plot(mean_tot,hgt,'b-');
 title({[char(NM(snum)),'.',char(ORB(snum)),' ',char(MN(snum)),...
     ' / ',char(DY(snum)),' / ',char(YR(snum)),' ',char(TM(snum)),' UTC'];...
     'Vertical Profile of Mean Reflectivity'});
legend('Convective','Stratiform','Other','Total');
ylabel('Mean Reflectivity (dBZ)');
hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            % Graphing PDF %
                            
hf=figure;
ha=axes;
set(ha,'Xlim',[0 70],'xtick',[0:5:70],'Ylim',[0 35],'Ytick',[0:5:35]);
hold on;
grid on;
plot(dBZbin,PDF_convect,'k-');
plot(dBZbin,PDF_strat,'c-');
plot(dBZbin,PDF_other,'r-');
plot(dBZbin,PDF_tot,'b-');
 title({[char(NM(snum)),'.',char(ORB(snum)),' ',char(MN(snum)),...%
     '/',char(DY(snum)),'/',char(YR(snum)),' ',char(TM(snum)),'UTC'];...
     'Probability Density Function Diagram'});
xlabel('Reflectivity (dBZ)');
ylabel('Frequency (%)');
legend('Convective','Stratiform','Other','Total');
hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                            % Graphing CDF %
                            
hf=figure;
ha=axes;
set(ha,'Xlim',[0 70],'xtick',[0:5:70],'Ylim',[0 100],'Ytick',[0:10:100]);
hold on;
grid on;
plot(dBZbin,convect_CDF,'k-');
plot(dBZbin,strat_CDF,'c-');
plot(dBZbin,other_CDF,'r-');
plot(dBZbin,tot_CDF,'b-');
title({[char(NM(snum)),'.',char(ORB(snum)),' ',char(MN(snum)),...
    '/',char(DY(snum)),'/',char(YR(snum)),' ',char(TM(snum)),'UTC'];...
    'Cumulative Frequency Diagram'});
xlabel('Reflectivity (dBZ)');
ylabel('Frequency (%)');
legend('Convective','Stratiform','Other','Total');
hold off;                           
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                       % Graphing Cumulative CFAD %
                    
run cbcustom_andy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                         % Graphing CFAD %
                    
run cbcustom_CFADandy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%