Modified_Right(f,m,N)

%
% Written by M. Harper Langston - 5/10/00
% [e-mail]
%
% Right-hand vector is modified according to whether you're using
% the 5,9, or 9-point modified scheme.  This is called bt Form_Right.m
% which is called by poisson.m  If you are using the 5 or 9-point scheme,
% the result is straightforward, merely throwing away the boundary values
% in f.
% If you are using the modified nine-point scheme (corresponding to N=10),
% then the right hand side has a five-point stencil applied to it as per 
% Iserles, page 127.
function [v] = Modified_Right(f,m,N)
% Form the blocks, S, T and Z
if N==5 | N==9
   f = f(2:m+1,2:m+1);
   v = f(:);
   % The modified 9-point scheme can be computed quickly using the Gamma
   % matrix found from before and then using fast sine transforms.
elseif N==10 % Modified step for right side.  See Iserles, page 127
   % Compute Gamma for the modified scheme.  I let N = 101 here.  
   % Here,S = [1/12 2/3 1/12], T = [1/12] in a five-point manner.  See Form_Gamma.m
   % for more info.
   [Gam] = Form_Gamma(m,101);
   % Remove boundaries of f since we will apply a five-point stencil to it
   fbot = f(1,1:end); fleft = f(1:end,1); ftop = f(end,1:end); fright = f(1:end,end);
   % After removing boundaries, make f one long vector
   f = f(2:m+1,2:m+1);   f = f(:);   len = length(f);
   % As in 5-point stencil, make a new vector of the boundary points
   fnew = fleft(2:length(fleft)-1);
	fnew(len-m+1:len) = fright(2:length(fright)-1);
	for k = 1:m
  		fnew((k-1)*m + 1) = fnew((k-1)*m + 1) + fbot(k+1);
     	fnew(m*k) = fnew(m*k) + ftop(k+1);
   end
c = Transform(m,f)';
% Having formed c, rearrange rows into columns, so have to convert c to
% matrix format.  Only takes O(m) time which is negligible
for j = 1:m
   cnew(1:m,j) = c((j-1)*m + 1:(j-1)*m + m);
end
c = cnew';   c = c(:);
% Gam is tridiagonal and stored as sparse, so Gam*c is cheap.
t = Gam*c;
% The resulting vector t must be rearranged back so columns are rows: O(m)
for j = 1:m
   tnew(1:m,j) = t((j-1)*m + 1:(j-1)*m + m);
end
t = tnew';   t = t(:);
% Call the DF sine Transform routine for speed to get the final result,u:
v = Transform(m,t)';
v = v + (1/12).*fnew;
else
   error('Must use 9 or 5 point scheme (N = 5 ot N=9 in poisson.m)')
end
%
% Written by M. Harper Langston - 5/10/00
% [e-mail]
%