用MATLAB编写正方体重磁正演并画出重力异常曲面的可执行程序

好的，下面是一个简单的正方体重磁正演的 MATLAB 代码，并画出重力异常曲面的程序，你可以根据自己的需求进行修改和完善。

% 正方体重磁正演程序
% 输入参数：
%   rho: 立方体密度，单位 kg/m^3
%   a:  立方体边长，单位 m
%   x0,y0,z0: 立方体中心坐标，单位 m
%   dx,dy,dz: 计算网格的间距，单位 m
%   Nx,Ny,Nz: 计算网格的数量
% 输出参数：
%   G: 重力异常，单位 mGal
%   B: 磁异常，单位 nT

% 计算网格坐标
[X,Y,Z] = meshgrid((0:Nx-1)*dx+x0-a/2, (0:Ny-1)*dy+y0-a/2, (0:Nz-1)*dz+z0-a/2);

% 计算欧拉距离
r = sqrt((X-x0).^2 + (Y-y0).^2 + (Z-z0).^2);

% 计算重力异常
G = 2*pi*G*rho*(a/2)^2*(log((r+a/2)./r)-1/2);

% 计算磁异常
m = 1; % 磁化强度，单位 A/m
sus = 0.01; % 磁化率，单位 SI
B = zeros(Nx,Ny,Nz); % 初始化磁异常矩阵
for i = 1:Nx
    for j = 1:Ny
        for k = 1:Nz
            r_vec = [X(i,j,k)-x0, Y(i,j,k)-y0, Z(i,j,k)-z0];
            B_vec = m/(4*pi)*((3*dot(r_vec,m)-norm(m)^2)/norm(r_vec)^5*r_vec + sus/norm(r_vec)^3*m);
            B(i,j,k) = norm(B_vec);
        end
    end
end

% 将结果转换为 mGal 和 nT
G = G*1e5;
B = B*1e9;

% 绘制重力异常曲面
figure;
surf(X(:,:,1), Y(:,:,1), G(:,:,1));
xlabel('x (m)');
ylabel('y (m)');
zlabel('g (mGal)');
title('重力异常曲面');
colorbar;

% 绘制磁异常曲面
figure;
surf(X(:,:,1), Y(:,:,1), B(:,:,1));
xlabel('x (m)');
ylabel('y (m)');
zlabel('B (nT)');
title('磁异常曲面');
colorbar;

你可以将以上代码保存为一个.m文件并运行，按照提示输入参数即可计算出相应的重力异常和磁异常，并绘制出相应的曲面图。注意，以上代码仅仅是一个简单的示例，实际应用中还需要考虑许多其他的因素，比如地球引力场、地形变化、地磁场变化等等。