% One-grid solver of the Lippmann-Schwinger equation
%
% (LS)  w = 1 - g * (V w), 
%
% where g is a Green function. Equation (LS) appears in scattering problems for the Helmholtz
% equation and in the computation of exponentially growing solutions of Faddeev to the Schrodinger 
% equation. 
%
% To use this solver it is needed to write Matlab routines for the Green function g and the 
% multiplicator function V. Names of these routines are given as string arguments. 
% Both functions should take as argument a matrix of complex evaluation points.
% It is also possible to pass one extra parameter as argument to these routines.
%
% This solver is based on GMRES. See Gennadi Vainikko 1997: Fast solvers of the Lippmann-Schwinger 
% equation, in "Direct and inverse problems of mathematical physics", Newark, DE,
% Int.Soc.Anal.Appl.Comput. 5, pp. 423-440, Kluwer acad.publ., Dordrecht 2000
%
% Arguments:
% gname      (string) name of the Green function 
% gparam     parameter passed as argument to the Green function (use [] for no parameter)
% mname      (string) name of the multiplicator function V 
% M          solution is given on a grid of size [2^M, 2^M] 
% s          (positive number) gives the square [-s,s-h]^2 filled by the grid where solution is given
% R          (positive number) Multiplicator function V is supported in D(0,R)
%
% Returns:
% w          (complex) solution of the Lippmann-Schwinger equation
% x1         x1-coordinates of the grid
% x2         x2-coordinates of the grid (example: mesh(x1,x2,real(w)) )
% h          grid parameter
%
% Calls to: GV_grids, GV_LS, gname, mname, gmres, fft2
%
% Samuli Siltanen April 2012

function [w, x1, x2, h] = GVLS_solve(gname, gparam, mname, M, s, R)

% Get the grid points and parameters for only one grid using the two-grid function GV_grids. 
[x1,x2,h,tmp,tmp,tmp] = GV_grids(M, M+1, s);

% Evaluate multiplicator function V
eval(['V = ', mname, '(x1 + i*x2);']);

% Constrct "incident wave" function f=1 as vertical vector
f = ones(size(x1));
f = f(:);

% construct initial guess for GMRES
iniguess = f;

% Evaluate Green function.
% Avoid singularity at origin by setting value there to zero
c         = 2^(M-1)+1; 
z         = (x1+i*x2);
z(c,c)    = 1;
fundevalcommand = ['fund = ', gname, '(z, gparam);'];
eval(fundevalcommand)
fund(c,c) = 0;               

% Smooth truncation of Green's function near the boundary
s  = abs(min(min(x1)));
ep = s-2*R; 
bigind       = abs(z)>=s;
fund(bigind) = 0;
medind       = (abs(z)<s) & (abs(z)>2*R);
fund(medind) = fund(medind).*(1-(abs(z(medind))-2*R)/ep);

% Compute discrete Fourier transform of Green's function
fund    = fftshift(fund);
fundfft = fft2(fund);

% Solve Lippmann-Schwinger equation using gmres
w = gmres('GV_LS', f, 10, 1e-6, 20, [], [], iniguess, fundfft, V, M, h);

% Reshape w to original form
w = reshape(w, 2^M, 2^M);


