Issues for large occultors

[rodluger]

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.

[rodluger]

The issue is for m >= 0 harmonics starting at l = 4 (though there seems to be a bit of instability for l = 3). This script reproduces the error. I'm going to investigate.

[rodluger]

Negative orders get unstable when the inclination isn't 90, so something about the symmetry was hiding it for those.

I figured out what is causing this. If I double the planet radius the instability mostly goes away. The oscillations are a numerical instability for large r and large b. I ran the plots above for the Sun occulting the Earth, so r = 100 and b is between 99 and 101.

The first idea I had was to factor k^2 as

k^2 = (1 - (b - r)) * (1 + (b - r)) / (4 * b * r)

to avoid error when b and r are large. I implemented this expression when r > 1 but this didn't fix it.

It turns out the issue is in the angle phi (and maybe in lam as well, but I'm not sure). When the occultor is very, very large, sin(phi) is always very close to 1. If I fix it at 1, the instability in the l = 1, m = 1 term (green in the plots below) goes away. The plot on the left is with the current code and the plot on the right is where I fix sin(phi) = 1:

The lmax = 4 curves also look much better, but there's still an instability, which I think is due to the linear limb darkening term (orange curve above), which I'll investigate next.

[rodluger]

FWIW, here's what I get using 24-digit precision (from the Boost library's MPFR wrapper).

[ericagol]

These look pretty good !I think there should be a way to do this with Float64. What does the brightening mean? Does the planet have negative flux in that model?

[rodluger]

Yeah, the map is negative near the limb. I set the Y{0,0} coefficient to 1 when I plot each of these curves, but for some spherical harmonics the absolute value of the specific intensity can actually exceed 2 * sqrt(PI) (the amplitude of Y{0,0}), leading to negative flux. If only there were an analytic way to compute whether a map goes negative! I'm still hopeful about Sturm's theorem...

[rodluger]

Also, a keen eye would notice a tiny bump at first contact for the m = 0 curve. This actually happens at all orders and degrees, including Y_{0,0}, regardless of occultor size, floating point precision, elliptic integral tolerance, and just about every other knob I can turn: It's negligible for all realistic cases, but it's still bugging me. It's not an issue with my out-of-transit flux, since a similar bump occurs just prior to complete occultation (when the flux should absolutely tend to zero).

[rodluger]

Ideas:

• Multiply by br32 in P() rather than in computeM() when ksq < 1
• Tabulate sin * cos in the I and H primitive integrals.
• Taylor expand sinphi and sinlam as follows:

[ericagol]

Is this the same issue you emailed me about? Is this resolved now?

[rodluger]

@ericagol Thanks for your approximation! That should fix the instability in the s2 term. There are still strong instabilities in the other terms, which I'm working on right now. I think I just found a fix, so I'll keep you posted!

[rodluger]

@ericagol FWIW the weird offset above was due to the value of PI I was using, so that's fixed!

[ericagol]

One time I used a radiation transfer code that used two slightly different values for the Thomson cross section in different parts of the code. It took a long time to find the bug which manifested as negative fluxes.

Eric Agol Astronomy Professor University of Washington

On Mar 31, 2018, at 10:49 AM, Rodrigo Luger notifications@github.com wrote:

@ericagol FWIW the weird offset above was due to the value of PI I was using, so that's fixed!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[rodluger]

Resolved using several Taylor expansions and re-parameterizations. However, for l > 8 things look pretty hopeless when the occultor is larger than the occulted. Discussion in paper.

[ericagol]

There should be a way to improve this

Eric Agol Astronomy Professor University of Washington

On Apr 2, 2018, at 4:56 PM, Rodrigo Luger notifications@github.com wrote:

Resolved using several Taylor expansions and re-parameterizations. However, for l > 8 things look pretty hopeless when the occultor is larger than the occulted. Discussion in paper.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[rodluger]

I'd love to chat when you get back. This is where we currently stand: I get sub-ppb precision at all Ylm degrees up to at least l = 10 for occultors that are smaller than the occulted body. When the radius ratio exceeds 1, I get sub-ppm precision for l < 5 and sub-ppt precision for l < 8. Even ppt precision is fine, since this is defined relative to the flux of the occulted body (i.e., the planet), which we likely won't be able to measure down to 1 ppt any time this century! But yes, there should be an elegant and analytic way around this. For Ylm terms that do not depend on z, I am currently approximating the occultor limb as a quartic and computing the line integral directly -- this works quite well. But for the terms with z dependence, I end up having to take the integral of functions that look like

