重力异常正演matlab程序,圆锥体完全均衡下重力异常正演 [MATLAB]

在完全均衡的模型下，若地表有一圆锥体(山峰等)，计算跨越山顶的截面上所得到的各种重力异常。

地壳密度 $kg\cdot m^{-3}$

上地幔密度 $g\cdot cm^{-3}$

地表地形圆锥体半径 (km)

地表地形圆锥体高度 (km)

计算莫霍面形变圆锥半径 (km)

计算莫霍面形变圆锥高度 (km)

地壳厚度 (km)

$2.8\times 10^{3}$

$3.5\times 10^{3}$

$2.0$

$2.0$

$2.0$

$8.0$

$30.0$

计算结果如下。横坐标单位：m，纵坐标单位：mGal

MATLAB代码如下：

% 生成地下体的布格重力异常

syms r a z x h;

f = r / sqrt((z + h)^2 + r^2 + x^2 - 2*r*x*cos(a))^3;

dense = - 700; G = 6.67E-11;

depth = 30000; subheight = 8000; height = 2000;

total_spots = 81;

total_anom = zeros(1, 81);

total_xvec = zeros(1, 81);

spots = 20; from = 50000; to = 2500; interval = -2500;

xvec = from:interval:to;

anom = zeros(1, spots);

for no = 1:spots

rad = from + (no - 1)*interval;

gap = 800;

N = subheight/gap;

anomaly = 0;

for n = 0:N-1

Radius = 2000 - (gap/4)*n;

fc = subs(f, [z, x, h], [depth + gap*n, rad, 0]);

func =