Scribd
Upload
191 views

Planet Orbits in Python

Planet Orbits in Python

Uploaded by

AlankarDutta
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd
191 views

Planet Orbits in Python

Planet Orbits in Python

Uploaded by

AlankarDutta
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as TXT, PDF, TXT or read online on Scribd
You are on page 1/ 7

"""

Acknowledgement: Converted from the C++ code


Author of the python code: Alankar Dutta & Ranajay Dutta
Pesidency University
Physics Department
Kolkata, West Bengal, India
presiuniv.ac.in
General Description:
====================
This module contains the most up-to-date physical and
astronomical constants in SI units.
"An Introduction to Modern Astrophysics", Appendix I
Bradley W. Carroll and Dale A. Ostlie
Addison Wesley, 2007
Weber State University
Ogden, UT
modastro@weber.edu
------------------------------------------------------------------"""

#The smallest non-zero number and the number of significant figures


tiny_sp
= 3.4e-38
tiny_dp
= 1.7e-308
tiny_qp
= tiny_dp
sig_fig_sp = 7
sig_fig_dp = 15
sig_fig_qp = sig_fig_dp
eps_sp
= 1E-6
eps_dp
= 1E-15
eps_qp
= eps_dp
#The largest number for given precision
biggest_sp = 3.4e38
biggest_dp = 1.7e308
biggest_qp = biggest_dp
biggest_i2 = 32767
biggest_i4 = 2147483647
biggest_i8 = 9223372036854775807
#Values
pi
two_pi
four_pi
four_pi_o3
pi_over_2
natural_e

related to pi and e
= 3.14159265358979323846264338327950
= 2*pi
= 4*pi
= four_pi/3
= pi/2
= 2.71828182845904523536028747135266

#Conversions for radians to degrees and degrees to radians


degrees_to_radians = pi/180
radians_to_degrees = 180/pi

#Physical constants
G
= 6.673e-11
c
= 2.99792458e08
mu_0
= four_pi*1e-07
epsilon_0 = 1/(mu_0*c*c)
e_C
= 1.602176462e-19
eV
= e_C
keV
= eV*1.0e3
MeV
= eV*1.0e6
GeV
= eV*1.0e9;
h
= 6.62606876e-34
hbar
= h/two_pi
k_B
= 1.3806503e-23
sigma
= 2*pow(pi,5)*pow(k_B,4)/(15*c*c*pow(h,3))
a_rad
= 4*sigma/c
a_rad_o3
= a_rad/3
four_ac_o3 = 4*a_rad_o3*c
m_e
= 9.10938188e-31
m_p
= 1.67262158e-27
m_n
= 1.67492716e-27
m_H
= 1.673532499e-27
u
= 1.66053873e-27
N_A
= 6.02214199e23
R_gas
= 8.314472
a_0
= four_pi*epsilon_0*hbar*hbar/(m_e*e_C*e_C)
R_infty
= m_e*pow(e_C,4)/(64*pow(pi,3)*epsilon_0*epsilon_0*pow(hbar,3)*c)
R_H
= m_p/(m_e + m_p)*R_infty
#Time constants
hr
= 3600
day
= 24*hr
J_yr
= 365.25*day
yr
= 3.15581450e7
T_yr
= 3.155692519e7
G_yr
= 3.1556952e7
#Astronomical length constants
= 1.4959787066e11
= 206264.806*AU
= c*J_yr
#Solar constants
M_Sun
= 1.9891e30
S_Sun
= 1.365e3
L_Sun
= four_pi*AU*AU*S_Sun
R_Sun
= 6.95508e8
Te_Sun
= pow((L_Sun/(four_pi*R_Sun*R_Sun*sigma)),0.25)
AU
pc
ly

#Solar magnitudes
Mbol_Sun
= 4.74
MU_Sun
= 5.67
MB_Sun
= 5.47
MV_Sun
= 4.82
Mbol_Sun_ap = -26.83
MU_Sun_ap = -25.91
MB_Sun_ap = -26.10
MV_Sun_ap = -26.75
BC_Sun
= -0.08
#Earth constants
M_Earth
= 5.9736e24

R_Earth

= 6.378136e6

#Unit Conversions
cm
= 1e-2
gram
= 1e-3
erg
= 1e-7
dyne
= 1e-5
esu
= 3.335640952e-10
statvolt
= 2.997924580e2
gauss
= 1e-4
angstrom
= 1e-10
jansky
= 1e-26

"""
Acknowledgement: Converted from the C++ code
Author of the python code: Alankar Dutta & Ranajay Dutta
Pesidency University
Physics Department
Kolkata, West Bengal, India
presiuniv.ac.in

General Description:
====================
Orbit computes the orbit of a small mass about a much larger mass,
or it can be considered as computing the motion of the reduced mass
about the center of mass.
"An Introduction to Modern Astrophysics", Appendix J
Bradley W. Carroll and Dale A. Ostlie
Second Edition, Addison Wesley, 2007
Weber State University
Ogden, UT
modastro@weber.edu
------------------------------------------------------------------"""
import Constants
import math
import matplotlib.pyplot as plt
from prettytable import PrettyTable

class orbit:

#This is the __str__() method to display information on orbit instances


def __str__(self):
message='This program computes the orbit of a small mass about a much la
rger mass,\n or it can be considered as computing the motion of the reduced mass
\n about the center of mass.'
return message

#This is the initializing method


