diff --git a/bin/RMT1D/MakeGif.m b/bin/RMT1D/MakeGif.m new file mode 100644 index 0000000..b691c1b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/bin/RMT1D/exmple/MT1Dexmple2.mlx b/bin/RMT1D/exmple/MT1Dexmple2.mlx new file mode 100644 index 0000000..3b7e986 Binary files /dev/null and b/bin/RMT1D/exmple/MT1Dexmple2.mlx differ diff --git a/bin/RMT1D/exmple/MT1Dexmple_ani_Pek.mlx b/bin/RMT1D/exmple/MT1Dexmple_ani_Pek.mlx new file mode 100644 index 0000000..ec5cb39 Binary files /dev/null and b/bin/RMT1D/exmple/MT1Dexmple_ani_Pek.mlx differ diff --git a/bin/RMT1D/exmple/RMT1Dexmple1.mlx b/bin/RMT1D/exmple/RMT1Dexmple1.mlx new file mode 100644 index 0000000..00eb4b9 Binary files /dev/null and b/bin/RMT1D/exmple/RMT1Dexmple1.mlx differ diff --git a/bin/RMT1D/exmple/RMT1Dexmple_Esp_sigz.mlx b/bin/RMT1D/exmple/RMT1Dexmple_Esp_sigz.mlx new file mode 100644 index 0000000..0e3db12 Binary files /dev/null and b/bin/RMT1D/exmple/RMT1Dexmple_Esp_sigz.mlx differ diff --git a/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp.mlx b/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp.mlx new file mode 100644 index 0000000..2956abe Binary files /dev/null and b/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp.mlx differ diff --git a/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_ChenHuang.mlx b/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_ChenHuang.mlx new file mode 100644 index 0000000..8dfafe6 Binary files /dev/null and b/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_ChenHuang.mlx differ diff --git a/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx b/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx new file mode 100644 index 0000000..b79bc2a Binary files /dev/null and b/bin/RMT1D/exmple/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx differ diff --git a/bin/RMT1D/exmple/phase_1000ohm_eps5_e4_7.fig b/bin/RMT1D/exmple/phase_1000ohm_eps5_e4_7.fig new file mode 100644 index 0000000..a89316d Binary files /dev/null and b/bin/RMT1D/exmple/phase_1000ohm_eps5_e4_7.fig differ diff --git a/bin/RMT1D/exmple/rho_1000ohm_eps5_e4_7.fig b/bin/RMT1D/exmple/rho_1000ohm_eps5_e4_7.fig new file mode 100644 index 0000000..648caa0 Binary files /dev/null and b/bin/RMT1D/exmple/rho_1000ohm_eps5_e4_7.fig differ diff --git a/bin/RMT1D/getEffSigma.m b/bin/RMT1D/getEffSigma.m new file mode 100644 index 0000000..a0041e0 --- /dev/null +++ b/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 diff --git a/bin/RMT1D/mt1DAnalyticFields.mlx b/bin/RMT1D/mt1DAnalyticFields.mlx new file mode 100644 index 0000000..07c7db2 Binary files /dev/null and b/bin/RMT1D/mt1DAnalyticFields.mlx differ diff --git a/bin/RMT1D/mt1DAnalyticFieldsEpsMu.mlx b/bin/RMT1D/mt1DAnalyticFieldsEpsMu.mlx new file mode 100644 index 0000000..09682a0 Binary files /dev/null and b/bin/RMT1D/mt1DAnalyticFieldsEpsMu.mlx differ diff --git a/bin/RMT1D/mt1DAniAnalyticField.m b/bin/RMT1D/mt1DAniAnalyticField.m new file mode 100644 index 0000000..8aba2de --- /dev/null +++ b/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 + + + + diff --git a/bin/RMT1D/mt1DAniAnalyticFieldEspMu.m b/bin/RMT1D/mt1DAniAnalyticFieldEspMu.m new file mode 100644 index 0000000..cca3ded --- /dev/null +++ b/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 + + + + diff --git a/bin/RMT1D/plot2d_Conduction_displacement_current_ratio.mlx b/bin/RMT1D/plot2d_Conduction_displacement_current_ratio.mlx new file mode 100644 index 0000000..193b999 Binary files /dev/null and b/bin/RMT1D/plot2d_Conduction_displacement_current_ratio.mlx differ diff --git a/bin/RMT1D/plot3d_Conduction_displacement_current_ratio.mlx b/bin/RMT1D/plot3d_Conduction_displacement_current_ratio.mlx new file mode 100644 index 0000000..a3fddbc Binary files /dev/null and b/bin/RMT1D/plot3d_Conduction_displacement_current_ratio.mlx differ