% Create initial variables

mindistance = 45;                  % Minimum distance in kilometers
maxdistance = 90;                  % Maximum distance in kilometers
elevationcut = 1;                  % Elevation cut

lowerz = 20;                        % Lower reflectivity range in dBz
upperz = 60;                        % Upper reflectivity range in dBz
rhomin = 0.98;                      % CC threshold
        
max = (upperz - lowerz) * 2 + 1;    % Another counter

zdraverage = zeros(1,max);          % Average ZDR Calculator
reflectivity = lowerz:0.5:upperz;   % Reflectivity range
rowcounter = 0;                     % ZDR Row counter
averagesum = zeros(1,max);          % 20-30 dBz ZDR Average

[cut, mindistance, maxdistance, elevationcut, lowerz, upperz, rhomin]...
    = inputchecker(cut, mindistance, maxdistance, elevationcut, lowerz,...
    upperz, rhomin);

% Apply mask based on CC (rho)

rho_mask(1:1192,1:360,1) = zeros(1192,360);
z_mask(1:1192,1:360,1) = zeros(1192,360);

rho_mask(mindistance*4:maxdistance*4,1:360,1) = ...
    (cut{1,elevationcut}.rho(mindistance*4:maxdistance*4,1:360) >= rhomin)...
    & (cut{1,elevationcut}.rho(mindistance*4:maxdistance*4,1:360) <= 1);

% Apply mask based on distance

z_mask(mindistance*4:maxdistance*4,1:360,1) = ...
    cut{1,elevationcut}.z(mindistance*4:maxdistance*4,1:360) >= ...
    lowerz & cut{1,elevationcut}.z(mindistance*4:maxdistance*4,1:360) ...
    <= upperz;

% Apply large mask

valid_mask(:,:,1) = and(z_mask,rho_mask);

% Create reflectivity array

z1 = valid_mask .* cut{1,elevationcut}.z(1:1192,1:360);

% Create ZDR array

zdr1 = valid_mask .* cut{1,elevationcut}.zdr(1:1192,1:360);

% Filter out zero values

outlier = z1 == 0;
z1(outlier) = NaN;

outlier = zdr1 == 0;
zdr1(outlier) = NaN;

% Calculate average

for z=lowerz:0.5:upperz
   
    zdrsum = 0;
    zdrcounter = 0;
    rowcounter = rowcounter + 1;
    
    for col=1:1:1192
        
        for row = 1:1:360
            
            if((z1(col,row) == z) && (~isnan(zdr1(col,row))))
                zdrsum = zdr1(col,row) + zdrsum;
                zdrcounter = zdrcounter + 1;
            end
            
        end
        
    end
    
    if(zdrcounter > 0)
        zdraverage(rowcounter) = zdrsum / zdrcounter;
    end
    
end

% Average of the averages

for i=1:1:max
    
    averagesum = zdraverage(i) + averagesum;
    
end

avgaverage = averagesum / max;

%directory - assumed the folder is previously selected in the GUI.
datest = sprintf('%s',cut{1,elevationcut}.Date);
datestrng = strrep(datest, '/','');
timest = sprintf('%s',cut{1,elevationcut}.Time);
timestr = strrep(timest,':','');
timestrng = timestr(1:6);
       
p1 = plot(reflectivity,zdraverage,'-r',reflectivity,avgaverage,'-k',...
    'LineWidth',2);

hold on

p2 = plot(z1,zdr1,'o','MarkerFaceColor','g','MarkerEdgeColor','k', ...
    'MarkerSize',5);

text((lowerz+upperz)/2-0.25,6,['ZDR Average ',num2str(avgaverage(1))]);
xlabel('Reflectivity (dBZ)');
ylabel('Differential Reflectivity (dB)');
title(['ZDR Average for reflectivity values of ',...
    num2str(lowerz), ' to ', num2str(upperz), ' dBZ']);

hold off

legend('ZDR Average','Cut ZDR Average');
grid on
grid(gca,'minor')