Some other questions:

[ericagol]

• [x] Do we need an expression near b=0? k^2 -> \infty, which shows up outside of elliptic integrals.
• [x] Code just in starry: paper focuses on spherical harmonics & math. Can compute arbitrary order limb darkening, but only need m=0 terms.
• [x] limbdark: rewriting Mandel & Agol linear term with stable formula. Write down m=0 relations.
• [x] Compare double-precision versions of Julia & python codes.
• [ ] Do we need handle special cases separately for n != 2 s_n? b+r = 1; b=r (or in vicinity)?

[rodluger]

@ericagol When b = 0, k^2 = infty, so my recursion formula for J_{u,v} doesn't work. The stable expression I derived near b = 0 comes from Taylor expanding the term in parentheses in the integrand of in terms of b. Each term in the series is then a constant times some power of s, so the integral can easily be written in terms of I_{u,v}. Here's a sketch of how I do this in Mathematica.

It would be nice to compute a general expression for the nth term in the series so I don't have to hard-code the coefficients. Mathematica returns a DifferenceRoot function that I don't quite understand.

Can you come up with an easy way to Taylor expand (1 - r^2 - b^2 - 2 b r s)^{3/2} about b = 0, assuming s is a constant?

[ericagol]

I'll think about it... this is going to be problematic when r ~ 1 since then the portion of the star near y=-1 will be unocculted.

[ericagol]

I think the solution is to use my transformed equations (as I did in section D.2.3 for k^2 <1). See issue #83 - the primitive integral is just a polynomial in \delta = (b-r)/(2r) for the \mu/2 even case, which should be stable to evaluate for small \delta. For the (\mu-1)/2 even case, I get the following expression:

where u=s_\varphi^2.

It's clear from the form of this expression that it should be stable as k^2 -> \infty. You can either expand the (1-u/k^2)^{3/2} term as a series expansion in k^2, or evaluate it in terms of elliptic integrals, and then do upward or downward iteration, whichever is stable. This approach should actually work for all k^2 >1.

\mathcal{P}(\mathbf{G}_n) = 2(2r)^{l-1}(1-(b-r)^2)^{3/2} \sum_{i=0}^{\tfrac{\mu-1}{4}+\tfrac{\nu-1}{2}}\mathcal{A}_{i,\tfrac{\mu-1}{4},\tfrac{\nu-1}{2}} \int_0^1 \frac{u^{i+\tfrac{\mu-1}{4}} \left(1+\frac{u}{k^2}\right)^{3/2}}{u^{1/2}(1-u)^{1/2}}du

[rodluger]

Awesome, I'll try this out. Thanks!

On Thu, Jun 7, 2018, 4:38 PM Eric Agol notifications@github.com wrote:

I think the solution is to use my transformed equations (as I did in section D.2.3 for k^2 <1). See issue #83 https://github.com/rodluger/starry/issues/83 - the primitive integral is just a polynomial in \delta = (b-r)/(2r) for the \mu/2 even case, which should be stable to evaluate for small \delta. For the (\mu-1)/2 even case, I get the following expression:

[image: image] https://user-images.githubusercontent.com/243664/41131475-d754b55e-6a70-11e8-91fa-31f5b493c57c.jpeg

It's clear from the form of this expression that it should be stable as k^2 -> \infty. You can either expand the (1-u/k^2)^{3/2} term as a series expansion in k^2, and do upward or downward iteration, whichever is stable. This approach should actually work for all k^2 >1.

\mathcal{P}(\mathbf{G}n) = 2(2r)^{l-1}(1-(b-r)^2)^{3/2} \sum{i=0}^{\tfrac{l-1}{2}}\mathcal{A}_{i,\tfrac{\mu-1}{4},\tfrac{\nu-1}{4}} \int_0^1 \frac{u^{i+\tfrac{\mu-1}{4}} \left(1+\frac{u}{k^2}\right)^{3/2}}{u^{1/2}(1-u)^{1/2}}

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/88#issuecomment-395599088, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK5LEWpT6JP6McV_YQb0zjkJBPNWYks5t6blogaJpZM4UbpBW .

[ericagol]

A more familiar looking form:

\mathcal{P}(\mathbf{G}_n) &= 2(2r)^{l-1}(1-(b-r)^2)^{3/2} \int_{-\pi/2}^{\pi/2} s_\varphi^{\tfrac{\mu-1}{2}} (1-s_\varphi^2)^{\tfrac{\mu-1}{4}} (\delta + s_\varphi^2)^{\tfrac{\nu-1}{2}}(1-k^{-2}s_\varphi^2)^{3/2} d\varphi\\ \mathcal{P}(\mathbf{G}_n) &= 2(2r)^{l-1}(1-(b-r)^2)^{3/2} \sum_{i=0}^{\tfrac{\mu-1}{4}+\tfrac{\nu-1}{2}}\mathcal{A}_{i,\tfrac{\mu-1}{4},\tfrac{\nu-1}{2}} \int_{-\pi/2}^{\pi/2} s_\varphi^{\tfrac{\mu-1}{2}+2i}\left(1-k^{-2}s_\varphi^2\right)^{3/2} d\varphi\\ \mathcal{P}(\mathbf{G}_n) &= 2(2r)^{l-1}(1-(b-r)^2)^{3/2} \sum_{i=0}^{\tfrac{\mu-1}{4}+\tfrac{\nu-1}{2}}\mathcal{A}_{i,\tfrac{\mu-1}{4},\tfrac{\nu-1}{2}} J_{(\mu-1)/4+i}\\ J_v & = \int_{-\pi/2}^{\pi/2} d\varphi s_\varphi^{2v} \left(1-k^{-2}s_\varphi^2\right)^{3/2}

[ericagol]

The only difference in this case is that J_v is evaluated from -\pi/2 to \pi/2 rather than -\kappa/2 to \kappa/2.

[rodluger]

Cool!

On Thu, Jun 7, 2018, 5:30 PM Eric Agol notifications@github.com wrote:

The only difference in this case is that J_v is evaluated from -\pi/2 to \pi/2 rather than -\kappa/2 to \kappa/2.

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/88#issuecomment-395607812, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK3mybY4HBG31A4rY4nJpxDmAEV0Zks5t6cWlgaJpZM4UbpBW .

[rodluger]

I get this for the expansion of J near k2 -> infty:

[rodluger]

@ericagol Ok, I implemented this for all cases and it works great! Right now, I have a constant small threshold for b below which I use the Taylor expansion. But I'm finding that the instability occurs at higher values of b as l increases. It looks like around l = 15 the instability becomes significant for b < 0.1. Your formula, on the other hand, seems to be stable at all l. We're going to need to do some more rigorous tests to check when we should switch between the two.

[ericagol]

@rodluger Great! Or perhaps just use the new formula everywhere?

[rodluger]

Yeah I think you finally convinced me :) Though I want to do the speed tests still. Most people will never exceed l = 3 or 4, so if my equations are faster, we should keep those in the code.

On Fri, Jun 8, 2018 at 12:23 PM Eric Agol notifications@github.com wrote:

@rodluger https://github.com/rodluger Great! Or perhaps just use the new formula everywhere?

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/88#issuecomment-395863613, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK2PVPCoMktc71fAHCi8rh8G9WO1eks5t6s8sgaJpZM4UbpBW .

[rodluger]

When you get a chance, can you write up the general solutions for ksq > 1 in the Tex document?

On Fri, Jun 8, 2018 at 12:24 PM Rodrigo Luger rodluger@gmail.com wrote:

Yeah I think you finally convinced me :) Though I want to do the speed tests still. Most people will never exceed l = 3 or 4, so if my equations are faster, we should keep those in the code.

On Fri, Jun 8, 2018 at 12:23 PM Eric Agol notifications@github.com wrote:

@rodluger https://github.com/rodluger Great! Or perhaps just use the new formula everywhere?

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/rodluger/starry/issues/88#issuecomment-395863613, or mute the thread https://github.com/notifications/unsubscribe-auth/AI5FK2PVPCoMktc71fAHCi8rh8G9WO1eks5t6s8sgaJpZM4UbpBW .

[ericagol]

Where: (3/2 n)

Incidentally,

\mathcal{J}_v & = \int_{-\pi/2}^{\pi/2} d\varphi \sum_{n=0}^\infty s_\varphi^{2n} (-1)^n \binom{3/2}{n} (k^2)^{-n}\\ & = \sum_{n=0}^\infty (-1)^n \binom{3/2}{n} (k^2)^{-n} \int_{-\pi/2}^{\pi/2} d\varphi s_\varphi^{2n}\\ & = \sum_{n=0}^\infty (-1)^n \binom{3/2}{n} (k^2)^{-n} \mathcal{I}_n

[ericagol]

@rodluger I've added these fairly hastily; I have yet to code them up to check for mistakes. I also added in the I_0, J_0, and J_1 terms for upward iteration (in both cases).

[rodluger]

@ericagol I coded up the ksq < 1 case and it seems to work great! Here's the stability plot for the Earth's secondary eclipse using your downward recursion (left) and my old approach (right): Your approach gives 3 orders of magnitude better precision. It's still not at the machine limit -- not sure why, but it could be because the actual flux calculation involves sums and differences of the terms in s. Even if s is computed to machine precision, the downstream operations can still cause small errors.

