rodluger/starry

Tools for mapping stars and planets.

Tools for mapping stars and exoplanets.

starry's Issues

Issues for large occultors

Earth secondary eclipse -- numerical instabilities at ingress/egress. Could be due to one of the instabilities in the s2 term or something new. I'm going to investigate.

Higher order limb darkening

This is easy. Figure out the general expression for the nth order limb darkening coefficient in terms of the spherical harmonic coefficients and allow users to provide an array of mu coefficients when instantiating a Star object.

Luminosity units don't make sense

I'm going to ditch the luminosity option in the Star object. This way users can just specify the luminosity of the planet relative to that of the star, which will make unit analysis less messy, especially when converting luminosity to flux and factoring in the telescope size and such.

I need to re-think how to compute the total luminosity for a limb-darkened star. I need to treat this case separately since the star isn't really described by a 2nd order spherical harmonic expansion; rather, it just appears to be from the observer's perspective (for instance, as the star rotates, a dark nightside doesn't come into view!). Note that because of this, the units of the light curve don't make any sense currently!

Numerical integrator needs work

The adaptive mesh technique doesn't converge to the right answer as the requested precision is increased indefinitely. We need stronger ways to ensure the mesh is refined near the occultor boundary.

Get autodiff working

Ultimately we want to use starry for inference, so having analytic derivatives of the light curve will be super useful.

Check my math

Check the equations in the paper for typos and errors and comb through the Mathematica scripts to ensure my proofs and derivations are correct.

simultaneous transits

Eventually I should use Green's theorem as in Pal (2012) to compute overlapping transit light curves, but for now the code should at least compute simultaneous-but-not-overlapping transits correctly. This line needs to be changed to make the flux calculation cumulative.

Flux integration

I still need to code up an exposure time option to numerically integrate the light curve over the exposure window.

Instabilities in the s2 term

The elliptic integrals K and Pi in the Mandel & Agol solution diverge when k --> 1. The result is still finite, since things cancel, but we should find a way to re-parametrize this.

I think some of the functions for the higher order terms might also diverge when k --> inf, but I need to investigate this further.

Equations (D37) & (D38)

In equation (D37), what conditions hold for p in the second case and q in the third case?

In equation (D38), d_4 simplifies to (3-p)*k^2

Compile speed-up?

The 12th order taylor expansion in taylor.h causes compile time to be several minutes. I could (should) reduce the order of the expansion, but are there any other tricks to get the compile time down?

Rodrigo doesn't understand consts in C++

Can someone check my code to see where I should and shouldn't precede declarations with const? Does it really matter, particularly if we're going to autodiff stuff in the future? I haven't at all been consistent with my usage.

Gradients are NaN when b = 0

>>> python -c "import starry; m = starry.grad.Map(); m[1,-1] = 1; print(m.flux(ro=0.1))"
[[ 0. nan nan nan nan nan nan nan]]

When b = 0, ksq is infinity and its derivatives are not defined. I tried overriding the computation of ksq for this case but wasn't able to fix the issue. This is an edge case, so not a huge deal, but would be nice to fix. We could easily bypass this by finding analytic closed-form expressions for the s vector when b = 0. Shouldn't be hard.

Gradient of elliptic Pi

@dfm Can you help me wrap my head around the gradient of the elliptic integral Pi on this line? The integral is a function of two parameters, so shouldn't there be two derivatives? This function currently only returns one -- which appears to be the sum of each of the derivatives. What's going on there?

Gradients broken when integrating flux over exposure time

I'm probably doing something dumb. These gradients don't look right (and there are small instabilities popping up in places):

Uninitialized variables...

...are probably in computesT() somewhere. I'm getting some weird values every now and then in the flux that are not consistently reproducible. I need to run valgrind to find the source.

Make the code photodynamical

Right now, the user inputs arrays of (x,y) positions for the occultor. This is great for some applications, but I want to allow users to specify orbital parameters and get a light curve out, either via a simple Keplerian solver or an N-body code. I coded all this up for planetplanet, so incorporating it into starry shouldn't be too hard.

Re-think how to integrate flux over exposure time

Integrated flux calculation is currently broken for individual bodies and when using starry.grad. This is a straightforward fix. I just need to sit down and re-code some of the stuff in orbital.h.

Rossiter-McLaughlin effect?

Does STARRY give a means of computing the Rossiter-McLaughlin effect?

The Doppler-shifted surface brightness could be decomposed as spherical harmonics.

Comments on the paper:

1). Are the spherical-harmonic angles \phi & \theta (which appears in section 2.1) different from the \phi and \theta which appear later in the paper (\theta is referenced just after equation 16, while \phi is discussed just after equation 24)? If so, it might be good to use some distinguishing marker to indicate that these are not the same variable.

2). Are the coordinates defined just after equation 2 left-handed? If so, this seems to contradict the description in Figure 1, which is right-handed, with the observer at +\infty.

3). Is the n subscript in equation (7) the same as the subscript n in equation 5?

4). In Figure 2, is it true that \phi is negative?

5). Can the Y_{lm} coefficients be allowed to vary with time? This might be useful, for instance, for modeling a spotted star with evolving spots.

6). Does first case in D32 require \mu to be even?

b = 1-r error in s_2

The computation of s_2 appears to have an outlier for r=0.1, b=0.9.

How do we ensure our surface maps are everywhere positive?

Given a set of spherical harmonic coefficients, how do we (easily | quickly | analytically) find out whether they lead to a map that is positive everywhere? This could be important for enforcing a physical prior when doing MCMC.

I think we need something analogous to Sturm's Theorem, but on the surface of the sphere in 2D. Here's some literature I found that may be helpful:

https://academic.oup.com/jlms/article-abstract/s2-19/3/451/849062?redirectedFrom=PDF
https://www.ppsloan.org/publications/shdering.pdf
http://adsbit.harvard.edu//full/1992SvA....36..220K/0000220.000.html

Better error handling

Most of the errors thrown in starry are sloppy and look like this, where I print to screen with std::cout and abort with exit(1). It would be much nicer to have a central logging system with error codes and such. I defined a few error codes here but this needs to be done more consistently throughout the code.

@dfm, is there a good C++ library for error handling/logging we could use?

Create documentation

Something like planetplanet's documentation: https://rodluger.github.io/planetplanet/

Rigorous test suite

I tested the code a ton, but we still need a rigorous set of benchmarks in the tests/ folder

Fix numerical instabilities as ksq --> inf.

I'll add some test cases soon.

Allow user to provide an input map as an array

Currently the user can provide an image or a healpix map, but it would be useful to allow the user to simply input a 2d numpy array of lat/lon coordinates. This will be trivial to implement.

Write introduction

Specifically, we need to compile references to all relevant papers in the literature.

Mapping of stellar pulsations in eclipsing binaries

As with many exoplanet ideas, there is cross-fertilization with stellar astrophysics. The starry formalism could perhaps be used in mapping pulsations of stars in eclipsing binaries during eclipses:

http://adsabs.harvard.edu/abs/1974ApJ...190..637N
http://cds.cern.ch/record/1010338/files/0701459.pdf
http://iopscience.iop.org/article/10.1086/491666
https://academic.oup.com/mnras/article/416/3/1601/959788
http://adsabs.harvard.edu/abs/2005ASPC..333..221B

Of course there is also the application of measuring starspots:

http://adsabs.harvard.edu/abs/1997MNRAS.287..556C
http://adsabs.harvard.edu/abs/1997MNRAS.287..567C

and doppler imaging:

http://adsabs.harvard.edu/abs/1987ApJ...321..496V

perhaps this could also be applied to polarimetric reconstruction of stars:

https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-2966.2010.16715.x

Move code to C++

The code is currently pretty slow since I've put zero effort into optimizing it. I plan to switch it over to C/C++ and use something like ctypes to interface with Python.

Compare to SPIDERMAN

Probably the easiest way to benchmark stuff.

Use new Mandel & Agol routines

Eric coded these up in Julia. Use these to compute the s_2 vector term. Note that we need to multiply by 2*pi/3 to get the right normalization.

ellk_bulirsch.txt
transit_quad.txt
ellec_bulirsch.txt
ellpic_bulirsch.txt

Application to Galilean moons

Probably in the next paper...
https://www.sciencedirect.com/science/article/pii/0019103576901858?via%3Dihub

Switch to templated Eigen types instead of homebuilt matrices

This will be easier to maintain and get higher performance. I'll see what I can do to take a stab at this over the weekend.

Mathematica scripts to compute matrices A1, A2, and A are wrong

The example matrices in the Appendix do not match up with what the C++ code outputs. The C++ code is correct and has been thoroughly vetted, so I need to go through the Mathematica scripts and find the bug!

Finish proofs

Some proofs still need to be added to the appendix.

r >= 1 and b=0 might cause an error in this line

If r >=1 and b=0, the source is completely occulted, so perhaps s2 isn't called. But, s2 is called with b=0 and r >=1, this line may cause a problem:

starry/starry/lld.h

Line 189 in 9476c4c

Lambda = T(-2. / 3.) * pow(1. - r2, 1.5);

Speed ups

When the occultor is not touching the limb of the planet, cos(lam) and cos(phi) are both zero and sin(lam) and sin(phi) are both unity. A lot of CPU cycles can be saved if we code this case up separately to avoid long recurrence relations whose answer is zero (or a constant).

phase curve v. rotation curve

