TableOfContents()

#=
$(RobustLocalResource("https://cfn-live-content-bucket-iop-org.s3.amazonaws.com/journals/1538-3881/153/5/224/revision1/ajaa6897f1_hr.jpg","../_assets/week3/ajaa6897f1_hr.jpg",:width=>"80%"))
Credit: 

$(RobustLocalResource("https://cfn-live-content-bucket-iop-org.s3.amazonaws.com/journals/0004-637X/743/2/200/revision1/apj408732f14_hr.jpg","../_assets/week3/apj408732f14_hr.jpg",:width=>"80%"))
Credit: 

$(RobustLocalResource("https://cfn-live-content-bucket-iop-org.s3.amazonaws.com/journals/1538-3881/153/5/224/revision1/ajaa6897t5_lr.gif","../_assets/week3/ajaa6897t5_lr.gif",:width=>"80%"))

=#
md"""
**A:** Yes, e.g., Kepler-19c 
See Fig. 1 of [Malavolta et al. (2017)](https://doi.org/10.3847/1538-3881/aa6897)

See Fig 14. of [Balllard et al. (2011)](https://doi.org/10.3847/1538-3881/aa6897)

See Table 6 of [Malavolta et al. (2017)](https://doi.org/10.3847/1538-3881/aa6897)
"""

let
  plt_in = plot(legend=:none)
  plt_out = plot(legend=:none)
  scatter!(plt_in,tt_in,ttvs_in)
  ylabel!(plt_in,"TTVᵢₙ (d)")
  scatter!(plt_out,tt_out,ttvs_out)
  ylabel!(plt_out,"TTVₒᵤₜ (d)")
  xlabel!(plt_out,"Time (d)")
  xlims!(plt_in,0,timespan)
  xlims!(plt_out,0,timespan)
  plot(plt_in, plt_out, layout=(2,1), link=:x, size=(1000,600), leftmargin=20px, bottommargin=18px, guidefontsize=18 )	
end

if false
md"""
### One transit
- Transit Depth: $\delta$ (dimensionless)
- Impact Parameter: $b$ (units of stellar radii)
- Ingress Duration: $\tau$
- Transit Durations:
   - Total duration: $t_{IV}-t_{I}$
   - Full-transit duration: $t_{III}-t_{II}$
   - Mathematically-convenient duration: $T$
   - Best-measured duration: $(t_1+t_2-t_3-t_4)/2$
"""
end

# $(Resource("https://www.aanda.org/articles/aa/full_html/2018/09/aa33085-18/aa33085-18-fig1.jpg", :width=>"80%"))
#Credit: 
md"""
See Fig 1. of [Rodenbeck et al. (2018)](https://doi.org/10.1051/0004-6361/201833085)
"""

#$(Resource("https://www.aanda.org/articles/aa/full_html/2018/09/aa33085-18/aa33085-18-fig2.jpg", :width=>"80%"))
#Credit:
md"""
See Fig. 2 of [Rodenbeck et al. (2018)](https://doi.org/10.1051/0004-6361/201833085)
"""

ChooseDisplayMode()

Full Width Mode       Present Mode

begin
    using PlutoUI, PlutoTeachingTools, HypertextLiteral
    using Plots, LaTeXStrings
    using Plots.PlotMeasures
end

nbsp = htl" "

 

mearth=3e-6

3.0e-6

ttvs_in, ttvs_out, tt_in, tt_out = test_ttv(6,ceil(Int64,(timespan-t0in)/Pin),ceil(Int64,(timespan-t0out)/Pout),[min*mearth,Pin,t0in,ein,ωin, mout*mearth,Pout,t0out,eout,ωout])

