代码收藏家 六自由度外弹道计算程序一步到位版（含气动参数计算）
声明：
本程序仅用于PYTHON编程学习交流用途。
相关计算公式参考《弹箭外弹道学 第二版》；关于六自由度外弹道计算的相关知识，推荐参考《弹箭外弹道学 第二版》《导弹飞行力学》《外弹道学》这几本书。
程序中使用的弹箭的六自由度刚体弹道方程如下所示（参考《弹箭外弹道学 第二版》）：
1. 程序介绍：
程序包含三个数据.csv文件用于输入计算数据并进行中转，此外包含一个SUPPORT文件夹用于存放气动参数计算软件（无界面版DATCOM），文件内容如下图所示：
其中：
001.csv用于存放六月自由度外弹道相关计算基础输入参数，如下图所示：
002.csv用于存放 DATCOM计算出的气动参数结果，如下图所示：
003.csv用于存放 DATCOM气动计算的输入参数，如下图所示：
程序计算过程：
用户根据外弹道计算需求对001.csv和003.csv中的相关参数进行设定，此后运行程序，执行六自由度外弹道计算过程。
计算结束后，程序自动绘制出基本的弹丸飞行轨迹曲线，如下图所示：
并同时在程序目录下输出一个名为result.csv的文件用于存放计算过程中的所有相关变量，以便后续的数据处理及绘图，result.csv文件中的内容如下图所示：
2. 程序代码实现：
2.1 基础内容设置
首先引入程序相关的库：
# 引入必要的库
#..........................................................................................
import math #数学计算库
from matplotlib import pyplot as plt #绘图库
import numpy as np #数值计算库
from scipy.integrate import solve_ivp #ODE求解库
from scipy.interpolate import interp1d, RectBivariateSpline #插值库
import pandas as pd #数据处理库
import chardet #编码检测库
import subprocess #子进程库
import os #操作系统库
import re #正则表达式库
#..........................................................................................并定义用于检测文件编码的函数：
# 定义函数，用于检测文件编码
def detect_encoding(file_path):
try:
with open(file_path, 'rb') as f:
rawdata = f.read()
result = chardet.detect(rawdata)
return result['encoding']
except Exception as e:
print(f"Error detecting encoding: {e}")
return None2.2 相关参数读取及气动参数计算
访问003.csv文件，读取用于进行气动计算的参数，并将其按照所需形式写入到DATCOM的运算输入文件中。
# 指定弹丸气动参数计算DATCOM软件的输入文件路径
csv_file_path_datcom_params = '003.csv'
# 尝试检测弹丸气动参数计算DATCOM软件的输入文件编码
detected_encoding = detect_encoding(csv_file_path_datcom_params)
if detected_encoding is None:
print("Failed to detect encoding. Exiting.")
exit()
else:
try:
# 读取CSV文件
csv_reader = pd.read_csv(csv_file_path_datcom_params, encoding=detected_encoding, header=0)
# 创建空列表来存储数据
first_column_data = []
second_column_data = []
third_column_data = []
# 遍历CSV文件的每一行
for index, row in csv_reader.iterrows():
first_column_data.append(row[0])
second_column_data.append(row[1])
third_column_data.append(row[2])
data = list(zip(first_column_data, second_column_data, third_column_data))
# 打印读取结果
#print(data)
# 赋值给变量
XCG = second_column_data[0]
LNOSE = second_column_data[1]
DNOSE = second_column_data[2]
BNOSE = second_column_data[3]
LCENTR = second_column_data[4]
DCENTR = second_column_data[5]
LAFT = second_column_data[6]
DAFT = second_column_data[7]
# 打印弹丸气动参数
print(f"XCG = {XCG}")
print(f"LNOSE = {LNOSE}")
print(f"DNOSE = {DNOSE}")
print(f"BNOSE = {BNOSE}")
print(f"LCENTR = {LCENTR}")
print(f"DCENTR = {DCENTR}")
print(f"LAFT = {LAFT}")
print(f"DAFT = {DAFT}")
except Exception as e:
print(f"Error processing CSV file 003.csv: {e}")
exit()
运行DATCOM程序并读取相关计算结果转存至002.csv中。
读取外弹道计算基础参数：
# 指定六自由度外弹道计算基础参数文件路径
basic_params_csv_file_path = '001.csv'
detected_encoding = detect_encoding(basic_params_csv_file_path)
if detected_encoding is None:
print("Failed to detect encoding. Exiting.")
else:
try:
# 读取CSV文件
csv_reader = pd.read_csv(basic_params_csv_file_path, encoding=detected_encoding, header=0)
# 创建空列表来存储数据
first_column_data = []
second_column_data = []
third_column_data = []
# 遍历CSV文件的每一行
for index, row in csv_reader.iterrows():
first_column_data.append(row[0])
second_column_data.append(row[1])
third_column_data.append(row[2])
data = list(zip(first_column_data, second_column_data, third_column_data))
# 打印读取结果
print(data)
except Exception as e:
print(f"Error processing CSV file 001.csv: {e}")读取外弹道计算气动参数：
# 指定六自由度外弹道气动基础参数文件路径
aero_params_csv_file_path = '002.csv'
detected_encoding = detect_encoding(aero_params_csv_file_path)
# 读取CSV文件
if detected_encoding is None:
print("Failed to detect encoding. Exiting.")
else:
try:
# 读取CSV文件
df = pd.read_csv(aero_params_csv_file_path, encoding=detected_encoding, header=0)
# 创建一个字典来存储二维数组
arrays = {}
# 遍历DataFrame，提取每个二维数组
current_array_name = None
current_array = []
for index, row in df.iterrows():
array_name = row[0]
if array_name != current_array_name: # 判断二维数组的名字是否改变
if current_array_name is not None: # 如果当前二维数组的名字不为空，则将上一个二维数组添加到字典中
arrays[current_array_name] = current_array # 添加到字典中
current_array_name = array_name # 更新当前二维数组的名字
current_array = [] # 重置当前二维数组
current_array.append(row[1:].tolist()) # 添加到当前二维数组中
# 添加最后一个数组
if current_array_name is not None:
arrays[current_array_name] = current_array
# 打印结果
for array_name, array in arrays.items():
print(f"{array_name}:")
for row in array:
print(row)
print()
except Exception as e:
print(f"Error processing CSV file 002.csv: {e}")
2.3 六自由度外弹道计算
调用ode求解器，使用RK45方法求解微分方程组。
# 调用ode求解器，使用RK45方法求解微分方程组
def wrap_sdofmotion(t, mp, global_params):
# 解包全局参数
X_cm, tau, a, C, A, m, d, l, S, v_0, theta_a_0, psi_2_0, phi_a_0, phi_2_0, gamma_0, omiga_xi_0, omiga_eta_0, omiga_zeta_0, x_0, y_0, z_0, w, a_W, g, rou, beta_D_1, beta_D_2, Lambda, delta_f, delta_N, delta_M, gamma_0_1, gamma_0_2, Omiga_E, a_N, ma1, Delta, CA, CN, X_CP, CNA, CYPA, Cnp1, Cnp3, Cnp5= global_params
return sdofmotion(t, mp, global_params)
# 设置初始条件和参数
t0 = 0.0
mp0 = np.array([v_0,theta_a_0,psi_2_0,omiga_xi_0,omiga_eta_0,omiga_zeta_0,phi_a_0,phi_2_0,gamma_0,x_0,y_0,z_0]) # 设置初始状态向量
t_end = t_estimation # 设置积分终止时间
global_params = [X_cm, tau, a, C, A, m, d, l, S, v_0, theta_a_0, psi_2_0, phi_a_0, phi_2_0, gamma_0,
omiga_xi_0, omiga_eta_0, omiga_zeta_0, x_0, y_0, z_0, w, a_W, g, rou, beta_D_1, beta_D_2,
Lambda, delta_f, delta_N, delta_M, gamma_0_1, gamma_0_2, Omiga_E, a_N,
ma1, Delta, CA, CN, X_CP, CNA, CYPA, Cnp1, Cnp3, Cnp5] # 设置全局参数列表
# 使用solve_ivp进行积分，指定方法为RK45并设置容差
sol = solve_ivp(wrap_sdofmotion, (t0, t_end), mp0, args=(global_params,), method='RK45', rtol=1e-7, atol=1e-7)
2.4 结果处理
将结果输出值result.csv文件，并绘制出基本外弹道飞行轨迹曲线。
# 处理结果
ts = sol.t # 时间序列
mps = sol.y.T # 状态向量序列
#..............................................................................................................................................
# 结果处理及导出
#..............................................................................................................................................
# 导出计算结果至外部csv文件中
df = pd.DataFrame(mps, columns=['v', 'theta_a', 'psi_2', 'omiga_xi', 'omiga_eta', 'omiga_zeta', 'phi_a', 'phi_2', 'gamma', 'x', 'y', 'z'])
df['t'] = ts
df.to_csv('result.csv', index=False)
print('The calculation result has been exported to the external csv file successfully!')
# 绘制结果图
fig, ax = plt.subplots(2, 2, figsize=(12, 10))
ax[0, 0].plot(ts, mps[:, 0], label='v')
ax[0, 0].set_xlabel('t')
ax[0, 0].set_ylabel('v')
ax[0, 0].legend()
ax[0, 0].set_title('Velocity')
ax[0, 1].plot(ts, mps[:, 9], label='X')
ax[0, 1].set_xlabel('t')
ax[0, 1].set_ylabel('X')
ax[0, 1].legend()
ax[0, 1].set_title('Position X')
ax[1, 0].plot(ts, mps[:, 10], label='Y')
ax[1, 0].set_xlabel('t')
ax[1, 0].set_ylabel('Y')
ax[1, 0].legend()
ax[1, 0].set_title('Position Y')
ax[1, 1].plot(ts, mps[:, 11], label='Z')
ax[1, 1].set_xlabel('t')
ax[1, 1].set_ylabel('Z')
ax[1, 1].legend()
ax[1, 1].set_title('Position Z')
fig.suptitle('Time Series Analysis of Motion Data')
plt.tight_layout() # 调整布局以防止标题重叠
plt.show()3. 程序打包
通过pyinstller将程序内容打包为可直接运行的exe文件，用户在设置好相关输入参数后，可通过直接运行exe文件进行外弹道解算。
4. 总结
借助PYTHON，本程序能够方便地将第三方气动参数计算软件（本例中为DATCOM无界面版）整合进程序中，简化了整个计算过程中的工作量（无需再从006.dat文件中找寻相关结果并进行复制粘贴和重新读取）。通过pyinstller将程序打包为可直接运行的exe文件，使得程序无需相关环境准备在不同设备上的流畅运行成为了可能。有关外弹道计算和相关DATCOM气动参数计算软件的资料及手册，以及相关编程过程中遇到的问题和经验总结，欢迎感兴趣的读者和作者沟通交流，qq：3383016937。
作者：曾三金同学