[rodluger]

Also, I don't see the I_0, J_0 and J_1 formulae -- did you add them to the document or just to your code?

[ericagol]

Just to the document, although I’m starting to add to code. I haven’t checked for accuracy yet.

[ericagol]

One thing I found that makes a big difference is being careful in computing k_c - you can see how I did it in the code. Could it also be that some precision is lost in the computation of Y_lm or the rotation matrices?

[ericagol]

I haven't computed or coded this up yet, but do you think we should be using the same approach for H_uv as we are now using for I_v? It should be straightforward to implement, I think, since this is just the same function as I_uv, but as a function of \lambda instead of \phi.

[ericagol]

I've coded up the equations for k^2 >1 (these are still in the julia code file sn_bigr.jl), which seems to work well, even near b=0! However, there are still significant errors for large l when b is near 0. I'm wondering whether this may be due to the computation of H_uv, for which I'm still using your old formalism.

[ericagol]

It turns out that the math for Q(G_n) is pretty similar to the \mu/2 even case for P(G_n). For Q(G_n) I find the expression:

where A_{i,u,v} needs to be evaluated at \delta = -1/2 and I_v needs to be evaluated with \beta instead of \kappa, or equivalently at k^2=((b+1)^2-r^2)/(4b) for r+b > 1, and for r+b < 1, \beta = \pi/2. It is interesting that for this integral the quantity ((b+1)^-r^2)/(4b) plays the role of k^2, and both quantities are >1 for b+r <1 and are less than one for b+r >1. So, I think we can just call the functions which evaluate A_{i,u,v} and I_v, but with \delta = -1/2 and k^2=((b+1)^2-r^2)/(4b) as input arguments.

\mathcal{Q}(\mathbf{G}_n) &= 2^{l+3}\int_{-\beta/2}^{\beta/2} s_\varphi^{\tfrac{\mu}{2}+2} (1-s_\varphi^2)^{\tfrac{\mu}{4}+1}(s_\varphi^2-1/2)^{\nu/2} d\varphi\\ &= 2^{l+3} \sum_{i=0}^{\tfrac{\mu}{4}+\tfrac{\nu}{2}+1} \mathcal{A}_{i,\tfrac{\mu}{4}+1,\tfrac{\nu}{2}} \int_{-\beta/2}^{\beta/2} (s_\varphi^2)^{i+\tfrac{\mu}{4}+1} d\varphi\\ &= 2^{l+3} \sum_{i=0}^{\tfrac{\mu}{4}+\tfrac{\nu}{2}+1} \mathcal{A}_{i,\tfrac{\mu}{4}+1,\tfrac{\nu}{2}} I_{i+\tfrac{\mu}{4}+1}\\ \beta/2 &= \sin^{-1} \sqrt{\frac{(b+1)^2-r^2}{4b}}

[ericagol]

@rodluger If you extend your Earth occultation calculation up to larger l, say 20, do you continue to achieve the same precision?

[ericagol]

@rodluger It's interesting how the precision drops to ~machine mid-ingress. This could be a clue of where to gain further precision.

[ericagol]

@rodluger I just noticed that H_uv is zero for odd values of v. So, we can add the extra condition that \nu/2 must be even, in addition to \mu/2 when computing H_uv.

[rodluger]

@ericagol Here's Earth with lmax=20: no signs of instabilities! This is so cool!

[rodluger]

@ericagol I think H_{u,v} = 0 for odd v only when the occultor isn't touching the limb, right? Also, looks like we won't need H_{u,v} any more if we can compute Q directly?

[ericagol]

@rodluger Oh, you are right; that was the case I was testing.

[rodluger]

@ericagol Your math for Q(G_n) checks out, but only if I use instead of the value you mentioned.

[rodluger]

I think it may be because you identified beta/2 with kappa instead of kappa/2...

[rodluger]

Never mind, I think there's a typo in the manuscript, page 46, where you say kappa = sin^-1(k). I think it should be kappa / 2 = sin^-1(k). Ignore my previous comments if that's the case.

[rodluger]

Yep, that was it: @ericagol can you fix that typo?

[rodluger]

@ericagol So, I tested it out and the trick of using A and I to compute Q works, but it can be extremely slow to converge because your re-defined k is typically much closer to 1 than it is to zero, so the terms in the series solution for I_vmax don't start getting small until after a few hundred iterations. But perhaps you were thinking of doing the upward iterations in this case? That might fix the issue.

[rodluger]

BTW I verified that the noise floor at 10^-10 is not due to the rotation matrices. It could be due to the evaluation of the Ylms. I'll keep digging.

[rodluger]

