编写长方体磁异常正演matlab程序

长方体磁异常正演是一个三维问题，需要考虑多个因素，包括长方体的位置、大小、磁化强度以及地球磁场的倾角和偏角等。下面是一个较为详细的长方体磁异常正演Matlab程序示例：

clc; clear all; close all;

% 定义长方体的位置、大小和磁化强度
x0 = 0; y0 = 0; z0 = 1000; % 长方体中心的位置
L = 100; W = 50; H = 20; % 长方体的长、宽和高
inc = 60; dec = 30; % 长方体的地磁倾角和偏角
mag = 0.1; % 长方体的磁化强度

% 定义计算网格
nx = 101; ny = 101; nz = 101; % 网格数量
x = linspace(-1000, 1000, nx); % x方向网格坐标
y = linspace(-1000, 1000, ny); % y方向网格坐标
z = linspace(0, 2000, nz); % z方向网格坐标
[X, Y, Z] = meshgrid(x, y, z);

% 计算磁异常
Bx = zeros(nz, ny, nx);
By = zeros(nz, ny, nx);
Bz = zeros(nz, ny, nx);
for i = 1:nx
    for j = 1:ny
        for k = 1:nz
            [bx, by, bz] = calc_mag(x0, y0, z0, L, W, H, inc, dec, mag, X(k,j,i), Y(k,j,i), Z(k,j,i));
            Bx(k,j,i) = bx;
            By(k,j,i) = by;
            Bz(k,j,i) = bz;
        end
    end
end

% 绘图
figure;
slice(X, Y, Z, sqrt(Bx.^2+By.^2+Bz.^2), [], [], [0 1000 2000]);
shading interp;
colorbar;
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('长方体磁异常');

% 计算磁场分量的函数
function [bx, by, bz] = calc_mag(x0, y0, z0, L, W, H, inc, dec, mag, x, y, z)
    mu = 4 * pi * 1e-7;
    R = sqrt((x-x0).^2 + (y-y0).^2 + (z-z0).^2);
    theta = atan2(y-y0, x-x0);
    phi = atan2(sqrt((x-x0).^2 + (y-y0).^2), z-z0);
    bx = 0; by = 0; bz = 0;
    for i = 1:2
        for j = 1:2
            for k = 1:2
                [bx1, by1, bz1] = calc_cube(L/2, W/2, H/2, R, theta, phi, inc, dec, mag, i, j, k);
                bx = bx + bx1;
                by = by + by1;
                bz = bz + bz1;
            end
        end
    end
end

% 计算立方体上每个面元的磁场分量的函数
function [bx, by, bz] = calc_cube(a, b, c, R, theta, phi, inc, dec, mag, i, j, k)
    mu = 4 * pi * 1e-7;
    if i == 1
        x1 = -a; x2 = a;
    else
        x1 = a; x2 = -a;
    end
    if j == 1
        y1 = -b; y2 = b;
    else
        y1 = b; y2 = -b;
    end
    if k == 1
        z1 = -c; z2 = c;
    else
        z1 = c; z2 = -c;
    end
    bx = 0; by = 0; bz = 0;
    for ii = 1:2
        for jj = 1:2
            for kk = 1:2
                x = x1 + (ii-1) * (x2-x1);
                y = y1 + (jj-1) * (y2-y1);
                z = z1 + (kk-1) * (z2-z1);
                R1 = sqrt((R.*sin(phi).*cos(theta)-x).^2 + (R.*sin(phi).*sin(theta)-y).^2 + (R.*cos(phi)-z).^2);
                theta1 = atan2(R.*sin(phi).*sin(theta)-y, R.*sin(phi).*cos(theta)-x);
                phi1 = atan2(sqrt((R.*sin(phi).*cos(theta)-x).^2 + (R.*sin(phi).*sin(theta)-y).^2), R.*cos(phi)-z);
                [bx1, by1, bz1] = calc_prism(R1, theta1, phi1, inc, dec, mag, a, b, c);
                bx = bx + bx1;
                by = by + by1;
                bz = bz + bz1;
            end
        end
    end
end

% 计算棱柱上每个面元的磁场分量的函数
function [bx, by, bz] = calc_prism(R, theta, phi, inc, dec, mag, a, b, c)
    mu = 4 * pi * 1e-7;
    bx = 0; by = 0; bz = 0;
    if R > 1e-6
        B = mag * (a*b*c) / (4*pi*R.^3);
        bx = bx + B * (cos(phi)*cos(theta)*cos(dec) + sin(phi)*sin(dec));
        by = by + B * (cos(phi)*sin(theta)*cos(dec) - sin(phi)*cos(dec));
        bz = bz + B * (sin(phi)*sin(inc) + cos(phi)*cos(inc)*cos(theta));
    end
end

该程序使用了三个嵌套的函数，分别用于计算长方体、立方体和棱柱的磁场分量。其中，`calc_mag`函数计算长方体在任意位置的磁场分量，`calc_cube`函数计算立方体上每个面元的磁场分量，`calc_prism`函数计算棱柱上每个面元的磁场分量。程序还包括了绘制三维图像的代码，可以直观地显示长方体的磁异常。