matlab球体重力异常正演

% 球体理论重力异常
clc;
clear;

% 球体参数
R=50; % 球体半径
p=1000; % 剩余密度 kg/m3
M=p*（4*pi*R^3）/3; % 剩余质量
G=6.67e-11; % 万有引力常数
D=100; % 球心深度

% 测点分布范围
dx=10; % X方向测点间距
dy=10; % Y方向测点间距
nx=41; % X方向测点数
ny=41; % Y方向测点数
xmin=-200; % X方向起点
ymin=-200; % Y方向起点
x=xmin:dx:（xmin+（nx-1）*dx）; % X方向范围
y=ymin:dy:（ymin+（ny-1）*dy）; % Y方向范围
[XY]=meshgrid（xy）; % 转化为排列

% 球体理论重力异常
data01=G*M*D./（（（x）.^2+D^2）.^（3/2））;
data01=data01*1e6;
data02=G*M*D./（（（X）.^2+Y.^2+D^2）.^（3/2））;
data02=data02*1e6;
figure（1）clf
subplot（331）plot（xdata01‘k‘xdata01‘b*‘）xlabel（‘X（m）‘）ylabel（‘重力异常 （g.u.）‘）grid on;
% subplot（312）contour（XYdata0220）xlabel（‘X（m）‘）ylabel（‘Y（m）‘）title（‘三度球体重力异常‘）
% axis equalaxis（[-200 200 -200 200]）colorbar;
subp