Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

map.flux() and sys.flux() give different results? #317

Open
ssagynbayeva opened this issue Jun 6, 2024 · 2 comments
Open

map.flux() and sys.flux() give different results? #317

ssagynbayeva opened this issue Jun 6, 2024 · 2 comments
Labels
bug Something isn't working

Comments

@ssagynbayeva
Copy link

ssagynbayeva commented Jun 6, 2024

Hi into the void!

I've been having this issue that map.flux() and map.design_matrix() are different from sys.flux() and sys.design_matrix(), when they should be the same (same map, same planet and star parameters), but the issue is only in in-transit flux.
image
Here's the zoomed-in transit flux:
image

Additionally, I calculated the transit duration given my parameters:
a = (1*np.square(truths['planet.porb']/365.25))**(1/3) * 215.03 # Solar radii;
b = 0
tdur = truths['planet.porb'] / np.pi * np.arcsin(np.sqrt((truths['planet.r']+1)**2 - b**2) / a)
which gives 0.08418049738331004
and it seems like map.flux() gives incorrect transit duration?

import starry
import matplotlib.pyplot as plt
import numpy as np

starry.config.quiet = True

true_obl = -30 

truths = {"planet.inc": 90,
          "planet.r": 0.1,
          "planet.t0": 0.5,
          "planet.porb": 1.,
          "star.inc": 90,
          "star.obl": true_obl,
          "star.prot": 5.,
          "sp.mu": 30,
          "sp.sigma": 1.,
          "u": [0.4, 0.2]}

map = starry.Map(15)
map.inc = truths['star.inc']
map.obl = truths['star.obl']

star = starry.Primary(map, r=1., m=1., prot=truths['star.prot']) 
planet = starry.Secondary(
    starry.Map(0,0),
    porb=truths['planet.porb'],
    t0=truths['planet.t0'],
    r=truths['planet.r'],
)

planet.inc = truths['planet.inc']

map.spot(contrast=0.9, radius=30, lon=40, lat=30)
map.spot(contrast=0.9, radius=20, lon=-10, lat=-30)

map.show()

sys = starry.System(star, planet)
t = np.arange(0, 5, 6 / 24 / 60)

X = sys.design_matrix(t).eval()

theta = (360 * t / truths["star.prot"]) % 360
xo, yo, zo = sys.position(t)
xo = xo.eval()[1]
yo = yo.eval()[1]
zo = zo.eval()[1]
A = map.design_matrix(
    theta=theta, xo=xo, yo=yo, zo=zo, ro=truths["planet.r"]
).eval()

plt.plot(X[:,:-1]-A)
plt.title('design matrix difference')
plt.show()

plt.plot(t, map.flux(theta=theta, xo=xo, yo=yo, zo=zo, ro=truths["planet.r"]).eval(), label='map.flux')
plt.plot(t, sys.flux(t, total=False)[0].eval(), label='sys.flux')
plt.legend()
plt.xlim(0.4,0.6)
plt.show()

a = (1*np.square(truths['planet.porb']/365.25))**(1/3) * 215.03 # Solar radii
b = 0
tdur = truths['planet.porb'] / np.pi * np.arcsin(np.sqrt((truths['planet.r']+1)**2 - b**2) / a)

print(tdur)

@ssagynbayeva ssagynbayeva added the bug Something isn't working label Jun 6, 2024
@ssagynbayeva
Copy link
Author

figured it out: map.design_matrix() does not consider planetary mass -- by default, the system assumes the planetary mass of 1, and caused the difference!

@ssagynbayeva
Copy link
Author

@lgrcia

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

1 participant