@ericagol Update on trying out your expression for Q. The series expression for I_vmax doesn't seem to converge near ingress for l = 20 and r = 110 (my Earth secondary eclipse example), so as I said above downward recursion doesn't work. I tried to do upward recursion, starting with I_0 = 2 asin(k), but that became unstable by the time I got to l = 20. I had to fall back to my two-index H_{u,v}, which is how I plotted the figure above. Can you check whether your implementation is stable in this case? It would be nice to do away with H (as well as lambda) entirely, but it's the only method I am able to get to work currently.

[ericagol]

@rodluger Yes, that's a typo; it should be \kappa = 2 \sin^{-1} k, which is written correctly just after D43, but incorrectly after D54; sorry about that! Fixing it now.

[ericagol]

@rodluger Upward recursion seems to work fine for computing H_v (for the particular set of parameters I'm using right now). I am getting larger errors than I would like for small b still, which I need to look into; I think this may be due to computing J_v.

[ericagol]

@rodluger Hmmm... I'm still getting errors larger than I would like, which seem to be due to J_v. I think the one way to solve this might be with a tridiagonal solver: compute the J_0 and J_{v_max} terms, and then iterate upwards and downwards, transforming the algebra so that you aren't multiplying by large factors \propto k^2 when k^2 >1 or k^{-2} when k^2 < 1.

[ericagol]

@rodluger We should map out where raising/lowering the v indices gives better precision. As I mentioned, I suspect raising will work better for k^2 ~ 1, and lowering in small and large k^2 limits.

[rodluger]

Yes, that would be super useful! I can probably code that up in the next few days. I'm going to keep trying to get a stable upward recursion for Q as above, but so far I've failed. I'll keep you posted. I haven't coded up the k2>1 case but will do that next.

[rodluger]

I didn't quite understand the tridiagonal method: do you mean we simply compute J_v from the top or bottom, whichever is closest?

[ericagol]

I mean the following - write the recursion as:

Then you have a system of equations going from J_0 to J_{v_max}. Compute these, and create a matrix of the coefficients -\alpha_v, 1, \beta_v, multiplying a vector from J_1 to J_{v_max - 1}, and on the right will be a vector with the first element \alpha_2 J_0 and the last element \beta_{v_max} J_{v_max}, with zeros in between. The first row of the matrix will just have 1 and -\beta_2 and the rest of the row zeros, while the last row will be zeros with -\alpha_{v_max} and 1 at the end. You can then solve this using the tri-diagonal solver which involves a forward sweep and backward substitution. I think this can be arranged to avoid underflow/overflow.

-\alpha_v J_{v-2} + J_{v-1} - \beta_v J_v = 0\\ \alpha_v = \frac{(2v-3)k^2}{2(v+1) + 2(v-1)k^2}\\ \beta_v = \frac{2v+3}{2(v+1)+2(v-1)k^2}

[ericagol]

In the k^2 < 1 case, I'm wondering if we just iterate with I_v/k^{1+2v} and J_v/k^{1+2v}, if that would be more stable?

RL: Probably!

[rodluger]

@ericagol Can you point me to the code you use to compute Q from A(-0.5) and I(beta)? I still can't get mine to converge for large l during eclipse ingress. Here's my Earth occultation plot using your new method and l = 20 (downward iteration for I on the left, upward on the right):

The plot on the left makes sense: during ingress, k^2=((b+1)^2-r^2)/(4b) is close to 1, so the series solution is very slow to converge. But the upward recursion (right) should work, and yet it goes unstable in almost exactly the same way. I'm wondering if I'm doing something wrong, or if there's some trick to how you defined kc or something...

EDIT: I'm noticing that kc is in fact quite small in this limit, so perhaps that's it?

[ericagol]

@rodluger It's in julia/sn_bigr.jl (although the code is no longer just for large r). I have two functions: Hv_raise! and Hv_lower. These both seem to give similar precisions.

[ericagol]

@rodluger The only difference with above is how you are computing Q(G_n)?

[rodluger]

@ericagol I think I found the issue. Look at (D55): If I re-write this to solve for I_v in terms of I_{v-1}, I get which for I_1 is equal to 0.5 * I_0 - k * kc. But in your code for I (here) and H (here), it looks like you have 0.5 * I_0 - kc / k. Is the error in (D55) or in your code (or in my understanding of it)?

[rodluger]

Evaluating I_1 analytically, I get which is consistent with my solution.

[ericagol]

@rodluger Thanks, good catch, and sorry about the mistake. Does this fix the computation of Q?

[rodluger]

@ericagol It doesn't; I coded up my version from the paper, so I didn't have that bug. I suspect that once you fix it you might see the instabilities.

[ericagol]

@rodluger But you could go back to the original version of computing Q from H_{uv}?