辐射实验数据处理一：通过有限点绘制曲面图

首先吐槽一下这种为了赶进度而在尚未学习核物理与辐射安全前就做实验的行为；
在实验六：环境伽马剂量监测实验中，要求绘制指定放射源的剂量分布图谱，但受制于实验场地，只能测量横纵两个相互垂直方向共十组数据。但实验同时要求我们使用软件绘制剂量分布的三维图像。显然已由的这些垂直的点放到Matlab里出来的是简单的平面堆叠。为此，需要找到一种方法绘制出优美的曲面出来。
当然，这和学习无关，和得分无关。纯粹是个人爱好。
首先，检查python和相关库，这里用的是python3.8 + matplotlib + scipy + numpy，都是很基础的库。不会pip方法的可以用Anaconda3一键安装。
接下来，需要对原始数据进行处理。

这些点在空间中的样子是这样的：

怎么画出来？代码如下：

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
#scatter就是描点命令，s是点的大小，c是点的颜色。需要复现的无脑改数据就完事了。
ax.scatter(65, 0, 43.2, s=20, c='r', depthshade=True)
ax.scatter(130, 0, 29.8, s=20, c='r', depthshade=True)
ax.scatter(195, 0, 24.05, s=20, c='r', depthshade=True)
ax.scatter(260, 0, 21.25, s=20, c='r', depthshade=True)
ax.scatter(325, 0, 20.75, s=20, c='r', depthshade=True)
ax.scatter( 0, 65,43.15, s=20, c='b', depthshade=True)
ax.scatter( 0, 130,30.5, s=20, c='b', depthshade=True)
ax.scatter( 0, 195,23.7, s=20, c='b', depthshade=True)
ax.scatter( 0,260, 21.6, s=20, c='b', depthshade=True)
ax.scatter( 0, 325,20.1, s=20, c='b', depthshade=True)

ax.set_xlabel('Parallel door direction')
ax.set_ylabel('Vertical door direction')
ax.set_zlabel('Radiation Dosage 10E-8Gy/h')

plt.show()

核心思路是：
辐射场的分布在空间内是均匀的，因此在xz和yz平面上的曲线符合该场在任意过z轴截面的曲线

沿着这个思路，我们需要找到一种方法，寻找出符合该平面上辐射随距离衰减的方式。通过阅读课本和实验手册，我们知道了这个衰减是与距离的平方呈反比的。为此，我打算直接采用最小二乘法拟合出过这些点的曲线。为了令这种拟合变得更加符合实验，我特定选取了如下拟合函数：
y=a1/(a2*(x*x)+a3*x+a4) +a5
由于换的MBP只有可怜的256g，为此不想在下载Oringin或者Matlab这类软件出拟合曲线了。Python也可以实现，毕竟他们都是图灵完备的。具体实现的方法参照如下代码：

import matplotlib.pyplot as plt
from scipy.optimize import leastsq
import  numpy as np


def Fun(p, x): # 定义拟合函数形式，这这里定义需要拟合的函数形式
    a1,a2,a3,a4,a5 = p #注意p是一个五维列表，最好用np.p的形式，会避免怪错
    y=a1/(a2*(x*x)+a3*x+a4) +a5
    return y

def error(p, x, y):  # 拟合残差
    return Fun(p, x) - y

def main():
    x = np.array([65.,130.,198.,260.,325.,1000.]) #输入你的实验数据，注意最后的1000是为了接近本底辐射
    y = np.array([43.2,30.,23.9,21.4,20.4,19.45]) #输入你的实验测得的辐射值
    p0 = [0.1,0.1,0.1,1,10]  # 拟合的初始参数设置
    para = leastsq(error, p0, args=(x, y))  # 进行拟合
    y_fitted = Fun(para[0], x)  # 画出拟合后的曲线

    plt.figure
    plt.plot(x, y, 'r', label='Original curve')
    plt.plot(x, y_fitted, '-b', label='Fitted curve')
    plt.legend()
    plt.show()
    print(para[0])

if __name__ == '__main__':
    main()

点击运行，不出意外的话会产生如下曲线：

是不是很接近？但需要注意的是这个曲线不是非常光滑，因为我懒得写了。
同样，也会给出拟合函数的五个未知参数大小： [-1.50903395e+09 -1.30671680e+04 1.38404720e+06 -9.74987756e+07 1.91459819e+01]这样，我就得到了一个漂亮的曲线： Z = (-1.50903395e+09) / ((-1.30671680e+04) * (X * X) + (1.38404720e+06) * X + (-9.74987756e+07)) + (1.91459819e+01);他是长成这样的：

这也就是我们要求的距离防护曲线。
最后一步就是绘制表面图了。这很简单，觉得难是因为大一的线性代数第二章空间几何没有好好听好好学。已知X-Z图这么绘制绕Z轴旋转的图？将X替换为根号下X的平方+Y的平方！（我讨厌web，我写不了公式）就像这个： Z = (-1.50903395e+09) / ((-1.30671680e+04) * (X * X+Y*Y) + (1.38404720e+06) * np.sqrt(X * X+Y*Y) + (-9.74987756e+07)) + (1.91459819e+01)
代码在这里：

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np

fig = plt.figure()
ax = plt.axes(projection='3d')
#这一片是把一开始的点描出来明白不？
ax.scatter(65, 0, 43.2, s=20, c='b', depthshade=True)
ax.scatter(130, 0, 29.8, s=20, c='b', depthshade=True)
ax.scatter(195, 0, 24.05, s=20, c='b', depthshade=True)
ax.scatter(260, 0, 21.25, s=20, c='b', depthshade=True)
ax.scatter(325, 0, 20.75, s=20, c='b', depthshade=True)
ax.scatter( 0, 65,43.15, s=20, c='y', depthshade=True)
ax.scatter( 0, 130,30.5, s=20, c='y', depthshade=True)
ax.scatter( 0, 195,23.7, s=20, c='y', depthshade=True)
ax.scatter( 0,260, 21.6, s=20, c='y', depthshade=True)
ax.scatter( 0, 325,20.1, s=20, c='y', depthshade=True)
# Make data.这里是XY的取值范围和步长
X = np.arange(0, 350, 0.25)
Y = np.arange(0, 350, 0.25)

X, Y = np.meshgrid(X, Y) #这里你Matlab肯定见过
Z = (-1.50903395e+09) / ((-1.30671680e+04) * (X * X+Y*Y) + (1.38404720e+06) * np.sqrt(X * X+Y*Y) + (-9.74987756e+07)) + (1.91459819e+01)#这是之前的函数


# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

# Customize the z axis.
ax.set_ylim(0, 350)#注意跳针范围
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_xlabel('Parallel door direction')#可以添加标签
ax.set_ylabel('Vertical door direction')
ax.set_zlabel('Radiation Dosage 10E-8Gy/h')
plt.show()

点击运行他应该不会报错，这样以来，就会得到这样的美丽图像：

美丽漂亮，还有被隐藏起来的瑕疵。