追加一个部分。说下原因，本人参加了一个数学建模大赛，其中用了本文部分内容作为模型。写那个论文的时候需要画图，买不起matlab，自己用python写了一个，就是画欧拉潘索情形的。自我感觉还不错，就发上来分享一下。偷懒就没写图形化窗口了，以后有时间可能会写，或者把它包装成库什么的。其他就不多说了，上源码。

import numpy as np
from matplotlib import pyplot as plt
from scipy import special,integrate
import random
#调整字体
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

def omega_t(t):
    #这里序号不一样！是定义的原因
    w1 = A[0]*special.ellipj(a*t+b,k*k)[1]
    w2 = A[1]*special.ellipj(a*t+b,k*k)[0]
    w3 = A[2]*special.ellipj(a*t+b,k*k)[2]
    return w1,w2,w3
def e3_ta(t):
    return np.sqrt((I[0]**2)*(omega_t(t)[0]**2)+(I[1]**2)*(omega_t(t)[1]**2))/((I[1])*(omega_t(t)[1]))
def e3_tb(t):
    return (I[0])*(omega_t(t)[0])/np.sqrt((I[0]**2)*(omega_t(t)[0]**2)+(I[1]**2)*(omega_t(t)[1]**2))
def e3_t(t):
    df_dt = (e3_tb(t+lim[3]/2)-e3_tb(t-lim[3]/2))/lim[3]
    return e3_ta(t)*df_dt
def e1_t(t):
    return L*(omega_t(t)[2]-e3_t(t))/(I[2]*omega_t(t)[2])
def eular_t(t):
    e2 = np.arccos(I[2]*omega_t(t)[2]/L)
    e3 = integrate.quad(e3_t,a=0,b=t)[0]+eular0[2]
    e1 = integrate.quad(e1_t,a=0,b=t)[0]+eular0[0]
    return e1,e2,e3

print('欢迎使用刚体自由定点转动（欧拉潘索情形）绘图器\n制作者为辰星引力航天工作室\n本绘图器将得出三主轴角速度及以角动量方向为极轴的欧拉角\n本绘图器得出欧拉角不含使得欧拉角处在定义域中的常数')
#初始化上、下限、作图步长、精确度、微扰
lim = [0,0,0,0,0]
while True:
    c = input('请输入时间上下限，作图步长，积分求导精确度，微扰。格式：t0,t,d1,d2,d3')
    d = c.split(',')
    if not len(d) == 5:
        print('请输入正确的格式')
        continue
    try:
        for i in range(5):
            lim[i] = float(d[i])
    except:
        print('请输入正确的格式')
        continue
    if lim[2] <= 0:
        print('作图步长必须大于零')
        continue
    if lim[3] <= 0:
        print('积分求导精确度必须大于零')
        continue
    if lim[4] == 0:
        print('微扰不可等于零')
        continue
    break

#初始化主轴惯量
I = [0,0,0]
#初始化三轴角速度
omega0 = [0,0,0]
#初始化欧拉角
eular0 = [0,0,0]

while True:
    c = input('请输入三个主轴惯量，和对应的初始角速度，初始欧拉角(章动角θ与初始角度无关，只是为了形式整齐，随便填个数就行)。格式：I1,I2,I3;O1,O2,O3;φ,θ,ψ')
    d = c.split(';')
    try:
        for i in range(len(d[0].split(','))):
            I[i]=float((d[0].split(','))[i])
        for i in range(len(d[1].split(','))):
            omega0[i]=float((d[1].split(','))[i])
        for i in range(len(d[2].split(','))):
            eular0[i]=float((d[2].split(','))[i])
    except:
        print('请输入正确的格式')
        continue
    if len(I)==3 and len(omega0) == 3 and len(eular0) == 3:
        for i in I:
            if i <= 0:
                print('主轴惯量应大于零')
                continue
        break
    print('请输入正确的格式')
#排序
if I[0]>I[1]:
    tr1 = I[0]
    I[0] = I[1]
    I[1] = tr1
    tr2 = omega0[0]
    omega0[0] = omega0[1]
    omega0[1] = tr2
if I[1]>I[2]:
    tr1 = I[1]
    I[1] = I[2]
    I[2] = tr1
    tr2 = omega0[1]
    omega0[1] = omega0[2]
    omega0[2] = tr2
