2

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!

Improve this question
5
  • 1
    ints in python are only limited by your memory Commented Mar 15, 2015 at 13:01
  • constants with a . are float Commented Mar 15, 2015 at 13:02
  • lose your fortran thinking to code in python Commented Mar 15, 2015 at 13:04
  • Well, from looking at your code, I don't see why, but for some reason your a(x,y) function is receiving a value of nan for either x or y. Commented Mar 15, 2015 at 13:11
  • And yes, you probably should define your arrays as floats. If you're using tile, it's as simple as specifying 0. instead of 0 as the repeated value: x = np.tile(0., N) Commented Mar 15, 2015 at 13:14

2 Answers 2

Reset to default
4

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)
Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

how to create it with float32??
@Zibri answered
You may run into problems if you do not have numpy, nor pip: w3schools.com/python/numpy/numpy_getting_started.asp Alternatives?
3

You need x = np.zeros(N), etc.: this declares the arrays as float arrays.

This is the standard way of putting zeros in an array (np.tile() is convenient for creating a tiling with a fixed array).

Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.