Planetary Eccentricity vs Eclipse Location

[exowanderer]

I think that STARRY might not be computing the eclipse location relative to eccentricity of the planetary orbit correctly. I may be attributing eccentricity or Omega incorrectly; please let me know if my code snippet below is correct.

Here is a plot of four different starry phase curves with different combinations of eccentricity and Omega (see legend)

This second figure below is identical to the one above; but zoomed in on the eclipse location, which has not changed; but was expected to.

The vertical dashed line shows the location of phase = 0.5. The shape of the eclipse bottom changes; it begins to look like a grazing eclipse at some versions of eccentricity and omega.

My understanding is that if the eccentricity is > 0.0 and the Omega ("The longitude of ascending node in degrees") is > 0.0, then the eclipse location should change (+/-) by up to several tenths in phase (0.3 - 0.7)

As described here:

def deltaphase_eclipse(ecc, omega):
    ''' Compute the phase location of the eclipse for a given orbit

        Parameters
        -----------
        ecc (float): eccentricity of the orbit; ranging from 0.0 - 1.0
        omega (float): argument of the ascending node in degrees; ranging 0-360

        Returns
        --------
        phase_shift (float): shift in phase for the center location of the eclipse
    '''
    phase_shift = 0.5*( 1 + (4. / np.pi) * ecc * np.cos(omega * 180/np.pi))

    return phase_shift

This script returns the following phase locations for the values of eccentricity and omega in the script below and figure above

print(deltaphase_eclipse(0.3,0)))
# 0.691

print(deltaphase_eclipse(0.3,90))
# 0.443

print(deltaphase_eclipse(0.3,180))
# 0.343

print(deltaphase_eclipse(0.3,270))
# 0.651

print(deltaphase_eclipse(0.3,360))
# 0.566

---

Code Snippet

from starry import kepler

from pylab import *; ion()
import numpy as np
try:
    import exomast_api # pip install git+https://github.com/exowanderer/exoMAST_API
except:
    !pip install git+https://github.com/exowanderer/exoMAST_API
    import exomast_api

def update_starry(system, times, model_params):
    ''' Update planet - system parameters '''
    planet = system.secondaries[0]

    edepth = model_params['edepth']
    day2night = model_params['day2night']
    phase_offset = model_params['phase_offset']

    planet.lambda0 = model_params['lambda0'] # Mean longitude in degrees at reference time
    planet.r = model_params['Rp_Rs'] # planetary radius in stellar radius
    planet.inc = model_params['inclination'] # orbital inclination 
    planet.a = model_params['a_Rs'] # orbital distance in stellar radius
    planet.prot = model_params['orbital_period'] # synchronous rotation
    planet.porb = model_params['orbital_period'] # synchronous rotation
    planet.tref = model_params['transit_time'] # MJD for transit center time

    planet.ecc = model_params['eccentricity'] # eccentricity of orbit
    planet.Omega = model_params['omega'] # argument of the ascending node

    max_day2night = np.sqrt(3)/2
    Y_1_0_base = day2night * max_day2night if edepth > 0 else 0

    # cos(x + x0) = Y_1_0*cos(x0) - Y_1n1*sin(x0)
    Y_1_0 = Y_1_0_base*cos(phase_offset/180*pi)
    Y_1n1 = -Y_1_0_base*sin(phase_offset/180*pi)
    Y_1p1 = 0.0

    planet[1,0] = Y_1_0
    planet[1, -1] = Y_1n1
    planet[1, 1] = Y_1p1
    planet.L = edepth / planet.flux()

    system.compute(times)

def plot_starry_phase_curve_phys(system, times, model_params, fig=None, ax=None, label='', figsize=(15,10)):

    if None in [fig, ax]:
        fig = plt.figure(figsize=figsize)
        ax = fig.add_subplot(111)

    update_starry(system, times, model_params)

    edepth = model_params['edepth']
    day2night = model_params['day2night']
    phase_offset = model_params['phase_offset']

    phase = (times - system.secondaries[0].tref) / system.secondaries[0].porb
    ax.plot(phase,system.lightcurve, label=label)

    # Universal
    if edepth > 0: ax.set_ylim(1-.1*edepth,1 + 1.1*edepth)

    ax.set_xlabel('Time [day]', fontsize=30)
    ax.set_ylabel('normalized flux', fontsize=30)

    for tick in ax.xaxis.get_major_ticks():
        tick.label.set_fontsize(20)

    for tick in ax.yaxis.get_major_ticks():
        tick.label.set_fontsize(20)

planet_info = exomast_api.exoMAST_API(planet_name='HD 189733 b')

lmax  = 1

phase = np.linspace(0, 1.0, 1000)
times = phase*planet_info.orbital_period + planet_info.transit_time

''' Instantiate Kepler STARRY model; taken from HD 189733b example'''
# Star
star = kepler.Primary()

# Planet
lambda0 = 90.0

planet = kepler.Secondary(lmax=lmax)
planet.lambda0 = lambda0 # Mean longitude in degrees at reference time
planet.r = planet_info.Rp_Rs # planetary radius in stellar radius
planet.L = 0.0 # flux from planet relative to star
planet.inc = planet_info.inclination # orbital inclination 
planet.a = planet_info.a_Rs # orbital distance in stellar radius
planet.prot = planet_info.orbital_period # synchronous rotation
planet.porb = planet_info.orbital_period # synchronous rotation
planet.tref = planet_info.transit_time # MJD for transit center time

