Question about returned timesteps/cadence

[nickvw97]

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

[hippke]

In the code, a day is simply defined to have 60*60*24 seconds. The params.cadences_per_day are then just an int based on that number. If you have a hard requirement to set an exact cadence length (in seconds), then you have to modify params.cadences_per_day so that it fits as good as possible. At the moment higher accuracy is not implemented. You would need to modify the Pandora code to do that, which should be relatively straightforward.