Search code examples
pythonarraysnumpyphysicsastronomy

Is there a way to define a float array in Python?


For my astronomy homework, I need to simulate the elliptical orbit of a planet around a sun. To do this, I need to use a for loop to repeatedly calculate the motion of the planet. However, every time I try to run the program, I get the following error:

RuntimeWarning: invalid value encountered in power
r=(x**2+y**2)**1.5
Traceback (most recent call last):
File "planetenstelsel3-4.py", line 25, in <module>
ax[i] = a(x[i],y[i])*x[i]
ValueError: cannot convert float NaN to integer

I've done some testing, and I think the problem lies in the fact that the values that are calculated are greater than what fits in an integer, and the arrays are defined as int arrays. So if there was a way do define them as float arrays, maybe it would work. Here is my code:

import numpy as np
import matplotlib.pyplot as plt

dt = 3600 #s
N = 5000
x = np.tile(0, N)
y = np.tile(0, N)
x[0] = 1.496e11 #m
y[0] = 0.0
vx = np.tile(0, N)
vy = np.tile(0, N)
vx[0] = 0.0
vy[0] = 28000 #m/s
ax = np.tile(0, N)
ay = np.tile(0, N)
m1 = 1.988e30 #kg
G = 6.67e-11 #Nm^2kg^-2

def a(x,y):
    r=(x**2+y**2)**1.5
    return (-G*m1)/r

for i in range (0,N):
    r = x[i],y[i]
    ax[i] = a(x[i],y[i])*x[i]
    ay[i] = a(x[i],y[i])*y[i]
    vx[i+1] = vx[i] + ax[i]*dt
    vy[i+1] = vy[i] + ay[i]*dt
    x[i+1] = x[i] + vx[i]*dt 
    y[i+1] = y[i] + vy[i]*dt    
plt.plot(x,y)
plt.show()

The first few lines are just some starting parameters.

Thanks for the help in advance!


Solution

  • When you are doing physics simulations you should definitely use floats for everything. 0 is an integer constant in Python, and thus np.tile creates integer arrays; use 0.0 as the argument to np.tile to do floating point arrays; or preferably use the np.zeros(N) instead:

    You can check the datatype of any array variable from its dtype member:

    >>> np.tile(0, 10).dtype
    dtype('int64')
    >>> np.tile(0.0, 10).dtype
    dtype('float64')
    >>> np.zeros(10)
    array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])
    >>> np.zeros(10).dtype
    dtype('float64')
    

    To get a zeroed array of float32 you'd need to give a float32 as the argument:

    >>> np.tile(np.float32(0), 10)
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)
    

    or, preferably, use zeros with a defined dtype:

    >>> np.zeros(10, dtype='float32')
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)