temfwd.m

matlab中常用电法各种装置的电磁正演计算

%
% b:         radius of Transmitter loop
% config :    structure defining the configuation (generated by cr1dmod) 
% layers :    structure defining layer-parameters (generated by cr1dmod)
% cparams:    structure defining calculation options (generated by compute)

omega = 2.*b./(1./layers(1).rho.*4e-7.*pi.*config.TxR.^2); 
if ~all([layers.m]==0)
    sigma =  1./Z_CR(omega, [layers.rho].', [layers.tau].',             ...
        [layers.c].', [layers.m].');
else
    sigma = 1./[layers.rho].'*ones(size(omega));
end

if strcmp(cparams.calc_type,'Quasi')
    % Quasi static Halfspace response assumes also nonmagnetic first layer
    % Here calculated using analytic expression
    disp('quasi+nonmag')
    layers(1).mu = 0; % ensures non-magnetic first layer
    TEM_handle = @TEM_kernel_quasi;
    z0 = i.*omega.*4.*pi.*1e-7;
%    k0_sq = 0;
    y = sigma + i.*8.85418782.*1e-12.*[layers.eps_r].'*omega;
    k_sq = -repmat(z0,length(layers),1).*y;
    Hz = 1./(k_sq(1,:).*config.TxR.^4).*...
        (-3 -(k_sq(1,:).*config.TxR.^2-3.*i.*k_sq(1,:).^0.5.*       ...
        config.TxR-3).*exp(-i.*k_sq(1,:).^0.5.*config.TxR));
    
elseif layers(1).mu == 0
    % Halfspace response non-magnetic first layer
    disp('nonmag')
    TEM_handle = @TEM_kernel_nonmag;
    z0 = i.*omega.*4.*pi.*1e-7;
    y0 = i.*omega.*8.85418782.*1e-12;
    k0_sq = -z0.*y0;
    y = sigma + i.*8.85418782.*1e-12.*[layers.eps_r].'*omega;
    k_sq = -repmat(z0,length(layers),1).*y;
    Hz = 1./((k_sq(1,:)-k0_sq).*config.TxR.^4).*(                   ...
        (k0_sq.*config.TxR.^2-3.*i.*k0_sq.^0.5.*config.TxR-3).*     ...
        exp(-i.*k0_sq.^0.5.*config.TxR) - (k_sq(1,:).*              ...
        config.TxR.^2-3.*i.*k_sq(1,:).^0.5.*config.TxR-3).*         ...
        exp(-i.*k_sq(1,:).^0.5.*config.TxR) );
    
else
    disp('full')    
    TEM_handle = @TEM_kernel_full;
    Hz = zeros(size(omega));
end

