立方体引起的引力异常计算＆画图

地球重力学需要我们计算立方体引起的重力异常，公式见《重力学与固体潮》。
这个程序取的z方向是竖直向下的，也就是说地面向下为正，地面向上为负

%定义一个立方体函数%
function [gravity]=draw_square(a,b,c,x0,y0,H,ph,z)
%长方体模型参数说明%
%a=2000;%长%b=200;%宽%c=100;%高%
%质心坐标x0,y0,z0 %H=1000立方体深埋深度;
%质心埋深H
%ph=2*10^3;%剩余密度
%z是测点的z坐标，比如如果在地面上就取0%
%采样区间%
x=(-40:2:40);
y=(-40:2:40);

%常数%
G=6.67e-11;

%计算异常%
[x1,y1]=meshgrid(x,y);   %生成计算用的网格线
r=(x0-x1).^2+(y0-y1).^2+(H-z).^2;
gravity=-G*ph.*(a.*log(r+b)+b.*log(r+a)-c.*atan(a*b./(r*c)))*10^5;%单位mGal

%画图%
%矩阵要和前面的那个矩阵取一样的%
x=(-40:2:40);
y=(-40:2:40);

[x1,y1]=meshgrid(x,y);%生成画图用的矩阵
figure(1)%图1
mesh(x1,y1,gravity);%三维
colorbar;
xlabel('x');
ylabel('y');
title('立方体异常');