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.

Sign up for GitHub

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

Jump to bottom

Enable modeling of starspots with both Fleck and Starry #627

Open
wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Enable modeling of starspots with both Fleck and Starry #627

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
0b4fb5b
initial commit of fleck/starry starspot models
erinmmay Feb 27, 2024
8aa829d
adding params to epf template, docs
erinmmay Feb 27, 2024
29c32a7
fixing things
erinmmay Feb 27, 2024
c50766f
fixing things part2
erinmmay Feb 27, 2024
6aedf18
fixing things part3
erinmmay Feb 27, 2024
9891a0c
adding fleck to yml, quick fix to fleck model
erinmmay Feb 27, 2024
2099338
sets fleck's slow option as the default for faster lc calculations
erinmmay May 1, 2024
4a32714
Merge branch 'main' into starspots
erinmmay May 1, 2024
040a769
addressing comments, updating fleck_tr to use new planetparams class
erinmmay May 25, 2024
5bca077
addressing formatting
erinmmay May 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
19 changes: 19 additions & 0 deletions demos/JWST/S5_fit_par_template.epf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,25 @@ u2 0.84 'free' 0 1 U
# osc_t0 0.3163 'fixed' 0 1 U
# osc_t1 0.3163 'free' 0.3 0.4 U

# ------------------
# ** Star Spot parameters **
Comment on lines +73 to +74
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please also add these changes to the docs/media/S5_fit_par_template.epf file which is shown on the https://eurekadocs.readthedocs.io/en/latest/ecf.html page

# each spot needs rad, lat, lon. see docs for details!
# for fleck: can only have one spot contrast (spotcon0), spotnpts is the number of temporal points to evaluate at, spotrad is relative to star
# for starry: assign one spotcon per spot (spotcon0, spotcon1, etc), spotres is the ydeg of the star map, spotrad is in degrees
# NOTE: fleck contrast = 1 - starry contrast
# ------------------
# spotstari 90 'fixed' 90 1 N
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While this parameter is set to fixed in the template, the demonstrated prior (should users choose to set the parameter to free) is 90 +/- 1 degree which seems excessively constraining. Wouldn't a prior of U(0,90) or something be more fitting if someone chose to fit for this parameter?

# spotrot 10 'fixed' 10 1 N
# spotnpts 100 'independent
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In your notes in ecf.rst, you mention that this should be ~200-500 for fleck - I recommend you stick with that for this template

# spotcon0 0.9 'free' 0 1 U
# spotrad0 0.2 'free' 0 360 U
# spotlat0 0 'free' -90 90 U
# spotlon0 0 'free' -180 180 U
# spotcon1 0.1 'free' 0 1 U
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably worth adding a comment at the end of this line to remind users that you should only add a spotcon1 parameter if you were using starry to model two spots since fleck only permits one spot contrast for all spots.

# spotrad1 0.2 'free' 0 1 U
# spotlat1 0 'free' -90 90 U
# spotlon1 0 'free' -180 180 U

# --------------------
# ** Systematic variables **
# Polynomial model variables (c0--c9 for 0th--3rd order polynomials in time); Fitting at least c0 is very strongly recommended!
Expand Down
2 changes: 1 addition & 1 deletion demos/JWST/S5_template.ecf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ rescale_err False # Rescale uncertainties to have reduced
fit_par ./S5_fit_par_template.epf # What fitting epf do you want to use?
verbose True # If True, more details will be printed about steps
fit_method [dynesty] #options are: lsq, emcee, dynesty (can list multiple types separated by commas)
run_myfuncs [batman_tr,polynomial] # options are: batman_tr, batman_ecl, sinusoid_pc, poet_tr, poet_ecl, poet_pc, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
run_myfuncs [batman_tr,polynomial] # options are: batman_tr, batman_ecl, fleck_tr, sinusoid_pc, poet_tr, poet_ecl, poet_pc, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas)
compute_ltt False # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect)
force_positivity False # Optional boolean for sinusoid_pc and poet_pc models. Set True to force positive phase variations.
num_planets 1 # Number of planets in the fit. Default is 1.
Expand Down
13 changes: 13 additions & 0 deletions docs/source/ecf.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1117,6 +1117,19 @@ This file describes the transit/eclipse and systematics parameters and their pri
- ``osc_t1`` - This offset determined the start phase of the sinusoidal function.

The time-dependent amplitude and period are defined as follows: ``amp = osc_amp * np.exp(-osc_amp_decay * (time - osc_t0))`` and ``per = osc_per * np.exp(-osc_per_decay * (time - osc_t0))``. The damped oscillator model is then defined as: ``osc = 1 + amp * np.sin(2 * np.pi * (time - osc_t1) / per)``. Note that ``osc[time < osc_t0] = 1``.

- Star Spot Parameters
- ``spotstari`` - The stellar inclination in degrees.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this parameter only matter if spotrot >> T_14? It would probably be helpful to add a small word of guidance to users on how they might set priors on this parameter or when they might choose to fix it to 90 degrees or something.

- ``spotrot`` - The stellar rotation rate in days. For fleck, only assign if you'd like to run in slow mode! (In slow mode the star rotates and spots move appropriately. Otherwise Eureka! will use fleck's fast mode which assumes the stellar rotation is >> transit time and spots are stationary)
- ``spotcon#`` - The spot contrast ratio. For fleck only assign one, for starry assign one per spot
- ``spotrad#`` - The spot radius. For fleck it is relative to the star, for starry it is in degrees
- ``spotlat#`` - The spot latitude.
- ``spotlon#`` - The spot longitude.
Fleck specific parameters:
- ``spotnpts`` - The number of temporal points to evalaute at. ~200-500 is good.
Starry specific parameters:
- ``spotnpts`` - The degree of spherical harmonics on the star (ydeg). ~30 is needed to appropriately model the spot.

- Systematics Parameters. Depends on the model specified in the Stage 5 ECF.
- ``c0--c9`` - Coefficients for 0th to 3rd order polynomials.

Expand Down
1 change: 1 addition & 0 deletions environment.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ dependencies:
- dynesty[version='>1.0'] # Lower limit needed for specific arguments
- emcee[version='>3.0.0'] # Lower limit needed for specific arguments
- flake8 # Needed for testing
- fleck
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please also add this to setup.cfg and environment_pymc3.yml

- george # Needed for GP
- h5py
- ipykernel
Expand Down
80 changes: 74 additions & 6 deletions src/eureka/S5_lightcurve_fitting/differentiable_models/StarryModel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,9 +132,43 @@ def setup(self):
p = temp.per*(24.*3600.)
Ms = ((2.*np.pi*a**(3./2.))/p)**2/const.G.value/const.M_sun.value

# Initialize star object
star = starry.Primary(starry.Map(udeg=self.udeg),
m=Ms, r=temp.Rs)
# check for spots and set parameters
if hasattr(self.parameters, 'spotrad0'):
spotrad = np.array([])
spotlat = np.array([])
spotlon = np.array([])
spotcon = np.array([])
for key in self.parameters.dict.keys():
if key.startswith('spot'):
if 'rad' in key:
spotrad = np.append(spotrad,
self.parameters.dict[key][0])
elif 'lat' in key:
spotlat = np.append(spotlat,
self.parameters.dict[key][0])
elif 'lon' in key:
spotlon = np.append(spotlon,
self.parameters.dict[key][0])
elif 'con' in key:
spotcon = np.append(spotcon,
self.parameters.dict[key][0])
elif 'rot' in key:
starrot = self.parameters.dict[key][0]
else:
spotres = self.parameters.dict[key][0]
Comment on lines +141 to +158
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code is going to break for mult_white or shared fits. I recommend you make use of the temp object I defined above which should only include the parameters for the current channel


# Initialize map object and add spots
map = starry.Map(ydeg=spotres, udeg=self.udeg)
for si in range(len(spotrad)):
map.spot(contrast=spotcon[si], radius=spotrad[si],
lat=spotlat[si], lon=spotlon[si])

# Initialize star object
star = starry.Primary(map, m=temp.Ms, r=temp.Rs, prot=starrot)
else:
# Initialize star object
star = starry.Primary(starry.Map(udeg=self.udeg),
m=temp.Ms, r=temp.Rs)

if hasattr(self.parameters, 'limb_dark'):
if self.parameters.limb_dark.value == 'kipping2013':
Expand Down Expand Up @@ -313,9 +347,43 @@ def update(self, newparams, **kwargs):
p = temp.per*(24.*3600.)
Ms = ((2.*np.pi*a**(3./2.))/p)**2/const.G.value/const.M_sun.value

