cnfft2D(varargin)

function w = cnfft2D(varargin)
% crankNicolson: uses Crank-Nicolson algorithm to approximate the solution
%   to the parabolic PDE (in 2-D):
%       u_{t}(x,t) - alpha^2 u_{xx}(x,t) = 0, 0<x<L, 0<t<T
%   subject to the boundary conditions
%       u(0,t) = u(L,t) = 0, 0<t<T
%   and the initial conditions
%       u(x,0) = f(x), 0<=x<=L
%
% arguments:
%   L (vector 1by2) - upper bound of spatial (x) variable
%       (Default L = [1,1])
%   T (scalar) - upper bound of time (t) variable
%       (Default T = 1)
%   alpha (scalar) - square root of coefficient of u_{xx} term
%       (Default alpha = 1)
%   m (vector 1by2) - number of discrete spatial intervals
%       (Default m = [100,100])
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (m+1 x m+1 x n) - approximation to u(x,t) at discrete space/time positions
%

% author: Troy J. Winkstern
% email: [e-mail]
% date: 30 Jan 2011

% parse input arguments
[L,T,alpha,m,n] = parseInputs(varargin{:});
L_x = L(1);
L_y = L(2);
m_x = m(1);
m_y = m(2);

% initialize h, k, lambda, and w
h_x = L_x/m_x;
h_y = L_y/m_y;
k = T./n;
lambda_x = (alpha.^2).*k./(h_x.^2);
lambda_y = (alpha.^2).*k./(h_y.^2);
w = zeros(m_y+1,m_x+1,n);

% initialize w (change this to sin(pix1/l)*sin(pix2/l)
for i=2:m_x
    for j=2:m_y
        w(j,i,1) = sin(pi.*h_x.*(i-1)./L_x).*sin(pi.*h_y.*(j-1)./L_y);
    end
end


for t=1:(n-1)
    y = zeros(m_y+1,m_x+1);
    
    %Computes second column of the matrix
    y(2,2) = (1-lambda_x.*lambda_y)*w(2,2,t) + ...
        (lambda_y/2)*w(3,2,t)+ ...
        (lambda_x/2)*w(2,3,t);
    
    for j = 3:(m_y-1)
        y(j,2) = (lambda_y/2)*w(j-1,2,t) + (1-lambda_x.*lambda_y)*w(j,2,t) + ...
            (lambda_y/2)*w(j+1,2,t)+ ...
            (lambda_x/2)*w(j,3,t);
    end
    
    y(m_y,2) = (lambda_y/2)*w(m_y-1,2,t) + (1-lambda_x.*lambda_y)*w(m_y,2,t) + ...
        (lambda_x/2)*w(m_y,3,t);
    
    %Computes the inner parts of the matrix
    
    for i = 3:(m_x-1)
        y(2,i) = (1-lambda_x.*lambda_y)*w(2,i,t) + ...
            (lambda_y/2)*w(3,i,t)+ (lambda_x/2)*w(2,i-1,t) + ...
            (lambda_x/2)*w(2,i+1,t);
        for j = 3:(m_y-1)
            y(j,i) = (lambda_y/2)*w(j-1,i,t) + (1-lambda_x.*lambda_y)*w(j,i,t) + ...
                (lambda_y/2)*w(j+1,i,t)+ (lambda_x/2)*w(j,i-1,t) + ...
                (lambda_x/2)*w(j,i+1,t);
        end
        y(m_y,i) = (lambda_y/2)*w(m_y-1,i,t) + (1-lambda_x.*lambda_y)*w(m_y,i,t) + ...
            (lambda_x/2)*w(m_y,i-1,t) + ...
            (lambda_x/2)*w(m_y,i+1,t);
    end
    
    %Computes the m_x column of the matrix
    
    y(2,m_x) = (1-lambda_x.*lambda_y)*w(2,m_x,t) + ...
        (lambda_y/2)*w(3,m_x,t)+ (lambda_x/2)*w(2,m_x-1,t);
    
    for j = 3:(m_y-1)
        y(j,m_x) = (lambda_y/2)*w(j-1,m_x,t) + (1-lambda_x.*lambda_y)*w(j,m_x,t) + ...
            (lambda_y/2)*w(j+1,m_x,t)+ (lambda_x/2)*w(j,m_x-1,t);
    end
    
    y(m_y,m_x) = (lambda_y/2)*w(m_y-1,m_x,t) + (1-lambda_x.*lambda_y)*w(m_y,m_x,t) + ...
        (lambda_x/2)*w(m_y,m_x-1,t);
    
   
   % went out to length 2m
   X = sqrt(2./m_y).*imag(fft(y, 2*m_y,1));
   X = X(1:m_y+1,:);
   
   Y = sqrt(2./m_x).*imag(fft(X, 2*m_x,2));
   Y = Y(:,1:m_x+1);
   
   % angles divided by m
   a_x = (pi./m_x).*(0:m_x);
   a_y = (pi./m_y).*(0:m_y).';
   
   % changed lambda to accomodate for 2-D
   c = 1 + (lambda_x.*lambda_y).*(1-cos(a_y))*(1-cos(a_x));

   % went out to length 2m
   Z = sqrt(2./m_y).*imag(fft(Y./c, 2*m_y,1));
   Z = Z(1:m_y+1,:);
   
   R = sqrt(2./m_x).*imag(fft(Z, 2*m_x,2));
   w(2:m_y,2:m_x,t+1) = R(2:m_y,2:m_x);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction parseInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,T,alpha,m,n] = parseInputs(varargin)

% check number of input arguments
nargs = length(varargin);
error(nargchk(0,5,nargs));

% get/set n
if nargs<5
    n = [];
else
    n = varargin{5};
end
if isempty(n)
    n = 100;
end

% check n has valid value
if (numel(n)>1) || (n<1) || ~isequal(round(n),n)
    error('n must be a positive integer.');
end

% get/set m
if nargs<4
    m = [];
else
    m = varargin{4};
end
if isempty(m)
    m = [100,100];
end

% check m has valid value
if (numel(m)~=2) || any(m<1) || any(~isequal(round(m),m))
    error('m must be a vector of positive integers.');
end

% get/set alpha
if nargs<3
    alpha = [];
else
    alpha = varargin{3};
end
if isempty(alpha)
    alpha = 1;
end

% check alpha has valid value
if numel(alpha)>1
    error('alpha must be a real scalar.');
end

% get/set T
if nargs<2
    T = [];
else
    T = varargin{2};
end
if isempty(T)
    T = 1;
end

% check T has valid value
if (numel(T)>1) || (T<=0)
    error('T must be a positive real scalar.');
end

% get/set L
if nargs<1
    L = [];
else
    L = varargin{1};
end
if isempty(L)
    L = [1,1];
end

% check L has valid value
if (numel(L)~=2) || any(L<=0)
    error('L must be a positive vector with 2 elements.');
end

% warning if L not integer
if any(~isequal(round(L),L))
    warning('L should be an vector of integers so that the inital condition takes on the value of 0 for t = 0 and t = T.');
end