fdemfwd.m

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

function out = FDEMfwd(varargin);

% out = FDEMfwd(config, layers, cparams);
% 
% FDEMfwd.m is part of the CR1Dmod forward modeling package, and contains
% the code used to calculate the response of a HCP configuration over
% 1D layered half space.
%
% Input:
% config :    structure defining the configuation (generated by cr1dmod) 
% layers :    structure defining layer-parameters (generated by cr1dmod)
% cparams:    structure defining calculation options (generated by compute)
%
% Output:
% out:        structure containing the model result in terms of
%                H:       Total measured field
%                H_prim:  Primary field
%
% Written by:
% Thomas Ingeman-Nielsen
% The Arctic Technology Center, BYG
% Technical University of Denmark
% Email: [e-mail]

cpu_t=cputime;

if ( nargin < 3 )
    disp('A call to fdemfwd must provide three input arguments!');
    return
end

% Set up structures of parameters...
ok = 0;
if ( isstruct(varargin{1})     )
    config = varargin{1};
    if ( isstruct(varargin{2}) )
        layers = varargin{2};
        if ( isstruct(varargin{3}) )
            cparams = varargin{3};
            ok = 1;
        end % if
    end % if
end % if

if ~ok
    disp('Input must be structures!');
    return
end % if

switch config.type
    case {'HCP FDEM (HLEM)'}
        out = Zdip_Hz(config, layers, cparams);  % Hz field of Z directed magnetic dipole 
    case {'VCP'}
        % not implemented yet
    case {'VCA'}
        % not implemented yet
end % switch

disp(['Time spent: ' num2str(cputime-cpu_t) ' s']);


% *************************************************************************
% * Code for Hz from a z-directed magnetic dipole                         *
% *************************************************************************
function [result] = Zdip_Hz(config, layers, cparams)

% [Hz] = Zdip_Hz(config, layers, cparams);
% Calculates the Hz from a z-directed magnetic dipole

omega = cparams.freq.*2.*pi;

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

z0 = i.*omega.*4.*pi.*1e-7;
y0 = i.*omega.*8.85418782.*1e-12;
k0_sq = -z0.*y0;

