# setKernel

Computes the Fourier coefficients of the kernel.

## Contents

- Syntax
- Description
- Example
- Input Arguments
- Output Arguments
- More About
- References
- See Also
- Code
- Code: subfunction: kernelConsis
- Code: subfunction: kernel2D
- Code: subfunction: kernel3D

## Syntax

    seti = setKernel(seti);
    seti = setKernel(seti,dispDepth);

## Description

`seti = setKernel(seti)` computes the Fourier coefficients of the kernel and stores them in `seti.kHat`. Note that the factor is not included in Fourier coefficients.

`seti = setKernel(seti,dispDepth)` does the same, but allows to control the depth of displayed messages.

## Example

    init;

    seti.dim = 2;
    seti.nCD = 256;
    seti.rCD = 0.2;
    seti.k = 250;
    seti.model = 'helmholtz2D';

    seti = setKernel(seti);

The computed Fourier coefficients of the kernel are in `seti.kHat`.

## Input Arguments

This arguments was set in rebis before setKernel is called:

- seti.dim : see setGrid.html.
- seti.nCD : see setGrid.html.
- seti.rCD : see setGrid.html.

Can be set by user (otherwise default is set):

- seti.k : Wave number, default: 250
- seti.model : Model of the problem, default: 'helmholtz2D'

Optional Input Argument

- dispDepth : Depth of displayed messages (0: no, 1 or greater: yes).

## Output Arguments

- seti.kHat : kernel, complex matrix of size nCD x nCD.

## More About

For theory see [1], [2] or Section 3.2 in [3].

Note that the factor is not included in Fourier coefficients in this code, but the factor is included in [3].

## References

- [1] Gennadi Vainikko. Fast Solvers of the Lippmann-Schwinger Equation. In Robert P. Gilbert, Joji Kajiwara, and Yongzhi S. Xu, editors, Direct and Inverse Problems of Mathematical Physics, pages 423–440. Springer, Boston, 2000. URL: https://doi.org/10.1007/978-1-4757-3214-6_25.
- [2] Thorsten Hohage. On the numerical solution of a three-dimensional inverse medium scattering problem. Inverse Problems, 17(6):1743, 2001. URL: https://doi.org/10.1088/0266-5611/17/6/314.
- [3] Florian Bürgel, Kamil S. Kazimierski, and Armin Lechleiter. A sparsity regularization and total variation based computational framework for the inverse medium problem in scattering. Journal of Computational Physics, 339:1-30, 2017. URL: https://doi.org/10.1016/j.jcp.2017.03.011.

## See Also

- setGeomSim.html

## Code

    function seti = setKernel(seti,varargin)

    if nargin == 2
        dispDepth = varargin{1};
    else
        dispDepth = 0;
    end

    seti = kernelConsis(seti,dispDepth);

    if strcmp(seti.model,'helmholtz2D') || strcmp(seti.model,'helmholtzHMode2D')
        seti.kHat = kernel2D(seti);
    elseif strcmp(seti.model,'helmholtz3D')
        seti.kHat = kernel3D(seti);
    else
        disp('setKernel.m: Error - Kernel coefficients of chosen model not implemented')
    end

    end

## Code: subfunction: kernelConsis

Checks consistency of kernel settings.

    function seti = kernelConsis(seti,dispDepth)

    seti = checkfield(seti,'k',250,dispDepth);

    if ~isfield(seti,'model')
        if seti.dim == 2
            seti.model = 'helmholtz2D';
        else
            seti.model = 'helmholtz3D';
        end
        setmessage(seti,'model',dispDepth);
    end

    if strcmp(seti.model,'helmholtz')
        if seti.dim == 2
            seti.model = 'helmholtz2D';
        else
            seti.model = 'helmholtz3D';
        end
        if dispDepth >= 1
            fprintf('   Adapted Helmholtz model to dimension %d.\n', seti.dim);
        end
    end

    if strcmp(seti.model,'helmholtz2D') && (seti.dim == 3)
        seti.dim = 2;
        if dispDepth >= 1
            disp('   Parameter "dim"=3 does not fit model "helmholtz2D". Set dim = 2.');
        end
    end

    if strcmp(seti.model,'helmholtz3D') && (seti.dim == 2)
        seti.dim = 3;
        if dispDepth >= 1
            disp('   Parameter "dim"=2 does not fit model "helmholtz3D". Set dim = 3.');
        end
    end

    end

## Code: subfunction: kernel2D

Computes the Fourier coefficients of the kernel in 2D case.

    function y = kernel2D(seti)

    s = ceil((seti.nCD-1)/2);
    j1 = (0:seti.nCD-1)-s;
    [j1,j2] = meshgrid(j1,j1);
    p = pi*sqrt(j1.^2+j2.^2); % p = pi*|j|

In this file we use notatio R = seti.rCD (this is the notation in Vainikko, see [1]). In [3] we used the notation 2 R = seti.rCD.

    kR = seti.k*seti.rCD;

    J  = @besselj; %besselj: Bessel function of first kind
    % besselh: Bessel function of third kind = Hankel function
    H0 = 1i*pi/2*besselh(0,1,kR); % using Hankel function of first kind and order 1
    H1 = 1i*pi/2*besselh(1,1,kR); % using Hankel function of first kind and order 0

The following two lines follow Vainikko (see [1]): * Note that we also deal with k != 1, while Vainikko presents only the case k = 1 on page 16. * Note the missing factor 1/(2R) in Fourier coefficients y = ... because Vainikko uses another basis function (with factor 1/(2R)).

    y = seti.rCD^2./(p.^2-kR^2).*(1 + H0*p.*J(1,p ) - H1*kR*J(0,p )); % for p != kR
    y(p == kR) = seti.rCD^2/2*(   H0*J(0,kR)    + H1*J(1,kR)); % can be proven by l'Hospital from the above formula

    s = s*[1 1]; % Do not delete this line!
    y = circshift(y,-s);

    end

## Code: subfunction: kernel3D

Computes the Fourier coefficients of the kernel in 3D case.

    function y = kernel3D(seti)

    kR = seti.k*seti.rCD;

    N2 = ceil((seti.nCD-1)/2);
    [JX,JY, JZ] = meshgrid((0:seti.nCD-1)-N2);

    % pi*|j|
    PIJ = pi*sqrt(JX.^2 + JY.^2 + JZ.^2);

    % The following computations of K follow [1] and [2],
    % but note the correction in case of PIJ == kR (as in Section 3.2 of [3]).

    K = kR^2./(PIJ.^2-kR^2) .* (1 - exp(kR*1.0i) * ( cos(PIJ) - 1.0i*kR./PIJ .* sin(PIJ) ) );
    %K(PIJ == kR) = -1.0i/sqrt(2) * kR *(1 - exp(kR*1.0i) / kR * sin(kR)); % wrong
    K(PIJ == kR) = 1/2*exp(1i*kR)*(kR*sin(kR) + 1i*(kR*cos(kR) - sin(kR))); % corrected (!) (see Section 3.2 in [3])
    K(PIJ == 0) = exp(kR*1.0i) * (1-kR*1.0i) - 1;

    % Exclude factor k^2
    K =  K/seti.k^2;

    % Note that circshift(K,-N2*[1,1,1]) is the same as fftshift(K).
    y = circshift(K,-N2*[1,1,1]);

    end