referenceData3D

Computes scattering data in 3D via series solution for a reference contrast q, that is a ball.

Contents

Syntax

[uScattRXRef, uScattROIRef] = referenceData3D(seti)

Description

[uScattRXRef, uScattROIRef] = referenceData3D(seti) computes the same in 3D as referenceData2D.html in 2D.

See referenceData2D.html for details. The differences in 3D are:

Example

init;
seti.dim = 3;
seti.incType = 'planeWave';
seti.incNb = 1;
seti.rBall = 0.015;
seti.qBall = 0.8;
seti.measType = 'farField'; % Note that 'nearField' is not available
seti.nCD = 64; % a high number will take a long computational time
seti = setGeomSim(seti); % set geometry and simulation
[uScattRXRef, uScattROIRef] = referenceData3D(seti);

See Also

Code

function [uScattRXRef, uScattROIRef] = referenceData3D(seti)

X = seti.grid(1,:);
Y = seti.grid(2,:);
Z = seti.grid(3,:);

[uScattRXRef,uScattRef] = measurementReference3D(seti.qBall, seti.rBall, seti.k, X, Y, Z, seti.measPnts, seti.measType);
uScattROIRef = uScattRef(seti.ROImask);

uScattROIRef = uScattROIRef(:);
uScattRXRef = uScattRXRef(:);
end

Code: subfunction measurementReference3D

Computation of reference solution in 3D and evaluation on mesh X, Y, Z.

function [uRXScatt, uScatt] = measurementReference3D(q, R, k, X, Y, Z, RX, type)
% k: wave number
% n = q^2 inside B(0,R)
%
% MODEL: \Delta u + k^2 (1+q) u = 0
%

uRXScatt = zeros(size(RX(1,:)));

% normalize RX
[~, NX] = size(RX);
RXdirs = zeros(size(RX)).';

for j=1:NX
    RXdirs(j,:) = RX(:,j) / norm(RX(:,j));
end

% Direction of incoming plane wave
dir = [1 0 0];

nMax=10;

norms = sqrt(X(:).^2+Y(:).^2+Z(:).^2);
directions = [X(:)./norms Y(:)./norms Z(:)./norms];
inside = norms < R;
outside = norms >= R;

zeroIdx = find(norms==0);
uScatt = zeros(size(norms));

for n=0:nMax
    A = [-spherHankel(n,k*R)  spherBessel(n,k*sqrt(1+q)*R)
        -dSpherHankel(n,k*R)  sqrt(1+q)*dSpherBessel(n,k*sqrt(1+q)*R)];

    detA = A(1,1)*A(2,2)-A(1,2)*A(2,1);
    AInv = 1/detA * [A(2,2), -A(1,2); -A(2,1), A(1,1)];

    rhs = [spherBessel(n,k.*R); dSpherBessel(n,k.*R)];
    t = AInv * rhs;

    coeff = zeros(2*n+1,2);
    for m=-n:n
        coeff(m+n+1,:) = (4*pi*1i^n)*conj(spherHarmonic(n, m, dir))*t;
    end

    % evaluation outside
    if ~all(outside==0)
        spherHar = zeros(length(norms(outside)), 1);

        for m=-n:n
            spherHar = spherHar + ...
                spherHarmonic(n,m,directions(outside,:))*coeff(m+n+1,1);
        end

        uScatt(outside) = uScatt(outside) + spherHankel(n,k*norms(outside)).*spherHar;
    end

    % evaluation inside (subtract incident wave)
    if ~all(inside==0)
        spherHarUi = zeros(length(norms(inside)), 1);
        spherHarInside = zeros(length(norms(inside)), 1);

        for m=-n:n
            tmp = spherHarmonic(n,m,directions(inside,:));

            spherHarInside = spherHarInside + tmp*coeff(m+n+1,2);
            spherHarUi = spherHarUi + tmp*4*pi*1.0i^n*conj(spherHarmonic(n,m,dir));
        end

        uScatt(inside) = uScatt(inside) + ...
            spherBessel(n,k*sqrt(1+q)*norms(inside)).*spherHarInside ...
            - spherBessel(n,k*norms(inside)).*spherHarUi;
    end

    if strcmp(type,'nearField')
        uRXScatt = 0;
        disp('uRXScatt was set to 0 in case of nearField in referenceData3D.m')
        % no near field implemented...
    elseif strcmp(type,'farField')
        for m=-n:n
            uRXScatt = uRXScatt + 1/k * 1/(1.0i^(n+1)) * spherHarmonic(n,m,RXdirs)' * coeff(m+n+1,1);
        end
    end
end

% not computable at 0, use average
uScatt(zeroIdx) = 0.5*(uScatt(zeroIdx+1)+uScatt(zeroIdx-1));
uScatt = reshape(uScatt, size(X));
end

Code: subfunction: spherHarmonic: Spherical harmonics

function y = spherHarmonic(n,m,x)
[phi,theta,~] = cart2sph(x(:,1),x(:,2),x(:,3));

leg = legendre(n,cos(theta+pi/2).');
legM = leg(abs(m)+1,:);

y = sqrt( ((2*n+1)/(4*pi))*(factorial(n-abs(m))/factorial(n+abs(m))) ) ...
    .* legM.' .* exp(1i*m*phi);
end

Code: subfunction: spherHankel: Spherical Hankel functions out of normal ones

function u = spherHankel(n,z)
u = sqrt(pi./(2.*z)).*besselh(n+1/2,z);
end

Code: subfunction: spherBessel: Spherical Bessel functions out of normal ones

function u = spherBessel(n,z)
u = sqrt(pi./(2.*z)).*besselj(n+1/2,z);
end

Code: subfunction: dSpherHankel: Derivative of spherical Hankel functions

function u = dSpherHankel(n,z)
u = -spherHankel(n+1,z) + (n./z).*spherHankel(n,z);
end

Code: subfunction: dSpherBessel: Derivative of spherical Bessel functions

function u = dSpherBessel(n,z)
u = -spherBessel(n+1,z) + (n./z).*spherBessel(n,z);
end