module binary use parameters, only: & verb_interact implicit none contains ! ======================================================================= ! Direct code ! ======================================================================= function direct(trel, y) result(yp) use, intrinsic :: & ieee_arithmetic, only: & ieee_value, ieee_quiet_nan use typedef, only: & int32, real64 use physconst, only: & GRAV, CLIGHT use stardata, only: & IMASS, IRADI, IAPSI, IINER, IINEZ, ITAUQ, & toffset, stars, & star_huntpol use vector, only: & cross_product, & cross_product_tensor_contract use deriv_data, only: & local_store, local_data use parameters, only: & cutoff, has_cutoff, & interact, has_interact, & epslim, omegalim use states, only: & status, status_stars, & STATUS_COLLIDE, STATUS_ESCAPE, STATUS_INTERACT use integrator_flags, only: & terminate_integration use flags_data, only: & use_sd, use_td, use_tf, use_pn, & use_gr, use_so, use_p2, use_ss, use_rr, & use_rxv, & use_ori, use_orq, use_ors, & use_pqm, use_pqq, use_cd use orientation, only: & rotation_w2A, & rotation_f2A, & rotation_dodt, & rotation_dqdt, & rotate_diag use collision, only: & collision_inspector implicit none real(kind=real64), intent(in) :: & trel real(kind=real64), intent(in), dimension(:) :: & y real(kind=real64), dimension(size(y, 1)) :: & yp real(kind=real64), parameter :: & CLIGHTm1 = 1.d0 / CLIGHT, & CLIGHTm2 = CLIGHTm1**2, & CLIGHTm4 = CLIGHTm2**2, & CLIGHTm5 = CLIGHTm4 * CLIGHTm1, & CLIGHTm6 = CLIGHTm5 * CLIGHTm1, & Gcm2 = GRAV * CLIGHTm2, & Gcm6 = Gcm2 * CLIGHTm4, & ONETHIRD = 1.d0 / 3.d0, & G3 = GRAV * 3.d0 integer(kind=int32), parameter :: & nstar = 2, & ndim = 3, & nquat = 4, & nedge = (nstar * (nstar - 1)) / 2 real(kind=real64), dimension(nstar), parameter :: & sc = (/1.d0, -1.d0/) integer(kind=int32), dimension(nstar), parameter :: & ir = (/2, 1/) ! local variables real(kind=real64), dimension(ndim) :: & r, v, re, & fpn, acc, w, & vx real(kind=real64), dimension(nquat, nstar) :: & orq12, roq real(kind=real64), dimension(ndim, nstar) :: & j12, ori12, w12, i12, & fsd12, ftd12, ftf12, f12, torque, rot, & wxr12, qr12, fpq12, tpq12, ttf12, tsd12 real(kind=real64), dimension(nstar) :: & m12, s12, k12, t12, & mr12, mi12, wn212, wre12, s5k12, mrat12 real(kind=real64), dimension(stars(1)%m, 2) :: & star12 integer(kind=int32) :: & i real(kind=real64) :: & t, f, g, h, b, & m, mu, rn, rnm1, rnm2, rnm3, rnm4, rnm8, & v2, phi, eps, an, nm1, vrn, eta, mm1, Gm, & ftgr logical, dimension(nstar) :: & lpqm12 logical, dimension(nedge) :: & collide ! higher-order terms real(kind=real64), dimension(ndim, nstar) :: & tgr12, jx12, tso12 real(kind=real64), dimension(ndim) :: & j, ddd, jo, & fso, f2pn, fss, frr, & x3sd !, lso real(kind=real64), dimension(nstar) :: & jr12 real(kind=real64) :: & j12dot, mum1, phi2, vr2, v4, vr4, phiv2, phivr2 ! quadrupole-quadrupole interaction real(kind=real64), dimension(nstar) :: & qrr12 real(kind=real64), dimension(ndim) :: & fqq, tqqq real(kind=real64), dimension(ndim, nstar) :: & q12, tqq, qqqr12 real(kind=real64), dimension(ndim, ndim, nstar) :: & qq12, A12 ! start of code t = trel + toffset do i=1,2 star12(:,i) = star_huntpol(t, i) end do m12(:) = star12(IMASS, :) s12(:) = star12(IRADI, :) i12(:,:) = star12(IINER, :) k12(:) = star12(IAPSI, :) t12(:) = star12(ITAUQ, :) ! use t < 0 as flag for containing -1/Q, tau otherwise instead r(:) = reshape(y(1:3), (/ndim/)) v(:) = reshape(y(4:6), (/ndim/)) j12(:,:) = reshape(y(7:12), (/ndim,nstar/)) if (use_ori) then if (use_orq) then orq12(:,:) = reshape(y(13:20), (/nquat,nstar/)) roq(:,:) = 0.d0 else ori12(:,:) = reshape(y(13:18), (/ndim,nstar/)) rot(:,:) = 0.d0 endif end if do i=1,nstar if (use_ori) then ! alternate approach if use i_x == i_y == i_z ! if (abs(i12(1,i)-i12(2,i)) + abs(i12(2,i)-i12(3,i)) > 1.d-12 * I12(3,i)) then ! here we assume i_x == i_y == 0 for spherical bodies ! and i_z == 0 if just point mass w/o spin if (i12(3,i) == 0.d0) then w12(:,i) = 0.d0 lpqm12(i) = .false. cycle endif if (i12(1,i) > 0.d0) then if (use_orq) then A12(:,:,i) = rotation_f2A(orq12(:,i)) else A12(:,:,i) = rotation_w2A(ori12(:,i)) endif w12(:,i) = matmul(A12(:,:,i), matmul(j12(:,i), A12(:,:,i)) / i12(:,i)) lpqm12(i) = .true. else w12(:,i) = j12(:,i) / i12(3,i) lpqm12(i) = .false. end if if (.not.(lpqm12(i) .or. use_ors)) & cycle if (use_orq) then roq(:,i) = rotation_dqdt(orq12(:,i), w12(:,i)) else rot(:,i) = rotation_dodt(ori12(:,i), w12(:,i)) endif ! TODO - modify rot(:,:) for GR/SR time delay (just because we can) if (.not.(lpqm12(i) .and. (use_pqm .or. use_pqq))) & cycle q12(:,i) = ONETHIRD * sum(i12(:,i)) - i12(:,i) qr12(:,i) = matmul(A12(:,:,i), q12(:,i) * matmul(r(:), A12(:,:,i))) qrr12(i) = sum(qr12(:,i) * r(:)) else if (i12(3,i) > 0.d0) then w12(:,i) = j12(:,i) / i12(3,i) else w12(:,i) = 0.d0 endif lpqm12(i) = .false. end if end do wn212 = sum(w12(:,:)**2, dim=1) mr12(:) = m12(2:1:-1) mi12(:) = 1.d0 / m12 mrat12(:) = mr12(:) * mi12(:) m = sum(m12(:)) mm1 = 1.d0 / m mu = product(m12) * mm1 rn = norm2(r) rnm1 = 1.d0 / rn rnm2 = rnm1**2 rnm3 = rnm2 * rnm1 rnm4 = rnm2**2 rnm8 = rnm4**2 re(:) = r(:) * rnm1 ! fail on interaction if (has_interact.and.(rn < interact(1,2))) then if (verb_interact) then write(6, "(' [direct] object interaction ',i1,' + ',i1)") 1, 2 endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = (/1, 2/) terminate_integration = .true. return endif if (use_ori.and.use_cd) then collide(:) = collision_inspector(i12, m12, s12, a12, r, nstar, nedge) else collide(1) = rn < sum(s12(:)) end if if (collide(1)) then if (verb_interact) then write(6, "(' [direct] object collision ',i1,' + ',i1)") 1, 2 endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/1, 2/) terminate_integration = .true. return endif ! fail on escape if (has_cutoff.and.(rn > cutoff(1))) then if (verb_interact) then write(6, "(' [direct] orbit escape')") endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/1, 0/) terminate_integration = .true. return endif do i=1,2 wre12(i) = dot_product(w12(:,i), re(:)) wxr12(:,i) = cross_product(w12(:,i), r(:)) end do v2 = sum(v(:)**2) Gm = GRAV * m phi = Gm * rnm1 ! check sign! w(:) = cross_product(re(:), v(:)) * rnm1 if (any(t12(:) < 0.d0)) then if (use_rxv) then nm1 = 1.d0 / max(norm2(w(:)), omegalim*sqrt(Gm*rnm3)) else eps = min(0.5d0 * v2 - phi, -epslim * phi) an = -0.5d0 * Gm / eps nm1 = sqrt(an**3 / Gm) end if endif s5k12(:) = s12(:)**5 * k12(:) vrn = dot_product(v(:), re(:)) eta = mu * mm1 f12(:,:) = 0.d0 torque(:,:) = 0.d0 ! permanent quadrupole-quadrupole interaction if (use_ori.and.use_pqq.and.all(lpqm12)) then mum1 = 1.d0 / mu g = 7.5d0 * GRAV * rnm3 * rnm4 do i=1, nstar qq12(:,:,i) = rotate_diag(q12(:,i),A12(:,:,i)) qqqr12(:,i) = matmul(qq12(:,:,i), qr12(:,ir(i))) enddo h = 7.d0 * rnm2 f = -31.5d0 * rnm4 * qrr12(1) * qrr12(2) & - sum(qq12(:,:,1) * qq12(:,:,2)) & + 14.d0 * rnm2 * dot_product(qr12(:,1), qr12(:,2)) ddd(:) = -2.d0 * (qqqr12(:,1) + qqqr12(:,2)) & + (h * qrr12(2)) * qr12(:,1) & + (h * qrr12(1)) * qr12(:,2) fqq(:) = (g * mum1) * (ddd(:) + f * r(:)) f12(:,1) = f12(:,1) + fqq(:) tqqq(:) = 2.d0 * g * cross_product(qr12(:,2), qr12(:,1)) & + 3.d0 * GRAV * rnm2 * rnm3 * & cross_product_tensor_contract(qq12(:,:,1), qq12(:,:,2)) do i = 1, nstar tqq(:,i) = sc(i) * tqqq(:) & - 2.d0 * g * cross_product(qqqr12(:,i), r(:)) & + g * h * qrr12(ir(i)) * cross_product(qr12(:,i), r(:)) enddo torque(:,:) = torque(:,:) + tqq(:,:) else fqq(:) = 0.d0 tqq(:,:) = 0.d0 endif ! permanent quadrupole-monopole interaction if (use_ori.and.use_pqm.and.any(lpqm12)) then mum1 = 1.d0 / mu g = G3 * rnm2 * rnm3 do i=1, nstar if (.not.lpqm12(i)) then fpq12(:,i) = 0.d0 tpq12(:,i) = 0.d0 cycle endif f = g * mr12(i) fpq12(:,i) = f * mum1 * ( & qr12(:,i) - 2.5d0 * rnm2 * qrr12(i) * r(:)) tpq12(:,i) = f * cross_product(qr12(:,i), r(:)) f12(:,i) = f12(:,i) + fpq12(:,i) torque(:,i) = torque(:,i) + tpq12(:,i) enddo else fpq12(:,:) = 0.d0 tpq12(:,:) = 0.d0 endif ! spin and tidal deformation, monopole interaction if (use_sd .or. use_td) then if (use_td) then f = -12.d0 * GRAV * rnm3 endif g = m * rnm4 do i=1,2 h = s5k12(i) * mi12(i) * g if (use_td) then ftd12(:,i) = (f * h * mr12(i)) * re(:) f12(:,i) = f12(:,i) + ftd12(:,i) else ftd12(:,i) = 0.d0 endif if (use_sd) then b = -2.d0 * h * wre12(i) fsd12(:,i) = & (h * (5.d0 * wre12(i)**2 - wn212(i))) * re(:) & + b * w12(:,i) tsd12(:,i) = (b * mu) * wxr12(:,i) torque(:,i) = torque(:,i) + tsd12(:,i) f12(:,i) = f12(:,i) + fsd12(:,i) else ftd12(:,i) = 0.d0 tsd12(:,i) = 0.d0 endif end do else ftd12(:,:) = 0.d0 fsd12(:,:) = 0.d0 tsd12(:,:) = 0.d0 endif if (use_td .and. use_tf .and. any(t12 /= 0.d0)) then f = 12.d0 * Gm * rnm8 vx(:) = (-2.d0 * vrn) * re(:) - v(:) do i=1,2 h = t12(i) if (h == 0.d0) then ftf12(:,i)=0.d0 cycle else if (h < 0.d0) then h = -0.5d0 * nm1 * h endif h = h * f * s5k12(i) * mrat12(i) ftf12(:,i) = h * (vx(:) + wxr12(:,i)) f12(:,i) = f12(:,i) + ftf12(:,i) b = h * mu * rn ttf12(:,i) = (b * wre12(i)) * r(:) & + (b * rn) * (w(:) - w12(:,i)) torque(:,i) = torque(:,i) + ttf12(:,i) end do else ftf12(:,:) = 0.d0 ttf12(:,i) = 0.d0 endif acc(:) = (- phi * rnm2) * r(:) + sum(f12(:,:), dim=2) if (use_pn.or.use_gr) then vr2 = vrn**2 fpn(:) = & phi * rnm1 * CLIGHTm2 * ( & v(:) * ( & ( 4.d0 - 2.d0 * eta) * vrn) & + re(:) * ( & + ( 4.d0 + 2.d0 * eta) * phi & + 1.5d0 * eta * vr2 & + (-1.d0 - 3.d0 * eta) * v2) & ) acc(:) = acc(:) + fpn(:) ftgr = CLIGHTm2 * (0.5d0 * v2 * (1.d0 - 3.d0 * eta) + (3.d0 + eta) * phi) else fpn(:) = 0.d0 endif ! higher order GR terms if (use_gr) then phi2 = phi**2 do i=1, 2 jr12(i) = dot_product(j12(:,3-i), re(:)) enddo if (use_so) then j(:) = sum(j12, dim=2) jo(:) = mu * w(:) * rn**2 jx12(:,1) = cross_product(j12(:,1), j12(:,1)) jx12(:,2) = -jx12(:,1) ddd(:) = (j12(:,2) * mi12(2) - j12(:,1) * mi12(1)) * (m12(1) - m12(2)) x3sd(:) = 3.d0 * j(:) + ddd(:) fso(:) = rnm3 * Gcm2 * ( & re(:) * ( 6.d0 * & dot_product(cross_product(re(:), v(:)), 2.d0 * j(:) + ddd(:))) & + cross_product(7.d0 * j(:) + 3.d0 * ddd(:) , v(:)) & + (3.d0 * vrn * cross_product(re(:), x3sd(:)))) do i=1,2 tso12(:,i) = rnm3 * Gcm2 * ( & cross_product(jo(:), j12(:,i)) * (2.d0 + 1.5d0 * mrat12(i)) & + jx12(:,i) & + cross_product(re(:), j12(:,i)) * (3.d0 * jr12(i))) end do ! TODO - add L_SO contribution to QD/TF ! lso(:) = eta * CLIGHTm2 * ( & ! + phi * double_cross_product(re(:), x3sd(:}, 1.d0) & ! - 0.5d0 * double_cross_product(v(:), j(:) + ddd(:), v2) ! forall (i=1:2) tgr12(:,i) = tgr12(:,i) + lso(:) * f12(:,i) ! TODO - end else fso(:)=0.d0 tso12(:,:)=0.d0 endif if (use_p2) then v4 = v2**2 vr4 = vr2**2 phiv2 = phi * v2 phivr2 = phi * vr2 f2pn(:) = phi * rnm1 * CLIGHTm4 * ( & re(:) * ( & + ( -9.d0 + eta * (-21.75d0) ) * phi2 & + ( eta * ( -3.d0 + eta * 4.d0 )) * v4 & + ( eta * ( -1.875d0 + eta * 5.625d0)) * vr4 & + ( eta * ( 4.5d0 + eta * (-6.d0 ))) * vr2 * v2 & + ( eta * ( 6.5d0 + eta * (-2.d0 ))) * phiv2 & + ( 2.d0 + eta * ( 25.d0 + eta * 2.d0 )) * phivr2) & + v(:) * (vrn * ( & + ( eta * ( 7.5d0 + eta * 2.d0 )) * v2 & + ( -2.d0 + eta * (-20.5d0 + eta * (-4.d0 ))) * phi & + ( eta * ( -4.5d0 + eta * (-3.d0 ))) * vr2))) ! This correction will need to be checked for more extreme systems ! before using it by default ! ! 2PN correction to QD, TF ftgr = ftgr + CLIGHTm4 * ( & + ( 0.375d0 + eta * (-2.625d0 + eta * 4.875d0 )) * v4 & + ( + eta * (-1.d0 + eta * (-2.5d0))) * phivr2 & + ( 3.5d0 + eta * (-5.d0 + eta * (-4.5d0))) * phiv2 & + ( 3.5d0 + eta * (-10.25d0 + eta )) * phi2) else f2pn(:) = 0.d0 end if if (use_ss) then mum1 = 1.d0 / mu j12dot = sum(product(j12, dim=2)) fss(:) = -3.d0 * rnm4 * mum1 * Gcm2 * ( & re(:) * (j12dot - 5.d0 * product(jr12)) & + j12(:,1) * jr12(1) & + j12(:,2) * jr12(2)) ! TODO - add L_SS contribution to QD/TF else fss(:)=0.d0 endif if (use_rr) then frr(:) = 1.6d0 * eta * phi2 * rnm1 * CLIGHTm5 * ( & re(:) * (vrn * (18.d0 * v2 + 1.5d0 * phi - 25.d0 * vr2)) & + v(:) * (2.d0 * phi + 15.d0 * vr2 - 6.d0 * v2)) else frr(:)=0.d0 endif acc(1:3) = acc(1:3) + fso(:) + f2pn(:) + fss(:) + frr(:) tgr12(:,:) = tso12(:,:) endif ! end higher-order GR terms if (use_pn) then torque(:,:) = torque(:,:) * (1.d0 + ftgr) endif if (use_gr) then torque(:,:) = torque(:,:) + tgr12(:,:) endif yp( 1: 3) = v(1:3) yp( 4: 6) = acc(1:3) yp( 7:12) = reshape(torque, (/6/)) if (use_ori) then if (use_orq) then yp(13:20) = reshape(roq, (/8/)) else yp(13:18) = reshape(rot, (/6/)) endif endif if (local_store) then local_data( 1:size(yp,1)) = yp(:) ! basic forces and torques local_data(21:26) = reshape(fsd12, (/6/)) local_data(27:32) = reshape(ftd12, (/6/)) local_data(33:38) = reshape(ftf12, (/6/)) local_data(39:44) = reshape(tsd12, (/6/)) local_data(45:50) = reshape(ttf12, (/6/)) local_data(51:53) = fpn(1:3) ! higher order GR terms local_data(54:56) = fso(1:3) local_data(57:59) = f2pn(1:3) local_data(60:62) = fss(1:3) local_data(63:65) = frr(1:3) local_data(66:71) = reshape(tgr12, (/6/)) ! permanent quadrupole momenent local_data(72:77) = reshape(fpq12, (/6/)) local_data(78:83) = reshape(tpq12, (/6/)) local_data(84:86) = fqq(1:3) local_data(87:92) = reshape(tqq, (/6/)) endif end function direct ! ======================================================================= ! Secular code ! ======================================================================= function secular(trel, y) result(yp) use, intrinsic :: & ieee_arithmetic, only: & ieee_value, ieee_quiet_nan use typedef, only: & int32, real64 use physconst, only: & GRAV, CLIGHT use stardata, only: & IMASS, IRADI, IAPSI, IINEZ, ITAUQ, & toffset, stars, & star_huntpol use vector, only: & cross_product use deriv_data, only: & local_store, local_data use parameters, only: & cutoff, has_cutoff, & interact, has_interact use states, only: & status, status_stars, & STATUS_COLLIDE, STATUS_ESCAPE, STATUS_INTERACT use integrator_flags, only: & terminate_integration use flags_data, only: & use_sd, use_td, use_tf, use_pn implicit none real(kind=real64), intent(in) :: & trel real(kind=real64), intent(in), dimension(:) :: & y real(kind=real64), dimension(size(y, 1)) :: & yp real(kind=real64), parameter :: & CLIGHTm2 = 1.d0 / CLIGHT**2 real(kind=real64), dimension(3) :: & e, h, & ee, he, qe, & gPN, edot, hdot real(kind=real64), dimension(3, 2) :: & j12, w12, & dQD12, gQD12, dTF12, gTF12, & g12, d12, torque real(kind=real64), dimension(2) :: & m12, s12, k12, i12, t12, & mr12, mi12, & we12, wh12, wq12, & s5kmi12, x1, wqh12 real(kind=real64) :: & m, mu, & t, Gm, en2, hn2, hn, en, & hnm1, hn3, & f, f2, f3, & c, c2, c3, c4, c5, & f_0, f_1, f_2, g_2, g_3, g_4, & hc2, s, s1 real(kind=real64), dimension(stars(1)%m, 2) :: & star12 integer(kind=int32) :: & i e(:) = y(1:3) h(:) = y(4:6) j12(:,:) = reshape(y(7:12), (/3,2/)) t = trel + toffset ! TODO - make star data type do i=1,2 star12(:,i) = star_huntpol(t, i) end do m12(:) = star12(IMASS, :) s12(:) = star12(IRADI, :) i12(:) = star12(IINEZ, :) k12(:) = star12(IAPSI, :) t12(:) = star12(ITAUQ, :) ! use q < 0 as flag for containing -1/Q, tau for q > 0 mr12(:) = m12(2:1:-1) mi12(:) = 1.d0 / m12(:) m = sum(m12) mu = product(m12) / m Gm = GRAV * m ! cm**3/s**2 do i=1, 2 w12(:,i) = j12(:,i) / i12(i) enddo en2 = sum(e(:)**2) hn2 = sum(h(:)**2) en = sqrt(en2) hn = sqrt(hn2) ! cm**2/s hnm1 = 1.d0 / hn hn3 = hn2 * hn ee = e(:) / en he = h(:) * hnm1 qe = cross_product(he(:), ee(:)) do concurrent (i=1:2) we12(i) = dot_product(w12(:,i), ee(:)) wh12(i) = dot_product(w12(:,i), he(:)) wq12(i) = dot_product(w12(:,i), qe(:)) end do f2 = 1.d0 - en2 f = sqrt(f2) f3 = f * f2 c = Gm * hnm1**2 ! 1/cm = 1/ln c2 = c**2 c3 = c2 * c c4 = c3 * c c5 = c4 * c s5kmi12(:) = s12(:)**5 * k12(:) * mi12(:) ! 1/a = c * f2 ! n = h * c2 * f3 ! fail on interaction: interact(1,2) > rp = (1-e)*a = l/(1+e) if (has_interact.and.(1.d0 < interact(1,2) * c * (1.d0 + en))) then if (verb_interact) then write(6, "(' [secular] object interaction ',i1,' + ',i1)") 1, 2 endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = (/1, 2/) terminate_integration = .true. return endif ! fail on merger: S1+S2 > rp = (1-e)*a = l/(1+e) if (1.d0 < sum(s12(:)) * c * (1.d0 + en)) then if (verb_interact) then write(6, "(' [secular] object collision')") endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/1, 2/) terminate_integration = .true. return endif ! fail escape: cutoff < a = l/(1-e**2) if (has_cutoff.and.(1.d0 > cutoff(1) * c * f2)) then if (verb_interact) then write(6, "(' [secular] orbit escape')") endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/1, 0/) terminate_integration = .true. return endif g12(:,:) = 0.d0 d12(:,:) = 0.d0 if (use_sd .or. use_td) then f_0 = 15.d0 + en2 * (22.5d0 + en2 * 1.875d0) wqh12 = wq12 * wh12 x1 = s5kmi12 * (m * c3 * f3) do i=1,2 dQD12(:,i) = x1(i)*( & qe(:)*( we12(i)*wh12(i)) & +ee(:)*(-wqh12(i))) gQD12(:,i) = x1(i)*(en*hnm1)*( & qe(:)*(wh12(i)**2-0.5d0*(we12(i)**2+wq12(i)**2) & +mr12(i)*(GRAV*c3*f_0)) & +he(:)*wqh12(i)) end do g12(:,:) = g12(:,:) + gQD12(:,:) d12(:,:) = d12(:,:) + dQD12(:,:) endif ! TODO - implement case for use_rxv if (use_tf .and. any(t12 /= 0.d0)) then f_1 = 9.d0 + en2 * (33.75d0 + en2 * (16.875d0 + en2 * 0.703125d0)) f_2 = 0.5d0 + en2 * ( 0.75d0 + en2 * 0.0625d0) ! g_1 same as f_2 g_2 = 0.5d0 + en2 * (2.25d0 + en2 * 0.3125d0) g_3 = 1.d0 + en2 * (7.5d0 + en2 * ( 5.625d0 + en2 * 0.3125d0 )) g_4 = 1.d0 + en2 * (3.0d0 + en2 * 0.375d0) hc2 = hn * c2 s1 = -6.d0 * c5 do i=1,2 s = t12(i) if (s == 0) then dTF12(:,i) = 0.d0 gTF12(:,i) = 0.d0 cycle else if (s > 0.d0) then s = -2.d0 * hc2 * f3 * s ! sadly, this is re-computed if bith Q > 0. endif s = s * s1 * s5kmi12(i) * mr12(i) dTF12(:,i) = s*hn*( & ee(:)*(f_2*we12(i)) & +qe(:)*(g_2*wq12(i)) & +he(:)*(g_4*wh12(i)-g_3*hc2)) gTF12(:,i) = s*en*( & ee(:)*(11.d0*f_2*wh12(i)-f_1*hc2) & -he(:)*(f_2*we12(i))) end do g12(:,:) = g12(:,:) + gTF12(:,:) d12(:,:) = d12(:,:) + dTF12(:,:) end if edot(:) = sum(g12, dim=2) hdot(:) = sum(d12, dim=2) ! TODO - implement torque correction and 2nd GR if (use_pn) then gPN(:) = qe(:)*(3.d0*en*hn3*f3*c4*CLIGHTm2) edot(:) = edot(:) + gPN(:) else gPN(:) = 0.d0 endif torque(:,:) = - mu * d12(:,:) yp(1:3) = edot(1:3) yp(4:6) = hdot(1:3) yp(7:12) = reshape(torque, (/6/)) if (local_store) then local_data(1:12) = yp(1:12) local_data(13:18) = reshape(dQD12, (/6/)) local_data(19:24) = reshape(gQD12, (/6/)) local_data(25:30) = reshape(dTF12, (/6/)) local_data(31:36) = reshape(gTF12, (/6/)) local_data(37:39) = gPN(1:3) endif end function secular end module binary