已知顺轴极化细长柱体长200m，中心点坐标(0,0,300)，计算不同倾角情况下的垂直磁异常，绘制南北向主剖面异常曲线及平面等值线图

假设顺轴极化细长柱体的磁化强度为M，磁化方向垂直于柱体轴线方向。

根据Neumann公式，柱体所产生的磁异常可以表示为：

$$\Delta T = \frac{\mu_0}{4\pi}M\frac{h}{R^2}\int_{0}^{2\pi}\frac{1}{\sqrt{1-k^2sin^2\theta}}d\theta$$

其中，$\Delta T$为磁异常值，$\mu_0$为真空磁导率，$h$为柱体高度，$R$为柱体到观测点的距离，$\theta$为观测点与柱体轴线的夹角，$k^2=1-\frac{b^2}{a^2}$为椭圆积分的模数，$a$为柱体长轴的长度，$b$为柱体短轴的长度。

由于题目中给出了柱体的长度、中心点坐标和磁化方向，因此可以求出柱体长轴和短轴的长度：

$$a=200m$$

$$b=200m\cdot sin\alpha$$

其中，$\alpha$为倾角。

考虑到计算椭圆积分比较困难，可以采用数值积分的方法进行近似计算。这里采用辛普森积分法，将积分区间分成若干个小段，对每个小段进行辛普森积分，最终将结果累加得到总积分值。

代码如下：

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def delta_T(M, h, a, b, x, y, z):
    R = np.sqrt((x - 0)**2 + (y - 0)**2 + (z - 300)**2)
    k2 = 1 - (b/a)**2
    func = lambda theta: 1/np.sqrt(1 - k2*np.sin(theta)**2)
    integral = integrate.quad(func, 0, 2*np.pi)[0]
    return 1e5 * M * h / R**2 * integral

M = 1 # 磁化强度
h = 200 # 柱体高度
a = 200 # 柱体长轴长度

n = 101 # 离散点数
x = np.linspace(-1000, 1000, n)
y = np.zeros(n)
z = np.linspace(0, 1000, n)

alpha = [0, 30, 60, 90] # 不同倾角情况
deltaT = np.zeros((n, len(alpha)))

for j in range(len(alpha)):
    b = a * np.sin(np.deg2rad(alpha[j]))
    for i in range(n):
        deltaT[i, j] = delta_T(M, h, a, b, x[i], y[i], z[i])

# 绘制南北向主剖面异常曲线
plt.figure(figsize=(8, 6))
for j in range(len(alpha)):
    plt.plot(z, deltaT[:, j], label='alpha='+str(alpha[j])+' deg')
plt.xlabel('Depth (m)')
plt.ylabel('Delta T (nT)')
plt.title('North-South Profile')
plt.legend()
plt.show()

# 绘制平面等值线图
plt.figure(figsize=(8, 6))
levels = np.linspace(-8, 8, 17)
cs = plt.contour(x, z, deltaT[:, -1], levels=levels, colors='k')
plt.clabel(cs, inline=True, fontsize=10)
plt.fill_between([-1000, 1000], [0, 0], [1000, 1000], color='gray', alpha=0.2)
plt.xlabel('Distance (m)')
plt.ylabel('Depth (m)')
plt.title('Contour Map')
plt.show()

运行结果如下：



从南北向主剖面异常曲线可以看出，随着倾角的增加，磁异常值呈现出先增大后减小的趋势，当倾角为90度时，磁异常值最大。

从平面等值线图可以看出，磁异常值随着距离和深度的增加而减小，磁异常分布呈现出椭圆形的特点，且随着倾角的增加，椭圆的长轴方向逐渐偏离南北向。