if length(layers)>1 || (layers(1).mu ~= 0 &&                            ...
        ~strcmp(cparams.calc_type,'Quasi'))
    switch cparams.hank_type
        case 'NHT'
            for k=1:max(size(omega))
                Hz(k) = Hz(k) + ...          % assume Tx is on the ground...
                    NJCST('J1', TEM_handle, config.TxR, cparams.Seg_tol, ...
                    cparams.NHT_tol, cparams.Max_seg, config.Txz,       ...
                    omega(k), [layers.eps_r].'.*8.85418782.*1e-12,      ...
                    ([layers.mu].'+1).*4e-7.*pi, sigma(:,k),            ...
                    [layers(1:end-1).thickness].')+1;
            end
        case 'FHT'
            for k=1:max(size(omega))
                Hz(k) = Hz(k) + ...          % assume Tx is on the ground...
                    FJCST('J14', config.TxR, -1, 0, 0,                ...
                    cparams.FHT_err, TEM_handle, config.Txz, omega(k),  ...
                    [layers.eps_r].'.*8.85418782.*1e-12,                ...
                    ([layers.mu].'+1).*4e-7.*pi, sigma(:,k),            ...
                    [layers(1:end-1).thickness].')+1;
            end
    end
end
Hz = Hz.*(2.*config.TxR.^2);


% *************************************************************************
% * Code for square transmitter loop and dipole reciever                  *
% *************************************************************************
function [Hz] = Hz_SqDip(b, config, layers, cparams);

% [Hz] = Hz_SqDip(b, config, layers, cparams)
% Calculates the Hz field for the central loop configuration
% at a speciffic frequency.
%
% b:          Equivalent radius of Transmitter loop b = sqrt(TxS^2/pi)
%                     where TxS is the side length of the square loop.
% config :    structure defining the configuation (generated by cr1dmod) 
% layers :    structure defining layer-parameters (generated by cr1dmod)
% cparams:    structure defining calculation options (generated by compute)

omega = 2.*b./(1./layers(1).rho.*4e-7.*pi.*config.TxR.^2);

if ~all([layers.m]==0)
    sigma =  1./Z_CR(omega, [layers.rho].', [layers.tau].',             ...
        [layers.c].', [layers.m].');
else
    sigma = 1./[layers.rho].'*ones(size(omega));
end

TEM_handle = @square_quad;

Txloop = [-config.TxS/2  -config.TxS/2;...
        -config.TxS/2   config.TxS/2;...
        config.TxS/2   config.TxS/2;...
        config.TxS/2  -config.TxS/2;...
        -config.TxS/2  -config.TxS/2];

S = sqrt(sum(diff(Txloop).^2,2));            % Calculate length of each transmitter loop element (side)

%h1 = waitbar(0,'Please wait...');

Hz = zeros(1,length(omega));

% We assume here the central loop configuration in which
% all four sides of a square loop contributes equally.

for k = 1:length(omega)
    %    for m = 1:4              % do calculations for the four sides 
    m = 1;                % This is true only for square loop in central loop config
    Hz(k) = Hz(k) + quad(TEM_handle, 0, S(m), cparams.Quad_tol, [], ...
        Txloop(m:m+1,:), S(m), omega(k),                            ...
        8.85418782.*1e-12.*[layers.eps_r].',                        ...
        ([layers.mu].'+1).*4e-7.*pi, sigma,                         ...   % 1./[layers.rho].'
        [layers(1:end-1).thickness].', cparams.FHT_err);
    %    end
    %    titlestr = ['Just did: \omega = ' num2str(round(omega(k)*1000)/1000)];
    %    waitbar(k/length(omega), h1, titlestr);
end

Hz = 4.* Hz;
% close(h1);        

% ---------------------------------------------------------------------        
function term = square_quad(s, dTx, S, omega, eps, mu, sigma, h, FHT_err);

% This function calculates the contribution from a linear wire segment to 
% the vertical magnetic field at a certain location.
% So far only implemented with the Fast Hankel Transform

TEM_handle = @TEM_kernel_full;

ds = repmat(diff(dTx)./S,length(s),1)...
    .*repmat(s', 1, 2) + repmat(dTx(1,:),length(s),1);  

r = sqrt(sum(ds.^2,2));
if sum(dTx(:,2))==0  % This is a y-directed integration
    if dTx(1,2)>dTx(2,2)  % negative integration direction
        y = dTx(1,1);
    else
        y = -dTx(1,1);
    end
else                 % This is an x-directed integration
    if dTx(1,1)>dTx(2,1)  % negative integration direction
        y = dTx(1,2);
    else
        y = -dTx(1,2);
    end
end

z = 0; % we assume source and receiver on the ground

for k = 1:length(r)
    term(k) = y./r(k).*(...
        FJCST('J14', r(k), -1, 0, 0, FHT_err, TEM_handle,             ...
        z, omega, eps, mu, sigma, h) + 1/(2.*r(k)^2)); % + r./(2.*(r.^2+z.^2).^(3/2)));
end

term = term./(4.*pi);

% *************************************************************************
% * Kernel functions for circular/square transmitter loop and dipole Rz   *
% *************************************************************************

function kern = TEM_kernel_full(lambda, varargin)

% kern = TEM_kernel_full(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((R_TE+1)*lambda/u0)*lambda
%
% This version is for the full solution!

% varargin{1} contains the height of the observation point and is not used
omega = varargin{2};
eps =   varargin{3};
mu =    varargin{4};
sigma = varargin{5};
h =     varargin{6};

lambda_sq = lambda.^2;
z0 = i.*omega.*4.*pi.*1e-7;
y0 = i.*omega.*8.85418782.*1e-12;
k0_sq = -z0.*y0;
u0 = sqrt(lambda_sq-k0_sq);
z = j.*omega.*mu;
y = sigma + i.*eps.*omega;
k_sq = -z.*y;
u = sqrt(repmat(lambda_sq,size(k_sq,1),1)-repmat(k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(z,1,size(lambda,2));
Y0 = u0./z0;

gamma = 0;
if ~isempty(h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(h,1,length(lambda)));
    phi = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)    
    for m = length(h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end
phi1 = (Y0-Y(1,:))./(Y0+Y(1,:));            %  as needed in the numerator
phi1mod = 1./(u0 + u(1,:).*z0./z(1,:));     %  as needed in the denominator

kern = (((gamma+1).*phi1mod)./(gamma.*phi1+1).*lambda-1./2).*lambda;   % = (R_TE+1)*lambda^2/u0


%--------------------------------------------------------------------------
function kern = TEM_kernel_nonmag(lambda, varargin)

% kern = TEM_kernel_nonmag(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((R_TE+1)*lambda/u0)*lambda
%
% This version is for the nonmagnetic first layer solution!

% varargin{1} contains the height of the observation point and is not used
omega = varargin{2};
eps =   varargin{3};
mu =    varargin{4};
sigma = varargin{5};
h =     varargin{6};

lambda_sq = lambda.^2;
z0 = i.*omega.*4e-7.*pi;
y0 = i.*omega.*8.85418782.*1e-12;
k0_sq = -z0.*y0;
u0 = sqrt(lambda_sq-k0_sq);
z = j.*omega.*mu;
y = sigma + i.*eps.*omega; 
k_sq = -z.*y;
u = sqrt(repmat(lambda_sq,size(k_sq,1),1)-repmat(k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(z,1,length(lambda));
Y0 = u0./z0;

gamma = 0;
if ~isempty(h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(h,1,length(lambda)));
    phi = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)    
    for m = length(h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

kern = ((2.*gamma.*u(1,:)./(gamma.*(k_sq(1,:)-k0_sq)+(u0+u(1,:)).^2)).* ...
    lambda-1./2).*lambda;   


%--------------------------------------------------------------------------
function kern = TEM_kernel_quasi(lambda, varargin)

% kern = TEM_kernel_quasi(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((R_TE+1)*lambda/u0)*lambda
%
% This version is for the quasi-static solution! k0=0

% varargin{1} contains the height of the observation point and is not used
omega = varargin{2};
eps =   varargin{3};
mu =    varargin{4};
sigma = varargin{5};
h =     varargin{6};

lambda_sq = lambda.^2;
z0 = i.*omega.*4.*pi.*1e-7;
z = j.*omega.*mu;
y = sigma + i.*eps.*omega;
k_sq = -z.*y;
u = sqrt(repmat(lambda_sq,size(k_sq,1),1)-repmat(k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(z,1,size(lambda,2));
Y0 = lambda./z0;

gamma = 0;
if ~isempty(h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(h,1,length(lambda)));
    phi = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)    
    for m = length(h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

phi1 = (Y0-Y(1,:))./(Y0+Y(1,:));    %  as needed in the numerator
phi1mod = (1./z0)./(Y0+Y(1,:));     %  as needed in the denominator

kern = ((2.*gamma.*u(1,:)./(gamma.*k_sq(1,:)+(lambda+u(1,:)).^2)).* ...
    lambda-1./2).*lambda;