From b59e1b03e10c642cb548fac40d2f74ed8d8f1108 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=88=98=E6=98=8E=E5=AE=8F?= Date: Sat, 3 Aug 2024 15:01:25 +0000 Subject: [PATCH] =?UTF-8?q?=E8=AE=A1=E7=AE=97MT=E5=93=8D=E5=BA=94?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: 刘明宏 --- MT1D/get_MTR.m | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 MT1D/get_MTR.m diff --git a/MT1D/get_MTR.m b/MT1D/get_MTR.m new file mode 100644 index 0000000..46cf374 --- /dev/null +++ b/MT1D/get_MTR.m @@ -0,0 +1,55 @@ +function [MTR] = get_MTR(Ex1,Ey1,Hx1,Hy1,Hz1,Ex2,Ey2,Hx2,Hy2,Hz2,freq) + +% Magnetotelluric response calculation +% Return: impedance, apparent resistivity, phase, tipper response + +mu0 = 4*pi*1e-7 ; % V·s/ (A·m) + + +omega = 2*pi*freq; +detH = (Hx1.*Hy2-Hx2.*Hy1); + +%% impedance +Zxx = (Ex1.*Hy2 - Ex2.*Hy1)./detH ; +Zxy = (Ex2.*Hx1 - Ex1.*Hx2)./detH ; +Zyx = (Ey1.*Hy2 - Ey2.*Hy1)./detH ; +Zyy = (Ey2.*Hx1 - Ey1.*Hx2)./detH ; + +MTR.Zxx = Zxx; +MTR.Zxy = Zxy; +MTR.Zyx = Zyx; +MTR.Zyy = Zyy; + +%% apparent resistivity +rhoXX = 1/(omega*mu0)*abs(Zxx).^2; +rhoXY = 1/(omega*mu0)*abs(Zxy).^2; +rhoYX = 1/(omega*mu0)*abs(Zyx).^2; +rhoYY = 1/(omega*mu0)*abs(Zyy).^2; + +MTR.rhoXX = rhoXX; +MTR.rhoXY = rhoXY; +MTR.rhoYX = rhoYX; +MTR.rhoYY = rhoYY; +%% phase +phaXX = -atan(imag(Zxx)./real(Zxx))*180/pi; +phaXY = -atan(imag(Zxy)./real(Zxy))*180/pi; +phaYX = -atan(imag(Zyx)./real(Zyx))*180/pi; +phaYY = -atan(imag(Zyy)./real(Zyy))*180/pi; + +% phaXX = -atan2(imag(Zxx),real(Zxx))*180/pi; +% phaXY = -atan2(imag(Zxy),real(Zxy))*180/pi; +% phaYX = -atan2(imag(Zyx),real(Zyx))*180/pi; +% phaYY = -atan2(imag(Zyy),real(Zyy))*180/pi; +MTR.phaXX = phaXX; +MTR.phaXY = phaXY; +MTR.phaYX = phaYX; +MTR.phaYY = phaYY; + +%% Tipper +Tzx = (Hz1.*Hy2 - Hz2.*Hy1)./detH ; +Tzy = (Hz2.*Hx1 - Hz1.*Hx2)./detH ; +MTR.Tzx = Tzx; +MTR.Tzy = Tzy; + +end +