rodluger / starry Goto Github PK
View Code? Open in Web Editor NEWTools for mapping stars and planets.
Home Page: https://starry.readthedocs.io
License: MIT License
Tools for mapping stars and planets.
Home Page: https://starry.readthedocs.io
License: MIT License
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.
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!
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.
Ultimately we want to use starry for inference, so having analytic derivatives of the light curve will be super useful.
Check the equations in the paper for typos and errors and comb through the Mathematica scripts to ensure my proofs and derivations are correct.
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.
I still need to code up an exposure time option to numerically integrate the light curve over the exposure window.
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.
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
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?
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.
>>> 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.
@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?
...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.
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.
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
.
Does STARRY give a means of computing the Rossiter-McLaughlin effect?
The Doppler-shifted surface brightness could be decomposed as spherical harmonics.
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?
The computation of s_2 appears to have an outlier for r=0.1, b=0.9.
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
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?
Something like planetplanet's documentation: https://rodluger.github.io/planetplanet/
I tested the code a ton, but we still need a rigorous set of benchmarks in the tests/ folder
I'll add some test cases soon.
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.
Specifically, we need to compile references to all relevant papers in the literature.
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
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.
Probably the easiest way to benchmark stuff.
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
Probably in the next paper...
https://www.sciencedirect.com/science/article/pii/0019103576901858?via%3Dihub
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.
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!
Some proofs still need to be added to the appendix.
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:
Line 189 in 9476c4c
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).
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?
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.
Currently using the GSL elliptic integral routines. The Bulirsch ones may or may not be faster, so we should do some speed tests!
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.
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
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.
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}
A sample problem where we try to recover the map of HD 189733b (for instance) would be nice to include.
@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.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.
>>>import starry
>>>m = starry.Map()
>>>m[1, 0] = 1
>>>m.flux(yo=2, ro=2)
nan
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.