for r = 1:length(config.Rspac)
    % Allocate memory for Hz
    result(r).H = zeros(size(omega));
    result(r).H_prim = zeros(size(omega));
    
    result(r).H_prim = 1./(4.*pi.*config.Rspac(r).^3).*                        ...
        exp(-i.*k0_sq.^0.5.*config.Rspac(r)).*(                                ...
        k0_sq.*config.Rspac(r).^2-i.*k0_sq.^0.5.*config.Rspac(r)-1);
    
    if strcmpi(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
        FDEM_handle = @Zdip_Hz_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;
        for k=1:max(size(omega))
            result(r).H(k) = 1/(2.*pi.*k_sq(1,k).*config.Rspac(r).^5).*           ...
                (9 -(9 + 9.*i.*k_sq(1,k).^0.5.*config.Rspac(r) -               ...
                4.*k_sq(1,k).*config.Rspac(r).^2-                              ...
                i.*k_sq(1,k).^1.5.*config.Rspac(r).^3)                         ...
                .*exp(-i.*k_sq(1,k).^0.5.*config.Rspac(r)));
        end % k
        result(r).H_prim = -ones(size(omega))./(4.*pi.*config.Rspac(r).^3);
        
    elseif layers(1).mu == 0
        % Halfspace response non-magnetic first layer
        % Here calculated using analytic expression    
        disp('nonmag')
        FDEM_handle = @Zdip_Hz_kernel_nonmag;
        y = sigma + i.*8.85418782.*1e-12.*[layers.eps_r].'*omega;
        k_sq = -repmat(z0,length(layers),1).*y;
        for k=1:max(size(omega))
            result(r).H(k) = 1/(2.*pi.*(k_sq(1,k)-k0_sq(k)).*config.Rspac(r).^5)  ...
                .*(                                                         ...
                (9 + 9.*i.*k0_sq(k).^0.5.*config.Rspac(r) -                    ...
                4.*k0_sq(k).*config.Rspac(r).^2-                               ...
                i.*k0_sq(k).^1.5.*config.Rspac(r).^3)                          ...
                .*exp(-i.*k0_sq(k).^0.5.*config.Rspac(r))                      ...
                -                                                           ...
                (9 + 9.*i.*k_sq(1,k).^0.5.*config.Rspac(r) -                   ...
                4.*k_sq(1,k).*config.Rspac(r).^2-                              ...
                i.*k_sq(1,k).^1.5.*config.Rspac(r).^3)                         ...
                .*exp(-i.*k_sq(1,k).^0.5.*config.Rspac(r))                     ...
                );
        end % k
    else
        disp('full')    
        FDEM_handle = @Zdip_Hz_kernel_full;
    end
    
    if length(layers) > 1 ||                                                ...
            ~strcmpi(cparams.calc_type,'Quasi')
        % If not homogeneous half-space, go through with hankel transform evaluation
        switch cparams.hank_type
            case 'NHT'
                for k=1:max(size(omega))
                    result(r).H(k) = result(r).H(k) + (1./(4.*pi)).*        ...
                        NJCST('J0', FDEM_handle, config.Rspac(r),           ...
                        cparams.Seg_tol, cparams.NHT_tol, cparams.Max_seg,  ...                             ...
                        omega(k), [layers.eps_r].'.*8.85418782.*1e-12,      ...
                        ([layers.mu].'+1).*4e-7.*pi, sigma(:,k),            ...
                        [layers(1:end-1).thickness].');
                end % k
            case 'FHT'
                for k=1:max(size(omega))
                    result(r).H(k) = result(r).H(k) + (1./(4.*pi)).*        ...
                        FJCST('J04', config.Rspac(r), -1, 0, 0,             ...
                        cparams.FHT_err, FDEM_handle,                       ...
                        omega(k), [layers.eps_r].'.*8.85418782.*1e-12,      ...
                        ([layers.mu].'+1).*4e-7.*pi, sigma(:,k),            ...
                        [layers(1:end-1).thickness].');
                end % k
        end % switch
        
    end % if
end % r

% *************************************************************************
% * Kernel functions for Hz from a z-directed magnetic dipole             *
% *************************************************************************

function kern = Zdip_Hz_kernel_full(lambda, varargin)

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

omega = varargin{1};
eps =   varargin{2};
mu =    varargin{3};
sigma = varargin{4};
h =     varargin{5};

lambda_sq = lambda.^2;
z0 = i.*omega.*4.*pi.*1e-7;
y0 = i.*omega.*8.854.*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 = (2./z0)./(Y0+Y(1,:));     %  as needed in the denominator

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


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

% kern = P_kernel(lambda, omega, eps, sigma, h)
%
% Kernelfunction:       (1+R_TE).*lambda^3./u0-homhalf
%
% This version is for the nonmagnetic ground case!

omega = varargin{1};
eps =   varargin{2};
mu =    varargin{3};
sigma = varargin{4};
h =     varargin{5};


lambda_sq = lambda.^2;
z0 = i.*omega.*4e-7.*pi;
y0 = i.*omega.*8.854.*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

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

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

% kern = P_kernel(lambda, omega, eps, sigma, h)
%
% Kernelfunction:       (1+R_TE).*lambda^3./u0-homhalf
%
% This version is for the quasistatic nonmagnetic ground case!

omega = varargin{1};
eps =   varargin{2};
mu =    varargin{3};
sigma = varargin{4};
h =     varargin{5};

lambda_sq = lambda.^2;
z0 = i.*omega.*4e-7.*pi;
% y0 = 0;
% k0_sq = -z0.*y0 = 0;
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

kern = (4.*gamma.*u(1,:).*lambda.^3)./                                  ...
    (gamma.*k_sq(1,:) + (lambda+u(1,:)).^2);                               % = layered-homhalf