x^n * (a + b x^2 + c x^4)^1.5

which look quite nasty. I couldn't find anything in the A&S book to help with these, but maybe you can come up with something!

P.S.: BTW, with quadruple precision, there are no numerical issues anywhere. It's just ~20x slower.

[ericagol]

Does this issue still need attention?

-Eric

On 4/3/18 10:14 AM, Rodrigo Luger wrote:

I'd love to chat when you get back. This is where we currently stand: image https://user-images.githubusercontent.com/9323819/38264128-f109ac8e-3726-11e8-892d-a84fc7f64c34.png I get sub-ppb precision at all |Ylm| degrees up to at least |l = 10| for occultors that are smaller than the occulted body. When the radius ratio exceeds 1, I get sub-ppm precision for |l < 5| and sub-ppt precision for |l < 8|. Even ppt precision is /fine/, since this is defined relative to the flux of the occulted body (i.e., the planet), which we likely won't be able to measure down to 1 ppt any time this century! But yes, there should be an elegant and analytic way around this. For Ylm terms that do not depend on |z|, I am currently approximating the occultor limb as a quartic and computing the line integral directly -- this works quite well. But for the terms with |z| dependence, I end up having to take the integral of functions that look like

|x^n * (a + b x^2 + c x^4)^1.5 |

which look quite nasty. I couldn't find anything in the A&S book to help with these, but maybe you can come up with something!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-378326959, or mute the thread https://github.com/notifications/unsubscribe-auth/AAO30CKU3qKJzm_9ycCQtPwcmCXNwdroks5tk630gaJpZM4S5eHD.

[ericagol]

For an occulter radius >> 1, it should be sufficient to approximate the edge of the occultor as being straight.

[rodluger]

For low degree spherical harmonics I think this is true. But as the degree increases, the flux becomes ever more sensitive to the curvature of the limb because the length scales on the surface get smaller. We could certainly suppress the numerical instabilities by doing this, but we would introduce error due to the approximation for large [image: l]. But it could certainly help -- I haven't tested it in detail, so we could still come out ahead with the approximation.

The other issue is that in order to do this, we'd have to take the line integral of Equation (34), whose general term looks like

[image: x^i y^j (1 - x^2 - y^2)^\frac{3}{2}]

The general case of the line integral is a sum of hypergeometric functions. For integer [image: i] and [image: j] they reduce to square roots and arc tangents. So there could be a nice recursive way to solve this, but I think it would take some serious work.

On Fri, Apr 20, 2018 at 9:57 PM Eric Agol notifications@github.com wrote:

For an occulter radius >> 1, it should be sufficient to approximate the edge of the occultor as being straight.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-383255632, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK-xvjMixzWUTv1C1tVNj_yYIsarVks5tqoP_gaJpZM4S5eHD .

[ericagol]

