function [ vort, lat_vort, lon_vort, div, lat_div, lon_div ] = vortcalc( lat, lon, u_wind, v_wind, i_ind, j_ind, edge_buffer )
% VORTCALC calculates vorticity and divergence using the circulation theorem
%
% Reference: Bourassa, M. A., and K. McBeth-Ford, 2010: Uncertainty in scatterometer-derived vorticity, 
% J. Atmos. Oceanic Technol. 27, 594 – 603. DOI: 10.1175/2009JTECHO689.1
%
% Based on original IDL code from Kelly McBeth-Ford, as updated by Heather Holbach (10/23/13;1/13/14)
% The IDL update includes a spline fit instead of a linear interpolation.
% Reference: Holbach, Heather M., Mark A. Bourassa, 2014: 
% The Effects of Gap-Wind-Induced Vorticity, the Monsoon Trough, and the ITCZ on East Pacific Tropical Cyclogenesis.
% Mon. Wea. Rev., 142, 1312–1325. doi: http://dx.doi.org/10.1175/MWR-D-13-00218.1
%
% This MATLAB translation is pretty literal, but some changes were made to
% better follow MATLAB conventions.  Note that MATLAB arrays are 1-based
% whereas IDL arrays are 0-based.  So, where an IDL index is also used as a
% count, separate alterations must be implemented in MATLAB regarding where
% it is used as an index vs where it is used as a count.
%
% This function is meant for computing the vorticity and divergence for a swath.
%
% Input variables
%  lat           latitude in a two dimensional array
%  lon           longitude in a two dimensional array
%  u_wind        Zonal component of the wind, in 2D array
%  v_wind        Meridional component of the wind, in 2D array   
%  rain_flag     quality indicator for wind data
%  NOTE: rain_flag wasn't used in IDL code and doesn't appear here
%  i_ind         change in index for across swath data. If used with data 
%                on a lat/lon grid, this would be the longitude index. Values
%                depend on the spatial scale used (storm size) and are set
%                in the subroutine ringsize.pro
%  j_ind         as for i_ind, but applied to along swath direction
%  edge_buffer   number of grid cells from the edge for which vorticity will
%                not be calculated. This variable is output from the
%                subroutine ringsize.pro
%  vort_thresh   threshold vorticity used to identify locations of interest.
%  NOTE: vort_thresh wasn't used in IDL code and doesn't appear here
%  miss          value used for missing data (e.g., -9999.)
%  NOTE: In MATLAB we're just using NaN for missing or bad values
%
% Output variables
%  vort          output vorticity in a two dimensional matrix
%  lat_vort      latitudes corresponding to the vorticity points
%  lon_vort      longitudes corresponding to the vorticity points
%  div          output divergence in a two dimensional matrix
%  lat_div      latitudes corresponding to the divergence points
%  lon_div      longitudes corresponding to the divergence 
%
% Other variables used within the subroutine
%  n_across      Number of data points across the swath
%  n_along       Number of data points along the swath
%  dist          Conversion factor from degrees to radians
%  n_ring        Number of line segments plus one in the shape used in the
%                calculation of circulation - IFF there are no flagged data
%                points.
%  pt_count      the actual number of segments plus one in the shape used 
%                in the calculation of circulation
%  u_wind_good   zonal wind component of points used to define the shape
%                used in the calculation of circulation
%  v_wind_good   meridional wind component of points used to define the shape
%                used in the calculation of circulation
%  x_loop        zonal distance from the center of mass of the points 
%                defining the perimeter of the shape. 
%  y_loop        meridional distance from the center of mass of the points 
%                defining the perimeter of the shape.
%  dx_loop       zonal change in position of adjacent points defining the 
%                perimeter of the shape. A flat earth assumption is applied,
%                the the mapping characteristics of the center of mass of 
%                the shape.
%  dy_loop       as for dx_loop, except applies to meridional changes in 
%                position
%  lat_sp        latitude for spline fit points
%  lon_sp        longitude for spline fit points
%  x_loop_sp     same as x_loop except for spline fit points
%  y_loop_sp     same as y_loop except for spline fit points
%  dx_loop_sp    same as dx_loop except for spline fit points
%  dy_loop_sp    same as dy_loop except for spline fit points
%  du_loop_sp    change in the zonal wind component of adjacent spline fit points
%                used for computing the circulation of individual segments.
%  dv_loop_sp    change in the meridional wind component of adjacent spline fit points.
%                used for computing the circulation of individual segments.
%  circ_vort     calculated circulation about one shape
%  circ_div      calculated circulation about one shape
%  area          calculated area for one shape

% SKIPPING THIS IDL CODE SINCE NOT NEEDED
%   flag to make sure that there is data in the u and v wind arrays
%   IF (N_ELEMENTS(size(u_wind)) EQ 5 AND N_ELEMENTS(size(v_wind)) EQ 5) THEN BEGIN

