%第二题，第1问程序，程序名：MT
function [rho_s,phase]=MT(rho,h)
% rho 为各层电阻率； h 为各层厚度
lp=length(rho);
u0=4.0e-7*pi;
freq=logspace(log10(1.0e-4),log10(1.0e4),200);
lf=length(freq);
for a=lf:-1:1;
for j=lp-1:-1:1;
S(lp)=1;
W(a)=2*pi*freq(a);
K(j)=sqrt(-1i*W(a)*u0./rho(j));
T(j)=(1-exp(-2*1i*K(j)*h(j)))/(1+exp(-2*1i*K(j)*h(j)));
B(j)=sqrt(rho(j+1))*S(j+1)+sqrt(rho(j))*T(j);
C(j)=sqrt(rho(j))+sqrt(rho(j+1))*S(j+1)*T(j);
S(j)=B(j)/C(j);
end
rho_s(a)=(abs(S(1))^2)*rho(1); %视电阻率
Z(1)=sqrt(1i*W(a)*u0*rho(1)); %波阻抗
phase(a)=atan(imag(Z(1)*S(1))/real(Z(1)*S(1)))*180/pi; %阻抗相位
end
subplot(2,1,1);
semilogx(freq,rho_s);
xlabel('f (Hz)');ylabel('\rho_s(\Omega\cdotm)');
title('二层介质模型（G型）视电阻率理论曲线');
subplot(2,1,2);
semilogx(freq,phase);%x freq
xlabel('f (Hz)');ylabel('phase/(°)');
title('二层介质模型（G型）波阻抗相位曲线');
end
