Matlab画长方体重力异常；三维图、剖面图、测线图

%剩余密度rd=2e14
%重力常数G =6.67e-11                        
%测量平面z =-10   
%长方体长、宽、高起终点a1=30;a2=70;b1=30;b2=40;c1=10;c2=30
%测线x最大值max_x=100
%测线y最大值max_y=100

%计算程序
function [draw_function]=Gravity_anomaly(G,z,rd,a1,a2,b1,b2,c1,c2,max_x,max_y)   
fid=fopen('重力异常.txt','w');
for x = 0:1:max_x
    for y = 0:1:max_y
        k = @(a,b,c)(c-z)./((a-x).^2 + (b-y).^2 + (c-z).^2).^(3/2)
        gravity = G.*rd.*integral3(k,a1,a2,b1,b2,c1,c2);
        
        disp(gravity)
        fprintf(fid,'%5.2f\n',gravity);
    end
end
fclose(fid);

%画图程序，可与计算程序分开运行
%clear;clc;close all                   %每个代码都会有的清空数据
gravity = importdata('重力异常.txt'); %读取数据，是个向量
gravity=reshape(gravity,max_x+1,max_y+1);     %改成预设大小的矩阵
[X, Y] = meshgrid(1:max_x+1,1:max_y+1);       %建立网格
gravity_profilex=1:1:(size(gravity,1))%剖面x数据
gravity_profiley=1:1:(size(gravity,2))%剖面y数据
x=gravity(:,50)%测线x=50数据
y=gravity(35,:)%测线y=35数据

figure(1)%画三维重力异常图
mesh(X,Y,gravity)
xlabel('x')
ylabel('y')
zlabel('value')
title('Gravity anomaly')
colorbar

figure(2)%画平面等值线图
%gravity=flipud(gravity)
imagesc(gravity)
xlabel('x')
ylabel('y')
zlabel('value')
title('Ga of Plane contour')
%set(gca,'XAxisLocation','top') 
set(gca,'YDir','normal')
colorbar

figure(3)%画剖面图
plot(gravity_profilex,gravity)
xlabel('y')
ylabel('value of line(x=1.2...)')
title('Ga of line(x=1.2...)')
figure(4)
gravity=gravity.'
plot(gravity_profiley,gravity)
xlabel('x')
ylabel('value of line(y=1.2...)')
title('Ga of line(y=1.2...)')

figure(5)%画经过长方体中心特殊测线图
plot(x,'LineWidth',2)
xlabel('x')
ylabel('value of line(y=35)')
title('Ga of survey line(y=50)')
figure(6)
plot(y,'LineWidth',2)
xlabel('y')
ylabel('value of line(x=50)')
title('Ga of survey line(x=50)')

Figure1

Figure2

Figure3

Figure4

Figure5

Figure6