if I[0]>I[1]:
    tr1 = I[0]
    I[0] = I[1]
    I[1] = tr1
    tr2 = omega0[0]
    omega0[0] = omega0[1]
    omega0[1] = tr2
#不想讨论了，直接加个小量防止出现除零错误
if I[1] == I[2]:
    I[2]+=lim[3]
#这里加微扰是为了模拟真实情况
for i in range(3):
    omega0[i] += lim[4]
#初始化动能和角动量平方
T = 0
L2 = 0
for i in range(3):
    T+=0.5*I[i]*(omega0[i]**2)
    L2+=(I[i]**2)*(omega0[i]**2)
L = np.sqrt(L2)

if (I[2]-I[1])*(L2-2*I[0]*T)<(I[1]-I[0])*(2*I[2]*T-L2):
    tr1 = I[2]
    I[2] = I[0]
    I[0] = tr1
    tr2 = omega0[2]
    omega0[2] = omega0[0]
    omega0[0] = tr2


#函数形如Ae(a*t+b,k)/e()为某个椭圆函数
k = np.sqrt(((I[1]-I[0])*(2*I[2]*T-L2))/((I[2]-I[1])*(L2-2*I[0]*T)))
a = np.sign(omega0[2])*np.sign(I[2]-I[0])*np.sqrt((I[2]-I[1])*(L2-2*I[0]*T)/(I[0]*I[1]*I[2]))
b = special.ellipkinc(np.arcsin(np.sqrt(I[1]*(I[2]-I[1])/(2*I[2]*T-L2))*omega0[1]),k*k)

A = []
A.append(np.sqrt((2*I[2]*T-L2)/(I[0]*(I[2]-I[0]))))
A.append(np.sqrt((2*I[2]*T-L2)/(I[1]*(I[2]-I[1]))))
A.append(np.sign(omega0[2])*np.sqrt((L2-2*I[0]*T)/(I[2]*(I[2]-I[0]))))

time = np.arange(lim[0],lim[1],lim[2])

omega = [[],[],[]]
omega = omega_t(time)

eular = [[],[],[]]
for i in time:
    for j in range(3):
        eular[j].append(eular_t(i)[j])
        
ti1='一般刚体自由定点转动三轴角速度随变化时间关系'
ti2='一般刚体自由定点转动欧拉角随时间变化关系'
#初始化画布对象，并调整窗口
fig = plt.figure(figsize=(12, 5), dpi=100)
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
#绘制三轴角速度随时间变化
ax1.set_title(ti1)
ax1.plot(time,omega[0])
ax1.plot(time,omega[1])
ax1.plot(time,omega[2])
ax1.set_xlabel('t')
ax1.set_ylabel('三轴角速度')
ax1.legend(['1轴角速度','2轴角速度','3轴角速度'])
ax1.axhline(0)
ax1.axvline(0)

#绘制欧拉角随时间变化
ax2.set_title(ti2)
ax2.plot(time,eular[0])
ax2.plot(time,eular[1])
ax2.plot(time,eular[2])
ax2.set_xlabel('t')
ax2.set_ylabel('欧拉角')
ax2.legend(['进动角','章动角','自转角'])
ax2.axhline(0)
ax2.axvline(0)

fig.savefig('刚体自由定点转动')
plt.show()

本程序为了体现变化情况，就舍去了欧拉角的定义域，还有，为了不讨论那么多分母会趋于零的极限，索性加个小量，不过呢，这样反而更好地看出微扰，而且由于直接加小量，章动角不可能严格为零，所以章动角很小的时候不会默认把进动算作自转。

放个效果图，记得之前说的网球拍定理吗？

这是只有1轴方向有角速度，其他方向是微扰（为了明显起见，微扰设得很大）的情况：

这是只有2轴方向有角速度，其他方向是微扰（为了明显起见，微扰设得很大）的情况：

这可以验证网球拍定理。

对了，说下，这个自用的，没考虑什么代码的效率，算起来很慢 几分钟才能出一张。

如果有发现程序出问题的，可以跟我说，不过我最近有些事，可能没法很快回……不过看看是不是库少了之类的先。