In [1]:
%matplotlib inline
import pymc as pm
import matplotlib.pyplot as plt
import numpy as np
In [2]:
# Define key parameters, Values from Zahnle 1998
vp = 13.070 # Orbital velocity of Jupiter, [km/s]
R_j = 71398.0 # jovian radii [km]
G = 6.67*10**(-11) # Universal Gravitational Constant [m3 kg−1 s−2]
moon = 'Ganymede'
v_orb = {'Io' : 17.3,
'Europa' : 13.7,
'Ganymede' : 10.9,
'Callisto' : 8.2} # Orbital velocity of Jupiter's Moons, [km/s]
R_sat = {'Io' : 1820,
'Europe' : 1570,
'Ganymede' : 2630,
'Callisto' : 2400} # Radius [km]
a_sat = {'Io' : 5.9*R_j,
'Europa' : 9.4*R_j,
'Ganymede' : 15.0*R_j,
'Callisto' : 26.4*R_j} # Semimajor axis of orbit [km]
g = {'Io' : 1.81,
'Europa' : 1.31,
'Ganymede' : 1.43,
'Callisto' : 1.25} # Surface Gravity [m/s2]
def v_esc(key,g=g,r=R_sat):
return np.sqrt(2*g[key]*R_sat[key]*1000)/1000 # Calculates surface velocity from the python dictionaries [km/s]
In [24]:
# Calculate distribution of encounter velocities
x1 = pm.Uniform('x1', lower=0.0, upper=1.0) # Uniformly distributed random number
def Vinf_eval(vp=vp,x1=x1):
return 0.3*vp*(1-0.45* np.log(1/(x1+0.2)-0.832))
def Uinf_eval(vp=vp,x1=x1,v_orb=v_orb[moon]):
return Vinf_eval(vp,x1)/v_orb
@pm.deterministic
def Vinf(vp=vp,x1=x1):
return Vinf_eval(vp,x1)
@pm.deterministic
def Uinf(vp=vp,x1=x1,v_orb=v_orb[moon]):
return Uinf_eval(vp,x1,v_orb)
In [25]:
#Create Model and
M = pm.Model([x1, Vinf])
mcmc = pm.MCMC(M)
In [26]:
#Generate Distribution
mcmc.sample(40000, 10000, 1)
Vinf_samples = mcmc.trace('Vinf')[:]
# Plot distribution of encounter velocities
nbins = 50
plt.figure(figsize=(12,6))
ax1 = plt.subplot(121)
ax1.hist(Vinf_samples,bins=nbins,normed=True)
# Plot culmulative distribution of encounter velocities
ax2 = plt.subplot(122)
ax2.hist(Vinf_samples,bins=nbins,cumulative=True,normed=True)
plt.show()
print "Min V_inf : %f" % Vinf_samples.min()
print "Max V_inf : %f" % Vinf_samples.max()
print "Median V_inf : %f" % np.median(Vinf_samples)
In [27]:
#Calculate Distribution of Inclinations
x2 = pm.Uniform('x2', lower=0.0, upper=1.0) # Uniformly distributed random number
def inclination(x2=x2):
return np.arccos(1-2*x2)
In [28]:
# Plot culmulative distribution of inclinations
# Kluge Fix to Populate a list for now so we can check as we go
inclination_samples = [inclination(x2.random()) for i in range(10000)]
values, base = np.histogram(inclination_samples,bins=50)
cumulative = np.cumsum(values)
plt.figure(figsize=(6,6))
ax1 = plt.subplot(111)
plt.plot(base[:-1]*180/np.pi, cumulative, '.')
plt.show()
In [29]:
x3 = pm.Uniform('x3', lower=0.0, upper=1.0) # Uniformly distributed random number
#WARNING, THERE IS AN ERROR IN ZAHNLE 2001
def eccentricity(Uinf=Uinf,x3=x3):
UU_inf = (Uinf)**2
return np.sqrt(1+x3*UU_inf*(UU_inf+2))
In [31]:
ecc_samples = [eccentricity(Uinf_eval(vp,x1.random()),x3.random()) for i in range(10000)]
In [32]:
plt.figure(figsize=(6,6))
ax1 = plt.subplot(111)
ax1 = plt.hist(ecc_samples,nbins)
In [49]:
def periapse(Uinf=Uinf,eccen='Nan'):
UU_inf = Uinf**2
return (eccen-1)/UU_inf
In [50]:
periapse_samples = [periapse(Uinf_eval(vp,x1.random()),eccentricity(Uinf_eval(vp,x1.random()),x3.random())) for i in range(10000)]
In [52]:
plt.figure(figsize=(6,6))
ax1 = plt.subplot(111)
ax1 = plt.hist(periapse_samples,nbins)
Impact probabilities and impact velocities are estimated using Opik’s formul as written for hyperbolic orbits.
We assume that the satellite’s orbit is circular. The normalized encounter velocity U ≡ venc/vorb is given as:
In [54]:
def U(e,q,i):
return np.sqrt(3-(1-e)/q-2*np.cos(i)*np.sqrt(q*(1+e)))
U_samples = [U() for i in range(10000)]
In [55]:
def Ux(e,q):
return np.sqrt(2-(1-e)/q-q(1+e))
In [56]:
def Uy(e,q,i):
return np.cos(i)*np.sqrt(q*(1+e))-1
In [57]:
def Uz(e,q,i):
return np.sin(i)*np.sqrt(q*(1-e))
In [58]:
def imact_probabilities(i,U,Ux,U_inf,R_sat,a_sat,v_orb,v_esc,R_j):
Pstar = ( R_sat/a_sat )**2 *( 1+(v_esc/(v_orb*U))**2 )*(U/Ux)/(np.pi*np.sin(i))
w = (2+ U_inf**2)/(2 + (R_j/a_sat)*U_inf**2 )
return w*Pstar
In [59]:
x4 = pm.Uniform('x4', lower=0.0, upper=1.0) # Uniformly distributed random number
def azimuth_eval(x4=x4):
return np.pi*(1-2*x4)
In [60]:
x5 = pm.Uniform('x5', lower=0.0, upper=1.0) # Uniformly distributed random number
def theta(x5=x5):
return np.arccos(np.sqrt(x5))
In [60]:
In [20]:
@pm.deterministic
def impact_event(x1=x1,x2=x2,x3=x3,x4=x4,x5=x5):
Uinf = Uinf_eval(x1)
return Uinf
In [21]:
mtest = pm.Model([x1,x2,x3,x4,x5,impact_event])
In [22]:
mcmc = pm.MCMC(mtest)
In [23]:
#Generate Distribution
mcmc.sample(40000, 10000, 1)
Vinf_samples = mcmc.trace('impact_event')[:]
In [ ]:
In [ ]:
In [ ]: