function [ result_quality, i_ind, j_ind, edge_buffer] = ringsize( storm_diameter, reverse_dir)
% RINGSIZE Determines the shape of the 'ring' about which circulation is calculated. 

% This function determines the shape of the 'ring' about which 
% circulation is calculated. The shape is a function of the size of the
% 'ring.'  
%
% Input variable
% storm_diameter  The size of the shape expressed in terms of the
%                 number of grid cells across the shape. That is, 
%                 the number of grid cells in the diameter. For example,
%                 if storm_diameter is 2, and the grid spacing is 25km,
%                 then the diameter is geophysical units is 50km.
% reverse_dir     If true, this creates indices in reverse order, so the 
%                 line integration follows the opposite direction.
%                 NOTE: MATLAB requires the indices to appear in the
%                 reverse order from that supplied in the IDL code
%
% Output variables
%  result_quality Set to zero if ringsize runs correctly.
%  i_ind          Change in index for across swath data. If used with data
%                 on a lat/lon grid, this would be the lonitude index. Values
%                 depend on the spatial scale used (storm_diameter).
%  j_ind          As for i_ind, but applies to along swath direction
%  edge_buffer    Number of grid cells from the edge for which vorticity will
%                 not be calculated. 
 
result_quality = 0;  % good output has result_quality = 0
storm_radius = 0.5 * storm_diameter;

switch(storm_diameter)
    case 0    
        fprintf('poor storm diameter: %d \n', storm_diameter);
        result_quality = -1;
        i_ind = [];
        j_ind = [];
        edge_buffer = 0;
        return

    case 1  % Setting vorticity ring size and shape for a 25 km diamter storm:
        i_ind = [0,1,1,0,0];
        j_ind = [0,0,1,1,0];
        edge_buffer = 0;

    case 2  % Setting vorticity ring size and shape for a 50km diameter storm 
        i_ind = [1,0,-1,0,1];
        j_ind = [0,1,0,-1,0];
        edge_buffer = 1;

    case 3  % Setting vorticity ring size and shape for a 75 km diameter storm:
        i_ind = [2,2,1,0,-1,-1,0,1,2];
        j_ind = [0,1,2,2,1,0,-1,-1,0];
        edge_buffer = 2;

    case 4  % Setting vorticity ring size and shape for a 100 km diameter storm:
        i_ind = [2,2,2,1,0,-1,-2,-2,-2,-1,0,1,2];
        j_ind = [-1,0,1,2,2,2,1,0,-1,-2,-2,-2,-1];
        edge_buffer = 2;

    case 5  % Setting vorticity ring size and shape for a 125km diameter storm 
        i_ind = [3,3,3,2,1,0,-1,-2,-2,-2,-2,-1,0,1,2,3,3];
        j_ind = [0,1,2,3,3,3,3,2,1,0,-1,-2,-2,-2,-2,-1,0];
        edge_buffer = 3;

    case 6 % Setting vorticity ring size and shape for a 150km diameter storm 
        i_ind = [3,3,2,1,0,-1,-2,-3,-3,-3,-2,-1,0,1,2,3,3];
        j_ind = [0,1,2,3,3,3,2,1,0,-1,-2,-3,-3,-3,-2,-1,0];
        edge_buffer = 3;

    case 7 % Setting vorticity ring size and shape for a 175km diameter storm 
        i_ind = [4,4,4,3,2,1,0,-1,-2,-3,-3,-3,-3,-2,-1,0,1,2,3,4,4];
        j_ind = [0,1,2,3,4,4,4,4,3,2,1,0,-1,-2,-3,-3,-3,-3,-2,-1,0];
        edge_buffer = 4;

    case 8  % Setting vorticity ring size and shape for a 200km diameter storm 
        i_ind = [4,4,4,3,2,1,0,-1,-2,-3,-4,-4,-4,-4,-4,-3,-2,-1,0,1,2,3,4,4,4];
        j_ind = [0,1,2,3,4,4,4,4,4,3,2,1,0,-1,-2,-3,-4,-4,-4,-4,-4,-3,-2,-1,0];
        edge_buffer = 4;

    case 9  % Setting vorticity ring size and shape for a 225km diameter storm 
        i_ind = [5,5,5,4,3,2,1,0,-1,-2,-3,-4,-4,-4,-4,-3,-2,-1,0,1,2,3,4,5,5];
        j_ind = [0,1,2,3,4,5,5,5,5,4,3,2,1,0,-1,-2,-3,-4,-4,-4,-4,-3,-2,-1,0];
        edge_buffer = 5;

    case 10  % Setting vorticity ring size and shape for a 250km diameter storm 
        i_ind = [5,5,5,4,3,2,1,0,-1,-2,-3,-4,-5,-5,-5,-5,-5,-4,-3,-2,-1,0,1,2,3,4,5,5,5];
        j_ind = [0,1,2,3,4,5,5,5,5,5,4,3,2,1,0,-1,-2,-3,-4,-5,-5,-5,-5,-5,-4,-3,-2,-1,0];
        edge_buffer = 5;

    otherwise % the following code figures out the ring characteristics
        max_pts = ceil(7 * storm_diameter / 2); % overestimate the number of points
                                                % used to define the shape
        i_ind = NaN( max_pts, 1 ); % slightly too large an array, initialized to NaN
        j_ind = i_ind; 
        edge_buffer = floor(( storm_diameter + 1 ) / 2);
        i_ind(1) = edge_buffer;
        j_ind(1) = 0;
        if ( 2 * edge_buffer == storm_diameter )
            i_cntr = 0.0;
            j_cntr = 0.0;
        else
            i_cntr = 0.5;
            j_cntr = 0.5;
        end

        for ind_count = 1: max_pts-1
            change_dir = [ j_cntr - j_ind(ind_count), i_ind(ind_count) - i_cntr ];
            change_dir = change_dir / sqrt(dot(change_dir, change_dir)) % unit vect
            % there are three loosely plausible points
            i_sign = change_dir(1) / abs( change_dir(1) );
            j_sign = change_dir(2) / abs( change_dir(2) );

            if ( change_dir(2) > 0.9238795325113) % sin(67.5 * !DTOR) ) 
                i_nxt_pt = i_ind( ind_count ) + [-1, 0, 1]; 
            elseif ( change_dir(2) > 0.3826834323651) && ... % sin(22.5 * !DTOR) 
                    (change_dir(2) <= 0.9238795325113) % sin(67.5 * !DTOR) )     
                i_nxt_pt = i_ind( ind_count ) + i_sign * [1, 1, 0];  
            elseif ( change_dir(2) > -0.3826834323651) && ... % -sin(22.5 * !DTOR)
                    (change_dir(2) <= 0.3826834323651) % sin(22.5 * !DTOR) )
                i_nxt_pt = i_ind( ind_count ) + i_sign * [1, 1, 1];  
            elseif ( change_dir(2) > -0.9238795325113) && ... % -sin(67.5 * !DTOR)
                    (change_dir(2) <= -0.3826834323651) % -sin(22.5 * !DTOR) )     
                i_nxt_pt = i_ind( ind_count ) + i_sign * [1, 1, 0];  
            else
                i_nxt_pt = i_ind( ind_count ) + [1, 0, -1]; 
            end

            if ( change_dir(1) > 0.9238795325113) % sin(67.5 * !DTOR) )
                j_nxt_pt = j_ind( ind_count ) + [-1, 0, 1];  
            elseif ( change_dir(1) > 0.3826834323651) && ... % sin(22.5 * !DTOR)
                    (change_dir(1) <= 0.9238795325113) % sin(67.5 * !DTOR) )  
                j_nxt_pt = j_ind( ind_count ) + j_sign * [0, 1, 1];  
            elseif ( change_dir(1) > -0.3826834323651) && ... % -sin(22.5 * !DTOR)
                    (change_dir(1) <= 0.3826834323651) % sin(22.5 * !DTOR) )
                j_nxt_pt = j_ind( ind_count ) + j_sign * [1, 1, 1]; 
            elseif ( change_dir(1) > -0.9238795325113) && ... % -sin(67.5 * !DTOR)
                    (change_dir(1) <= -0.3826834323651) % -sin(22.5 * !DTOR) )
                j_nxt_pt = j_ind( ind_count ) + j_sign * [0, 1, 1];
            else
                j_nxt_pt = j_ind( ind_count ) + [1, 0, -1];
            end

            % Choose the point that is closest to the correct storm radius 
            % away from the center location (IDL code had small bug here -- 
            % it missed subtracting the center point.  No real problem for large 
            % diameters but this fixed code works even for small diameters.)
            test_range = ( (i_nxt_pt - i_cntr).^2 + (j_nxt_pt - j_cntr).^2 - (storm_radius^2*[1 1 1]) ).^2;
            test_range_min = min(test_range);
            for i=1:3
                if(test_range(i) == test_range_min)
                    best_index = i; % index for best point (where test_range is at minimum)
                    break
                end
            end

            i_ind( ind_count + 1 ) = i_nxt_pt( best_index );
            j_ind( ind_count + 1 ) = j_nxt_pt( best_index ); 

        end

        %    truncate the arrays to the actual number of points (plus one) used to
        %    define the shape about which the circulation will be calcuated. Note that
        %    the last point is set equal to the first point.
        
        % IDL: i_test = where( i_ind EQ edge_buffer AND j_ind EQ 0 )
        % ALTERNATE MATLAB:  use logical vector times index vector
        best_index = max((1:max_pts)' .* ((i_ind == edge_buffer) .* (j_ind == 0)));
        i_ind = i_ind(1:best_index);
        j_ind = j_ind(1:best_index);

end % end of switch

if reverse_dir
    i_ind = flip(i_ind);
    j_ind = flip(j_ind);
end

% fprintf('Done in ringsize! Thank you!');

end