emgsafwd.m

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

                        end
                end
            end                    
        else
            % full solution
            %disp('full solution');
            switch cparams.hank_type
                case {'NHT'}
                    for k = 1:length(cparams.R)
                        kern(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                            NJCST('J0', @P_kernel_full, cparams.R(k),   ...
                            cparams.Seg_tol, cparams.NHT_tol,           ...
                            cparams.Max_seg, omega, lparams);
                    end
                case {'FHT'}
                    for k = 1:length(cparams.R)
                        kern(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                            FJCST('J04', cparams.R(k), -1, 0, 1,        ...
                            cparams.FHT_err, @P_kernel_full, omega,     ...
                            lparams);
                    end
            end
        end
    case {'TD'}
        disp('Not implemented yet!');
end % switch

strct = spline(cparams.R,kern);

strct.R = cparams.R;
strct.kern = kern;
strct.estimate = [];

if ~isempty(Rkern)
%    Rkern(end+1) = strct;    % This line is used to save the splines for
%    later review
    Rkern = strct;
else
    Rkern = strct;    
end
if debug
    disp('End of GSA driver');
end
% -------------------------------------------------------------------------
function kern = eval_Rspline(y, x, Pwire, Cwire, T, S, omega,           ...
    lparams, cparams);

global Rkern

P = repmat(diff(Pwire)./T,length(y),1)...
    .*repmat(y', 1, 3) + repmat(Pwire(1,:),length(y),1);  % Pwire-element coordinates (will be a vector)
C = diff(Cwire)./S.*x + Cwire(1,:);                       % Cwire-element coordinates (will be a scalar)

r = sqrt(sum((P-repmat(C,length(y),1)).^2,2));            % distances are calculated

kern = ppval(r, Rkern(end));
%Rkern(end).estimate = [Rkern(end).estimate; r]; % use this line to save
%r-values of estimation points.

% -------------------------------------------------------------------------
function P = P_fun(y, x, Pwire, Cwire, T, S, omega, lparams, cparams);

global estimate

P = repmat(diff(Pwire)./T,length(y),1)...
    .*repmat(y', 1, 3) + repmat(Pwire(1,:),length(y),1);  % Pwire-element coordinates (will be a vector)
C = diff(Cwire)./S.*x + Cwire(1,:);                       % Cwire-element coordinates (will be a scalar)

r = sqrt(sum((P-repmat(C,length(y),1)).^2,2));            % distances are calculated

% cparams.calc_type = 'Quasi';

if strcmp(cparams.calc_type,'Quasi')
    % Quasi static Halfspace response assumes also nonmagnetic first layer
    %disp('quasi+nonmag')
    P = 2./(lparams.k_sq(1).*r.^3).*((i.*lparams.k_sq(1).^0.5.*r+1).*   ...
        exp(-i.*lparams.k_sq(1).^0.5.*r) - 1);
    if ~isempty(lparams.h)
        switch cparams.hank_type
            case {'FHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err, ...
                        @P_kernel_quasinonmag, omega, lparams);
                end
            case {'NHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        NJCST('J0', @P_kernel_quasinonmag, r(k),         ...
                        cparams.Seg_tol, cparams.NHT_tol,               ...
                        cparams.Max_seg, omega, lparams);
                end
        end
    end
elseif cparams.nonmag
    % nonmagnetic layered halfspace
    %disp('nonmag');
    P = 2./((lparams.k_sq(1)-lparams.k0_sq).*r.^3).*(                   ...
        (i.*lparams.k_sq(1).^0.5.*r+1).*exp(-i.*lparams.k_sq(1).^0.5.*r)...
        - (i.*lparams.k0_sq.^0.5.*r+1).*exp(-i.*lparams.k0_sq.^0.5.*r));
    if ~isempty(lparams.h)
        switch cparams.hank_type
            case {'FHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err, ...
                        @P_kernel_full, omega, lparams);
                end
            case {'NHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        NJCST('J0', @P_kernel_nonmag, r(k),              ...
                        cparams.Seg_tol, cparams.NHT_tol,               ...
                        cparams.Max_seg, omega, lparams);
                end
        end
    end
else
    % full solution
    %disp('fullsol');
    switch cparams.hank_type
        case {'FHT'}
            for k = 1:length(r)
                P(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                    FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err,     ...
                    @P_kernel_full, omega, lparams);
            end
        case {'NHT'}
            for k = 1:length(r)
                P(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                    NJCST('J0', @P_kernel_full, r(k), cparams.Seg_tol,   ...
                    cparams.NHT_tol, cparams.Max_seg, omega, lparams);
            end
    end
end

%estimate = [estimate; r, P];

% -------------------------------------------------------------------------
function Q = Q_fun(r, omega, lparams, cparams)

if strcmp(cparams.calc_type,'Quasi')
    % Quasi static Halfspace response assumes also nonmagnetic first layer
    %disp('quasi+nonmag')
    Q = -lparams.z0./(2.*pi.*lparams.k_sq(1,:).*r);
    if ~isempty(lparams.h)
        switch cparams.hank_type
            case {'FHT'}
                for k = 1:length(r)
                    Q(k) = -1./4./pi.*... 
                        (FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err,        ...
                        @Q_kernel_quasinonmag, omega, lparams)) + Q(k);
                end
            case {'NHT'}
                for k = 1:length(r)
                    Q(k) = -1./4./pi.*... 
                        (NJCST('J0', @Q_kernel_quasinonmag, r(k),                ...
                        cparams.Seg_tol, cparams.NHT_tol, cparams.Max_seg,      ...
                        omega, lparams)) + Q(k);
                end
        end
    end
else
    % full solution
    %disp('fullsol');
    switch cparams.hank_type
        case {'FHT'}
            for k = 1:length(r)
                Q(k) = 1./4./pi.*... 
                    (FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err,            ...
                    @Q_kernel_full, omega, lparams));
            end
        case {'NHT'}
            for k = 1:length(r)
                Q(k) = 1./4./pi.*... 
                    (NJCST('J0', @Q_kernel_full, r(k), cparams.Seg_tol,          ...
                    cparams.NHT_tol, cparams.Max_seg, omega, lparams));
            end
    end
end


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

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

omega =   varargin{1};
lparams = varargin{2};

lambda_sq = lambda.^2;
u0 = sqrt(lambda_sq-lparams.k0_sq);
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Y0 = u0./lparams.z0;

gamma = 0;
if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.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(lparams.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./lparams.z0)./(Y0+Y(1,:));     %  as needed in the denominator

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


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

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

omega =   varargin{1};
lparams = varargin{2};

lambda_sq = lambda.^2;
u0 = sqrt(lambda_sq-lparams.k0_sq);
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Y0 = u0./lparams.z0;

gamma = 0;
if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.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(lparams.h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

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

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

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

omega =   varargin{1};
lparams = varargin{2};

lambda_sq = lambda.^2;
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Y0 = lambda./lparams.z0;

gamma = 0;
if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.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(lparams.h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

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

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

% kern = Q_kernel(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((1-R_TM)*u0/y0 - (1+r_TE)*z0/u0)/lambda

omega =   varargin{1};  
lparams = varargin{2};

lambda_sq = lambda.^2;
u0 = sqrt(lambda_sq-lparams.k0_sq);
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Z = u./repmat(lparams.y,1,size(lambda,2));
Y0 = u0./lparams.z0;
Z0 = u0./lparams.y0;

gamma_TE = 0;
gamma_TM = 0;

if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi_TE = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    phi_TM = (Z(1:end-1,:)-Z(2:end,:))./(Z(1:end-1,:)+Z(2:end,:)); %  = (Zn-Zn+1)/(Zn+Zn+1)
    for m = length(lparams.h):-1:1
        gamma_TE = expuh2(m,:).*(gamma_TE+phi_TE(m,:))./(gamma_TE.*phi_TE(m,:)+1);
        gamma_TM = expuh2(m,:).*(gamma_TM+phi_TM(m,:))./(gamma_TM.*phi_TM(m,:)+1);    
    end
end

phi1_TE = (Y0-Y(1,:))./(Y0+Y(1,:)); % as needed in the numerator
phi1mod_TE = 2./(Y0+Y(1,:));        % as needed in the denominator

phi1_TM = (Z0-Z(1,:))./(Z0+Z(1,:)); % as needed in the numerator
phi1mod_TM = -2.*Z(1,:).*Z0./(Z0+Z(1,:));        % as needed in the denominator

kern_TE = ((gamma_TE+1).*phi1mod_TE)./(gamma_TE.*phi1_TE + 1); 
kern_TM = ((gamma_TM-1).*phi1mod_TM)./(gamma_TM.*phi1_TM + 1);

kern = (kern_TM-kern_TE)./lambda;   


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

% kern = Q_kernel(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((1-R_TM)*u0/y0 - (1+r_TE)*z0/u0)/lambda

omega =   varargin{1};  
lparams = varargin{2};

%lambda = lambda.*induct;
lambda_sq = lambda.^2;
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Z = u./repmat(lparams.y,1,size(lambda,2));
Y0 = lambda./lparams.z0;

gamma_TE = 0;
gamma_TM = 0;

if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi_TE = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    phi_TM = (Z(1:end-1,:)-Z(2:end,:))./(Z(1:end-1,:)+Z(2:end,:)); %  = (Zn-Zn+1)/(Zn+Zn+1)
    for m = length(lparams.h):-1:1
        gamma_TE = expuh2(m,:).*(gamma_TE+phi_TE(m,:))./(gamma_TE.*phi_TE(m,:)+1);
        gamma_TM = expuh2(m,:).*(gamma_TM+phi_TM(m,:))./(gamma_TM.*phi_TM(m,:)+1);    
    end
end

kern_TE = gamma_TE.*lparams.z0./(gamma_TE.*lparams.k_sq(1,:) +          ...
    (lambda + u(1,:)).^2);
kern_TM = gamma_TM./((gamma_TM+1).*lparams.y(1,:));

kern = (kern_TM+kern_TE).*4.*u(1,:)./lambda;