磁法正演matlab实现

clear;x=-50:1:50;b=1;h=20;B1=90;I0s=90;Is=90;Ms=400;r1=B1-I0s;B2=135;r2=B2-I0s;E3=90+r1-I0s;E2=90+r2-I0s;Hax1=2*Ms*b*sin(B1)*(-h*sin(r1)-x.*cos(r1))./(x.^2+h^2);Za1=2*Ms*b*sin(B1)*(h*cos(r1)-x.*sin(r1))./(x.^2+h^2);T1=2*Ms*b*sin(B1)*(h*cos(E3)-x.*sin(E3))./(x.^2+h^2);Hax2=2*Ms*b*sin(B2)*(-h*sin(r2)-x.*cos(r2))./(x.^2+h^2);Za2=2*Ms*b*sin(B2)*(h*cos(r2)-x.*sin(r2))./(x.^2+h^2);T2=2*Ms*b*sin(B2)*(h*cos(E2)-x.*sin(E2))./(x.^2+h^2);figure(1);subplot(2,1,1);plot(x,Hax1(:),'b*',x,Za1(:),'r-',x,T1(:),'x');legend('Hax1','Za1','T1');title('垂直磁化直立无限延深板状体磁异常剖面特征：B=90Is=90');subplot(2,1,2);plot(x,Hax2(:),'*',x,Za2(:),'+',x,T2(:),'-');legend('Hax2','Za2','T2');title('垂直磁化倾斜无限延深板状体磁异常剖面特征：B=135Is=90');ans='OK!';