📄 dcgsafwd.m
字号:
function [result] = dcgsafwd(varargin);
% [result] = dcgsafwd(config, layers, cparams);
%
% dcgsafwd.m is part of the CR1Dmod forward modeling package, and contains
% the code used to calculate the DC 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)
%
% 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;
if nargin < 3
disp('Not enough input parameters');
return
else
config = varargin{1};
layers = varargin{2};
cparams = varargin{3};
end
h = [layers.thickness]; h(end) = [];
config.r = [];
% setup config structure
switch config.type
case {'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) = 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 P1 a P2
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];
result(k,m).Aspac = config(k,m).Aspac(k);
result(k,m).Nspac = config(k,m).Nspac(m);
if result(k,m).Aspac == 0 || result(k,m).Nspac == 0
result(k,m).skip = 1;
else
result(k,m).skip = 0;
config(k,m).r = sqrt(sum(([ ... % calculate distance between electrodes
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));
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);
end % if
end % m
end % k
case {'Wenner'}
for k = 1:length(config(1,1).Aspac)
config(k) = config(1);
config(k).C1 = [-1.5.*config(k).Aspac(1) 0 0]; % transmitter electrode C1 P1 P2 C2
config(k).C2 = [1.5.*config(k).Aspac(1) 0 0]; % transmitter electrode o------o------o------o
config(k).P1 = [-0.5.*config(k).Aspac(1) 0 0]; % receiver electrode a a a
config(k).P2 = [0.5.*config(k).Aspac(1) 0 0]; % receiver electrode
result(k).Aspac = config(k).Aspac(k);
if result(k).Aspac == 0
result(k).skip = 1;
else
result(k).skip = 0;
config(k).r = sqrt(sum(([ ... % calculate distance between electrodes
config(k).P2-config(k).C2; ...
config(k).P2-config(k).C1; ...
config(k).P1-config(k).C2; ...
config(k).P1-config(k).C1]).^2,2));
result(k).G_factor = 2*pi*(1/config(k).r(1)- ...
1/config(k).r(2)-1/config(k).r(3)+ ...
1/config(k).r(4))^(-1);
end % if
end % k
case {'Schlumberger'}
for k = 1:length(config(1,1).OA)
for m = 1:length(config(1,1).OM)
config(k,m) = config(1,1);
config(k,m).C1 = [-config(k,m).OA(k) 0 0]; % transmitter electrode C1 P1 P2 C2
config(k,m).C2 = [config(k,m).OA(k) 0 0]; % transmitter electrode o--------o---o--------o
config(k,m).P1 = [-config(k,m).OM(m) 0 0]; % receiver electrode | OA |
config(k,m).P2 = [config(k,m).OM(m) 0 0]; % receiver electrode OM| |
result(k,m).OA = config(k,m).OA(k);
result(k,m).OM = config(k,m).OM(m);
if result(k,m).OA == result(k,m).OM
result(k,m).skip = 1;
else
result(k,m).skip = 0;
config(k,m).r = sqrt(sum(([ ... % calculate distance between electrodes
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));
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);
end % if
end % m
end % k
otherwise
config.r = sqrt(sum(([ ... % calculate distance between electrodes
config.P2-config.C2; ...
config.P2-config.C1; ...
config.P1-config.C2; ...
config.P1-config.C1]).^2,2));
if any(config.r) == 0
result.skip = 1;
else
result.skip = 0;
result.G_factor = 2*pi*(1/config.r(1)- ...
1/config.r(2)-1/config.r(3)+ ...
1/config.r(4))^(-1);
result.Z = [];
end
end
r = unique([config.r]);
if isfinite(layers(1).thickness)
% Do the calculations
Q = Q_DC(r, [layers.rho], h, cparams);
% Distribute results
for k = 1:prod(size(config))
if ~result(k).skip
result(k).Z = Q(r==config(k).r(1)) - ...
Q(r==config(k).r(2)) - ...
Q(r==config(k).r(3)) + ...
Q(r==config(k).r(4));
else
result(k).Z = NaN;
result(k).G_factor = 1;
end % if
end % k
else
% Distribute results
for k = 1:prod(size(config))
if ~result(k).skip
result(k).Z = layers(1).rho./result(k).G_factor;
else
result(k).Z = NaN;
result(k).G_factor = 1;
end % if
end % k
end % if
if length(config) == 1
disp(['Resistivity: ' num2str(result.Z.*result.G_factor) ' Ohmm']);
end
disp(['Time spent: ' num2str(cputime-cpu_t) ' s']);
% ******************************************************
% * Subfunctions *
% ******************************************************
% ------------------------------------------------------
function Q = Q_DC(r, rho, h, cparams)
for k = 1:length(r)
if strcmpi(cparams.hank_type,'NHT')
K(k) = NJCST('J0', @K_DC, r(k), cparams.Seg_tol, cparams.NHT_tol,...
cparams.Max_seg, 1./rho, h);
else
K(k) = FJCST('J02', r(k), -1, 0, 1, cparams.FHT_err, @K_DC, ...
1./rho, h);
end % if
end % k
Q = rho(1)./(2.*pi).*(1./r.' + K);
% ------------------------------------------------------
function Out = K_DC(lambda, sigma, d);
K = ones(size(d,2)+1,length(lambda));
for k = length(d):-1:1
K(k,:) = (K(k+1,:)+sigma(k+1)/sigma(k).*tanh(lambda.*d(k)))./ ...
(sigma(k+1)/sigma(k)+K(k+1,:).*tanh(lambda.*d(k)));
end
Out = K(1,:)-1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -