function y = BirmanSchwinger(x,type,beta,lambda,L)

% BirmanSchwinger.m - Application of the NLS BS operator Kminus and Kplus in R^3
%
% In :  x : input vector N^3-by-1 which will be reshaped into a N-by-N-by-N array
%       type : 'pluss' for Kplus, 'minus' for Kminus
%       beta : exponent in soliton's equation is 2*beta+1
%       lambda : parameter of Kplus and Kminus
%       L : truncation in x to a cube of sidelength L
%
% Out : y = K(x)
%
% By Laurent Demanet, 2005

global soliton;
global Lkernel;

N3 = length(x);
N = round((length(x))^(1/3));
if N3 == N^3,
    x = reshape(x,N,N,N);
else
    error('The length of x should be the cube of some integer');
end;
if nargin < 2, type = 'pluss'; end;
if nargin < 3, beta = 1; end;
if nargin < 4, lambda = 0; end;
if nargin < 5, L = 15; end;

center = floor(N/2 + 1);

if type == 'pluss',
    y = (sqrt(2*beta+1)*soliton.^(beta)).*x;
else
    y = (soliton.^(beta)).*x;
end
y = fftshift(fftn(ifftshift(y)));

% Convolution kernel as fundamental solution of the Laplace or Helmholtz equation

Klambda = (L/N)^3 * Lkernel;
Klambda(center,center,center) = L^2/N^2 * 1/(2*pi);
Klambda = real(fftshift(fftn(ifftshift(Klambda))));
        % Morally this is Klambda = 1./(4*pi^2*absXIsquared/L^2 + lambda^2)
        % except we get rid of the undersampling problem near the origin.

y = y.*Klambda;
clear Klambda;
y = fftshift(ifftn(ifftshift(y)));
y = (sqrt(2*beta+1)*soliton.^(beta)).*real(y);
y = reshape(y,N^3,1,1);
