function U = FastPoissonSolve(F,hx,hy)
% F is n1-by-n2
% hx and hy are positive real scalars
% U is n1-by-n2 such that
% (1/hy^2)*T(n1)*U + U*T(n2)*(1/hy^2) = F (1)
%

[n1,n2] = size(F);

% Refer to the subfunctions T(n), S(n), and d(n) in Show121Props
% Substitute
% T(n1) = S(n1)*D(n1)*S(n1)*(2/(n1+1))
% and
% T(n2) = S(n2)*D(n2)*S(n2)*(2/(n2+1))
%
% into (1) and simplify. This gives
%
% (1/hy^2)D(n1)*Utilde + Utilde*D(n2)(1/hx^2) = Ftilde (2)
% where
% Utilde = S(n1)*U*S(n2)
% and
% Ftilde = S(n1)*(F*S(n2)) = S(n1)*(S(n2)*F')'

% Compute the diagonals of D(n1) and D(n2)

dn1 = d(n1);
dn2 = d(n2);

% Compare (k,j) entries in (2):
%
% (dn1(k)/hy^2)*Utilde(k,j) + Utilde(k,j)*(dn2(j)/hx^2) = Ftilde(k,j)
%
% Thus, Utilde(k,j) = Ftilde(k,j)/(dn1(k)/hy^2 + dn2(j)/hx^2);

% We overwrite F with Utilde...

F = FastSine(FastSine(F')'); % Ftilde

% The vectorized version of

%for k=1:n1
% for j=1:n2
% F(k,j) = F(k,j)/( dn1(k)/hy^2 + dn2(j)/hx^2 );
% end
%end

% is this

F = F./((dn1/hy^2)*ones(1,n2)+ones(n1,1)*(dn2/hx^2)');

% Since Utilde = S(n1)*U*S(n2) it follows that
%
% U = (4/((n1+1)*(n2+1)))*S(n1)*Utilde*S(n2)

U = (4/((n1+1)*(n2+1)))*FastSine(FastSine(F')');

function dn = d(n)
% A column n-vector with the property that inv(S(n))*T(n)*S(n) = diag(dn)
dn = 4*sin((1:n)'*pi/(2*n+2)).^2;

function Y = FastSine(X)
% X is n1-by-n2
% Y = S(n1)*X where S(n1) i sthe sine transform matrix.
[n1,n2] = size(X);
Y = fft([zeros(1,n2);X;zeros(1,n2);-X(n1:-1:1,:)]);
Y = imag(Y(2:1+n1,:))/(-2);