zhengyan.rar_Matlab 磁化_磁化 matlab

clc close all % 测点分布范围 dx=5; % X方向测点间距 dy=5; % Y方向测点间距 nx=81; % X方向测点数 ny=81; % Y方向测点数 xmin=-200; % X方向起点 ymin=-200; % Y方向起点 x=xmin:dx:(xmin+(nx-1)*dx); % X方向范围 y=ymin:dy:(ymin+(ny-1)*dy); % Y方向范围 [X,Y]=meshgrid(x,y); % 转化为排列 % 球体参数 a=0; %剖面磁方位角 r=20; % 球体半径 m v=4*pi*r^3; u=4*pi*10^(-7); %磁导率 M=0.7; %磁化强度 A/m m=M*v; %磁矩 R=30; % 球体埋深 m I=pi/3;%有效磁化倾角is % I=0:pi/6:pi/2; %有效磁化倾角is % Za=(u*m*((2*R.^2-(X-50).^2-Y.^2)*sin(i)-3*R*(X-50).*cos(i)*cos(a)-3*R*Y.*cos(i)*sin(a)))./(4*pi*((X-50).^2+Y.^2+R.^2).^(5/2)); hax=(u*m*((2*X.^2-Y.^2-R.^2)*cos(I)*cos(a)-3*R*X.*sin(I)+3*X.*Y.*cos(I)*sin(a)))./(4*pi*(X.^2+Y.^2+R.^2).^(5/2)); hay=(u*m*((2*Y.^2-X.^2-R.^2)*cos(I)*sin(a)-3*R*Y.*sin(I)+3*X.*Y.*cos(I)*cos(a)))./(4*pi*(X.^2+Y.^2+R.^2).^(5/2)); za=(u*m*((2*R.^2-X.^2-Y.^2)*sin(I)-3*R*X.*cos(I)*cos(a)-3*R*Y*cos(I)*sin(a)))./(4*pi*(X.^2+Y.^2+R.^2).^(5/2)); dt=hax.*cos(I)*cos(a)+hay.*cos(I)*sin(a)+za.*sin(I);%总磁场强度异常 %za对应的平面等值线图、曲面图及主剖面异常图； figure(1); subplot(121); mesh(X,Y,za); title('za曲面图'); subplot(122); contour(za,25) title('za平面等值线图'); grid on figure(5) % 改变I时za主剖面异常图 subplot(221); for j=1:4; I=0:pi/6:pi/2; %有效磁化倾角is za(:,:,j)=(u*m*((2*R.^2-X.^2-Y.^2)*sin(I(j))-3*R*X.*cos(I(j))*cos(a)-3*R*Y*cos(I(j))*sin(a)))./(4*pi*(X.^2+Y.^2+R.^2).^(5/2)); za1(j,:)=za(41,:,j); end plot(x,za1) grid on; legend('I=0','I=\pi/6','I=\pi/3','I=\pi/2'); %hax对应的平面等值线图、曲面图及主剖面异常图； title('改变I时za主剖面异常图 a=0 r=20 R=30'); xlabel('za主剖面方向距离x'); ylabel('za值大小'); grid on % 改变a时za主剖面异常图 subplot(222); for j=1:4; a=0:pi/6:pi/2; za(:,:,j)=(u*m*((2*R.^2-X.^2-Y.^2)*sin(I(3))-3*R*X.*cos(I(3))*cos(a(j))-3*R*Y*cos(I(3))*sin(a(j))))./(4*pi*(X.^2+Y.^2+R.^2).^(5/2)); za2(j,:)=za(41,:,j); end plot(x,za2) grid on; legend('a=0','a=\pi/6','a=\pi/3','a=\pi/2'); title('改变a时za主剖面异常图 I=pi/3 r=20 R=30'); xlabel('za主剖面方向距离x'); ylabel('za值大小'); grid on % 改变r时za主剖面异常图 subplot(223); m1=[0,0,0]; v1=[0,0,0]; for j=1:3; r=[10,15,20]; v1(j)=4*pi*r(j)^3; m1(j)=M*v1(j); za(:,:,j)=(u*m1(j)*((2*R.^2-X.^2-Y.^2)*sin(I(3))-3*R*X.*cos(I(3))*cos(a(1))-3*R*Y*cos(I(3))*sin(a(1))))./(4*pi*(X.^2+Y.^2+R.^2).^(5/2)); za3(j,:)=za(41,:,j); end plot(x,za3) grid on; legend('r=10','r=15','a=20'); title('改变r时za主剖面异常图 I=pi/3 a=0 R=30'); xlabel('za主剖面方向距离x'); ylabel('za值大小'); grid on % 改变R时za主剖面异常图 subplot(224); for j=1:3; R=[30,40,50]; za(:,:,j)=(u*m*((2*R(j).^2-X.^2-Y.^2)*sin(I(3))-3*R(j)*X.*cos(I(3))*cos(a(1))-3*R(j)*Y*cos(I(3))*sin(a(1))))./(4*pi*(X.^2+Y.^2+R(j).^2).^(5/2)); za4(j,:)=za(41,:,j); end plot(x,za4) grid on; legend('R=30','R=40','R=50'); title('改变R时za主剖面异常图 a=0 I=pi/3 r=20'); xlabel('za主剖面方向距离x'); ylabel('za值大小'); grid on %hax对应的平面等值线图、曲面图及主剖面异常图； figure(2); subplot(1,3,1); mesh(X,Y,hax); title('hax曲面图'); subplot(1,3,2); contour(hax,25) title('hax平面等值线图'); grid on subplot(1,3,3); hax1=hax(41,:); plot(x,hax1) title('hax主剖面异常图'); xlabel('hax主剖面方向距离x'); ylabel('hax值大小'); grid on %hay对应的平面等值线图、曲面图及主剖面异常图； figure(3); subplot(1,3,1); mesh(X,Y,hay); title('hay曲面图'); subplot(1,3,2); contour(hay,25) title('hay平面等值线图'); grid on subplot(1,3,3); hay1=hay(41,:); plot(x,hay1) title('hay主剖面异常图'); xlabel('hay主剖面方向距离x'); ylabel('hay值大小'); grid on %dt对应的平面等值线图、曲面图及主剖面异常图； figure(4); subplot(1,3,1); mesh(X,Y,dt); title('dt曲面图'); subplot(1,3,2); contour(dt,25) title('dt平面等值线图'); grid on subplot(1,3,3); dt1=dt(41,:); plot(x,dt1) title('dt主剖面异常图'); xlabel('dt主剖面方向距离x'); ylabel('dt值大小'); grid on