I may have figured out a computational way forward on this (although I haven't fully worked it out).

For M_{0,0} at small k, there is a cancellation in leading orders in k between the two terms, causing a k^4 dependence.

I can show where this arises by transforming to:

where m = k^2. Since the leading (constant) terms in dE/dm and dK/dm have opposite signs, these cancel, causing the first term to be of order m^2. Since the last term is also proportional to m^2, the entire equation is proportional to m^2=k^4.

A similar approach could be taken with the other three first terms (I have yet to work these out), and then the expansion can be carried out after applying recursion to obtain the higher order terms.

[rodluger]

Awesome, thanks! I am currently Taylor expanding M for large occultors to 10th order, and I think we might be doing the same thing. But I hard code the expansions as a large table -- your way could help automate this at arbitrary order. Check out my Mathematica file here:

https://github.com/rodluger/starry/blob/master/tex/notebooks/computeMTaylor.pdf?raw=true

The expansion helps but doesn't fix the issue completely: see Figure 12. The peaks in the curves in the lower panel are where I switch to the Taylor version. Going to higher order in k^2 doesn't seem to bring the peak error down since I think once we get to 10th order the expressions get dominated by floating point error. Note that they also diverge at very large radius, albeit at a slower rate than the original expressions.

On Sun, Apr 22, 2018, 12:38 AM Eric Agol notifications@github.com wrote:

I may have figured out a computational way forward on this (although I haven't fully worked it out).

For M_{0,0} at small k, there is a cancellation between the two terms, causing a k^4 dependence.

I can show where this arises by transforming to:

[image: untitled] https://user-images.githubusercontent.com/243664/39091175-c10a9670-45a3-11e8-98c6-8cdb17b9cd75.jpg

where m = k^2. Since the leading (constant) terms in dE/dm and dK/dm have opposite signs, these cancel, causing the first term to be of order m^2. Since the last term is also proportional to m^2, the entire equation is proportional to m^2=k^4.

A similar approach can be taken with the other four first terms, and then the expansion can be carried out after iterating.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-383352861, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK1WjibyPkLIUEXmsVj_QoOhZWSR6ks5tq_sZgaJpZM4S5eHD .

[ericagol]

I figured out the series expansion for the M_{p,q} coefficients, so there's now no need to hard-code the coefficients and compute them to 10th order, I think.

I've re-expressed the M_{p,q} terms in terms of the function \varepsilon_3, which has a leading term proportional to k^4:

\varepsilon_3 &= -m(dE/dm + (1-m)dK/dm)\ &= m^2\frac{1}{8} \int0^{\pi/2} \frac{\sin^2{2\phi}}{(1-m\sin^2{\phi})^{3/2}} d\phi\ &= m^2 \frac{\pi}{4} \sum{n=1}^\infty \frac{[(2n-1)!!]^2 m^{n-1}}{4^n (n+1)!(n-1)!}

where m=k^2.  The series converges quickly for k < 0.1.

Then, the M_{p,q} terms for p=0,2 q=0,2 are given by:

{\cal M}_{0,0} &= [(8-16m)\varepsilon3 + 4m^2K(m)]/3\ {\cal M}{0,2} &= [(8-28m-12m^2)\varepsilon3 + (22-6m)m^2K(m)]/15\ {\cal M}{2,0} &= [(32-52m+12m^2)\varepsilon3 + (-2+6m)m^2K(m)]/15\ {\cal M}{2,2} &= [(32-76m+36m^2-24m^3)\varepsilon_3 + (-2+30m-12m^2)m^2K(m)]/105

I haven't set up the machinery to compute the light curves with these coefficients, so I don't yet know if there is still inaccuracy for large p/q when computing the light curves.

Attached is some Julia code for computing these coefficients.  I tried to apply the recursion relations separately to the terms multiplying \varepsilon_3 and K(m), but this didn't seem to help the precision.

-Eric

Awesome, thanks! I am currently Taylor expanding M for large occultors to 10th order, and I think we might be doing the same thing. But I hard code the expansions as a large table -- your way could help automate this at arbitrary order. Check out my Mathematica file here:

https://github.com/rodluger/starry/blob/master/tex/notebooks/computeMTaylor.pdf?raw=true

The expansion helps but doesn't fix the issue completely: see Figure 12. The peaks in the curves in the lower panel are where I switch to the Taylor version. Going to higher order in k^2 doesn't seem to bring the peak error down since I think once we get to 10th order the expressions get dominated by floating point error. Note that they also diverge at very large radius, albeit at a slower rate than the original expressions.

On Sun, Apr 22, 2018, 12:38 AM Eric Agol notifications@github.com wrote:

I may have figured out a computational way forward on this (although I haven't fully worked it out).

For M_{0,0} at small k, there is a cancellation between the two terms, causing a k^4 dependence.

I can show where this arises by transforming to:

[image: untitled]

https://user-images.githubusercontent.com/243664/39091175-c10a9670-45a3-11e8-98c6-8cdb17b9cd75.jpg

where m = k^2. Since the leading (constant) terms in dE/dm and dK/dm have opposite signs, these cancel, causing the first term to be of order m^2. Since the last term is also proportional to m^2, the entire equation is proportional to m^2=k^4.

A similar approach can be taken with the other four first terms, and then the expansion can be carried out after iterating.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub

https://github.com/rodluger/starry/issues/23#issuecomment-383352861, or mute the thread

https://github.com/notifications/unsubscribe-auth/AI5FK1WjibyPkLIUEXmsVj_QoOhZWSR6ks5tq_sZgaJpZM4S5eHD .

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-383376628, or mute the thread https://github.com/notifications/unsubscribe-auth/AAO30IDjhnZ0YVHqFtDhS0SppBUSfNyzks5trHHegaJpZM4S5eHD.

Computes the M_{p,q} coefficients defined in (D37) for

the r0 >> 1 limit in which k2 << 1 and r0-1 < b < r0+1

(this is a huge occultor).

using PyPlot include("ellk_bulirsch.jl") include("ellec_bulirsch.jl")

function eps3_of(k2::Real) kk = ellk_bulirsch(k2) eps3 = kk-ellec_bulirsch(k2)-0.5k2kk return eps3,kk end

function eps3_series(k2::Real)

Series evaluation of eps3 for k2 < 1:

@assert(k2 < 1) tol = eps(k2) eps3 = zero(k2)

Compute n=1 term

term = pi/32 eps3 += term n=2 while abs(term) > tolabs(eps3) term = (2n-1)^2k2/4/(n^2-1) eps3 += term n +=1 end return eps3*k2^2 end

Test these

k2 = logspace(-4,-1,10) for i=1:10 eps3_x,kk = eps3_of(big(k2[i])) eps3_a =eps3_series(k2[i]) println(convert(Float64,eps3_x)," ",eps3_a," ",eps3_a/convert(Float64,eps3_x)-1.0) end

function mpq_init(k2::Real)

Initialize the M_{p,q} computation:

Alternate form:

if k2 < 0.1

eps3 = eps3_series(k2)

kk = ellk_bulirsch(k2)

else

eps3,kk=eps3_of(k2)

end

m00 = (8.-16.k2)/3.eps3 + 4.k2^2kk/3

m02 = (8.-28.k2-12.k2^2)/15.eps3 + (22.-6.k2)k2^2kk/15

m20 = (32.-52.k2+12.k2^2)/15.eps3 + (-2.+6.k2)k2^2kk/15

m22 = (32.-76.k2+36.k2^2-24.k2^3)/105.eps3 + (-2.+30.k2-12.k2^2)k2^2kk/105.

m001 = (8.-16.k2)/3. # eps3 term m002 = 4.k2^2/3 # K(m) term m021 = (8.-28.k2-12.k2^2)/15. m022 = (22.-6.k2)k2^2/15 m201 = (32.-52.k2+12.k2^2)/15. m202 = (-2.+6.k2)k2^2/15 m221 = (32.-76.k2+36.k2^2-24.k2^3)/105. m222 = (-2.+30.k2-12.k2^2)k2^2/105. return m001,m002,m021,m022,m201,m202,m221,m222 end

r0 = big(100.0)

r0 = 100.0 nb = 1000 b = linspace(r0-1,r0+1,nb)

k2 = (1.-(r0-b).^2)./(4.b.r0) k2_big = big.(k2) pmax = 10; qmax=10 mpq = zeros(Float64,nb,pmax+1,qmax+1) mpq1 = zeros(Float64,nb,pmax+1,qmax+1) mpq2 = zeros(Float64,nb,pmax+1,qmax+1) mpq_big = zeros(BigFloat,nb,pmax+1,qmax+1) mpq1_big = zeros(BigFloat,nb,pmax+1,qmax+1) mpq2_big = zeros(BigFloat,nb,pmax+1,qmax+1) m001 = 0.; m021=0.; m201=0.; m221=0. m002 = 0.; m022=0.; m202=0.; m222=0. m001_big = 0.; m021_big=0.; m201_big=0.; m221_big=0. m002_big = 0.; m022_big=0.; m202_big=0.; m222_big=0. eps3 = zeros(Float64,nb); kk=zeros(Float64,nb) eps3_big = zeros(BigFloat,nb); kk_big=zeros(BigFloat,nb) for i=1:nb

Alternate form:

if k2[i] < 0.1 eps3[i] = eps3_series(k2[i]) kk[i] = ellk_bulirsch(k2[i]) else eps3[i],kk[i]=eps3_of(k2[i]) end m001,m002,m021,m022,m201,m202,m221,m222 = mpq_init(k2[i]) mpq[i,1,1]=m001eps3[i]+m002kk[i]; mpq[i,1,3]=m021eps3[i]+m022kk[i]; mpq[i,3,1]=m201eps3[i]+m202kk[i]; mpq[i,3,3]=m221eps3[i]+m222kk[i] mpq1[i,1,1]=m001; mpq1[i,1,3]=m021; mpq1[i,3,1]=m201; mpq1[i,3,3]=m221 mpq2[i,1,1]=m002; mpq2[i,1,3]=m022; mpq2[i,3,1]=m202; mpq2[i,3,3]=m222 if k2_big[i] < 0.1 eps3_big[i] = eps3_series(k2_big[i]) kk_big[i] = ellk_bulirsch(k2_big[i]) else eps3_big[i],kk_big[i] =eps3_of(k2_big[i]) end m001_big,m002_big,m021_big,m022_big,m201_big,m202_big,m221_big,m222_big = mpq_init(k2_big[i]) mpq_big[i,1,1]=m001_bigeps3_big[i]+m002kk_big[i]; mpq_big[i,1,3]=m021_bigeps3_big[i]+m022_bigkk_big[i] mpq1_big[i,1,1]=m001_big; mpq1_big[i,1,3]=m021_big mpq2_big[i,1,1]=m002_big; mpq2_big[i,1,3]=m022_big mpq_big[i,3,1]=m201_bigeps3_big[i]+m202kk_big[i]; mpq_big[i,3,3]=m221_bigeps3_big[i]+m222_bigkk_big[i] mpq1_big[i,3,1]=m201_big; mpq1_big[i,3,3]=m221_big mpq2_big[i,3,1]=m202_big; mpq2_big[i,3,3]=m222_big end

m00 = (8.0-12.0k2)/3.eps1+(-8+16k2)/3.eps2

m02 = (8.-24.k2)/15.eps1 +(-8+28k2+12k2.^2)/15.*eps2

m20 = (32-36k2)/15.eps1+(-32+52k2-12k2.^2)/15.*eps2

m22 = (32-60k2+12k2.^2)/105.eps1+(-32+76k2-36k2.^2+24k2.^3)/105.*eps2

println("M_{0,0}: ",m00," ",m00_alt)

clf()

Now, use recursion relation to compute M_{p,q}:

for p=4:2:pmax for q = 0:2:2 d3 = 2p+q-(p+q-2)(1.-k2) d4 = (3-p)k2 mpq1[:,p+1,q+1] = (d3.mpq1[:,p-1,q+1]+d4.mpq1[:,p-3,q+1])/(p+q+3) mpq2[:,p+1,q+1] = (d3.mpq2[:,p-1,q+1]+d4.mpq2[:,p-3,q+1])/(p+q+3) mpq[:,p+1,q+1] = mpq1[:,p+1,q+1].eps3+mpq2[:,p+1,q+1].kk d3_big = 2p+q-(p+q-2)(1.-k2_big) d4_big = (3-p)k2_big mpq1_big[:,p+1,q+1] = (d3_big.mpq1_big[:,p-1,q+1]+d4_big.mpq1_big[:,p-3,q+1])/(p+q+3) mpq2_big[:,p+1,q+1] = (d3_big.mpq2_big[:,p-1,q+1]+d4_big.mpq2_big[:,p-3,q+1])/(p+q+3) mpq_big[:,p+1,q+1] = mpq1_big[:,p+1,q+1].eps3_big+mpq2_big[:,p+1,q+1].kk_big println("p ",p," q ",q," max diff: ",maximum(abs.(mpq[:,p+1,q+1]-convert(Array{Float64,1},mpq_big[:,p+1,q+1])))) clf() plot(b,mpq[:,p+1,q+1]) plot(b,convert(Array{Float64,1},mpq_big[:,p+1,q+1])) read(STDIN,Char) end for q=4:2:qmax d1 = q+2+(p+q-2)(1.-k2) d2 = (3-q)(1.-k2) mpq1[:,p+1,q+1]=(d1.mpq1[:,p+1,q-1]+d2.mpq1[:,p+1,q-3])/(p+q+3) mpq2[:,p+1,q+1]=(d1.mpq2[:,p+1,q-1]+d2.mpq2[:,p+1,q-3])/(p+q+3) mpq[:,p+1,q+1]=mpq1[:,p+1,q+1].eps3+mpq2[:,p+1,q+1].kk d1_big = q+2+(p+q-2)(1.-k2_big) d2_big = (3-q)(1.-k2_big) mpq1_big[:,p+1,q+1]=(d1.mpq1_big[:,p+1,q-1]+d2_big.mpq1_big[:,p+1,q-3])/(p+q+3) mpq2_big[:,p+1,q+1]=(d1.mpq2_big[:,p+1,q-1]+d2_big.mpq2_big[:,p+1,q-3])/(p+q+3) mpq_big[:,p+1,q+1]=mpq1_big[:,p+1,q+1].eps3_big+mpq2_big[:,p+1,q+1].kk_big println("p ",p," q ",q," max diff: ",maximum(abs.(mpq[:,p+1,q+1]-convert(Array{Float64,1},mpq_big[:,p+1,q+1])))) clf() plot(b,mpq[:,p+1,q+1]) plot(b,convert(Array{Float64,1},mpq_big[:,p+1,q+1])) read(STDIN,Char) end end

[rodluger]

Fantastic, thanks! I'll look into this tomorrow.

[rodluger]

@ericagol what happened to the double factorial on this line:

  term *= (2*n-1)^2*k2/4/(n^2-1)

of function eps3_series in your code?

[ericagol]

@ericagol https://github.com/ericagol what happened to the double factorial on this line:

|term = (2n-1)^2*k2/4/(n^2-1) |

of |function eps3_series| in your code?

It is supposed to be computed recursively with the (2*n-1)^2 term (which should yield the [(2n-1)!!]^2 term under recursion).

[rodluger]

Gotcha. Had to look up the double factorial!

On Tue, Apr 24, 2018 at 6:04 PM Eric Agol notifications@github.com wrote:

@ericagol https://github.com/ericagol what happened to the double factorial on this line:

|term = (2n-1)^2*k2/4/(n^2-1) |

of |function eps3_series| in your code?

It is supposed to be computed recursively with the (2*n-1)^2 term (which should yield the [(2n-1)!!]^2 term under recursion).

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-384079327, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK6GgQPVU7rw90bX9LvRhEEH8FKFCks5tr5NIgaJpZM4S5eHD .

[rodluger]

I tabulated all the half-factorials, so I'm computing it as I'm going to do some stability tests with these new expressions.

[rodluger]

@ericagol I got this to work, and it does indeed help with the stability. But it's not quite as good as my earlier brute-force 10th order Taylor expansion. Here are the diagnostic plots:

Original version, no Taylor expansion:

My 10th order expansion:

Eric's varepsilon3 expansion:

The curves are not monotonic because I treat terms with even l + m and odd l + m differently. Here's what they look like for even terms (top) and odd terms (bottom) separately:

For even terms, I can easily approximate the occultor limb as a quartic function and compute the line integral analytically. This leads to convergence at large occultor radius in the top plot. The peak in the curves occurs where I switch over to the quartic expression (there's a different sweet spot for each spherical harmonic degree, which I figured out manually). If we really wanted to push the error down below ppm, I could probably figure out a sextic or

For odd terms, I use your varepsilon expansion. It pushes the instability out to larger occultor radius, but the expressions still diverge as the radius increases. If we could find a limiting value of these expressions at very large r, that might fix this issue. Ideally we would approximate the occultor limb as a line or a parabola, but I wasn't able to find an analytic solution to the line integral -- it looks even nastier than the expressions we are currently evaluating to compute M().

Any ideas?

P.S.: I created the stability branch for us to play around with this stuff. I coded up your varepsilon3 approach in that branch.

[rodluger]

We clearly have some work to do still.

[rodluger]

After a lot of hacking on this, I'm mostly happy with the workarounds Eric and I came up with. I came up with better stability plots that explore the error on the s terms in b-r-l space. They run pretty quickly, so they're always available here. Basically, our relative errors are sub-ppm everywhere for lmax <= 8 and for 1e-5 <= r <= 300, which is all the cases anyone should ever need (for the next several decades at least!)

Note that the fractional errors are larger, but that's because the flux can get very small in places. We ultimately don't care whether our model gives 2e-16 when the true answer is 1e-16, so I think that's a bad metric.

We can obviously improve on these in the future. A relatively straightforward thing to do is expand the occultor limb with a sixth- or eighth-order polynomial (it's currently a quartic) for the even terms. For the odd terms, I'm still stumped, but we can probably come up with something if we need to.

I'm closing this issue for now.

[ericagol]

I'm re-opening with a couple of questions:

1). Did transforming from the occulter from radius r to radius 1 make a difference?

2). Can you also make some plots with the derivatives as well?

3). Also, to transition from fractional to absolute errors, you can use asinh errors (i.e. luptitudes, where the transition (f_0) occurs for depth of transit for a uniform source.

[rodluger]

@ericagol

1. The transformation from r --> 1, b --> b / r, and 1 --> 1 / r didn't seem to make any difference as far as I could tell. I didn't test this exhaustively, though.
2. I'll work on derivative plots.
3. Interesting. I'll take a look at that!

[rodluger]

What actually helped was switching between the original expressions and the Taylor versions depending on the values of r, l, and b. I sat down yesterday and mapped out the errors in this 3-parameter space and found much better transitions between the expressions. In principle I could use the derivative information as well.

[ericagol]

I tried using the generalized elliptic integrals from Bulirsch (1969) to compete the M_{p,q} integrals, but the results are very unstable.

[ericagol]

I computed M{p,q} coefficients in polynomials in m= k^2 in BigFloat precision (256 bit), and find that all coefficients are zero up to m^(p/2+1)! So, part of the problem for large l (lowercase L) is that the Taylor series expansion to a fixed order gets less precise. Also, there may be some cancellation in computing J{uv}.

The nice thing is that these polynomial coefficients only need to be computed once - which can be done at high precision, and converted to Floats and stored for a particular problem.

[rodluger]

That's a good point. I noticed this in my coefficients but didn't think to take the expansion to higher order for larger p. Do you think this would help with the error?

On Thu, May 24, 2018, 8:05 AM Eric Agol notifications@github.com wrote:

I computed M_{p,q} coefficients in polynomials in m= k^2 in BigFloat precision (256 bit), and find that all coefficients are zero up to m^(p/2+1)! So, part of the problem for large l (lowercase L) is that the Taylor series expansion to a fixed order gets less precise.

The nice thing is that these polynomial coefficients only need to be computed once - which can be done at high precision, and converted to Floats and stored.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-391748738, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK7KKtY7csQ-8sc6q2nOiNdncj7iSks5t1swugaJpZM4S5eHD .

[ericagol]

Yes, I think so. I'm computing M_{p,q} with the polynomial expansion formula, and getting exact agreement between floating point and BigFloat (for the floating point case, I just set the coefficients up to (k^2)^(p/2+1) to zero). There is a natural transition point at k^2 = 1/2 for which values of r > \sqrt{1-b^2} we can transition to the series expansion (for k^2 below this value, the transition becomes more complicated, although we could still do this as well if k^2=1/2 requires too many terms in the expansion).

[ericagol]

There is probably a fundamental reason that all of these cancellations occur, but I haven't divined it yet.

[rodluger]

Wow, this is awesome. Thanks for digging into this!

On Thu, May 24, 2018 at 9:02 AM Eric Agol notifications@github.com wrote:

There is probably a fundamental reason that all of these cancellations occur, but I haven't divined it yet.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-391769143, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK-36sGiIgbysOOk0GSn3wOpAwDzOks5t1tmfgaJpZM4S5eHD .

[ericagol]

I've got the odd terms computing to ~machine precision now (I think)!

However, the even terms still have large errors - I need to look into what's causing this problem.

[rodluger]

Fantastic!!! The even terms don't depend on M_{p,q}. Those are the ones I performed a quartic expansion https://github.com/rodluger/starry/blob/a4ef7b0744013a127d27314d87a337ed1d3f5377/starry/taylor.h#L203 of the limb for.

On Thu, May 24, 2018 at 12:45 PM Eric Agol notifications@github.com wrote:

I've got the odd terms computing to ~machine precision now (I think)!

However, the even terms still have large errors - I need to look into what's causing this problem.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-391836341, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK5vDgm0Ly6FVUZj5jv9_52Hvu62Sks5t1w3agaJpZM4S5eHD .

[ericagol]

I think the cancellation is coming from when \phi = -\pi/2 + \epsilon, for small \epsilon.

In this case, -2\cos{\phi} ~ -2\epsilon, while 2\phi+\pi ~ 2\epsilon, so all of the even and odd terms have nearly equal values, but opposite signs! There should be a way to solve this problem without having to carry out the approximate quartic limb expansion; I'll look into this.

[rodluger]

Great, thanks!

On Thu, May 24, 2018 at 1:24 PM Eric Agol notifications@github.com wrote:

I think the cancellation is coming from when \phi = -\pi/2 + \epsilon, for small \epsilon.

In this case, -2\cos{\phi} ~ -2\epsilon, while 2\phi+\pi ~ 2\epsilon, so all of the even and odd terms have nearly equal values, but opposite signs! There should be a way to solve this problem without having to carry out the approximate quartic limb expansion; I'll look into this.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/23#issuecomment-391847539, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK_4hw5ewp7gXJM61bCDyYlpKqOu7ks5t1xbrgaJpZM4S5eHD .

[rodluger]

@ericagol I increased the order of my Taylor expansion of M_{p,q} from 12 to 16 but saw no change in the error of the s2 term for large occultor radius. I'm still getting errors that creep up to 1 ppm for l = 8. What's the maximum difference between your s2 terms computed using BigFloat and the Taylor expansion for l = 8?

[ericagol]

What r value? 100?

RL: Yeah.

[ericagol]

I may have figured out part of the issue with K_uv, caused by alternating signs in I_uv. For r>>1, let \phi = -\pi/2 + \kappa, where \kappa << 1. Then, the integral definition for I_uv gives:

to first order in \kappa, where \phi^\prime = \phi-3\pi/2. I tried simply transforming from \phi to \kappa, which helps the rounding error a bit, but not completely.

I{u,v} = (-1)^v \int{-\kappa}^\kappa d\phi^\prime (\sin{\phi^\prime})^u (\cos{\phi^\prime})^v \approx (-1)^v \frac{2}{u+1} \kappa^{u+1}

So, for fixed u, the v terms alternate sign, which causes several cancelling in the K_uv sum since b ~ r. In practice, we should keep higher order terms. The even u terms have \kappa & \sin{\kappa} dependence, while the odd terms have \sin{\kappa} dependence. We can keep track of the coefficients of these terms, which have the most significant cancellation, and then the subsequent terms give the final values of K_uv.

[ericagol]

By the way: I_uv(\phi) = H_uv(\phi) (they are the same integral, save for the difference in being a function of \phi or \lambda). H_uv doesn't seem to have this rounding error problem since there is no sum over terms, and \lambda is closer to \pi/2, not -\pi/2, for large \r.

[ericagol]

@rodluger Regarding your question above, for r = 100, the fractional M_pq values are being computed to machine precision compared with the BigFloat version. I'm still having trouble computing with theK_uv terms due to the alternating signs of I_uv, so the so I'm getting huge errors at large l values for the even mu, nu terms. The odd terms seem to be computed to machine precision (although I just noticed a bug whereby some odd terms give huge, unphysical fluxes: mu=5, nu=11 is going up to 2x10^7, while mu=1,nu=15 is fine).

[rodluger]

Interesting. I'll have to run some tests, but I thought my M_{p,q} values were good to machine precision, but somewhere upstream the summations lead to cancellation errors. Are you saying you get machine precision in the actual s2 values for odd mu simply by Taylor expanding M? Have you done anything else to J or P?

[ericagol]

@rodluger I take that back - I'm getting agreement between the double precision and BigFloat computations, but I think both are wrong. I'm also still getting numerical noise for the even values. Not sure what is going on or how to fix it yet.

[rodluger]

Ah, too bad. Keep me posted, but don't stress too much over this! The code is stable for all major, practical applications, and I'm happy submitting it even if we don't find a general solution to these instabilities. But it would be nice, of course, to make it work to arbitraryl!

[ericagol]

I'm getting there... still debugging code. I'll revisit next week.