function [filterout,vrncout,ampspec,ampspecrat] = get_wavenumber(varin,numwave)

% Usage: [filterout,vrncout,ampspec,ampspecrat] = get_wavenumber(varin,numwave);
%
% Written by Sylvie Lorsolo and Altug Aksoy, August 2010
%
% This function performs a wavenumber analysis, through FFT, on a given 2-d field
% in polar coordinates. Use get_polar_var to obtain the field. Wavenumber 0
% corresponds to the axisymmetric (i.e., azimuthal mean) structure as a function
% of radial distance from the center.
%
% Inputs:
% varin:      2-d variable to be analyzed converted to polar coordinates
%               -> radial direction in first dimension, azimuthal direction
%                  second dimension (i.e., azimuth in columns)
% numwave:    Maximum wavenumber to be analyzed and returned
%
% Output:
% filterout:  Filtered out data for each wavenumber (starting from 0), in cell format.
%             Last element will contain all of >= numwave+1 accumulated.
% vrncout:    Variance explained by each wavenumber component.
% ampspec:    Amplitude spectrum as a function of radial distance and wavenumber.
%             Last element will contain all of >= numwave+1 accumulated.
% ampspecrat: Ratio of amplitude spectrum relative to total amplitude as a
%             function of radial distance and wavenumber.


filterout  = [];
vrncout    = [];
ampspec    = [];
ampspecrat = [];

if (nargin ~= 2)
    warning('Following error occurred:')
    disp(['This function expects 2 input arguments as follows:                        ';...
          '    [filterout,vrncout,ampspec,ampspecrat] = get_wavenumber(varin,numwave);'])
    return
end

% Actual index of last wavenumber
n = numwave + 1;

% Get size in radial and azimuth directions
[nrad,naz] = size(varin);

% Number of FFT components in the azimuth direction
N = ceil((naz+1)/2);

% Detrend field
%varin = varin - mean(mean(varin(~isnan(varin))));

% Compute FFT in azimuth
F = (fft(varin')).';

% Obtain spectrum
%   - Requested part, 1 to numwave
ampspec        = abs(F(:,1:n))/naz;
%   - Sum of rest
ampspec(:,n+1) = sum(abs(F(:,n+1:N)),2)/naz;

% Obtain percent of total amplitude spectrum
ampspecrat     = ampspec ./ repmat(sum(ampspec,2),[1 (n+1)]);

% Compute individual filtered-out fields
for i = 1:(n+1)
    FF = F;
    if (i == (n+1))
        % Special case of all wavenumbers greater than or equal to (n+1)
        % accumulated
        if (i > 2)
            FF(:,1:i-1)           = 0;
            FF(:,end-(i-1)+2:end) = 0;
        else
            FF(:,1:i-1)           = 0;
        end
    else
        % Retain current wavenumber amplitude and set all else to zero
        if (i > 1)
            % Wavenumbers 1 and greater - need to treat two sides of the
            % spectrum separately
            FF(:,1:i-1)           = 0;
            FF(:,i+1:end-(i-1))   = 0;
            FF(:,end-(i-1)+2:end) = 0;
        else
            % Wavenumber 0 - Just keep first
            FF(:,i+1:end)         = 0;
        end
    end
    % Compute inverse FFT to get to filtered field at current wavenumber
    filterout{i} = real(ifft(FF.').');
    % Compute variance explained
    cc = corrcoef_nonan(filterout{i},varin);
    vrncout(i) = cc.^2;
end