(ttv_in = [0.0, -0.000261540136363074, -0.0005231189372340493, -0.000783587050550149, -0.0010441513005326468, -0.0013025381454352494, -0.0015610193352999582, -0.001816343708617378, -0.002071662878941905, -0.002322975206513561  …  0.0043464216838199005, 0.004500483531700491, 0.004780498534725676, 0.004927160223065881, 0.005196166342563159, 0.005334757796921942, 0.005591809972496657, 0.005721691738533963, 0.005965889353426056, 0.0060864546301831685], ttv_out = [0.0, 0.0003297990549240407, 0.0006582266761751396, 0.0009839195417966154, 0.0013055303235018022, 0.0016217351304131118, 0.0019312403462597835, 0.002232788745955712, 0.002525164839579767, 0.0028071994546337627  …  -0.0012772267961582285, -0.001587332426770693, -0.0018913490070818655, -0.002188104194005995, -0.002476450922305607, -0.0027552718540606288, -0.0030234837972269142, -0.003280042025930563, -0.003523944436163147, -0.003754235480328827], tt_in = [0.0, 9.999738459863638, 19.999476881062765, 29.99921641294945, 39.99895584869947, 49.99869746185456, 59.9984389806647, 69.99818365629139, 79.99792833712105, 89.99767702479349  …  1180.0043464216837, 1190.0045004835317, 1200.0047804985347, 1210.0049271602231, 1220.0051961663426, 1230.0053347577968, 1240.0055918099724, 1250.0057216917385, 1260.0059658893533, 1270.0060864546301], tt_out = [0.0, 20.100329799054926, 40.200658226676175, 60.30098391954179, 80.40130553032351, 100.50162173513041, 120.60193124034626, 140.70223278874593, 160.80252516483958, 180.90280719945463  …  1085.3987227732039, 1105.4984126675731, 1125.5981086509928, 1145.6978118958061, 1165.7975235490776, 1185.897244728146, 1205.9969765162027, 1226.096719957974, 1246.1964760555638, 1266.2962457645197])

