Matlab画任意形体的重力异常三维图、剖面或平面图

%剩余密度rd=2e14
%重力常数G =6.67e-11                        
%测量平面z =-10   
%测线x最大值max_x=100
%测线y最大值max_y=100

%画由长方体组成的各种不规则形体程序
clear;clc

%画一个小正方体
x=10;y=10;z=10;%角点
l=10;w=10;h=10;%长，宽，高
[a,b,c]=meshgrid([0 1]);
p=alphaShape(l*a(:)-(l-x),w*b(:)-(w-y),h*c(:)-(0-z));
plot(p,'edgecolor','none')
xlabel('x');ylabel('y');zlabel('z');
camlight
grid on;

%画长棱柱1
for x = 5:2:95
    for y = 5:2:95
        if x+y==60
            z = 30
            l=5;w=5;h=5
            hold on;
            [a,b,c]=meshgrid([0 1]);
            p=alphaShape(l*a(:)-(l-x),w*b(:)-(w-y),h*c(:)-(0-z));
            plot(p,'edgecolor','none')
            xlabel('x');ylabel('y');zlabel('z');
        end
    end
end

%画长棱柱2
for x = 5:1:95
    for y = 5:1:95
        if x+y==60
            z = 20
            l=10;w=10;h=10
            hold on;
            [a,b,c]=meshgrid([0 1]);
            p=alphaShape(l*a(:)-(l-x),w*b(:)-(w-y),h*c(:)-(0-z));
            plot(p,'edgecolor','none')
            xlabel('x');ylabel('y');zlabel('z');
        end
    end
end

%画球
hold on;
[m1,n1,k1]=sphere(20)
for i =1:1:20
    for j = 1:1:20
        x=m1(i,j)*20+70
        y=n1(i,j)*20+70
        z=k1(i,j)*20+70
        l=8;w=8;h=8
        hold on;
        [a,b,c]=meshgrid([0 1]);
        p=alphaShape(l*a(:)-(l-x),w*b(:)-(w-y),h*c(:)-(0-z));
        plot(p,'edgecolor','none')
        xlabel('x');ylabel('y');zlabel('z');
    end
end

%画随机图形
for x = 5:10:50
    for y = 50:10:95
        z = 40
        l = rand(1,1)*5
        w = rand(1,1)*5
        h = rand(1,1)*5 %长宽高随机
        hold on;
        [a,b,c]=meshgrid([0 1]);
        p=alphaShape(l*a(:)-(l-x),w*b(:)-(w-y),h*c(:)-(0-z));
        plot(p,'edgecolor','none')
        xlabel('x');ylabel('y');zlabel('z');
    end
end

%画圆
[m2,n2,k2]=sphere(20)
for i =1:1:20
    for j = 1:1:20
        if (k2(i,j)*20+40<42) &&(k2(i,j)*20+40>38)
            x=m2(i,j)*20+30
            y=n2(i,j)*20+60
            z=k2(i,j)*20+40
            l=6;w=6;h=3
            hold on;
            [a,b,c]=meshgrid([0 1]);
            p=alphaShape(l*a(:)-(l-x),w*b(:)-(w-y),h*c(:)-(0-z));
            plot(p,'edgecolor','none')
            xlabel('x');ylabel('y');zlabel('z');
        end
    end
end
axis([0,100,0,100,0,100]);

%计算程序
 
rd=2e14
G =6.67e-11                        
z =-10   
max_x=100
max_y=100
fid=fopen('重力异常.txt','w');
for x = 0:5:max_x
    for y = 0:5:max_y
        k = @(a,b,c)(c-z)./((a-x).^2 + (b-y).^2 + (c-z).^2).^(3/2)
        
        %小正方体
        smallcube_value = 0
        smallcube_value = integral3(k,0,10,0,10,10,20);
        
        %球体
        sphere_value=0;
        for i =1:1:20
            for j = 1:1:20
                sphere_value =  sphere_value + integral3(k,m1(i,j)*20+70-8,m1(i,j)*20+70,n1(i,j)*20+70-8,n1(i,j)*20+70,k1(i,j)*20+70,k1(i,j)*20+70+8) 
            end
        end
        
        %圆形体
        circular_value=0;
        [m2,n2,k2]=sphere(20)
        for i =1:1:20
            for j = 1:1:20
                if (k2(i,j)*20+40<42) &&(k2(i,j)*20+40>38)
                    circular_value =  circular_value + integral3(k,m2(i,j)*20+30-6,m2(i,j)*20+30,n2(i,j)*20+60-6,n2(i,j)*20+60,k2(i,j)*20+40,k2(i,j)*20+40+3) 
                end
            end
        end
 
        %长棱柱1
        prism1_value=0
        for i = 5:2:95
            for j = 5:2:95
                if i+j==60
                    k3 = 30
                    prism1_value = prism1_value + integral3(k,i-5,i,j-5,j,k3,k3+5);
                end
            end
        end
                    
        %长棱柱2
        prism2_value=0
        for i = 5:1:95
            for j = 5:1:95
                if i+j==60
                    k2 = 20
                    prism2_value = prism2_value + integral3(k,i-10,i,j-10,j,k2,k2+10);
                end
            end
        end
        
        eventual_value = smallcube_value + sphere_value + circular_value + prism1_value + prism2_value;
        gravity = G.*rd.*eventual_value;
        disp(gravity)
        fprintf(fid,'%5.2f\n',gravity);
    end
end
fclose(fid);
 
%画图程序，可与计算程序分开运行
%clear;clc;close all                   %每个代码都会有的清空数据
gravity = importdata('重力异常.txt'); %读取数据，是个向量
gravity=reshape(gravity,21,21);     %改成预设大小的矩阵
[X, Y] = meshgrid(1:21,1:21);       %建立网格
gravity_profilex=1:1:(size(gravity,1))%剖面x数据
gravity_profiley=1:1:(size(gravity,2))%剖面y数据
x=gravity(:,4)%测线x=20数据
y=gravity(14,:)%测线y=70数据
 
figure(2)%画三维重力异常图
mesh(X,Y,gravity)
xlabel('x')
ylabel('y')
zlabel('value')
title('Gravity anomaly')
colorbar
 
figure(3)%画平面等值线图
%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(4)%画剖面图
plot(gravity_profilex,gravity)
xlabel('y')
ylabel('value of line(x=1.2...)')
title('Ga of line(x=1.2...)')
figure(5)
gravity=gravity.'
plot(gravity_profiley,gravity)
xlabel('x')
ylabel('value of line(y=1.2...)')
title('Ga of line(y=1.2...)')
 
figure(6)%画经过长方体中心特殊测线图
plot(x,'LineWidth',2)
xlabel('x')
ylabel('value of line(y=35)')
title('Ga of survey line(y=50)')
figure(7)
plot(y,'LineWidth',2)
xlabel('y')
ylabel('value of line(x=50)')
title('Ga of survey line(x=50)')
%As randomgraph_value is lost,if you are interested in this, try to add it.

Figure1 异常体在地下分布情况

Z轴默认垂直指向地下，z值越小则越浅，越大则越深；比如球体最深，方块最浅。 

Figure2 三维重力异常

Figure3 重力异常平面等值线

 Figure4  重力异常y方向剖面

 

 Figure5  重力异常x方向剖面

  Figure6  测线y=14

   Figure7  测线x=4

 

 如果把网格设置为100*100，那么运行时间大约需要800分钟。