%   i_bad = WHERE( u_wind LT miss + 1 OR v_wind LT miss + 1, n_bad )
%   IF ( n_bad GT 0 ) THEN BEGIN
%     u_wind (i_bad) = miss  ; replace bad values with the missing value
%     v_wind (i_bad) = miss  ; replace bad values with the missing value
%   ENDIF
%
%   Equivalent MATLAB code, but not needed since we already have NaNs for bad or missing values
%     bad_values = logical((u_wind < miss+1) + isnan(u_wind) + (v_wind < miss+1) + isnan(v_wind));
%     u_wind(bad_values) = NaN; % miss; % put in missing value wherever bad values were found
%     v_wind(bad_values) = NaN; % miss;

% find the dimensions of the swath
n_across = size(lat, 1);
n_along = size(lat, 2);

n_ring = size(i_ind, 2); % number of elements in i_ind

% allocate arrays
u_loop = NaN(n_ring-1, 1);
v_loop = NaN(n_ring-1, 1);

dx_loop = NaN(n_ring-1, 1);
dy_loop = NaN(n_ring-1, 1);

x_loop = NaN(n_ring, 1);
y_loop = NaN(n_ring, 1);

vort = NaN(n_across, n_along); % initialize with missing value
lat_vort = NaN(n_across, n_along);
lon_vort = NaN(n_across, n_along);

div = NaN(n_across, n_along);
lat_div = NaN(n_across, n_along);
lon_div = NaN(n_across, n_along);

lat_good = NaN(n_ring, 1);
lon_good = NaN(n_ring, 1);

u_wind_good = NaN(n_ring, 1);
v_wind_good = NaN(n_ring, 1);

%{
  edge_buffer, i_ind and j_ind are determined based on the size of the 
  vorticity 'ring'. edge_buffer prevents problems at the edge of a swath 
  and i_ind/j_ind represent the size and shape of the 'rings'.
%}

% pass in edge buffer variable set in case statements
% ensuring that there is no missing data in the small 'ring':

for i = edge_buffer + 1 : n_across - 1 - edge_buffer
    for j = edge_buffer + 1 : n_along - 1 - edge_buffer
        
        pt_count = 0;
        
        % IDL code sets array values above to NaN, but we already did this
        
        % Find all of the points along the 'ring' that have no missing values
        for k = 1 : n_ring - 1
            if ~isnan( u_wind( i + i_ind(k), j + j_ind(k) )) && ...
             ~isnan( v_wind( i + i_ind(k), j + j_ind(k) )) && ...
             ~isnan( lat( i + i_ind(k), j + j_ind(k) )) && ...
             ~isnan( lon( i + i_ind(k), j + j_ind(k) )) 

                lat_good(pt_count + 1) = lat( i + i_ind(k), j + j_ind(k) );
                lon_good(pt_count + 1) = lon( i + i_ind(k), j + j_ind(k) );
                u_wind_good(pt_count + 1) = u_wind( i + i_ind(k), j + j_ind(k) );
                v_wind_good(pt_count + 1) = v_wind( i + i_ind(k), j + j_ind(k) );
                
                % Correct longitudes to account for crossing the prime meridian
                % adjust this for different applications
                % CHECK: DB I don't think we want to do this for our use in MATLAB
                %{
                if lon_good(pt_count + 1) >= 0.0 && lon_good(pt_count + 1) < 60.0 
                    lon_good(pt_count + 1) = lon_good(pt_count + 1) + 360.0;
                end
                %}
                pt_count = pt_count + 1;
            end
        end
        
        % pt_count now equals the number of elements in the "good" arrays. 
        
        % Ensure that the starting point equals the finishing point of the 'ring':
        lat_good(pt_count + 1) = lat_good(1);
        lon_good(pt_count + 1) = lon_good(1);
        u_wind_good(pt_count + 1) = u_wind_good(1);
        v_wind_good(pt_count + 1) = v_wind_good(1);

        dist = 111300.0; % one degree of latitude in meters
        
        if max(lon_good) - min(lon_good) > 60
            for idx=(1:pt_count+1)
                if lon_good(idx) < 180
                    lon_good(idx) = lon_good(idx) + 360;
                end
            end
        end
        
        % Calculating dx and dy in meters
        % Calculating area, vorticity, and divergence only when the number of points is greater than 2

        if ( pt_count >= 0.75 * n_ring && pt_count == 4 ) || pt_count >= 0.8 * n_ring 
 
            lat_avg = mean(lat_good(1:pt_count)); 
            lon_avg = mean(lon_good(1:pt_count));

            for k = 1 : pt_count
                x_loop(k) = (lon_good(k)-lon_avg) * dist * cos( deg2rad(0.5 * (lat_good(k) + lat_avg)));
                y_loop(k) = (lat_good(k) - lat_avg) * dist;
            end
            
            x_loop(pt_count + 1) = x_loop(1);
            y_loop(pt_count + 1) = y_loop(1);
            
            % Calculating area
            if pt_count > 3
                area = x_loop(1) * ( y_loop(2) - y_loop(1) ) + x_loop(1) * ...
                 ( y_loop(1) - y_loop( pt_count ) );
                for m = 2 : pt_count
                  area = area + x_loop(m) * ( y_loop(m + 1) - y_loop(m - 1) );
                end
            else
                area = ( x_loop(2) - x_loop(1) ) * ...
                    ( y_loop(3) - y_loop(1) ) - ( x_loop(3) - x_loop(1) ) * ...
                    ( y_loop(2) - y_loop(1) );
            end
            
            area_final = abs(0.5 * area);
            
            circ_vort = 0.0;    % initialize circulation for vorticity
            circ_div = 0.0;     % initialize circulation for divergence
            
            if area_final > 0.0

                % Added spline fit for more accurate vorticity calculation
                % Applies a spline fit to obtain u_wind and v_wind in between data points along the 'ring'

                [u_loop, v_loop, nn] = spline_db(u_wind_good(1:pt_count+1),v_wind_good(1:pt_count+1));
                % nn is an array containing the number of points used for each segment in the spline fit

                % Create arrays for spline fit points (lat, lon, x and y in meters, dx and dy in meters, du and dv)					
                total_fit_points = sum(nn);   % IDL: TOTAL(nn, /INTEGER)
                lat_sp = NaN(total_fit_points + 1, 1);
                lon_sp = NaN(total_fit_points + 1, 1);
                x_loop_sp = NaN(total_fit_points + 1, 1);
                y_loop_sp = NaN(total_fit_points + 1, 1);
                dx_loop_sp = NaN(total_fit_points, 1);
                dy_loop_sp = NaN(total_fit_points, 1);
                du_loop_sp = NaN(total_fit_points, 1);
                dv_loop_sp = NaN(total_fit_points, 1);

                lat_sp(1) = lat_good(1);
                lon_sp(1) = lon_good(1);

                for k = 1 : pt_count
                    % compute lat and lon for spline points in between known points
                    if k == 1		
                        for h = 2 : nn(k)+1
                            lat_sp(h) = lat_good(k)+(h-1)*((lat_good(k+1)-lat_good(k))/nn(k));
                            lon_sp(h) = lon_good(k)+(h-1)*((lon_good(k+1)-lon_good(k))/nn(k));
                        end
                    elseif k > 1
                        for h = sum(nn(1:k-1)) : sum(nn(1:k))
                            h_idx = h + 1;
                            lat_sp(h_idx) = lat_good(k)+(h-sum(nn(1:k-1)))*((lat_good(k+1)-lat_good(k))/nn(k));
                            lon_sp(h_idx) = lon_good(k)+(h-sum(nn(1:k-1)))*((lon_good(k+1)-lon_good(k))/nn(k));
                        end
                    end
                end

                % compute the zonal and meridional distances from the center of mass for all spline fit points
                for cnt1 = 1 : total_fit_points + 1
                    x_loop_sp(cnt1) = ( lon_sp(cnt1) - lon_avg ) * dist * ...
                        cos( deg2rad( 0.5 * ( lat_sp(cnt1) + lat_avg ) ));  
                    y_loop_sp(cnt1) = ( lat_sp(cnt1) - lat_avg ) * dist;
                end

                % compute the change in the zonal and meridional distances (dx and dy) and 
                % wind components (du and dv) for adjacent spline fit points
                for cnt2 = 1 : total_fit_points
                    dx_loop_sp(cnt2) = x_loop_sp(cnt2+1) - x_loop_sp(cnt2); 
                    dy_loop_sp(cnt2) = y_loop_sp(cnt2+1) - y_loop_sp(cnt2);
                    du_loop_sp(cnt2) = u_loop(cnt2+1) + u_loop(cnt2);   
                    dv_loop_sp(cnt2) = v_loop(cnt2+1) + v_loop(cnt2);
                    % Compute and sum circulation about the entire 'ring'
                    % Circulation is computed for each segment
                    circ_vort = circ_vort + (0.5*( du_loop_sp(cnt2) * dx_loop_sp(cnt2) + dv_loop_sp(cnt2) * dy_loop_sp(cnt2) ));
                    circ_div = circ_div + (0.5*( -dv_loop_sp(cnt2) * dx_loop_sp(cnt2) + du_loop_sp(cnt2) * dy_loop_sp(cnt2) ));
                end
			
                % End of addition for spline fit
                
                % Compute vorticity and divergence
                vort(i,j) = circ_vort/area_final;
                lat_vort(i,j) = lat_avg;
                lon_vort(i,j) = lon_avg;
                div(i,j) = circ_div/area_final;
                lat_div(i,j) = lat_avg;
                lon_div(i,j) = lon_avg;

            else
                % In case area_final equals 0.0, vort and div cannot be calculated
                % CHECK: DB NOT REALLY NEEDED since we initiallize to NaN
                vort(i,j) = NaN;
                lat_vort(i,j) = NaN;
                lon_vort(i,j) = NaN;
                div(i,j) = NaN;
                lat_div(i,j) = NaN;
                lon_div(i,j) = NaN;
            end

        else
            % In case there are too few points, area cannot be calculated:
            % CHECK: DB NOT REALLY NEEDED since we initiallize to NaN
            vort(i,j) = NaN;
            lat_vort(i,j) = NaN;
            lon_vort(i,j) = NaN;
            div(i,j) = NaN;
            lat_div(i,j) = NaN;
            lon_div(i,j) = NaN;
        end
    end
end

end % end of function