用matlab编写一段代码：地球物理勘探中长方体模型的重磁正演，分别得到重力和磁异常的三维图

好的，我可以为你提供一段MATLAB代码实现地球物理勘探中长方体模型的重磁正演，并分别得到重力和磁异常的三维图。请注意，这段代码仅供参考，需要根据具体需求进行修改和调整。

% 地球物理勘探中长方体模型的重磁正演
% 作者：CSDN开发的ChitGPT

% 定义长方体模型参数
% 长度、宽度、深度
L = 5000;
W = 5000;
D = 500;

% 每个网格的大小
dx = 50;
dy = 50;
dz = 50;

% 定义观测点
% 观测点的x、y、z坐标
xobs = [-2500:dx:2500];
yobs = [-2500:dy:2500];
zobs = 100;

% 定义地球磁场强度
B0 = 50000;

% 定义磁化率和密度
k = 0.05;
rho = 2.67;

% 初始化重力和磁异常
grav = zeros(length(xobs), length(yobs));
mag = zeros(length(xobs), length(yobs));

% 计算重力和磁异常
for i = 1:length(xobs)
    for j = 1:length(yobs)
        % 计算重力和磁异常的贡献
        for ii = 1:floor(L/dx)
            for jj = 1:floor(W/dy)
                for kk = 1:floor(D/dz)
                    % 计算长方体的中心点坐标
                    xc = (ii-0.5)*dx - L/2;
                    yc = (jj-0.5)*dy - W/2;
                    zc = (kk-0.5)*dz;
                    
                    % 计算观测点到长方体中心点的距离
                    r = sqrt((xobs(i)-xc)^2 + (yobs(j)-yc)^2 + (zobs-zc)^2);
                    
                    % 计算重力和磁异常的贡献
                    grav(i,j) = grav(i,j) + rho*dx*dy*dz*9.81*(zobs-zc)/r^3;
                    mag(i,j) = mag(i,j) + B0*k*dx*dy*dz*(zobs-zc)/r^3;
                end
            end
        end
    end
end

% 绘制三维图
[X,Y] = meshgrid(xobs, yobs);
surf(X,Y,grav);
title('重力异常三维图');
xlabel('x (m)');
ylabel('y (m)');
zlabel('重力异常 (mGal)');

figure;
surf(X,Y,mag);
title('磁异常三维图');
xlabel('x (m)');
ylabel('y (m)');
zlabel('磁异常 (nT)');

请注意，这段代码并没有进行错误处理和边界处理，需要根据具体需求进行修改和调整。