planet.ecc = planet_info.eccentricity # eccentricity of orbit
planet.Omega = planet_info.omega # argument of the ascending node

# Instantiate the system
system = kepler.System(star, planet)

# Specific plottings
fpfs = 1000/1e6

### Run Plots

fig = plt.figure(figsize=(20,6))
ax = fig.add_subplot(111)

''' Instantiate Kepler STARRY model; taken from HD 189733b example'''
lambda0 = 90.0

model_params = {}
model_params['lambda0'] = lambda0 # Mean longitude in degrees at reference time
model_params['Rp_Rs'] = planet_info.Rp_Rs # planetary radius in stellar radius
model_params['inclination'] = planet_info.inclination # orbital inclination 
model_params['a_Rs'] = planet_info.a_Rs # orbital distance in stellar radius
model_params['orbital_period'] = planet_info.orbital_period # synchronous rotation
model_params['transit_time'] = planet_info.transit_time # MJD for transit center time

model_params['eccentricity'] = planet_info.eccentricity # eccentricity of orbit
model_params['omega'] = planet_info.omega # argument of the ascending node

model_params['edepth'] = fpfs
model_params['day2night'] = 1.0
model_params['phase_offset'] = 0

model_params['eccentricity'] = 0.0
model_params['omega'] = 0.0
label='eccentricity={}; omega={}'.format(model_params['eccentricity'], model_params['omega'])
plot_starry_phase_curve_phys(system, times, model_params, fig=fig, ax=ax, label=label)

model_params['eccentricity'] = 0.3
model_params['omega'] = 0.0
label='eccentricity={}; omega={}'.format(model_params['eccentricity'], model_params['omega'])

plot_starry_phase_curve_phys(system, times, model_params, fig=fig, ax=ax, label=label)

model_params['eccentricity'] = 0.3
model_params['omega'] = 90.0
label='eccentricity={}; omega={}'.format(model_params['eccentricity'], model_params['omega'])

plot_starry_phase_curve_phys(system, times, model_params, fig=fig, ax=ax, label=label)

model_params['omega'] = 180.0
label='eccentricity={}; omega={}'.format(model_params['eccentricity'], model_params['omega'])

plot_starry_phase_curve_phys(system, times, model_params, fig=fig, ax=ax, label=label)

model_params['omega'] = 270.0
label='eccentricity={}; omega={}'.format(model_params['eccentricity'], model_params['omega'])

plot_starry_phase_curve_phys(system, times, model_params, fig=fig, ax=ax, label=label)

ax.legend(loc=0, fontsize=15)

plt.axvline(0.5, ls='--');
plt.axhline(1.0, ls='--');

[rodluger]

@exowanderer I think you want to change w (the longitude of pericenter, which rotates the ellipse of the orbit about the orbital plane) instead of Omega (the argument of ascending node, which rotates the orbit on the plane of the sky, assuming it is viewed ~edge on). Changing Omega should have no observational effect for a single-planet system, since it just rotates everything on the plane of the sky.

[exowanderer]

Excellent. I didn't see that parameter!

I tried this out, keeping planet.Omega = 90; but varying planet.w between [0,90,180,270].

The eclipse location definitely varies with respect to eccentricity and omega; but the transit center location also moves. I can imagine why this is true; but I have not seen this behaviour before. Note that I did not modify planet.tref for any of these plots.

Could you confirm that the follow plots make sense when modifying omega (as per the legend).

[exowanderer]

Because the transit center is changing with planet.w, is it possible that the observer is being rotated around the star, instead of the planetary orbit?

[@kevin218 and I had to draw orbits on the white board to figure that out]

[rodluger]

@exowanderer @kevin218 The behavior in your plot is correct. The issue is in your interpretation of tref, which is not the time of transit t0. They only coincide for a circular orbit. Instead, tref is the time at which lambda0 is specified. Take a look at the last section of this tutorial. If you want to hold the time of transit fixed, you'll need to solve Kepler's equation to figure out what lambda0, the initial mean longitude, is at the time of transit.

I could code up a utility to do this for you, if you think that would be useful. Just let me know (or open an issue).

[exowanderer]

That makes sense. Also, this is the third time that the answer was "well, [Jonathan's] interpretation is at fault". Thank you for your help.

Would it be sufficient to marginalize over lambda0 along with eccentricity? I would hope that the statistics work out fine; and that I can avoid inflating the error bars. I'm sure solving Kepler's equations is more direct; but marginalizing is trivial.

[rodluger]

@exowanderer That should work as long as you're not conditioning on a given observed value of the transit time, t0. You'll need to solve Kepler's equation to figure out what t0 is at every step in the chain if you want to (say) weight the likelihood by some prior on t0. I'm going to re-open this issue, since I should implement a get_t0() function (or something like that) anyway.

[exowanderer]

I assume that you are offering to develop the get_t0 as a means to translate from lambda0 and tref into transit_time by solving Kepler's equations; but I wonder if this is not the same as marginalizing over t0 itself.

In the similar behaviour to our wrapping of the day2night contrast and phase_offset from before, would it be possible to input the transit time and eccentricity, instead of the reference time?

I don't mean to up end your existing infrastructure; but instead to give the user an option of one or the other. I assume that this would take the same amount of time as developing a get_t0 method; but would allow observers to introduce more observation driven parameters (similar to above).

[rodluger]

@exowanderer Thanks, that's a better way to think about it. I'll implement the general conversion between t0 and tref, lambda0 so the user can provide either one or the other.