[matlab]斜阶梯式长条引起的重力异常计算

画的是《重力学与固体潮》那本书的斜阶梯的重力异常

function [gravity]=draw_gravity_dixing(h1,h2,alpha,xigma,x0,z,direction)
%h1=100;%上表面
%h2=200;%下表面
%x0=0;%移动原点位置，向右为正
%xigma=3*10^3;%剩余密度
pi=3.1415926535
%alpha=alpha*pi;%倾角

%采样区间%
x=(-5000:50:5000);
y=(-5000:50:5000);
%z=-100;%测点高度
G=6.67e-11;
[x2,y2]=meshgrid(x,y);   %生成网线节点矩阵
if direction==1
    gravity=G*xigma*(pi*(h2-h1)+2*h2.*atan(((x2-x0)+h2*cot(alpha))/h2)-2*h1.*atan(((x2-x0)+h1*cot(alpha))/h1)+(x2-x0)*sin(alpha)*sin(alpha).*log(((h2+(x2-x0)*sin(alpha)*cos(alpha)).^2+(x2-x0).^2*sin(alpha)^4)./(((h1+(x2-x0)*sin(alpha)*cos(alpha)).^2+(x2-x0).^2*sin(alpha)^4)))-2.*(x2-x0)*sin(alpha)*cos(alpha).*atan(((x2-x0)*(h2-h1)*sin(alpha)*sin(alpha))./((x2-x0).^2*sin(alpha)^2+(h1+h2).*(x2-x0)*sin(alpha)*cos(alpha)+h1*h2)));
end
if direction==2
    gravity=G*xigma*(pi*(h2-h1)+2*h2.*atan(((x0-x2)+h2*cot(alpha))/h2)-2*h1.*