# Initialize star object
star = starry.Primary(starry.Map(udeg=self.udeg),
m=Ms, r=temp.Rs)
# check for spots and set parameters
if hasattr(self.parameters, 'spotrad0'):
spotrad = np.array([])
spotlat = np.array([])
spotlon = np.array([])
spotcon = np.array([])
for key in self.parameters.dict.keys():
if key.startswith('spot'):
if 'rad' in key:
spotrad = np.append(spotrad,
self.parameters.dict[key][0])
elif 'lat' in key:
spotlat = np.append(spotlat,
self.parameters.dict[key][0])
elif 'lon' in key:
spotlon = np.append(spotlon,
self.parameters.dict[key][0])
elif 'con' in key:
spotcon = np.append(spotcon,
self.parameters.dict[key][0])
elif 'rot' in key:
starrot = self.parameters.dict[key][0]
else:
spotres = self.parameters.dict[key][0]
Comment on lines +356 to +373
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See my comment above about the use of the temp object to make sure this code works for shared fits and mult_white fits


# Initialize map object and add spots
map = starry.Map(ydeg=spotres, udeg=self.udeg)
for si in range(len(spotrad)):
map.spot(contrast=spotcon[si], radius=spotrad[si],
lat=spotlat[si], lon=spotlon[si])

# Initialize star object
star = starry.Primary(map, m=temp.Ms, r=temp.Rs, prot=starrot)
else:
# Initialize star object
star = starry.Primary(starry.Map(udeg=self.udeg),
m=temp.Ms, r=temp.Rs)

if hasattr(self.parameters, 'limb_dark'):
if self.parameters.limb_dark.value == 'kipping2013':
Expand Down
16 changes: 16 additions & 0 deletions src/eureka/S5_lightcurve_fitting/fitters.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,10 @@ def callback(theta):
if meta.isplots_S5 >= 1:
plots.plot_fit(lc, model, meta, fitter=calling_function)

# Plot star spots
if 'fleck_tr' in meta.run_myfuncs and meta.isplots_S5:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if 'fleck_tr' in meta.run_myfuncs and meta.isplots_S5 >= 3:

plots.plot_fleck_star(lc, model, meta, fitter=calling_function)
Comment on lines +165 to +167
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we add a starry spot plotting function (which I think we should), then we'd need to call that in each of these edited places


# Plot GP fit + components
if model.GP and meta.isplots_S5 >= 1:
plots.plot_GP_components(lc, model, meta, fitter=calling_function)
Expand Down Expand Up @@ -434,6 +438,10 @@ def emceefitter(lc, model, meta, log, **kwargs):
# Plot fit
if meta.isplots_S5 >= 1:
plots.plot_fit(lc, model, meta, fitter='emcee')

# Plot star spots
if 'fleck_tr' in meta.run_myfuncs and meta.isplots_S5 >= 3:
plots.plot_fleck_star(lc, model, meta, fitter='emcee')

# Plot GP fit + components
if model.GP and meta.isplots_S5 >= 1:
Expand Down Expand Up @@ -900,6 +908,10 @@ def dynestyfitter(lc, model, meta, log, **kwargs):
if meta.isplots_S5 >= 1:
plots.plot_fit(lc, model, meta, fitter='dynesty')

# Plot star spots
if 'fleck_tr' in meta.run_myfuncs and meta.isplots_S5 >=3:
plots.plot_fleck_star(lc, model, meta, fitter='dynesty')

# Plot GP fit + components
if model.GP and meta.isplots_S5 >= 1:
plots.plot_GP_components(lc, model, meta, fitter='dynesty')
Expand Down Expand Up @@ -1012,6 +1024,10 @@ def lmfitter(lc, model, meta, log, **kwargs):
if meta.isplots_S5 >= 1:
plots.plot_fit(lc, model, meta, fitter='lmfitter')

# Plot star spots
if 'fleck_tr' in meta.run_myfuncs and meta.isplots_S5>=3:
plots.plot_fleck_star(lc, model, meta, fitter='lmfitter')

# Plot GP fit + components
if model.GP and meta.isplots_S5 >= 1:
plots.plot_GP_components(lc, model, meta, fitter='lmfitter')
Expand Down