Going back to basics: accelerating exoplanet transit modelling using Taylor-series expansion of the orbital motion
MMNRAS , 1–6 (2020) Preprint 22 September 2020 Compiled using MNRAS L A TEX style file v3.0
Going back to basics: accelerating exoplanet transit modelling usingTaylor-series expansion of the orbital motion
H. Parviainen , (cid:63) , J. Korth Instituto de Astrofísica de Canarias (IAC), E-38200 La Laguna, Tenerife, Spain Dept. Astrofísica, Universidad de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain Rheinisches Institut für Umweltforschung, Abteilung Planetenforschung an der Universität zu Köln, Universität zu Köln, Aachenerstraçe 209, 50931 Köln, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
A significant fraction of an exoplanet transit model evaluation time is spent calculating projected distances between the planet andits host star. This is a relatively fast operation for a circular orbit, but slower for an eccentric one. However, because the planet’sposition and its time derivatives are constant for any specific point in orbital phase, the projected distance can be calculatedrapidly and accurately in the vicinity of the transit by expanding the planet’s x and y positions in the sky plane into a Taylor seriesat mid-transit. Calculating the projected distance for an elliptical orbit using the four first time derivatives of the position vector(velocity, acceleration, jerk, and snap) is ∼
100 times faster than calculating it using the Newton’s method, and also significantlyfaster than calculating z for a circular orbit because the approach does not use numerically expensive trigonometric functions.The speed gain in the projected distance calculation leads to 2-25 times faster transit model evaluation speed, depending on thetransit model complexity and orbital eccentricity. Calculation of the four position derivatives using numerical differentiationtakes ∼ µ s with a modern laptop and needs to be done only once for a given orbit, and the maximum error the approximationintroduces to a transit light curve is below 1 ppm for the major part of the physically plausible orbital parameter space. Key words:
Methods: numerical – Techniques: photometric – Planets and satellites
An exoplanet transit model aims to reproduce the photometric signalcaused by a planet crossing over the limb-darkened disk of its host star(Mandel & Agol 2002; Seager & Mallen-Ornelas 2003; Winn 2010).Evaluation of the transit model can generally be divided into twoparts: a) calculation of the projected planet-star centre distance, z ;and b) calculation of the flux decrement caused by a planet occludinga part of the stellar disk visible to the observer.The main focus in transit model development has been on thesecond part, but the calculation of projected distances can actuallytake a significant fraction of the total model evaluation time. Thestandard approach for calculating z for a single point in time requiresseveral ( ≈
6) trigonometric function calls and solving the Kepler’sequation numerically. While worrying about the computational costof using trigonometric functions might seem frivolous, z needs tobe calculated at least once for each photometric data point whenevaluating an exoplanet transit model, and multiple times if the modelneeds to be supersampled (such as for Kepler and
TESS long cadencelight curves, Kipping 2010). Further, it is already common to havephotometric data sets of tens or hundreds of thousands of data points(such as a four-year
Kepler light curve), and a transit light curveanalysis consisting of a posterior optimisation and Markov ChainMonte Carlo (MCMC) sampling steps can require the model to beevaluated a large number ( ∼ ) of times over all the data points. (cid:63) E-mail: [email protected]
Thus, while calculating z using the standard approaches is a trivialmatter for small data sets, speeding up the calculation has a poten-tial to yield significant real-life performance gains when modellingmodern data sets. While accuracy is more important than speed fora scientific code, a speed increase without any significant sacrificesin accuracy gives freedom for exploratory analyses and experimen-tation, which can lead to new interesting discoveries, or, at least,increase the reliability of our analyses.In this short paper we show how a very simple change in thecomputation of z can lead to a significant speed-up of a transit modelwithout sacrificing model accuracy. The approach is based on high-school level mathematics (Taylor series expansion, a tool that hasbeen used in astronomy and astrophysics for centuries, especially inthe research of eclipsing binaries) and has been tested thoroughly. Theapproach still requires the ability to calculate the eccentric anomalyto a high precision in order to calculate the position derivativesusing numerical differentiation, but this needs to be done only for asmall number of points in time (seven in our implementation) for asingle Keplerian orbit, rather than calculating it for each datapointseparately.We provide an example Python implementation of the methodin Appendix A, and the approach has been adopted as the main z computation method in the PyTransit (Parviainen 2015) transitmodelling package. https://github.com/hpparvi/PyTransit © 2020 The Authors a r X i v : . [ a s t r o - ph . E P ] S e p H. Parviainen et al. b T b T c zy x R Figure 1.
Orbit of a transiting short-period exoplanet on an eccentric orbit.The figure shows the projected planet-star centre distance ( z ), the impactparameter ( b ), the stellar radius ( R (cid:63) ), transit centre time ( T c ), and the timeof minimum projected distance ( T b ). The two latter are equal for a circularorbit, but generally differ slightly for an eccentric orbit. Vertical dotted lineshows x = y = Calculation of the projected planet-star separation ( z , see Fig. 1) asa function of time is a necessary step for exoplanet transit modelevaluation. The standard approach for calculating z for a genericeccentric orbit requires the calculation of the eccentric anomaly fromthe mean anomaly, which requires us to solve Kepler’s equation, forwhich no closed-form solutions exist. Thus, the Kepler’s equationneeds to be solved using numerical methods, such as iteration or theNewton’s method. After the eccentric anomaly has been solved, thecomputation of z still requires six trigonometric function calls, whichare relatively expensive operations.Could there be a way to calculate z without the need to solveKepler’s equation or use trigonometric functions? Planet’s x and y positions in the sky-plane draw smooth and well-behaved curves asa function of time, as shown in Fig. 2. The x position is a monotoni-cally increasing function of time near the transit, and the y position isa smooth unimodal function with a single minimum near the transit(this for non-zero impact parameter since y is constant for b = l ( t ) = (cid:213) n = l ( n ) ( t ) n ! ( t − t ) n , (1)where l is the position (either x or y ), t is the point around whichthe Taylor series is expanded, l ( n ) is the nth derivative of l evaluatedat point t , and n ! is the factorial of n . Mid-transit time where x = t (although other possibilities exists, such asthe time of minimum projected distance or the time of minimum y position), after which we only need to select n to ensure a sufficientaccuracy so that the approximation does not affect the transit modelin any significant fashion. After testing the accuracy of different n (see discussion about accuracy later in Sect. 4), we chose to use thefour first time derivatives of position: velocity, acceleration, jerk, andsnap.The first step is to calculate the planet’s position in the sky-planeat mid-transit time and its four time derivatives. For a circular orbit,the position at mid-transit is [ , − b ] , where b is the impact parameter.However, for an eccentric orbit the y position differs from − b .We use a seven-point central finite difference method (Fornberg1988) to calculate the derivatives. This requires us to calculate thepositions at seven uniformly spaced times centred around the mid-transit time. The calculation of these locations requires using a stan- dard accurate method for evaluating Keplerian orbits, but this needsto be done only once for a given orbit.The velocity, acceleration, jerk, and snap vector component d (either x or y ) can be computed given the points (cid:174) l i = (cid:174) l ( t − ih ) where i = [− , − , − , , , , ] and h is the time step, as v d = − d − + d − − d − + d − d + d h , (2) a d = d − − d − + d − − d + d − d + d h , (3) j d = d − − d − + d − − d + d − d h , (4) s d = − d − + d − − d − + d − d + d − d h . (5)After the derivatives have been calculated, the projected distance canbe computed for time t c by first calculating the time difference to thenearest transit centre, t , E = (cid:22) t c − t + . pp (cid:23) , (6) t = t c − ( t + E p ) , (7)where E is the epoch, t the mid-transit time, (cid:98)(cid:99) denotes the flooroperation, and p the orbital period, and then evaluating the Taylorseries at t as (cid:174) l = (cid:174) l + (cid:174) v t + (cid:174) at + (cid:174) jt + (cid:174) st , (8) z = |(cid:174) l | (9)where (cid:174) l is the position vector at mid-transit, and (cid:174) v , (cid:174) a , (cid:174) j , and (cid:174) s arethe velocity, acceleration, jerk, and snap vectors, respectively.The approximation requires that the planet’s position is evaluatedat seven points in time for an orbit that does not evolve in time (thatis, the orbital parameters do not evolve in time). However, if the orbitis perturbed by external forces, such as other massive bodies in amultiplanet system, the series terms need to be calculated separatelyfor each transit. This leads to a photodynamical model where theterms are calculated using a set of positions calculated with an n-bodyintegrator, as done in PyTTV by Korth et al. (2020, in preparation).The derivatives known, the computation of z requires only mul-tiplications, summation, and a single square-root operation. Giventhe simplicity of the approximation, we provide an example Pythonimplementation in Appendix A. The real-world improvement in the transit model evaluation speeddepends on how heavy the transit shape model is relative to the z calculation method (that is, how large fraction of the transit modelexecution time is spent on computing the orbit). For the transit modelassuming quadratic stellar limb darkening by Mandel & Agol (2002),the speed gain is between 6 (eccentric orbit calculated using theNewton’s method) and 2 (circular orbit), that is, the model is 6 timesfaster to evaluate for an eccentric orbit when z is calculated using aTaylor series expansion rather than Newton’s method. For the mostsimple transit shape model that assumes uniform stellar disk, thespeed gain is between 24 (eccentric orbit) and 2 (circular orbit). In We group the expressions slightly differently for the actual implemen-tation to reduce sensitivity to floating point round-off errors, as shown inAppendix A.MNRAS000
Orbit of a transiting short-period exoplanet on an eccentric orbit.The figure shows the projected planet-star centre distance ( z ), the impactparameter ( b ), the stellar radius ( R (cid:63) ), transit centre time ( T c ), and the timeof minimum projected distance ( T b ). The two latter are equal for a circularorbit, but generally differ slightly for an eccentric orbit. Vertical dotted lineshows x = y = Calculation of the projected planet-star separation ( z , see Fig. 1) asa function of time is a necessary step for exoplanet transit modelevaluation. The standard approach for calculating z for a genericeccentric orbit requires the calculation of the eccentric anomaly fromthe mean anomaly, which requires us to solve Kepler’s equation, forwhich no closed-form solutions exist. Thus, the Kepler’s equationneeds to be solved using numerical methods, such as iteration or theNewton’s method. After the eccentric anomaly has been solved, thecomputation of z still requires six trigonometric function calls, whichare relatively expensive operations.Could there be a way to calculate z without the need to solveKepler’s equation or use trigonometric functions? Planet’s x and y positions in the sky-plane draw smooth and well-behaved curves asa function of time, as shown in Fig. 2. The x position is a monotoni-cally increasing function of time near the transit, and the y position isa smooth unimodal function with a single minimum near the transit(this for non-zero impact parameter since y is constant for b = l ( t ) = (cid:213) n = l ( n ) ( t ) n ! ( t − t ) n , (1)where l is the position (either x or y ), t is the point around whichthe Taylor series is expanded, l ( n ) is the nth derivative of l evaluatedat point t , and n ! is the factorial of n . Mid-transit time where x = t (although other possibilities exists, such asthe time of minimum projected distance or the time of minimum y position), after which we only need to select n to ensure a sufficientaccuracy so that the approximation does not affect the transit modelin any significant fashion. After testing the accuracy of different n (see discussion about accuracy later in Sect. 4), we chose to use thefour first time derivatives of position: velocity, acceleration, jerk, andsnap.The first step is to calculate the planet’s position in the sky-planeat mid-transit time and its four time derivatives. For a circular orbit,the position at mid-transit is [ , − b ] , where b is the impact parameter.However, for an eccentric orbit the y position differs from − b .We use a seven-point central finite difference method (Fornberg1988) to calculate the derivatives. This requires us to calculate thepositions at seven uniformly spaced times centred around the mid-transit time. The calculation of these locations requires using a stan- dard accurate method for evaluating Keplerian orbits, but this needsto be done only once for a given orbit.The velocity, acceleration, jerk, and snap vector component d (either x or y ) can be computed given the points (cid:174) l i = (cid:174) l ( t − ih ) where i = [− , − , − , , , , ] and h is the time step, as v d = − d − + d − − d − + d − d + d h , (2) a d = d − − d − + d − − d + d − d + d h , (3) j d = d − − d − + d − − d + d − d h , (4) s d = − d − + d − − d − + d − d + d − d h . (5)After the derivatives have been calculated, the projected distance canbe computed for time t c by first calculating the time difference to thenearest transit centre, t , E = (cid:22) t c − t + . pp (cid:23) , (6) t = t c − ( t + E p ) , (7)where E is the epoch, t the mid-transit time, (cid:98)(cid:99) denotes the flooroperation, and p the orbital period, and then evaluating the Taylorseries at t as (cid:174) l = (cid:174) l + (cid:174) v t + (cid:174) at + (cid:174) jt + (cid:174) st , (8) z = |(cid:174) l | (9)where (cid:174) l is the position vector at mid-transit, and (cid:174) v , (cid:174) a , (cid:174) j , and (cid:174) s arethe velocity, acceleration, jerk, and snap vectors, respectively.The approximation requires that the planet’s position is evaluatedat seven points in time for an orbit that does not evolve in time (thatis, the orbital parameters do not evolve in time). However, if the orbitis perturbed by external forces, such as other massive bodies in amultiplanet system, the series terms need to be calculated separatelyfor each transit. This leads to a photodynamical model where theterms are calculated using a set of positions calculated with an n-bodyintegrator, as done in PyTTV by Korth et al. (2020, in preparation).The derivatives known, the computation of z requires only mul-tiplications, summation, and a single square-root operation. Giventhe simplicity of the approximation, we provide an example Pythonimplementation in Appendix A. The real-world improvement in the transit model evaluation speeddepends on how heavy the transit shape model is relative to the z calculation method (that is, how large fraction of the transit modelexecution time is spent on computing the orbit). For the transit modelassuming quadratic stellar limb darkening by Mandel & Agol (2002),the speed gain is between 6 (eccentric orbit calculated using theNewton’s method) and 2 (circular orbit), that is, the model is 6 timesfaster to evaluate for an eccentric orbit when z is calculated using aTaylor series expansion rather than Newton’s method. For the mostsimple transit shape model that assumes uniform stellar disk, thespeed gain is between 24 (eccentric orbit) and 2 (circular orbit). In We group the expressions slightly differently for the actual implemen-tation to reduce sensitivity to floating point round-off errors, as shown inAppendix A.MNRAS000 , 1–6 (2020) ccelerating exoplanet transit modelling using Taylor series X [ R ] e = 0.00, d min = 4.0 R e = 0.25, d min = 3.0 R e = 0.50, d min = 2.0 R Y [ R ] Time [d] Z [ R ] Time [d]
Time [d]
Figure 2.
The exact (solid black line) and approximate (dashed black line) sky-plane x and y values and the projected distance z for three short-period orbits withdifferent eccentricities. The vertical lines mark the beginning and the end of a transit. The orbits have a common period (1 d), semi-major axis (4 R (cid:63) ), impactparameter (0.5), and argument of periastron (0). The minimum planet-star separation, d min , tells the separation between the planet and the star at periastron. both cases, the minimum speed gain is around 2 (that is, the modelis at least twice as fast to calculate). While the Taylor series approximation of z is significantly computa-tionally faster than the other approaches for calculating z , its practicalusability depends on the error caused to the exoplanet transit model.The accuracy of the approximation depends on the three-dimensionalcurvature of the orbit at the mid transit time, what again depends onthe semi-major axis, eccentricity, and argument of periastron.Figure 2 shows the actual and approximated x and y coordinatesand the projected distance z for three increasingly eccentric short-period orbits with orbital period, p , of 1 d, scaled semi-major axis, a ,of 4 R (cid:63) , and impact parameter, b , of 0.5. Figure 3 shows the maximumabsolute errors in a transit light curve caused by the approximationfor a circular orbit as a function of the planet-star separation (that is,the semi-major axis) at mid-transit (upper panel) and orbital period(lower panel). The orbits correspond to three planets with radiusratio of 0.15, 0.1, and 0.05 orbiting a star with a stellar density, ρ (cid:63) ,of 1 . − with an impact parameter of 0.5. The figure focuses onthe ultra-short-period and short-period regime because the error isbelow 1 ppm for semi-major axes larger than 5 R (cid:63) . Considering thecurrently known transiting exoplanets, the maximum absolute errorintroduced by the approximation would be ∼
10 ppm, staying below1 ppm for all but the most extreme ultra-short-period planets.Figure 4 shows the maximum absolute errors in transit light curvescaused by the approximation for four sets of orbital parameters as afunction of increasing eccentricity. The scenarios are: a) p = . a = . R (cid:63) , b) p = a = R (cid:63) , c) p =
15 d and a = R (cid:63) ,and d) p =
30 d and a = . R (cid:63) , and the eccentricities cover therange of eccentricities for known planets with periods less or equal to the scenario period. The maximum error for most of the physicallyplausible orbits is below 1 ppm, and still below 10 ppm for the veryeccentric orbits. A planet’s normalised planet-star centre distance near a transit (ora secondary eclipse) can be calculated using the planet’s sky-planeposition at mid-transit time and its four first time derivatives for thewhole physically plausible orbital parameter space without sacrific-ing transit model accuracy. The approach is ∼
100 times faster tocalculate than an approach using Newton’s method to solve the Ke-pler’s equation, and yields a 2-24 gain in transit model evaluationspeed. Further, since the approach is based on expanding the sky-plane position, the position can be used directly with transit modelsthat break the radial symmetry, such as the gravity-darkened transitmodel for rapidly rotating stars by Barnes (2009). A gravity-darkenedmodel utilising the approach to compute the ( x , y ) position has beenadded to a coming PyTransit version (v2.4), but here the speed gainover the standard approach is relatively small due to computationalcost of the transit model itself.As clear from Figs. 3 and 4, the errors introduced by the approxi-mation into the transit model are negligible. The absolute maximumerror is below 1 ppm for all but the shortest orbital periods and highesteccentricities, and generally below 10 ppm for any currently knownplanets.We could also expand z directly into a Taylor series instead of thesky-plane x and y positions. However, the projected distance has arelatively sharp minimum (compared to the behaviour of x and y positions), and the time of the minimum does not necessarily matchour mid-transit time for which x =
0. Thus, expanding z would MNRAS , 1–6 (2020)
H. Parviainen et al.
Planet-star separation [R ] M a x e rr o r [ pp m ] k=0.15k=0.10k=0.05 Period [d] M a x e rr o r [ pp m ] Figure 3.
Maximum absolute error to a transit light curve introduced bythe approximation for a circular orbit. The maximum error is shown for threeplanet sizes as a function of the planet-star separation (upper panel) and period(lower panel) assuming a stellar density of 1 . − , impact parameter of0.5, and quadratic limb darkening with coefficients ( u = . , v = . ) . eccentricity M a x e rr o r [ pp m ] a b c d Figure 4.
Maximum absolute error to a transit light curve introduced bythe approximation for a circular orbit for four different scenarios averagedover 2 × samples in argument of periastron and impact parameter. Thescenarios a, b, c, and d are described in Sect. 4 require one to first find the minimum z time and then include higher-order derivatives into the series. This increases complexity of theimplementation and would also likely reduce numerical stability, sowe decided to prefer the approach described here.The approach naturally works when modelling transits (oreclipses) only, and the full Keplerian orbit needs to be evaluatedwhen modelling phase curves. However, even then it may be bene-ficial to calculate the projected distances for the transit model usingthe Taylor series approach, especially if the transit model needs to besupersampled.The planet-star contact points (beginning of ingress, T , end ofingress, T , beginning of egress, T , and end of egress T , Winn 2010)are easy to compute numerically. Calculation of a single point takes ≈
500 ns, and the calculation of different durations ( T , T , T , and T ) takes between 1-2 µ s. We do not include the code to calculate thecontact points here, but make it available from PyTransit repositoryin GitHub. PyTransit also uses the T and T points to create atransit bounding box in time that is used to ensure we do not wastetime evaluating the model over the out-of-transit points.The centre time for the series expansion, t , affects the accuracy. Itcould be beneficial to choose t to match the time where y is minimum(so that y velocity is zero), or the time of minimum z . It could also bepossible to choose a different t for the x and y expansion. However,both approaches would require more computation to solve thoselocations than just choosing the mid-transit time, and are probablynot worth the work considering that the current approach alreadyreaches an accuracy that has basically no effect on the transit lightcurve model.The speed gains discussed in Sect. 3 depend significantly on theoverall implementation of the whole transit model. The examplesin this study consider light curves where most of the points are intransit (that is, most of the out-of-transit data has been removed).Having a light curve with a small fraction of in-transit points (suchas when modelling a full Kepler or TESS light curve directly) willsignificantly increase the speed gain unless the transit model is smartenough to skip the out-of-transit points.The final effect on the evaluation speed also depends on the otherparts of the posterior computation, such as the noise model. The gainwill be smaller when the posterior computation time is dominated bythe noise model evaluation (such as when using brute-force GaussianProcesses), and greatest in an analysis with a computationally cheapnoise model and a large number of data points.The approach has been adopted as the main z computation methodin the PyTransit exoplanet transit modelling package by Parviainen(2015). However, considering the simplicity of the approach, webelieve it can be useful to everyone developing exoplanet transitmodels and modelling frameworks independent of the programminglanguage used. Thus, the approach can be easily added to othercommonly used transit modelling packages, such as EXOFAST byEastman et al. (2013), batman by Kreidberg (2015), ellc by Maxted(2016), or TLCM by Csizmadia (2020). ACKNOWLEDGEMENTS
We thank E. Agol for his valuable comments that substantially helpedto improve the manuscript. HP acknowledges financial support fromthe Agencia Estatal de InvestigaciÃşn del Ministerio de Ciencia,InnovaciÃşn y Universidades (MICIU) and UniÃşn Europea FondosFEDER (EU FEDER) funds through the project PGC2018-098153-B-C31. JK acknowledges support by DFG grants PA525/19-1 and
MNRAS000
MNRAS000 , 1–6 (2020) ccelerating exoplanet transit modelling using Taylor series PA525/18-1 within the DFG Schwerpunkt SPP 1992, Exploring theDiversity of Extrasolar Planets.
DATA AVAILABILITY
There are no new data associated with this article.
REFERENCES
Barnes J., 2009, ApJ, 705, 683Csizmadia S., 2020, MNRAS, 496, 4442Eastman J., Gaudi B. S., Agol E., 2013, Publ. Astron. Soc. Pacific, 125, 83Fornberg B., 1988, Math. Comput., 51, 699Kipping D. M., 2010, MNRAS, pp no–noKreidberg L., 2015, ArXiv, 1507.08285Mandel K., Agol E., 2002, ApJ, 580, L171Maxted P. F., 2016, A&A, 591, 1Parviainen H., 2015, MNRAS, 450, 3233Seager S., Mallen-Ornelas G., 2003, ApJ, 585, 1038Winn J. N., 2010, in Seager S., ed., , EXOPLANETS. University of ArizonaPress, Tucson, AZ, Chapt. Transits a ( arXiv:1001.2010 ), http://arxiv.org/abs/1001.2010 APPENDIX A: EXAMPLE IMPLEMENTATION
Here we show an example numba -accelerated
Python implemen-tation of the approach used by the
PyTransit transit modellingpackage. First, a method to calculate the sky-plane x and y deriva-tives from numba import njit from numpy import (arctan2 , cos , sin , sqrt , floor ,mod , pi)@njit def ta_newton_s(t, t0, p, e, w):offset = arctan2(sqrt(1.0-e**2)*sin(0.5*pi-w),e + cos(0.5*pi - w))offset -= e*sin(offset)ma = mod(2*pi*(t-(t0-offset*p/2*pi))/p, 2*pi)ea = maerr = 0.05k = 0 while abs (err) > 1e-8 and k < 1000:err = ea - e*sin(ea) - maea = ea - err/(1.0-e*cos(ea))k += 1sta = sqrt(1.0-e**2)*sin(ea)/(1.0-e*cos(ea))cta = (cos(ea)-e)/(1.0-e*cos(ea)) return arctan2(sta , cta)@njit def xyeo(t, t0, p, e, w, ae, ci):f = ta_newton_s(t, t0, p, e, w)r = ae / (1.+ e*cos(f))x = -r * cos(w + f)y = -r * sin(w + f) * ci return x, y@njit def vajs_from_paiew(t0, p, a, i, e, w): """Planet velocity , acceleration , jerk , andsnap at mid -transit in [R_star / day]""" dt = 2e-2 ae = a*(1.-e**2)ci = cos(i)x0, y0 = xyeo(t0 -3*dt, t0, p, e, w, ae, ci)x1, y1 = xyeo(t0 -2*dt, t0, p, e, w, ae, ci)x2, y2 = xyeo(t0 -1*dt, t0, p, e, w, ae, ci)x3, y3 = xyeo(t0 , t0, p, e, w, ae, ci)x4, y4 = xyeo(t0+1*dt, t0, p, e, w, ae, ci)x5, y5 = xyeo(t0+2*dt, t0, p, e, w, ae, ci)x6, y6 = xyeo(t0+3*dt, t0, p, e, w, ae, ci) a, b, c = 1/60, 9/60, 45/60vx=(a*(x6-x0)+b*(x1-x5)+c*(x4-x2))/dtvy=(a*(y6-y0)+b*(y1-y5)+c*(y4-y2))/dt a, b, c, d = 1/90, 3/20, 3/2, 49/18ax=(a*(x0+x6)-b*(x1+x5)+c*(x2+x4)-d*x3)/dt**2ay=(a*(y0+y6)-b*(y1+y5)+c*(y2+y4)-d*y3)/dt**2 a, b, c = 1/8, 1, 13/8jx=(a*(x0-x6)+b*(x5-x1)+c*(x2-x4))/dt**3jy=(a*(y0-y6)+b*(y5-y1)+c*(y2-y4))/dt**3 a, b, c, d = 1/6, 2, 13/2, 28/3sx=(-a*(x0+x6)+b*(x1+x5)-c*(x2+x4)+d*x3)/dt**4sy=(-a*(y0+y6)+b*(y1+y5)-c*(y2+y4)+d*y3)/dt**4 return y3, vx, vy, ax, ay, jx, jy, sx, sy
Here xyeo calculates the x and y positions at given times usingthe Newton’s method to calculate the true anomaly ( ta_newton_s ).Now, the projected distance can be calculated using the derivativesas a Taylor series @njit(fastmath=True) def z_taylor(tc, t0, p, y0,vx, vy, ax, ay,jx, jy, sx, sy): """Projected planet -star distance using aTaylor series expansion.""" epoch = floor((tc - t0 + 0.5*p) / p)t = tc - (t0 + epoch * p)t2 = t*t
MNRAS , 1–6 (2020)
H. Parviainen et al. t3 = t*t2t4 = t*t3px = vx*t + ax*t2/2 + jx*t3/6 + sx*t4/24py = y0 + vy*t + ay*t2/2 + jy*t3/6 + sy*t4/24 return sqrt(px**2 + py**2)
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000