I think what I'm technically computing are "rotational light curves", not phase curves, since the phase of a planet refers specifically to the illumination pattern due to the star, which doesn't matter for emission. Should I change the terminology throughout the paper?

Rename `u` to `axis`

I'm reserving u to refer to the limb darkening coefficient vector, so the rotation axis will now be called axis. I'm going to update all code, scripts, and notebooks, as well as the texfile to reflect this.

Switch to Bulirsch routines?

Currently using the GSL elliptic integral routines. The Bulirsch ones may or may not be faster, so we should do some speed tests!

ellk_bulirsch.txt
ellec_bulirsch.txt
ellpic_bulirsch.txt

Speed up the limb darkening calculations

Currently, the calculation for a quadratically limb-darkened star is about 6 times slower than batman:

Batman currently uses a polynomial approximation to the elliptic integrals of the first and second kind, so that could be accounting for most of that difference. Another factor is that I'm currently computing all terms in the solution vector up to lmax=2. We only need the m = 0 terms for limb-darkened stars. I tried to prevent computation of m!=0 terms when limb darkening is applied but got weird results. Stay tuned.

Deriving surface brightness of the Earth from EPOXI moon occultation

Figure 6 reminded me that it would be cool to see what could be inferred from the color-dependent EPOXI observation of the Earth transited by the Moon:

https://www.nasa.gov/topics/solarsystem/features/epoxi_transit.html

Here's a four-color version of the movie:

https://www.dropbox.com/s/yykr3n91zqmf0oq/Earth4_4colors-1202.mov?dl=0

The one challenge, of course, is that the Earth's emission is due to the convolution of the albedo pattern with the illumination profile of the Sun (and perhaps a scattering function?).

-Eric

One more instability in the s2 term

The solution for the s2 term (linear limb darkening) is unstable as k2 --> 1, which occurs when b + r --> 1. Here's a plot of the flux in the vicinity of the instability. It's present on both sides of the boundary, so we should probably find a Taylor expansion to apply there.

This script reproduces the issue.

Problem computing M_{p,q} and J_{u,v} when b=0 (and solution).

For b=0, k^2 = \infty, so M_{p,q} and J_{u,v} become problematic to compute.

However, for b=0,

J_{u,v} = \int_{\pi-\phi}^{2\pi+\phi} c_\varphi^u s_\varphi^v (1-r^2)^{3/2} d\varphi

= (1-r^2)^{3/2} I_{u,v}

Application to real data

A sample problem where we try to recover the map of HD 189733b (for instance) would be nice to include.

Instabilities!

@ericagol: I created a branch to keep track of the instabilities in the code. Whenever anyone commits changes, add [stability] somewhere to the commit message to automatically generate the diagnostic plots on travis. These will be force-pushed to the master-stability branch (or whatever-stability, where whatever is the current branch name).

The plots are generated for every spherical harmonic, currently up to 4th order. They are plots in the b-r plane of the fractional and relative error compared to calculations done using 128-bit precision. The subplots below the main plots are zoomed-in regions about the usual suspects: b = r, b = r - 1 and b = r + 1.

The attached plot is for the Y_{1,0} spherical harmonic (the s2 linear limb darkening term). As you can see, we've managed to fix most of the instabilities except

• b = r for 0.5 < r < 2. That's described in this issue.
• b = r - 1 for r > 1. This corresponds to the infinitesimal time before totality, where only a tiny sliver of the body is visible and the total flux is very close to zero. Since the relative error is small (see right panel), I'm not worried about this instability. We shouldn't really care if our calculations give 1e-10 when the true flux is 2e-10, for instance. That's a large (100%) fractional error, but completely unmeasurable.

Yet another instability in s2

This one happens when b ~ r and k^2 < 1 on this line because the elliptic integral is discontinuous at b = r. I suspect the numerical integration scheme is choking near that value.

When b = r, I can evaluate vPI in the limit:

but it would be good to have a reparametrization in the vicinity.

NaNs when b = r > 1 for s2 term

>>>import starry
>>>m = starry.Map()
>>>m[1, 0] = 1
>>>m.flux(yo=2, ro=2)
nan

s_n for m=-1; and possibly for other Y_lm values.

I just noticed that s_n = 0 whenever m=-1.

This makes sense since in this case there is odd symmetry about the x-axis,
so the flux below the axis (y < 0) cancels with the flux of opposite sign
above the axis (y > 0).

We can either skip computing these, or perhaps use this as a check that the code is working fine?

This symmetry can also be used to check that m=1 is being computed correctly if
the occultor is placed on the y-axis (x=0). A rotation needs to be made, but then
the resulting flux should be zero if the only non-zero Y_lm has m=1.

In addition, some other symmetries might be exploited as a check. For b=0, some Y_lm values alternate between positive and negative flux around the origin, so these should give zero flux, I think.