function test_ttv(jmax::Integer,n1::Integer,n2::Integer,data::Vector)
    @assert(jmax>=1)  # Should there be a larger minimum?
    @assert(n1>2)
    @assert(n2>2)
    @assert(length(data)==10)
    # Performs a test of the transit_times.jl routine
    # Set up planets planar-planet types for the inner and outer planets:
    p1=TTVFaster.Planet_plane_hk(data[1],data[2],data[3],data[4],data[ 5])
    p2=TTVFaster.Planet_plane_hk(data[6],data[7],data[8],data[9],data[10])
    time1 = collect(p1.trans0 .+ range(0, step=1, length=n1) .* p1.period)
    time2 = collect(p2.trans0 .+ range(0, step=1, length=n2) .* p2.period)
    alpha0=(p1.period/p2.period)^(2//3)
    # Initialize the computation of the Laplace coefficients:
    b0=TTVFaster.LaplaceCoefficients.initialize(jmax+1,alpha0)
    # Define arrays to hold the TTVs:
    ttv_el_type = eltype(data) == Float64 ? Float64 : Number
    ttv1=Array{ttv_el_type}(undef,n1)
    ttv2=Array{ttv_el_type}(undef,n2)
    # Define arrays to hold the TTV coefficients and Laplace coefficients:
    f1=Array{Float64}(undef,jmax+2,5)
    f2=Array{Float64}(undef,jmax+2,5)
    b=Array{Float64}(undef,jmax+2,3)
    TTVFaster.compute_ttv!(jmax,p1,p2,time1,time2,ttv1,ttv2,f1,f2,b,alpha0,b0)
    time1 .+= ttv1
    time2 .+= ttv2
    return (;ttv_in=ttv1,ttv_out=ttv2, tt_in=time1, tt_out=time2)
end

test_ttv (generic function with 1 method)

# WARNING:  This code is very old (written for Julia v0.3) and should not be used as an example for writing your own code
# Computes transit timing variations to linear order in eccentricity.  
# Please cite Agol & Deck (2015) if you make use of this in published research.
# https://github.com/ericagol/TTVFaster
module TTVFaster

export Planet_plane_hk, compute_ttv!

struct Planet_plane
# Parameters of a planet in a plane-parallel system
  # Mass ratio of the planet to the star:
  mass_ratio :: Float64
  # Initial time of transit:
  period   :: Float64
  trans0   :: Float64
  eccen    :: Float64
# longitude of periastron measured from line of sight, in radians:
  omega    :: Float64
end


struct Planet_plane_hk{T<:Number} # Parameters of a planet in a plane-parallel system
  # Mass ratio of the planet to the star:
  mass_ratio :: T
  # Initial time of transit:
  period   :: T
  trans0   :: T
  # e times cos or sin of longitude of periastron measured from line of sight, in radians:
  ecosw    :: T
  esinw    :: T
end

 
module LaplaceCoefficients

export laplace_coefficients_initialize
export laplace_wisdom

"""
#/* compute Laplace coefficients and Leverrier derivative/
#          j
#     j   d     i
#    a   ---   b (a)
#          j    s
#        da
#
#   by series summation */
#/* Code due to Jack Wisdom */
"""
function laplace_wisdom(s::Rational,i::Integer,j::Integer,a::Number)
# function laplace_wisdom,s,i,j,a  IDL
# double laplace(double s, int i, int j, double a);  c

##define LAPLACE_EPS 1.0e-12
LAPLACE_EPS = convert(eltype(a),1.0e-12)

#if (i lt 0) then i = -i
i=abs(i)

if j <= i     #/* compute first term in sum */
  factor4 = one(a)
  for k=0:j-1
    factor4 *= i - k
  end
  lap_coef_sum = factor4
  q0 = 0
else
  q0 = fld(j + 1 - i,2)
  lap_coef_sum = zero(a)
  factor4 = one(a)
end

#  /* compute factors for terms in lap_coef_sum */

factor1 = s
factor2 = s + i
factor3 = i + 1
for q=1:q0-1       #/* no contribution for q = 0 */
  factor1 *= s + q
  factor2 *= s + i + q
  factor3 *= i+q+1
end
if q0 > 1
  q=q0
else
  q=1
end
#println(j+1-i,q0)
term = a*a * factor1 * factor2 / (factor3 * q)

#  /* sum series */

while  (term*factor4) > LAPLACE_EPS
  factor4 = one(a)
  for k=0:j-1
    factor4 *= (2*q + i - k)
  end
  lap_coef_sum += term * factor4
  factor1 += 1
  factor2 += 1
  factor3 += 1
  q = q+1
  term *= a*a * factor1 * factor2 / (factor3 * q)
end

#  /* fix coefficient */

for k=0:i-1
  lap_coef_sum *= (s+k)/(k+1)
end

apower = (q0 <= 0) ?  i : 2*q0 + i - 2
lap_coef_sum *= 2 * a^apower
# Return the Laplace Coefficient:
return lap_coef_sum
end

"""
# This computes the Laplace coefficients via recursion.
"""
function initialize(jmax::Integer,alpha::Number)
  nmax=7
  b0=Array{eltype(alpha)}(undef,nmax,jmax+1) # Array to hold the coefficients
  # Compute the highest two Laplace coefficients using Wisdom's series approach:
  for j=0:jmax
    for i=0:nmax-1
      b0[i+1,j+1]=laplace_wisdom(1//2,j,i,alpha)/alpha^i
    end
  end
  return b0
end

end # module

laplace_coefficients_initialize(jmax::Integer,alpha::Number) = LaplaceCoefficients.initialize(jmax,alpha)

# Computes TTV coefficients for first-order eccentricity
# solution from Agol & Deck (2015).  Please cite this paper
# if you make use of this in your research.

u(gamma::T,c1::T,c2::T) where { T<:Number} = ((3+gamma*gamma)*c1+2*gamma*c2)/(gamma*gamma*(1-gamma*gamma))
# m=+/-1
v(z::T,d1::T,d2::T,m::Integer) where T<:Number = ((m*(1-z*z)+6*z)*d1+(2+z*z)*d2)/(z*(1-z*z)*(z+m)*(z+2*m))

function ttv_succinct!(jmax::Integer,alpha::Number,f1::Array{Float64,2},f2::Array{Float64,2},b::Array{Float64,2},alpha0::Number,b0::Array{Float64,2})

# See simple_solution.pdf 7/16/2015

# Fourth-order Taylor expansion approximation of Laplace coefficients:
dalpha = alpha-alpha0
for i=0:2
  for j=0:jmax
    b[j+1,i+1]=b0[i+1,j+1]+dalpha*(b0[i+2,j+1]+0.5*dalpha*(b0[i+3,j+1]+dalpha/3.0*(b0[i+4,j+1]+dalpha*0.25*b0[i+5,j+1])))
  end
end

sqrtalpha = sqrt(alpha)

# Loop over j:
@inbounds for j=0:jmax
  # \delta_{j1} (this is indirect coefficient which is only needed for j=1)
  dj1 = j==1 ? 1.0 : 0.0

  # Compute dimensionless frequencies (equation 30):
  beta = j*(1-alpha*sqrtalpha)
  kappa =  beta /  (alpha*sqrtalpha)

  # Compute disturbing function coefficients (equation 31):
  A_j00 = b[j+1,1]
  A_j10 =  alpha* b[j+1,2]
  A_j01 = -(A_j10 + A_j00)
  A_j20 =  alpha*alpha * b[j+1,3]
  A_j11 = -(2*A_j10 + A_j20)
  A_j02 = 2*A_j00 + 4*A_j10 + A_j20
  jd=convert(eltype(alpha),j)
  # Inner planet coefficients, in order k=0,-1,1,-2,2 (see Table 1):
  if j >=2
    f1[j+1,1]=alpha*u(beta          ,jd*(    A_j00-alpha*dj1),A_j10-alpha*dj1)
    f1[j+1,2]=alpha*u(beta-1.0      ,jd*(-jd*A_j00-0.5*A_j10+1.5*alpha*dj1),-jd*A_j10-0.5*A_j20+alpha*dj1)
    f1[j+1,3]=alpha*u(beta+1.0      ,jd*( jd*A_j00-0.5*A_j10-0.5*alpha*dj1), jd*A_j10-0.5*A_j20-alpha*dj1)
    f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
    f1[j+1,5]=alpha*u(beta+alpha*sqrtalpha,jd*(-jd*A_j00-0.5*A_j01),-jd*A_j10-0.5*A_j11)
  else
    if j==0
      f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
    else
      f1[j+1,1]=alpha*u(beta          ,jd*(    A_j00-alpha*dj1),A_j10-alpha*dj1)
      f1[j+1,2]=alpha*u(beta-1.0      ,jd*(-jd*A_j00-0.5*A_j10+1.5*alpha*dj1),-jd*A_j10-0.5*A_j20+alpha*dj1)
      f1[j+1,3]=alpha*u(beta+1.0      ,jd*( jd*A_j00-0.5*A_j10-0.5*alpha*dj1), jd*A_j10-0.5*A_j20-alpha*dj1)
      f1[j+1,4]=alpha*u(beta-alpha*sqrtalpha,jd*( jd*A_j00-0.5*A_j01-2.0*alpha*dj1), jd*A_j10-0.5*A_j11-2.0*alpha*dj1)
    end
  end
  # Add in the k=\pm 1 coefficients (note that d1 & d2 are the same as c1 & c2 for k=0):
  if j >= 1
    f1[j+1,2]=f1[j+1,2]+alpha*v(beta,jd*(A_j00-alpha*dj1),A_j10-alpha*dj1,-1)
    f1[j+1,3]=f1[j+1,3]+alpha*v(beta,jd*(A_j00-alpha*dj1),A_j10-alpha*dj1, 1)
  end

# Now for the outer planet:
# Outer planet coefficients, in order k=0,-2,2,-1,1 (see Table 1):
  one_over_alpha_squared = 1/(alpha*alpha)
  if j >= 2
    f2[j+1,1]=u(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared)
    f2[j+1,2]=u(kappa-1,-jd*(jd*A_j00-0.5*A_j01-0.5*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j02-dj1*one_over_alpha_squared)
    f2[j+1,3]=u(kappa+1,-jd*(-jd*A_j00-0.5*A_j01+1.5*dj1*one_over_alpha_squared),-jd*A_j01-0.5*A_j02+dj1*one_over_alpha_squared)
    f2[j+1,4]=u(kappa-1/(alpha*sqrtalpha),-jd*(-jd*A_j00-0.5*A_j10),-jd*A_j01-0.5*A_j11)
  else
    if j == 1
      f2[j+1,1]=u(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared)
      f2[j+1,2]=u(kappa-1,-jd*(jd*A_j00-0.5*A_j01-0.5*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j02-dj1*one_over_alpha_squared)
      f2[j+1,3]=u(kappa+1,-jd*(-jd*A_j00-0.5*A_j01+1.5*dj1*one_over_alpha_squared),-jd*A_j01-0.5*A_j02+dj1*one_over_alpha_squared)
    end
  end
  f2[j+1,5]=u(kappa+1/(alpha*sqrtalpha),-jd*(jd*A_j00-0.5*A_j10-2.0*dj1*one_over_alpha_squared),jd*A_j01-0.5*A_j11-2.0*dj1*one_over_alpha_squared)
# Add in the k=\pm 2 coefficients (note that d1 & d2 are the same as c1 & c2 for k=0):
  if j >= 1
    f2[j+1,2]=f2[j+1,2]+v(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared,-1)
    f2[j+1,3]=f2[j+1,3]+v(kappa,-jd*(A_j00-dj1*one_over_alpha_squared),A_j01-dj1*one_over_alpha_squared, 1)
  end
# That's it!
end
return 
end

"""
#/* compute Laplace coefficients and Leverrier derivative/
#          j
#     j   d     i
#    a   ---   b (a)
#          j    s
#        da
#
#   by series summation */
#/* Code due to Jack Wisdom */
"""
function laplace_wisdom(s::Rational,i::Integer,j::Integer,a::Number)
# function laplace_wisdom,s,i,j,a  IDL
# double laplace(double s, int i, int j, double a);  c

##define LAPLACE_EPS 1.0e-12
LAPLACE_EPS = convert(eltype(a),1.0e-12)

#if (i lt 0) then i = -i
i=abs(i)

if j <= i     #/* compute first term in sum */
  factor4 = one(a)
  for k=0:j-1
    factor4 *= i - k
  end
  lap_coef_sum = factor4
  q0 = 0
else
  q0 = fld(j + 1 - i,2)
  lap_coef_sum = zero(a)
  factor4 = one(a)
end

#  /* compute factors for terms in lap_coef_sum */

factor1 = s
factor2 = s + i
factor3 = i + 1
for q=1:q0-1       #/* no contribution for q = 0 */
  factor1 *= s + q
  factor2 *= s + i + q
  factor3 *= i+q+1
end
if q0 > 1
  q=q0
else
  q=1
end
#println(j+1-i,q0)
term = a*a * factor1 * factor2 / (factor3 * q)

#  /* sum series */

while  (term*factor4) > LAPLACE_EPS
  factor4 = one(a)
  for k=0:j-1
    factor4 *= (2*q + i - k)
  end
  lap_coef_sum += term * factor4
  factor1 += 1
  factor2 += 1
  factor3 += 1
  q = q+1
  term *= a*a * factor1 * factor2 / (factor3 * q)
end

#  /* fix coefficient */

for k=0:i-1
  lap_coef_sum *= (s+k)/(k+1)
end

apower = (q0 <= 0) ?  i : 2*q0 + i - 2
lap_coef_sum *= 2 * a^apower
# Return the Laplace Coefficient:
return lap_coef_sum
end




"""
# Computes transit-timing variations to linear order in
# eccentricity for non-resonant, plane-parallel planets.
# Input:
#   jmax:  Maximum j over which to sum the TTV calculation for both planets
#     p1:  Planet type for inner planet
#     p2:  Planet type for outer planet
#  time1:  Transit times for inner planet
#  time2:  Transit times for outer planet
# alpha0:  Initial alpha for Taylor expansion of coefficients
#     b0:  Laplace coefficients and derivatives for use in Taylor expansion
# Output:
#   ttv1: TTVs of the inner planet
#   ttv2: TTVs of the outer planet
#     f1: TTV coefficients for inner planet
#     f2: TTV coefficients for outer planet
#      b: Laplace coefficients (& derivatives) for outer planet
"""
function compute_ttv!(jmax::Integer,p1::Planet_plane_hk,p2::Planet_plane_hk,time1::Vector,time2::Vector,ttv1::Vector,ttv2::Vector,f1::Array,f2::Array,b::Array,alpha0::Number,b0::Array)

# Compute the semi-major axis ratio of the planets:
# println(p1.period,p2.period)
alpha = (p1.period/p2.period)^(2//3)  # Julia supports rational numbers!
# Number of times:
ntime1 = length(time1)
ntime2 = length(time2)
# Compute the coefficients:
ttv_succinct!(jmax+1,alpha,f1,f2,b,alpha0,b0)  # I need to compute coefficients one higher than jmax
# Compute TTVs for inner planet (equation 33):
# Compute since of \pomegas:
e1 = sqrt(p1.esinw*p1.esinw+p1.ecosw*p1.ecosw)
e2 = sqrt(p2.esinw*p2.esinw+p2.ecosw*p2.ecosw)
sin1om=p1.esinw/e1 
sin2om=p2.esinw/e2
cos1om=p1.ecosw/e1 
cos2om=p2.ecosw/e2
# Compute mean motions:
n1=2pi/p1.period
n2=2pi/p2.period
# Compute initial longitudes:
lam10=-n1*p1.trans0 + 2*p1.esinw # 2*p1.eccen*sin1om
lam20=-n2*p2.trans0 + 2*p2.esinw # 2*p2.eccen*sin2om
@inbounds  for i=1:ntime1
# Compute the longitudes of the planets at times of transit of planet 1 (equation 49):
  lam11 = n1*time1[i]+lam10
  lam21 = n2*time1[i]+lam20
  psi1  = lam11-lam21 # Compute difference in longitudes at times of transit of planet 1
  sinpsi1=sin(psi1)
  cospsi1=cos(psi1)
  sinlam11 = sin(lam11)
  coslam11 = cos(lam11) 
  sinlam1om1=sinlam11*cos1om-coslam11*sin1om
  coslam1om1=coslam11*cos1om+sinlam11*sin1om
  sinlam1om2=sinlam11*cos2om-coslam11*sin2om
  coslam1om2=coslam11*cos2om+sinlam11*sin2om
  ttv1[i]=zero(p1.period) #0.0
  sinjm1psi1=zero(p1.period) #0.0
  cosjm1psi1=one(p1.period) #1.0
# Sum over j:
  for j=1:jmax
    sinjpsi1=sinjm1psi1*cospsi1+cosjm1psi1*sinpsi1
    cosjpsi1=cosjm1psi1*cospsi1-sinjm1psi1*sinpsi1
    ttv1[i] += f1[j+1,1]*sinjpsi1
    ttv1[i] += f1[j+1,2]*e1*(sinjpsi1*coslam1om1-cosjpsi1*sinlam1om1)
    ttv1[i] += f1[j+1,3]*e1*(sinjpsi1*coslam1om1+cosjpsi1*sinlam1om1)
    ttv1[i] += f1[j  ,4]*e2*(sinjpsi1*coslam1om2-cosjpsi1*sinlam1om2)
    ttv1[i] += f1[j+2,5]*e2*(sinjpsi1*coslam1om2+cosjpsi1*sinlam1om2)
    sinjm1psi1=sinjpsi1
    cosjm1psi1=cosjpsi1
  end
# Multiply by period and mass ratio, and divide by 2*Pi:
  ttv1[i] = ttv1[i]*p1.period*p2.mass_ratio/(2pi)
end
# Compute TTVs for outer planet (equation 33):
@inbounds for i=1:ntime2
# Compute the longitudes of the planets at times of transit of planet 2:
  lam12 = n1*time2[i]+lam10
  lam22 = n2*time2[i]+lam20
  sinlam22 = sin(lam22)
  coslam22 = cos(lam22) 
  psi2  = lam12-lam22 # Compute difference in longitudes at times of transit of planet 2
  sinpsi2=sin(psi2)
  cospsi2=cos(psi2)
  sinlam2om1=sinlam22*cos1om-coslam22*sin1om
  coslam2om1=coslam22*cos1om+sinlam22*sin1om
  sinlam2om2=sinlam22*cos2om-coslam22*sin2om
  coslam2om2=coslam22*cos2om+sinlam22*sin2om
  ttv2[i]=zero(p2.period) #0.0
  sinjm1psi2=zero(p2.period) #0.0
  cosjm1psi2=one(p2.period) #1.0
# Sum over j:
  for j=1:jmax
    sinjpsi2=sinjm1psi2*cospsi2+cosjm1psi2*sinpsi2
    cosjpsi2=cosjm1psi2*cospsi2-sinjm1psi2*sinpsi2
    ttv2[i] += f2[j+1,1]*sinjpsi2
    ttv2[i] += f2[j+1,2]*e2*(sinjpsi2*coslam2om2-cosjpsi2*sinlam2om2)
    ttv2[i] += f2[j+1,3]*e2*(sinjpsi2*coslam2om2+cosjpsi2*sinlam2om2)
    ttv2[i] += f2[j+2,4]*e1*(sinjpsi2*coslam2om1-cosjpsi2*sinlam2om1)
    ttv2[i] += f2[j  ,5]*e1*(sinjpsi2*coslam2om1+cosjpsi2*sinlam2om1)
    sinjm1psi2=sinjpsi2
    cosjm1psi2=cosjpsi2
  end
# Multiply by period and mass ratio, and divide by 2*Pi:
  ttv2[i] = ttv2[i]*p2.period*p1.mass_ratio/(2pi)
end
# Finished!
return
end  # compute_ttv!

end # module

Main.var"workspace#3".TTVFaster

To run this tutorial locally, download this file and open it with Pluto.jl.

To run this tutorial locally, download this file and open it with Pluto.jl.

To run this tutorial locally, download this file and open it with Pluto.jl.

To run this tutorial locally, download this file and open it with Pluto.jl.