def __init__(self):
"""
From Constants module the following constants will be used:
Constants.G
Constants.AU
Constants.M_Sun
Constants.pi
Constants.two_pi
Constants.yr
Constants.radians_to_degrees
Constants.eps_dp
"""
#Write introductory information for user
print('
Orbit computes the orbit of a small mass about a much larger
mass.\n')
print('
print('
print('
print('
print('

Details of the code are described in: \n')


An Introduction to Modern Astrophysics \n')
Bradley W. Carroll and Dale A. Ostlie \n')
Addison Wesley \n
')
copyright 2007\n\n
')

#Trying to open a new Orbit.txt file or appending to a previous one


#Its just a check to see whether a file can be written on orbital infos
if not the program terminates
try:
orbit_file = open('Orbit.txt', 'a')
orbit_file.write("\n")
except:
try:
orbit_file = open('Orbit.txt', 'w')
orbit_file.write("\n")

except:
print('Unable to open Orbit.txt --- Terminating computation')
self.dummy=input('Enter any character to quit...')
exit()
orbit_file.close()

#Get input from user


def user_input(self):
self.Mstrsun
star (in solar masses): '))
self.aAU
he orbit (in AU):
'))
self.ecc
y:
'))
self.user_input_tuple

= float(input(' Enter the mass of the parent


= float(input(' Enter the semimajor axis of t
= float(input(' Enter the orbital eccentricit
= (self.Mstrsun, self.aAU, self.ecc)

#Convert user entered values to conventional SI units


Mstar
= self.user_input_tuple[0]*Constants.M_Sun
a
= self.user_input_tuple[1]*Constants.AU
self.user_input_tuple_conv = ( Mstar, a, self.ecc)
del self.user_input_tuple
return self.user_input_tuple_conv

#Calculate the orbital period in seconds using Kepler's Third Law


def orbital_period(self,inp_tuple):
P = ((4*(Constants.pi**2)*(inp_tuple[1]**3))/(Constants.G*inp_tuple[0]))
**(0.5)
#Convert the orbital period to years and print the result
print('\n The period of this orbit is ' , (P/Constants.yr), ' yr')
return P
#Enter the number of time steps and the time interval to be printed
def step_time_input(self):
n=int(input('\n\n'+
'Please enter the number of time steps to be calculated and the \n'+
'frequency with which you want time steps printed. \n\n'+
'Note that taking too large a time step during the computation \n'+
'will produce inaccurate results. \n\n'+
'Enter the number of time steps desired for the computation: '))
n+=1 #increment to include t=0 (initial) point
kmax=int(input('\n\n'
'How often do you want time steps to be printed? \n'
'
1 = every time step \n'
'
2 = every second time step \n'
'
3 = every third time step \n'
'
etc. \n\n'
'Frequency: '))
step_tuple = (n, kmax)
return step_tuple
#Print header information for output file

def ouput_file(self,inp_tuple,time):
orbit_file = open('Orbit.txt', 'a')
orbit_file.write('
Elliptical Orbit Data\n\n')
orbit_file.write('
Mstar = ')
orbit_file.write(str(inp_tuple[0]))
orbit_file.write(' Msun\n')
orbit_file.write('
a
= ')
orbit_file.write(str((inp_tuple[1]/Constants.AU)))
orbit_file.write(' AU\n')
orbit_file.write('
revolution time period
= ')
orbit_file.write(str((time/Constants.yr)))
orbit_file.write(' yr\n')
orbit_file.write('
eccentricity
= ')
orbit_file.write(str(inp_tuple[2]))
orbit_file.write('\n\n')

#Start main time step loop


def time_step_loop(self, data_tuple, step_tuple,time_period):
orbit_file = open('Orbit.txt', 'a') #open the file Orbit.txt for printin
g
#Initialize print counter, angle, elapsed time, and time step
counter = 1
#print counter
theta = 0
#angle from direction to per
ihelion (radians)
t
= 0
#elapsed time (s)
dt
= time_period/(step_tuple[0]-1)
#time step (s)
delta = 2*Constants.eps_dp
#allowable error at end of p
eriod
#Start of loop
plotx = []
#collects x co ordinates for plotting
ploty = []
#collects y co ordinates for plotting
r
= 0.0
#Variables inside the loop
LoM
= 0.0
#Variables inside the loop
dtheta = 0.0
#Variables inside the loop
x = PrettyTable(['t (yr)', 'r (AU)', 'theta (deg)', 'x (AU)', 'y (AU)'])
#prepare data Table Headers
x.align['t (yr)'] = 'l'
x.padding_width = 1
while (((theta - Constants.two_pi) < (dtheta/2)) and ((t - time_period)
< dt/2)):
#Calculate the distance from the principal focus using Kepler's Firs
t Law
r = data_tuple[1]*((1 - data_tuple[2]**2)/(1 + (data_tuple[2]*(math.
cos(theta)))))
#If time to print, convert to cartesian coordinates. Be sure to pri
nt last point also
if (counter==1):

x.add_row([(t/Constants.yr),(r/Constants.AU), (theta*Constants.r
adians_to_degrees), ((r*math.cos(theta))/Constants.AU),((r*math.sin(theta))/Cons
tants.AU)])
plotx.append((r*math.cos(theta))/Constants.AU)

ploty.append((r*math.sin(theta))/Constants.AU)

x co ordinates
y co ordinates
orbit_file.write('\n')
#Prepare for the next time step: Update the elapsed time.
t += dt
#Calculate the angular momentum per unit mass, L/m
LoM
= (Constants.G*data_tuple[0]*data_tuple[1]*(1 - data_tuple[
2]**2))**(0.5)
#Compute the next value for theta
dtheta = (LoM/(r**2))*dt
theta += dtheta
#Reset the print counter if necessary
counter += 1
if (counter > step_tuple[1] or (theta - Constants.two_pi)/Constants.
two_pi > delta or (t - time_period)/time_period > delta):
counter = 1
plt.plot(plotx,ploty)
orbit_file.write(str(x))
#plots the Orbit curve
print( '\nOrbit Graph successfully Generated!\n\n' )
print( '\n\nThe computation is finished and the data are in Orbit.txt\n\
n' )
dummy = 'dumb'
dummy = input ( 'Enter any character and hit <enter> to exit: ')

You might also like