时间戳

首页 > 解决方案 > Python上的潮汐力计算

python - Python上的潮汐力计算

问题描述

我目前正在开发一个绘制木星月球系统轨道的程序，部分代码用于计算轨道上的潮汐力。我的目标是能够将不同的潮汐力值绘制为距离的函数。但是，我似乎无法让它工作，即它甚至没有打印出潮汐力值。到目前为止，这是我的代码：

#This function requires user to input mass, semi-major axis, and the eccentricity of the satellite
def JovianSystem(m, A, e, Vi): 
    tmin = 0
    tmax = 0.0005
    dt = 0.00000001 
    M_jup = 1.898e27
    mass_ratio = m/M_jup #Mass ratio of satellite and Jupiter 
    pi = 3.14 
    G = 4*pi**2 
    alpha = 0.001
    d = []

    #Initialize Arrays for time, x and y positions, x and y velocities, and theta 
    t = np.arange(tmin, tmax, dt)
    x = np.zeros(len(t))
    y = np.zeros(len(t))
    Vx = np.zeros(len(t))
    Vy = np.zeros(len(t))
    theta = np.zeros(len(t))
    r = np.zeros(len(t))
    thetap = [0,]

    # Specify initial conditions
    x[0] = A*(1+e)                                                # initialize position of satellite in AU
    y[0] = 0
    Vx[0] = 0                                                    # initialize velocity of satellite in AU/yr
    Vy[0] = Vi* np.pi * np.sqrt((1+e)/(A*(1-e))*(1+mass_ratio))
    theta[0] = 0                                                 # initialize starting angle of satellite in Radians

    # Calculating trajectory
    for i in range(1, len(t), 1):
        theta[i] = np.arctan2(y[i-1],x[i-1])
        r = A*(1-e**2)/(1+e*np.cos(theta[i]))

        Vx[i] = Vx[i-1] - ((G*x[i-1])/(r**3)+((G*x[i-1]*alpha)/r**5))*dt
        Vy[i] = Vy[i-1] - ((G*y[i-1])/(r**3)+((G*y[i-1]*alpha)/r**5))*dt

        x[i] = x[i-1] + Vx[i]*dt
        y[i] = y[i-1] + Vy[i]*dt

        #thetap.append(theta)
        d.append(r)

    for i in range(1, len(t), 1): 
        F_tidal[i] = (2*G*M_jup*m*d)/x[i]
    print(F_tidal)
    return [x,y, F_tidal]

#x, y, thetap = JovianSystem(m,    A,     e,    Vi): 
x, y, F_tidal = JovianSystem(893.2e20, 0.0028, 0.0041, 22) 

a = 0.0028
e = 0.0041
r= sqrt(x**2+y**2)

aphelion = np.max(r)
perihelion = np.min(r)

daphelion = a*(1+e)
dperihelion = a*(1-e)

print('calculated aphelion=', aphelion)
print('calculated perihelion=', perihelion)
print('accepted aphelion=', daphelion)
print('accepted perihelion=', dperihelion)

plt.plot(x,y)
plt.title('Orbit of Io' )
plt.plot(0,0, marker='o', color = 'orange', markersize = 20)
plt.plot(0.0028,0, marker='o', color='y', markersize=10)
plt.xlabel('X (AU)')
plt.ylabel('Y (AU)')
plt.show()

这是我得到的错误：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-38-14530b6d92c9> in <module>
      1 #x, y, thetap = JovianSystem(m,    A,     e,    Vi):
----> 2 x, y, F_tidal = JovianSystem(893.2e20, 0.0028, 0.0041, 22)
      3 
      4 a = 0.0028
      5 e = 0.0041

<ipython-input-37-b88c29847a34> in JovianSystem(m, A, e, Vi)
     43 
     44    for i in range(1, len(t), 1):
---> 45        F_tidal[i] = (2*G*M_jup*m*d)/x[i]
     46    print(F_tidal)
     47    return [x,y, F_tidal]

TypeError: can't multiply sequence by non-int of type 'float'

任何帮助都感激不尽！我只是想了解我能做些什么来解决这个问题并获得潮汐力值。

标签： python

解决方案

您的问题在第 45 行，看起来您有一些列表d，该列表的长度由 确定，t并且您正在尝试针对列表添加多个浮点值(2*G*M_jup*m*d)。这可能不是您想要做的。

例如：

> d = [1 ,2, 3]
[1, 2, 3]
> d * 3
[1, 2, 3, 1, 2, 3, 1, 2, 3]

首先回答 ifF_tidal[i]应该是一个列表、单个浮点值还是什么，对于 and 是一样d的x。根据这一点，您可能希望对dor进行矩阵运算x，在这种情况下，我会指出您使用Numpy，特别是 numpy 数组而不是d变量列表。

另一方面，如果您想从中提取单个数据点d，您只需要使用索引调用来指示它，即d[i]迭代其值（在这种情况下，您的问题只是一种类型或缺少语法）。即你的第 45 行应该是：

F_tidal[i] = (2*G*M_jup*m*d[i])/x[i]

推荐阅读