hippke / pandora Goto Github PK
View Code? Open in Web Editor NEWA fast transit-fitting algorithm to search for transits of exomoons
License: GNU General Public License v3.0
A fast transit-fitting algorithm to search for transits of exomoons
License: GNU General Public License v3.0
Currently we use a flat prior [0,1] of 10^24 to 10^28 kg for the planetary mass.
Example: 0.01 ==> 0.01 * 10^28 = 10^26 kg
Example: 0.5 ==> 0.5*10^28 kg
Should it be logg'ed instead? So that 0.5 --> 10^26 kg?
Hi,
Thanks for this fun and useful software!
I was testing the example.py
and got the following error simply copying the code into a ipython shell:
In [5]: model = pandora.moon_model(params)
---------------------------------------------------------------------------
TypingError Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 model = pandora.moon_model(params)
File ~/.cache/pypoetry/virtualenvs/platonium-nqdyOqtx-py3.8/lib/python3.8/site-packages/pandoramoon/pandora.py:114, in moon_model.__init__(self, params)
112 self.numerical_grid = params.numerical_grid
113 self.time = params.time
--> 114 self.cache = create_occult_cache(self.u1, self.u2, dim=300)
File ~/.cache/pypoetry/virtualenvs/platonium-nqdyOqtx-py3.8/lib/python3.8/site-packages/numba/core/dispatcher.py:415, in _DispatcherBase._compile_for_args(self, *args, **kws)
409 msg = str(e).rstrip() + (
410 "\n\nThis error may have been caused by the following argument(s):\n%s\n"
411 % "\n".join("- argument %d: %s" % (i, err)
412 for i, err in failed_args))
413 e.patch_message(msg)
--> 415 error_rewrite(e, 'typing')
416 except errors.UnsupportedError as e:
417 # Something unsupported is present in the user code, add help info
418 error_rewrite(e, 'unsupported_error')
File ~/.cache/pypoetry/virtualenvs/platonium-nqdyOqtx-py3.8/lib/python3.8/site-packages/numba/core/dispatcher.py:358, in _DispatcherBase._compile_for_args.<locals>.error_rewrite(e, issue_type)
356 raise e
357 else:
--> 358 reraise(type(e), e, None)
File ~/.cache/pypoetry/virtualenvs/platonium-nqdyOqtx-py3.8/lib/python3.8/site-packages/numba/core/utils.py:80, in reraise(tp, value, tb)
78 value = tp()
79 if value.__traceback__ is not tb:
---> 80 raise value.with_traceback(tb)
81 raise value
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<built-in function empty>) found for signature:
>>> empty(UniTuple(int64 x 2), dtype=Literal[str](float32))
There are 2 candidate implementations:
- Of which 2 did not match due to:
Overload of function 'empty': File: numba/core/typing/npydecl.py: Line 504.
With argument(s): '(UniTuple(int64 x 2), dtype=unicode_type)':
No match.
During: resolving callee type: Function(<built-in function empty>)
During: typing of call at /lhome/nicholas/.cache/pypoetry/virtualenvs/platonium-nqdyOqtx-py3.8/lib/python3.8/site-packages/pandoramoon/occult.py (15)
File "../../.cache/pypoetry/virtualenvs/platonium-nqdyOqtx-py3.8/lib/python3.8/site-packages/pandoramoon/occult.py", line 15:
def create_occult_cache(u1, u2, dim):
<source elided>
ks = np.linspace(0.001, k_max, dim)
fs = np.empty((dim, dim), dtype="float32")
^
I'm using numba 0.51.2
and pandoramoon 1.0
.
Cheers,
Nicholas
Toshio Fukushima, "Fast computation of a general complete elliptic integral of third kind by half and double argument transformations" (Journal of Computational and Applied Mathematics, Volume 253, 2013).
https://www.sciencedirect.com/science/article/pii/S037704271300201X#!
Fixed: Bug was here:
Planetary transit duration at b=0 equals the width of the star
Formally correct would be: (R_star+r_planet) for the transit duration T1-T4
Here, however, we need T2-T3 (points where center of planet is on stellar limb)
https://www.paulanthonywilson.com/exoplanets/exoplanet-detection-techniques/the-exoplanet-transit-method/
New:
tdur_p = per_planet / pi * arcsin(sqrt(R_star ** 2) / a_planet)
Verification with compare_to_batman.py
at very high time resolution. All differences are <1e-7 now.
Take care of constant values like 2*
@jit(nopython=True, error_model='numpy')
def frac_diff3(x, y):
out = np.empty_like(x)
dt = x.dtype # Cast the constant using the dtype of the input
for i in range(x.shape[0]):
# Could also use np.float32(2) to always use same type, regardless of input
out[i] = dt.type(2) * (x[i] - y[i]) / (x[i] + y[i])
return out
https://nbviewer.org/github/numba/numba-examples/blob/master/notebooks/simd.ipynb
When plotting the planet-only model in the tutorial, it says
# Convert limb darkening back from q to u type:
u1, u2 = ld_invert(mll_planet_only_q1, mll_planet_only_q2)
This is the wrong way round and makes for an ugly result plot.
The correct way to transform LDs in q form to the usual (for plotting) u form is:
u1, u2 = ld_convert(mll_planet_only_q1, mll_planet_only_q2)
Make animation of transit a function of the main class
Use matplotlib ffmpeg directly, instead of writing PNGs and calling imagemagick afterwards
Currently many units in [km]. Change to R_P/R_{\odot}?
About 30% of the entire Pandora runtime goes into calculating an arctan and a tan:
def ellipse_pos(a, per, tau, Omega, i, time):
"""2D x-y Kepler solver WITHOUT eccentricity, WITHOUT mass"""
Q = 2 * arctan(tan((pi * (time - tau * per) / per))) # <== This is the offender
V = sin(Q) * cos((i / 180 * pi))
x = (cos(Omega / 180 * pi) * cos(Q) - sin(Omega / 180 * pi) * V) * a
y = (sin(Omega / 180 * pi) * cos(Q) + cos(Omega / 180 * pi) * V) * a
return x, y
Nestle has best speed (efficient samples per time) in some online comparisons
At the moment: Numerical estimate
In Luger Appendix C an analytical version is mentioned
Dear Mr. Hippke and Mr. de Leon,
I am currently utilizing your code as part of my master's thesis, which focuses on detecting exomoons using PLATO.
During my work, I encountered an issue related to the observation cadence and the timestep in the output data. Specifically, when configuring the system for 3456 observations per day, which equates to a 25 second cadence, the resulting timestep in the data is approximately 25.1 seconds.
Although I recognize that this problem may be caused by an oversight on my part, I have not found any way to fix this problem. Following my supervisor's suggestion, I am reaching out to provide feedback through this channel. I will present the code below, where I have made the problem as clear as possible.
I appreciate your time and any guidance you can offer on this matter.
With kind regards,
Nick
import numpy as np
import pandas as pd
import pandoramoon as pandora
import matplotlib.pyplot as plt
# Gravitational constant
G = 6.6743e-11 # m^3 kg^-1 s-2
#################
# Kepplers third law #
#################
def KEP_a_from_P(P, Mp, Ms): # p denotes primary, s denotes secondary
M_tot = Mp + Ms
a3 = P**2 * G * M_tot / (4 * np.pi**2)
return a3**(1/3)
def KEP_P_from_a(a, Mp, Ms): # p denotes primary, s denotes secondary
M_tot = Mp + Ms
P2 = a**3 * ( G * M_tot / (4 * np.pi**2) )**(-1)
return np.sqrt(P2)
#################
# System parameters #
#################
# Radii
R_star = 696342000 # Radius of the star in meters (Sun)
R_planet = 69911000 # Radius of the planet in meters (Jupiter)
R_moon = 6371000 # Radius of the moon in meters (Earth)
# Masses
M_star = 1.989e30 # Mass of the star in kg (Sun)
M_planet = 1.899e27 # Mass of the planet in kg (Jupiter)
M_moon = 5.972e24 # Mass of the moon in kg (Earth)
# Period of planet around star in days
P_sp_days = 4
# ... in seconds
P_sp_secs = 4 * 24 * 60 * 60
# Semi-major axis of the planet in meters
a_sp_meters = KEP_a_from_P(P_sp_secs, M_star, M_planet)
# ... in stellar radii
a_sp_Rstar = a_sp_meters / R_star
# Semi-major axis of the moon in meters
a_pm_meters = 128000000
# Period of the moon in days
P_pm_days = KEP_P_from_a(a_pm_meters, M_planet, M_moon) / (24 * 60 * 60)
# transit duration
def transit_duration(P, R_star, R_planet, a):
ff = P / np.pi
Rtot = R_star + R_planet
T = ff * np.arcsin(Rtot / a)
return T
#########
# Pandora #
#########
# Call Pandora and get model with these parameters
params = pandora.model_params()
params.R_star = R_star # [m]
params.u1 = 0.4089 # Default value
params.u2 = 0.2556 # Default value
# Planet parameters
params.per_bary = P_sp_days # Period of the planet [days]
params.a_bary = a_sp_Rstar # SMA planet [R_star]
params.r_planet = R_planet / R_star # Planetary radius [R_star]
params.b_bary = 0 # [0..1]
params.t0_bary = 1 # [days]
params.t0_bary_offset = 0 # [days]
params.M_planet = M_planet # [kg]
params.w_bary = 0 # [deg]
params.ecc_bary = 0 # [0..1]
# Moon parameters
params.r_moon = R_moon / R_star # [R_star]
params.per_moon = P_pm_days # [days]
params.tau_moon = 0 # [0..1]
params.Omega_moon = 0 # [0..180]
params.i_moon = 90 # [0..180]
params.e_moon = 0 # [0..1]
params.w_moon = 0 # [deg]
params.M_moon = M_moon
# Other model parameters
params.cadences_per_day = 3456 # [int] cadence of 25 seconds (24 * 60 * 60) / 25
params.epoch_distance = 4 # [days]
params.epochs = 183 # 2* 365.25 / 4 # [int]
params.supersampling_factor = 1 # [int] default value
params.hill_sphere_threshold = 1.1 # default value
params.epoch_duration = transit_duration(P_sp_days, R_star, R_planet, a_sp_meters)
# Obtain time grid
time = pandora.time(params).grid()
# Define model
model = pandora.moon_model(params)
# Evaluate model for each point in time grid
flux_total, flux_planet, flux_moon = model.light_curve(time)
# Get coordinates
xp, yp, xm, ym = model.coordinates(time)
#########
# Problem #
#########
print(f'First time point in days: {time[0]}')
print(f'Second time point in days: {time[1]}')
print(f'Time step in days: {time[1] - time[0]}')
print(f'Time step in seconds: {(time[1] - time[0]) * 24 * 60 * 60}')
print(f'Time step in seconds: {(time[2] - time[1]) * 24 * 60 * 60}')
The result:
First time point in days: 0.933785251160733
Second time point in days: 0.9340756667258174
Time step in days: 0.00029041556508446753
Time step in seconds: 25.091904823297995
Time step in seconds: 25.091904823307587
range (iterable (ndim,))
customizable (e.g., full range 0..180deg, not just the peak) see documentationshould give 10% speed gain
Currently stellar mass is not in the model
planetary semimajor axis and planetary period are free parameters
better to have a mass (density?) prior?
Speed
Accuracy vs. double precision
Priors are fixed parameter ranges
Some combinations of these are unphysical
Currently, no check is made and the model is always calculated
We can introduce sensible boundaries of moons e.g. outside of Hill sphere (like TLS: not search for things that can't exist)
Then, Pandora can return a model for those unphysical situations so that their LogLikelihood is very bad (e.g., np.ones for the moon)
Difficult but possible due to post-processing with Pumkin occultations
Current implementation:
transit_threshold_x = 3 * (a_moon / R_star) + 2 * r_planet + 2 * r_moon
if transit_threshold_x < 2:
transit_threshold_x = 2
OK for 90° inclination through testing.
Too loose for other configurations --> Can be made faster
Näherungsverfahren? Welcher Fehler?
Exakte Lösung? Wie Koordinatensystemtransformation?
e.g. Lightkurve quantity is converted automatically, but then performance is 10x slower
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.