📄 emgsafwd.m
字号:
function [result] = emgsafwd(varargin);
global debug
% [result] = emgsafwd(config, layers, cparams);
%
% emgsafwd.m is part of the CR1Dmod forward modeling package, and contains
% the code used to calculate the response of electrode arrays on the surface
% of a 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)
% name: optional name to display in the waitbar window
%
% Output:
% result: structure containing the calculated responses as well
% as information on the configuration, such as the
% geometric factor.
%
% Written by:
% Thomas Ingeman-Nielsen
% The Arctic Technology Center, BYG
% Technical University of Denmark
% Email: tin@byg.dtu.dk
cpu_t=cputime;
global Rkern estimate
Rkern = [];
estimate = [];
if isappdata(0,'debug') && getappdata(0,'debug')
dbstop if all error
debug = 1;
else
debug = 0;
end
if ( nargin < 3 )
disp('A call to emgsafwd must provide three input arguments!');
return
end
if ( nargin == 4 ) && ~isempty(varargin{4})
modelname = [varargin{4} ': '];
else
modelname = [];
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
end
end
if ~ok
disp('Non-structure inputs have not been implemented yet!');
return
end
% Set up vectors of layer parameters
eps0 = 8.85418782.*1e-12;
mu0 = 4e-7.*pi;
h = [layers.thickness]'; h(end) = []; % only include layers above halfspace
omega = cparams.freq.*2.*pi;
sigma = 1./Z_CR(omega, [layers.rho].', [layers.tau].', [layers.c].', ...
[layers.m].');
if layers(1).mu == 0
cparams.nonmag = 1;
else
cparams.nonmag = 0;
end
h1 = waitbar(0,[modelname 'Please wait...']);
% setup config structure for dipole-dipole calculations
if strcmpi(config.type, 'Dipole-Dipole')
for k = 1:length(config(1,1).Aspac)
for m = 1:length(config(1,1).Nspac)
config(k,m) = config(1,1);
config(k,m).C1 = [-config(k,m).Aspac(k)-...
0.5*config(k,m).Aspac(k)*config(k,m).Nspac(m) 0 0]; % transmitter dipole C1 a C2 P2 a P1
config(k,m).C2 = [config(k,m).C1(1)+...
config(k,m).Aspac(1) 0 0]; % transmitter dipole o------o o------o
config(k,m).P1 = [config(k,m).C2(1)+... % receiver dipole Tx n*a Rx
config(k,m).Nspac(m)*config(k,m).Aspac(k) 0 0];
config(k,m).P2 = [config(k,m).P1(1)+... % receiver dipole
config(k,m).Aspac(k) 0 0];
config(k,m).Cwire = [config(k,m).C1;config(k,m).C2];
% Calculate length of each wire element
config(k,m).S = sqrt(sum(diff(config(k,m).Cwire).^2,2));
config(k,m).Pwire = [config(k,m).P1;config(k,m).P2];
% Calculate length of each wire element
config(k,m).T = sqrt(sum(diff(config(k,m).Pwire).^2,2));
config(k,m).cosines = 1;
config(k,m).R = [];
for s = 1:length(config(k,m).S)
for t = 1:length(config(k,m).T)
config(k,m).R = [config(k,m).R; sqrt(sum( ...
(config(k,m).Cwire(s:s+1,:)- ...
config(k,m).Pwire(t:t+1,:)).^2,2))];
config(k,m).R = [config(k,m).R; sqrt(sum( ...
(config(k,m).Cwire(s:s+1,:)- ...
config(k,m).Pwire([t+1,t],:)).^2,2))];
end
end
config(k,m).Rmin = min(config(k,m).R);
config(k,m).Rmax = max(config(k,m).R);
% calculate distance between electrodes
config(k,m).r = sqrt(sum(([...
config(k,m).P2-config(k,m).C2;...
config(k,m).P2-config(k,m).C1;...
config(k,m).P1-config(k,m).C2;...
config(k,m).P1-config(k,m).C1]).^2,2));
% prepare result structure
result(k,m).Aspac = config(k,m).Aspac(k);
result(k,m).Nspac = config(k,m).Nspac(m);
% result(k,m).G_factor = 2*pi*(...
% 1/config(k,m).r(1)- ...
% 1/config(k,m).r(2)- ...
% 1/config(k,m).r(3)+ ...
% 1/config(k,m).r(4))^(-1);
result(k,m).G_factor = 2*pi*(...
1/config(k,m).r(4)- ...
1/config(k,m).r(3)- ...
1/config(k,m).r(2)+ ...
1/config(k,m).r(1))^(-1);
result(k,m).Z = omega'; %allocate space for Z
end % m
end % k
cparams.Rmin = min(min([config.Rmin]));
cparams.Rmax = max(max([config.Rmax]));
else % if it is a general surface array
config.Cwire = [config.C1;config.Cwire;config.C2];
config.S = sqrt(sum(diff(config.Cwire).^2,2)); % Calculate length of each wire element
config.Pwire = [config.P1;config.Pwire;config.P2];
config.T = sqrt(sum(diff(config.Pwire).^2,2)); % Calculate length of each wire element
config.cosines = (diff(config.Cwire)*diff(config.Pwire).') ...
./(abs(config.S)*abs(config.T).'); % calculate cosines of angles between wire elements
config.R = [];
for s = 1:length(config.S)
for t = 1:length(config.T)
if config.cosines(s,t) ~=0
config.R = [config.R; sqrt(sum((config.Cwire(s:s+1,:)...
-config.Pwire(t:t+1,:)).^2,2))];
config.R = [config.R; sqrt(sum((config.Cwire(s:s+1,:)...
-config.Pwire([t+1,t],:)).^2,2))];
end % if
end % t
end % s
cparams.Rmin = min(config.R);
cparams.Rmax = max(config.R);
% calculate distance between electrodes
config.r = sqrt(sum(([config.P2-config.C2; config.P2-config.C1; ...
config.P1-config.C2; config.P1-config.C1]).^2,2));
% prepare result structure
result.G_factor = 2*pi*(1/config.r(1)-1/config.r(2)-1/config.r(3)+ ...
1/config.r(4))^(-1)
result.Z = omega'; %allocate space for Z
end % if
if ( ~cparams.Rspline )
% In the case no spline interpolation is wanted
% return handle to FD calculation routine
estimate = [];
switch cparams.domain
case {'FD'}
P_handle = @P_fun;
Q_handle = @Q_fun;
case {'TD'}
disp('Not implemented yet!');
end % switch
elseif ~isempty(cparams.Rmin) % If there is coupling
% If spline interpolation is wanted
% setup some parameters and return handle to spline evaluation
Rlim = (floor(log2([cparams.Rmin cparams.Rmax]).*cparams.Rsp_NDEC+100)-[102 97])./cparams.Rsp_NDEC;
cparams.R = 2.^(Rlim(1):1./cparams.Rsp_NDEC:Rlim(2));
if isappdata(0, 'debug') && getappdata(0, 'debug')
disp('Evaluation distances for the creation of the spline function:');
disp(cparams.R);
assignin('base','R_values',cparams.R);
end
P_handle = @eval_Rspline;
Q_handle = @Q_fun;
else
Q_handle = @Q_fun;
end
lparams.h = h;
original_Quad_tol = cparams.Quad_tol;
lo = length(omega); % number of frequencies
lc = prod(size(config)); % number of configurations
lt = lo.*lc; % total number of calculations
for k = 1:lo % cycle over frequencies
lparams.z0 = i.*omega(k).*mu0;
lparams.y0 = i.*omega(k).*eps0;
lparams.k0_sq = -lparams.z0.*lparams.y0;
lparams.y = sigma(:,k) + i.*eps0.*[layers.eps_r].'*omega(k);
lparams.z = i.*omega(k).*mu0.*([layers.mu]'+1);
if strcmp(cparams.calc_type,'Quasi')
% Quasi static Halfspace response assumes also nonmagnetic first layer
lparams.z(1,:) = lparams.z0(1,:);
% ensures non-magnetic first layer
end
lparams.k_sq = -lparams.z.*lparams.y;
if ( cparams.Rspline && isfield(cparams, 'R'))
gsa_driver(cparams, omega(k), lparams);
end
for m = 1:lc % cycle over configurations
r_un = unique(config(m).r);
drawnow;
P_term = 0;
if ~all(isequal(config(m).cosines, 0))
% tolerance for each segment integration should add up to Quad_tol
cparams.Quad_tol = original_Quad_tol/ ...
sum(sum(config(m).cosines&1));
total_num = length(config(m).S).*length(config(m).T);
for s = 1:length(config(m).S)
for t = 1:length(config(m).T)
if config(m).cosines(s,t) ~=0
if debug
disp(['Integration pass: ' ...
num2str((s-1)*length(config(m).T)+t) ...
' of ' num2str(total_num)]);
end
P_term = P_term + ...
j.*omega(k).*4.*pi.*1e-7./4./pi.* ... % removed from the P-function and added here!
dblquad(P_handle, 0, config(m).T(t), 0, ...
config(m).S(s), cparams.Quad_tol, @quad, ...
config(m).Pwire(t:t+1,:), ...
config(m).Cwire(s:s+1,:), ...
config(m).T(t), config(m).S(s), ...
omega(k), lparams, cparams).* ...
config(m).cosines(s,t);
end
end
end
end
if debug
disp('Evaluating Q-function');
end
Q = feval(Q_handle, r_un, omega(k), lparams, cparams);
result(m).Z(k) = (P_term + ...
Q(config(m).r(1)==r_un)-Q(config(m).r(2)==r_un)- ...
Q(config(m).r(3)==r_un)+Q(config(m).r(4)==r_un));
titlestr = [modelname 'Just did: \omega = ' num2str(round(omega(k) ...
*1000)/1000) ' Config: ' num2str(m)]; % 'Hz, Aspac = ' ...
% num2str(Aspac(m)) 'm, Nspac = ' sprintf('%4.2f',Nspac(n))];
waitbar(((k-1)*lc+m)/lt, h1, titlestr);
end % m (config)
end % k (omega)
close(h1);
disp(['Time spent: ' num2str(cputime-cpu_t) ' s']);
%****************************************************
% Subfunctions
%****************************************************
% -------------------------------------------------------------------------
function gsa_driver(cparams, omega, lparams);
global Rkern estimate debug
if debug
disp(['GSA driver: ' num2str(omega)])
end
switch cparams.domain
case {'FD'}
if strcmp(cparams.calc_type,'Quasi')
% Quasi static Halfspace response assumes also nonmagnetic first layer
%disp('quasi+nonmag')
r = cparams.R;
kern = 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 {'NHT'}
for k = 1:length(r)
kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
NJCST('J0', @P_kernel_quasinonmag, ...
cparams.R(k), cparams.Seg_tol, ...
cparams.NHT_tol, cparams.Max_seg, ...
omega, lparams);
end
case {'FHT'}
for k = 1:length(r)
kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
FJCST('J04', cparams.R(k), -1, 0, 1, ...
cparams.FHT_err, @P_kernel_quasinonmag, ...
omega, lparams);
end
end
end
elseif cparams.nonmag
% nonmagnetic layered halfspace
%disp('nonmag');
r = cparams.R;
kern = 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 {'NHT'}
for k = 1:length(r)
kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
NJCST('J0', @P_kernel_nonmag, ...
cparams.R(k), cparams.Seg_tol, ...
cparams.NHT_tol, cparams.Max_seg, ...
omega, lparams);
end
case {'FHT'}
for k = 1:length(r)
kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
FJCST('J04', cparams.R(k), -1, 0, 1, ...
cparams.FHT_err, @P_kernel_nonmag, ...
omega, lparams);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -