openSource/em3d-mt
1
0
派生 0
镜像自地址 https://gitee.com/sduem/em3d-mt.git 已同步 2025-08-03 03:16:52 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

rmt1d

浏览代码 
Signed-off-by: 刘明宏 <liuminghong@mail.sdu.edu.cn>
这个提交包含在：
刘明宏
2025-04-29 07:05:27 +00:00
提交者 Gitee
父节点 4fdab7f7c5
当前提交 fefd9de858
共有 17 个文件被更改，包括 425 次插入 和 0 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

10
bin/RMT1D/MakeGif.m 普通文件
查看文件

@@ -0,0 +1,10 @@
function MakeGif(filename,index)
f = getframe(gcf);
imind = frame2im(f);
[imind,cm] = rgb2ind(imind,256);
if index==1
imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',1);
else
imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.5);
end
end

二进制
bin/RMT1D/exmple/MT1Dexmple2.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/MT1Dexmple_ani_Pek.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/RMT1Dexmple1.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/RMT1Dexmple_Esp_sigz.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/RMT1Dexmple_ani_Esp.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_ChenHuang.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/phase_1000ohm_eps5_e4_7.fig 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/exmple/rho_1000ohm_eps5_e4_7.fig 普通文件
查看文件

二进制文件未显示。

20
bin/RMT1D/getEffSigma.m 普通文件
查看文件

@@ -0,0 +1,20 @@
function A=getEffSigma(sigma,esp,freq)
%computes effective azimuthal anisotropic (see Josef Pek et al., 2002). conductivity
if nargin == 1
A.xx = sigma(:,1) - (sigma(:,5).^2) ./ sigma(:,3);
A.yy = sigma(:,2) - (sigma(:,6).^2) ./ sigma(:,3);
A.xy = sigma(:,4) - (sigma(:,5) .* sigma(:,6)) ./ sigma(:,3);
elseif nargin == 3
omega = 2 * pi * freq;
iom = 1i*omega;
y1 = sigma(:,1)-iom.*esp(:,1);
y2 = sigma(:,2)-iom.*esp(:,2);
y3 = sigma(:,3)-iom.*esp(:,3);
y4 = sigma(:,4)-iom.*esp(:,4);
y5 = sigma(:,5)-iom.*esp(:,5);
y6 = sigma(:,6)-iom.*esp(:,6);
A.xx = y1 - y5.^2 ./y3;
A.yy = y2 - y6.^2 ./y3;
A.xy = y4 - (y5.*y6) ./y3;
end

二进制
bin/RMT1D/mt1DAnalyticFields.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/mt1DAnalyticFieldsEpsMu.mlx 普通文件
查看文件

二进制文件未显示。

193
bin/RMT1D/mt1DAniAnalyticField.m 普通文件
查看文件

@@ -0,0 +1,193 @@
function [Ex, Ey, Ez, Hx, Hy]=mt1DAniAnalyticField(freqs,sigma,zNode,eTop)
%input:
% fre:are the computational frequencies
% sigma includes the conductivity and Euler angle of the calculation element
% znode is the node coordinate
% eTop is Calculation mode: XY,eTop=[1;0]
% YX,eTop=[0;1]
%output:
%Ex is the boundary electric field along the x direction
%Ey is the boundary electric field along the y direction
%Ez is the boundary electric field along the z direction
%Hx is the boundary magnetic field along the x direction
%Hy is the boundary magnetic field along the y direction
if size(sigma,1)~=length(zNode)-1
fprintf('error')
end
h = diff(zNode);
% assume halfspace at the bottom of model
sigma = [sigma; sigma(end:end, :)];
% effective azimuthal anisotropic conductivity
A= getEffSigma(sigma);% get Axx, Axy, Ayx
nFreq = length(freqs);
nLayer = size(sigma, 1);
Ex = zeros(nLayer, nFreq);
Ey = zeros(nLayer, nFreq);
Ez = zeros(nLayer-1, nFreq);
Hx = zeros(nLayer, nFreq);
Hy = zeros(nLayer, nFreq);
sh = @(x) 0.5*(exp(x)-exp(-x)); % sinh()
ch = @(x) 0.5*(exp(x)+exp(-x)); % cosh()
MU0 = 4*pi*1e-7;
for iFreq=1:nFreq
omega = 2 * pi * freqs(iFreq);
iom = 1i*omega*MU0;
Sprod = zeros(4, 4, nLayer);
Sprod(:,:,end) = eye(4);
% loop over layers to get the accumulative product of S (i.e. Sprod).
for j = nLayer:-1:1
ada = A.xx(j) + A.yy(j);
add = A.xx(j) - A.yy(j);
if add>=0
dA12 = sqrt( add^2 + 4*A.xy(j)^2 );
elseif add<0
dA12 = -sqrt( add^2 + 4*A.xy(j)^2 );
end
A1 = 0.5 * (ada + dA12);
A2 = 0.5 * (ada - dA12);
kp = sqrt(-iom) * sqrt(A1);
km = sqrt(-iom) * sqrt(A2);
% xip, xim, Qp, and Qm for each layer are saved for later use.
xip(j) = -kp/iom;
xim(j) = -km/iom;
if A.xy(j) == 0
QpNeg(j) = 0;
QmNegRec(j) = 0;
else
QpNeg(j) = 2 * A.xy(j) / (add+dA12); % -Qp
QmNegRec(j) = 0.5 * (add-dA12) / A.xy(j); % -1/Qm
end
QpDerQm = QpNeg(j) * QmNegRec(j); % Q_p/Q_m
dq = 1-QpDerQm; % 1-Q_p/Q_m
% M of the basement
if j==nLayer
% 1st and 3rd columns won't be used, so set them to zeros
M = [ 0 1 0 QmNegRec(j); ...
0 QpNeg(j) 0 1; ...
0 -xip(j)*QpNeg(j) 0 -xim(j); ...
0 xip(j) 0 xim(j)*QmNegRec(j) ];
continue;
end
sp = sh(kp*h(j));
sm = sh(km*h(j));
cp = ch(kp*h(j));
cm = ch(km*h(j));
S(1,1) = cp-QpDerQm*cm;
S(1,2) = -QmNegRec(j)*(cp-cm);
S(1,3) = QmNegRec(j)*(sp/xip(j)-sm/xim(j));
S(1,4) = sp/xip(j)-QpDerQm*sm/xim(j);
S(2,1) = QpNeg(j)*(cp-cm);
S(2,2) = -QpDerQm*cp+cm;
S(2,3) = QpDerQm*sp/xip(j)-sm/xim(j);
S(2,4) = QpNeg(j)*(sp/xip(j)-sm/xim(j));
S(3,1) = -QpNeg(j)*(xip(j)*sp-xim(j)*sm);
S(3,2) = QpDerQm*xip(j)*sp-xim(j)*sm;
S(3,3) = -QpDerQm*cp+cm;
S(3,4) = -QpNeg(j)*(cp-cm);
S(4,1) = xip(j)*sp-QpDerQm*xim(j)*sm;
S(4,2) = -QmNegRec(j)*(xip(j)*sp-xim(j)*sm);
S(4,3) = QmNegRec(j)*(cp-cm);
S(4,4) = cp-QpDerQm*cm;
S = S/dq;
% disp(['Layer#: ', num2str(j)]);
% rcond(S)
Sprod(:,:,j) = S * Sprod(:,:,j+1);
end % nLayer
SM = Sprod(:,:,1) * M;
% assume Ex0 & Ey0 are known, compute Cm and Dm of the bottom layer
Ex0 = eTop(1);
Ey0 = eTop(2);
G = [SM(1,2) SM(1,4); SM(2,2) SM(2,4)];
% CD = inv(G) * [Ex0; Ey0];
CD = zeros(2,1);
CD(1) = ( G(2,2)*Ex0 - G(1,2)*Ey0)/det(G);
CD(2) = (-G(2,1)*Ex0 + G(1,1)*Ey0)/det(G);
% % Alternatively, one can assume Hx0 & Hy0 are known
% Hx0 = 1.;
% Hy0 = 1.;
% G = [SM(3,2) SM(3,4); SM(4,2) SM(4,4)];
% CD = inv(G) * [Hx0; Hy0];
% the full C of the bottom layer
Cbot = [0; CD(1); 0; CD(2)];
MC = M * Cbot;
for j=1:nLayer
F = Sprod(:,:,j) * MC;
Ex(j, iFreq) = F(1);
Ey(j, iFreq) = F(2);
Hx(j, iFreq) = F(3);
Hy(j, iFreq) = F(4);
if j>1
kp = -iom * xip(j-1);
km = -iom * xim(j-1);
QpDerQm = QpNeg(j-1)*QmNegRec(j-1);
dq = 1-QpDerQm;
dz = 0.5 * h(j-1);
sp = sh(kp*dz);
sm = sh(km*dz);
cp = ch(kp*dz);
cm = ch(km*dz);
S(1,1) = cp-QpDerQm*cm;
S(1,2) = -QmNegRec(j-1)*(cp-cm);
S(1,3) = QmNegRec(j-1)*(sp/xip(j-1)-sm/xim(j-1));
S(1,4) = sp/xip(j-1)-QpDerQm*sm/xim(j-1);
S(2,1) = QpNeg(j-1)*(cp-cm);
S(2,2) = -QpDerQm*cp+cm;
S(2,3) = QpDerQm*sp/xip(j-1)-sm/xim(j-1);
S(2,4) = QpNeg(j-1)*(sp/xip(j-1)-sm/xim(j-1));
S(3,1) = -QpNeg(j-1)*(xip(j-1)*sp-xim(j-1)*sm);
S(3,2) = QpDerQm*xip(j-1)*sp-xim(j-1)*sm;
S(3,3) = -QpDerQm*cp+cm;
S(3,4) = -QpNeg(j-1)*(cp-cm);
S(4,1) = xip(j-1)*sp-QpDerQm*xim(j-1)*sm;
S(4,2) = -QmNegRec(j-1)*(xip(j-1)*sp-xim(j-1)*sm);
S(4,3) = QmNegRec(j-1)*(cp-cm);
S(4,4) = cp-QpDerQm*cm;
S = S/dq;
% computes E-fields at layer center
SF = S * F;
Exc = SF(1);
Eyc = SF(2);
Ez(j-1, iFreq) = -(sigma(j-1, 5)*Exc + sigma(j-1, 6)*Eyc)/sigma(j-1, 3);
end
end
end

202
bin/RMT1D/mt1DAniAnalyticFieldEspMu.m 普通文件
查看文件

@@ -0,0 +1,202 @@
function [Ex, Ey, Ez, Hx, Hy]=mt1DAniAnalyticFieldEspMu(freqs,sigma,epsr,mur,zNode,eTop)
%input:
% fre:are the computational frequencies
% sigma includes the conductivity and Euler angle of the calculation element
% znode is the node coordinate
% eTop is Calculation mode: XY,eTop=[1;0]
% YX,eTop=[0;1]
%output:
%Ex is the boundary electric field along the x direction
%Ey is the boundary electric field along the y direction
%Ez is the boundary electric field along the z direction
%Hx is the boundary magnetic field along the x direction
%Hy is the boundary magnetic field along the y direction
if size(sigma,1)~=length(zNode)-1
fprintf('error')
end
h = diff(zNode);
MU0 = 4*pi*1e-7;
MU = MU0 * mur;
EPS0 = 8.854187817620389 * 1e-12 ;
EPS = epsr*EPS0;
% assume halfspace at the bottom of model
sigma = [sigma; sigma(end:end, :)];
EPS = [ EPS; EPS(end:end, :)];
MU = [ MU; MU(end:end, :)];
% effective azimuthal anisotropic conductivity
A = getEffSigma(sigma,EPS,freqs);% get Axx, Axy, Ayx
nFreq = length(freqs);
nLayer = size(sigma, 1);
Ex = zeros(nLayer, nFreq);
Ey = zeros(nLayer, nFreq);
% Ez = zeros(nLayer-1, nFreq);
Ez = zeros(nLayer-1, nFreq);
Hx = zeros(nLayer, nFreq);
Hy = zeros(nLayer, nFreq);
sh = @(x) 0.5*(exp(x)-exp(-x)); % sinh()
ch = @(x) 0.5*(exp(x)+exp(-x)); % cosh()
for iFreq=1:nFreq
omega = 2 * pi * freqs(iFreq);
iom = 1i*omega*MU;
Sprod = zeros(4, 4, nLayer);
Sprod(:,:,end) = eye(4);
% loop over layers to get the accumulative product of S (i.e. Sprod).
for j = nLayer:-1:1
ada = A.xx(j) + A.yy(j);
add = A.xx(j) - A.yy(j);
if add>=0
dA12 = sqrt( add^2 + 4*A.xy(j)^2 );
elseif add<0
dA12 = -sqrt( add^2 + 4*A.xy(j)^2 );
end
A1 = 0.5 * (ada + dA12);
A2 = 0.5 * (ada - dA12);
kp = sqrt(-iom(j)) * sqrt(A1);
km = sqrt(-iom(j)) * sqrt(A2);
% xip, xim, Qp, and Qm for each layer are saved for later use.
xip(j) = -kp/iom(j);
xim(j) = -km/iom(j);
if A.xy(j) == 0
QpNeg(j) = 0;
QmNegRec(j) = 0;
else
QpNeg(j) = 2 * A.xy(j) / (add+dA12); % -Qp
QmNegRec(j) = 0.5 * (add-dA12) / A.xy(j); % -1/Qm
end
QpDerQm = QpNeg(j) * QmNegRec(j); % Q_p/Q_m
dq = 1-QpDerQm; % 1-Q_p/Q_m
% M of the basement
if j==nLayer
% 1st and 3rd columns won't be used, so set them to zeros
M = [ 0 1 0 QmNegRec(j); ...
0 QpNeg(j) 0 1; ...
0 -xip(j)*QpNeg(j) 0 -xim(j); ...
0 xip(j) 0 xim(j)*QmNegRec(j) ];
continue;
end
sp = sh(kp*h(j));
sm = sh(km*h(j));
cp = ch(kp*h(j));
cm = ch(km*h(j));
S(1,1) = cp-QpDerQm*cm;
S(1,2) = -QmNegRec(j)*(cp-cm);
S(1,3) = QmNegRec(j)*(sp/xip(j)-sm/xim(j));
S(1,4) = sp/xip(j)-QpDerQm*sm/xim(j);
S(2,1) = QpNeg(j)*(cp-cm);
S(2,2) = -QpDerQm*cp+cm;
S(2,3) = QpDerQm*sp/xip(j)-sm/xim(j);
S(2,4) = QpNeg(j)*(sp/xip(j)-sm/xim(j));
S(3,1) = -QpNeg(j)*(xip(j)*sp-xim(j)*sm);
S(3,2) = QpDerQm*xip(j)*sp-xim(j)*sm;
S(3,3) = -QpDerQm*cp+cm;
S(3,4) = -QpNeg(j)*(cp-cm);
S(4,1) = xip(j)*sp-QpDerQm*xim(j)*sm;
S(4,2) = -QmNegRec(j)*(xip(j)*sp-xim(j)*sm);
S(4,3) = QmNegRec(j)*(cp-cm);
S(4,4) = cp-QpDerQm*cm;
S = S/dq;
% disp(['Layer#: ', num2str(j)]);
% rcond(S)
Sprod(:,:,j) = S * Sprod(:,:,j+1);
end % nLayer
SM = Sprod(:,:,1) * M;
% assume Ex0 & Ey0 are known, compute Cm and Dm of the bottom layer
Ex0 = eTop(1);
Ey0 = eTop(2);
G = [SM(1,2) SM(1,4); SM(2,2) SM(2,4)];
% CD = inv(G) * [Ex0; Ey0];
CD = zeros(2,1);
CD(1) = ( G(2,2)*Ex0 - G(1,2)*Ey0)/det(G);
CD(2) = (-G(2,1)*Ex0 + G(1,1)*Ey0)/det(G);
% % Alternatively, one can assume Hx0 & Hy0 are known
% Hx0 = 1.;
% Hy0 = 1.;
% G = [SM(3,2) SM(3,4); SM(4,2) SM(4,4)];
% CD = inv(G) * [Hx0; Hy0];
% the full C of the bottom layer
Cbot = [0; CD(1); 0; CD(2)];
MC = M * Cbot;
for j=1:nLayer
F = Sprod(:,:,j) * MC;
Ex(j, iFreq) = F(1);
Ey(j, iFreq) = F(2);
Hx(j, iFreq) = F(3);
Hy(j, iFreq) = F(4);
if j>1
kp = -iom(j-1) * xip(j-1);
km = -iom(j-1) * xim(j-1);
QpDerQm = QpNeg(j-1)*QmNegRec(j-1);
dq = 1-QpDerQm;
dz = 0.5 * h(j-1);
sp = sh(kp*dz);
sm = sh(km*dz);
cp = ch(kp*dz);
cm = ch(km*dz);
S(1,1) = cp-QpDerQm*cm;
S(1,2) = -QmNegRec(j-1)*(cp-cm);
S(1,3) = QmNegRec(j-1)*(sp/xip(j-1)-sm/xim(j-1));
S(1,4) = sp/xip(j-1)-QpDerQm*sm/xim(j-1);
S(2,1) = QpNeg(j-1)*(cp-cm);
S(2,2) = -QpDerQm*cp+cm;
S(2,3) = QpDerQm*sp/xip(j-1)-sm/xim(j-1);
S(2,4) = QpNeg(j-1)*(sp/xip(j-1)-sm/xim(j-1));
S(3,1) = -QpNeg(j-1)*(xip(j-1)*sp-xim(j-1)*sm);
S(3,2) = QpDerQm*xip(j-1)*sp-xim(j-1)*sm;
S(3,3) = -QpDerQm*cp+cm;
S(3,4) = -QpNeg(j-1)*(cp-cm);
S(4,1) = xip(j-1)*sp-QpDerQm*xim(j-1)*sm;
S(4,2) = -QmNegRec(j-1)*(xip(j-1)*sp-xim(j-1)*sm);
S(4,3) = QmNegRec(j-1)*(cp-cm);
S(4,4) = cp-QpDerQm*cm;
S = S/dq;
% computes E-fields at layer center
SF = S * F;
Exc = SF(1);
Eyc = SF(2);
Ez(j-1, iFreq) = -(sigma(j-1, 5)*Exc + sigma(j-1, 6)*Eyc)/sigma(j-1, 3);
end
end
Ez(j, iFreq) = Ez(j-1);
end

二进制
bin/RMT1D/plot2d_Conduction_displacement_current_ratio.mlx 普通文件
查看文件

二进制文件未显示。

二进制
bin/RMT1D/plot3d_Conduction_displacement_current_ratio.mlx 普通文件
查看文件

二进制文件未显示。