#Calculate the distance between two satellites that are initially
#aligned with the "sun" M, given M, r1, r2, and the elapsed time t.
#Get the necessary constants and functions from visual, which also imports numpy
from visual import *
#gravitational constant in mks units
G = 6.67300e-11 # m3 kg-1 s-2
#get r1
r1 = input('Enter radius r1 in metres and press return. ')
#get r2
r2 = input('Enter radius r2 in metres and press return. ')
#get M
M = input('Enter mass M in kilogramsand press return. ')
#get t
t = input('Enter elapsed time t in days and press return. ')
#T1 in seconds
T1 = ((4*pi**2)*r1**3/(G*M))**0.5
#T2 in seconds
T2 = ((4*pi**2)*r2**3/(G*M))**0.5
#Convert t in days to t in seconds
t_seconds = t*24.0*60.0*60.0
#Calculate angular displacement
delta_theta = 2*pi*t_seconds*(1/T1-1/T2)
#delta theta could be positive or negative; it doesn't matter for this calculation.
#(Why doesn't it matter?)
delta_r = (r1**2+r2**2-2*r1*r2*cos(delta_theta))**0.5
seconds_in_year = 365.0*24.0*60.0*60.0
#For output purposes, calculate periods in years.
T1_years = T1/seconds_in_year
T2_years = T2/seconds_in_year
#For output purposes, calculate angular separation in degrees.
#Use the "%" operator to make the angle fall between 0 and 360 degrees.
delta_theta_degrees = delta_theta*180.0/pi
delta_theta_degrees = delta_theta_degrees%360
#print out results
print 'You set M = ', M, 'kg, r1 = ', r1, 'm, r2 = ', r2, 'm, t = ',t,' days.'
print 'Orbital period 1 is ', T1, 'years. Orbital period 2 is ', T2, 'years.'
print 'Angular separation is ', delta_theta_degrees, 'degrees.'
print 'The distance between the satellites is ', delta_r,' m.'
sun = sphere(pos=(0,0,0),radius = 0.2, color=color.yellow)
sat1 = sphere(radius = 0.02, color = color.red)
sat2 = sphere(radius = 0.02, color = color.blue)
trail1 = curve(color = color.red)
trail2 = curve(color = color.blue)
dr = curve(color = color.white)
drlabel = label()
theta1 = 0
theta2 = 0
T1_anim = 1.0
T2_anim = T2/T1*T1_anim
t = 0
dt = 0.01
theta1 = 2*pi*t/T1_anim
theta2 = 2*pi*t/T2_anim
r1_anim = 1
r2_anim = r2/r1*r1_anim
while (1==1):
delta_theta = theta2-theta1
delta_r = (r1**2+r2**2-2*r1*r2*cos(delta_theta))**0.5
pos1 = array((r1_anim*cos(theta1), r1_anim*sin(theta1), 0))
pos2 = array((r2_anim*cos(theta2), r2_anim*sin(theta2), 0))
sat1.pos = pos1
sat2.pos = pos2
trail1.append(pos=sat1.pos)
trail2.append(pos=sat2.pos)
#wait before showing the line joining the two
if (t>0.3):
dr.pos = (pos1,pos2)
pos_label = 0.5*(pos1+pos2)
drlabel.pos = pos_label
drlabel.text = str(float(delta_r/1000.0))+'km'
t = t + dt
theta1 = 2*pi*t/T1_anim
theta2 = 2*pi*t/T2_anim
rate(10)