function [rho_a,phase]=MT1D_FWD(rho,h,T)
mu=(4e-7)*pi;
%T=logspace(-3,4,40);
%T=logspace(-3,3,20);
k=zeros(size(rho,2),size(T,2));
for N=1:size(rho,2)
k(N,:)=sqrt(-i*2*pi*mu./(T.*rho(N))); %k矩阵存储各个层的复波数；mu后的点表示这里是点除运算，用mu去除以数组或矩阵里每一个元素;T.*表示用T的每一个元素乘以rho
%k(N,:)=-i*2*pi*mu./(T.*rho(N));
end
m=size(rho,2);
Z=-(i*mu*2*pi)./(T.*k(m,:)); %底层阻抗
for n=m-1:-1:1
A=-(i*2*pi*mu)./(T.*k(n,:)); %A相当于书中公式的Z0m,B相当于公式中的： e^-2kmhm
B=exp(-2*k(n,:)*h(n));
Z=A.*(A.*(1-B)+Z.*(1+B))./(A.*(1+B)+Z.*(1-B)); %书中求Zm的公式按照上面两行定义的A和B代入化简后就是这个形式；
end
rho_a=(T./(mu*2*pi)).*(abs(Z).^2);
rho_a=roundn(rho_a,-4);
%phase=-atan(imag(Z)./real(Z)).*180/pi;
%phase=angle(Z).*180/pi
%rho_1=T.*((Z.^2)./(-i*mu*2*pi)) %全视电阻率
%phase=angle(rho_1).*180/pi
phase=angle(Z).*180/pi
phase=roundn(phase,-4);
