正方体重磁异常matlab正演程序
时间: 2023-07-29 10:05:44 浏览: 47
以下是一个简单的正方体重磁异常正演的 Matlab 程序:
```matlab
% 定义正方体的参数
L = 1000; % 正方体边长,单位:米
rho = 3000; % 正方体密度,单位:千克/立方米
sus = 0.01; % 正方体磁化率,无单位
inc = 30; % 正方体倾角,单位:度
dec = 45; % 正方体偏角,单位:度
% 定义计算区域
Nx = 100; % x 方向采样点数
Ny = 100; % y 方向采样点数
Nz = 100; % z 方向采样点数
dx = L/Nx; % x 方向采样间距,单位:米
dy = L/Ny; % y 方向采样间距,单位:米
dz = L/Nz; % z 方向采样间距,单位:米
x = linspace(-L/2,L/2,Nx); % x 方向采样位置,单位:米
y = linspace(-L/2,L/2,Ny); % y 方向采样位置,单位:米
z = linspace(0,L,Nz); % z 方向采样位置,单位:米
% 计算正方体磁场
[Bx,By,Bz] = cube_mag_field(x,y,z,L,rho,sus,inc,dec);
% 计算正方体重磁异常
[Bxg,Byg,Bzg] = cube_grav_field(x,y,z,L,rho,sus);
dGz = gravity_anomaly(x,y,z,Bxg,Byg,Bzg);
% 绘制重磁异常剖面图
figure;
imagesc(x/1000,y/1000,dGz); % 将 x,y 轴单位转化为千米
xlabel('x (km)');
ylabel('y (km)');
colorbar;
```
其中,`cube_mag_field` 和 `cube_grav_field` 分别是计算正方体磁场和重力场的函数。以下是这两个函数的代码:
```matlab
function [Bx,By,Bz] = cube_mag_field(x,y,z,L,rho,sus,inc,dec)
% 计算正方体磁场
% 输入参数:
% x,y,z: 采样位置,单位:米
% L: 正方体边长,单位:米
% rho: 正方体密度,单位:千克/立方米
% sus: 正方体磁化率,无单位
% inc: 正方体倾角,单位:度
% dec: 正方体偏角,单位:度
% 输出参数:
% Bx,By,Bz: 正方体磁场分量,单位:纳特
% 计算正方体中心位置坐标
xc = mean(x);
yc = mean(y);
zc = mean(z);
% 计算正方体顶点坐标
xv = xc + [-L/2,-L/2,-L/2,-L/2,L/2,L/2,L/2,L/2];
yv = yc + [-L/2,-L/2,L/2,L/2,-L/2,-L/2,L/2,L/2];
zv = zc + [-L/2,L/2,-L/2,L/2,-L/2,L/2,-L/2,L/2];
% 计算正方体中心到采样点的距离和方向
r = sqrt((x-xc).^2 + (y-yc).^2 + (z-zc).^2);
theta = acos((z-zc)./r);
phi = atan2(y-yc,x-xc);
% 计算正方体磁场
[Bx,By,Bz] = magnet(x,y,z,xv,yv,zv,rho,sus,inc,dec,theta,phi);
end
function [Bx,By,Bz] = cube_grav_field(x,y,z,L,rho,sus)
% 计算正方体重力场
% 输入参数:
% x,y,z: 采样位置,单位:米
% L: 正方体边长,单位:米
% rho: 正方体密度,单位:千克/立方米
% sus: 正方体磁化率,无单位
% 输出参数:
% Bx,By,Bz: 正方体重力场分量,单位:mGal
% 计算正方体中心位置坐标
xc = mean(x);
yc = mean(y);
zc = mean(z);
% 计算正方体顶点坐标
xv = xc + [-L/2,-L/2,-L/2,-L/2,L/2,L/2,L/2,L/2];
yv = yc + [-L/2,-L/2,L/2,L/2,-L/2,-L/2,L/2,L/2];
zv = zc + [-L/2,L/2,-L/2,L/2,-L/2,L/2,-L/2,L/2];
% 计算正方体重力场
[Bx,By,Bz] = gravity(x,y,z,xv,yv,zv,rho,sus);
end
```
其中,`magnet` 和 `gravity` 函数是计算磁场和重力场的核心函数,这里不再赘述,可以参考其他资料。