磁异常的分析MATLAB

写出球体的磁异常，分析典型设置下磁异常的变化特征，注重分析磁倾角变化、磁化率变化、球体埋深以及球体半径变化的磁场的特征

我可以给出一些参考代码，但需要先说明的是，球体的磁异常特征与球体本身的磁化率、半径、埋深等因素有关，因此代码需要先定义这些参数，再进行模拟。

参考代码：

假设球体半径为r，埋深为depth，球体磁化率为k，地磁场磁倾角为delta，水平方向磁场为bx，竖直方向磁场为bz，我们可以先定义这些参数：

r = 1; %球体半径 depth = 100; %球体埋深 k = 0.1; %球体磁化率 delta = 30; %磁倾角为30度 bx = 20000; %水平方向磁场为20000nT bz = 50000; %竖直方向磁场为50000nT

接下来，我们可以编写函数来计算球体在某一点产生的磁场：

function [bxp,bzp] = calcMagField(xp,zp) re = 6371.2; % 地球半径 mu = 4pi1e-7; % 真空磁导率 inc = 0; % 磁倾角为0 dec = 90; % 极轴方向 rp = sqrt(xp.^2+zp.^2); latp = asin(zp./rp); lonp = atan2(xp,zp); [bxh,bzh] = wrldmagm([latp180/pi,lonp180/pi,rp/1000,dec,inc],'2015'); bxp = bxhcos(deltapi/180) - bzhsin(deltapi/180); bzp = bxhsin(deltapi/180) + bzhcos(deltapi/180); bxp = bxp + (mu/4/3.14)(3k(rp^3)/((rp^2+r^2-2rprcos(latp))^2))... (bzp+(r^3/(2rp^3))((3(xp/rp)bxp-3(zp/rp)bzp))); bzp = bzp + (mu/4/3.14)(3k(rp^3)/((rp^2+r^2-2rprcos(latp))^2))... (bxp+(r^3/(2rp^3))((3(zp/rp)bzp-3(xp/rp)bxp))); end

其中，我们调用了Matlab内置的函数wrldmagm来计算地球磁场，然后再加上球体产生的磁场。值得注意的是，这里计算磁场的单位是nT。

最后，我们可以定义一组坐标点，计算球体在这些点上的磁场：

x = -500:50:500; z = -500:50:500; [X,Z] = meshgrid(x,z); Bxp = zeros(size(X)); Bzp = zeros(size(Z)); for i=1:numel(X) [Bxp(i),Bzp(i)] = calcMagField(X(i),Z(i)+depth); end

这里我们用meshgrid生成了一组坐标点，然后通过循环计算每个点上的磁场。

最后，我们可以用Matlab自带的subplot和contourf函数绘制磁场的等值线图：

figure; subplot(121); contourf(X,Z,Bzp,30); colorbar; xlabel('x (m)'); ylabel('z (m)'); title('Vertical Component (nT)'); subplot(122); contourf(X,Z,Bxp,30); colorbar; xlabel('x (m)'); ylabel('z (m)'); title('Horizontal Component (nT)');

运行结果：