%?球體理論重力異常
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