I am using python package pytransit to model exoplanet transits. Here is my current code:
from pytransit import UniformModel
import numpy as np
def calc_params(period, rp_over_rs):
period_s = period*3600
G = 6.674*10**(-11)
M = 0.6 * 1.989*10**30
R_s = 0.01 * 6.957*10**8
R_p = rp_over_rs * R_s
a = ((period_s**2*G*M)/(4*np.pi**2))**(1/3)
theta = (R_s + R_p)/a
return a/R_s, theta
model_times = np.arange(0, 100, 1)
tm = UniformModel()
tm.set_data(model_times)
period = np.arange(3, 30, 0.25)
for p in period:
a_over_rs, theta = calc_params(p, 1)
inclination = np.arange(np.pi/2-theta, np.pi/2+theta, 0.0001)
min_fluxes = []
for i in inclination:
fluxes = tm.evaluate_ps(k=1, t0=1, p=p, a=a_over_rs, i=i)
min_fluxes.append(np.min(fluxes))
I want to put minimum transit depth np.min(fluxes) as a function of orbital period and inclination angle on a 3D plot. The problems I am having come from:
The range of inclination angles
90 +/- thetadepends on the semi-major axis -- which is itself dependent on orbital period -- so the number of values in the arraymin_fluxeschanges based on the period.The function that produces the list of fluxes in the transit
tm.evaluate_pscannot take numpy lists as input, so doing something like in this post usingmeshgridseems impossible.Of course, using a colormap on a 2D plot like in this post seems plausible, but still begs the question of how to ensure a given minimum flux is tied back to the period and inclination that produced it. I thought that using a 2D numpy array could work (i.e
array[period][inclination] = min_fluxbut numpy arrays cannot be indexed by decimals like I have in my range of inclinations.
For me this is a complicated problem, both in terms of how to approach it and what tools I should use to solve it, and I am currently very stuck on this. If anyone could offer advice or help me out, I would greatly appreciate it.