module binary use types, only: & real64, int32 implicit none integer(kind=int32), parameter :: & VECLEN = 12 contains ! ======================================================================= ! Direct code ! ======================================================================= function direct(trel, y) result(yp) use, intrinsic :: & ieee_arithmetic, only: & ieee_value, ieee_quiet_nan use types, only: & int32, real64 use physconst, only: & GRAV, CLIGHT use stardata, only: & IMASS, IRADI, IAPSI, IINER, IQUAL, & toffset, stars, & star_huntpol use vec, only: & cross_product use deriv_data, only: & local_store, local_data use parameters, only: & cutoff, has_cutoff use states, only: & status, status_stars, & STATUS_OK, STATUS_COLLIDE, STATUS_ESCAPE use flags_data, only: & use_sd, use_td, use_tf, use_pn, & use_gr, use_so, use_p2, use_ss, use_rr, & use_rxv 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 real(kind=real64), dimension(3) :: & r, v, re, & fpn, acc, w, & vx real(kind=real64), dimension(3, 2) :: & j12, w12, & fsd12, ftd12, ftf12, f12, torque, & wxr12 real(kind=real64), dimension(2) :: & m12, s12, k12, i12, q12, & mr12, mi12, wn212, wre12, s5k12, qi12, 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 ! higher-order terms real(kind=real64), dimension(3, 2) :: & tgr12, jx12, tso12 real(kind=real64), dimension(3) :: & j, ddd, jo, & fso, f2pn, fss, frr, & x3sd !, lso real(kind=real64), dimension(2) :: & jr12 real(kind=real64) :: & j12dot, mum1, phi2, vr2, v4, vr4, phiv2, phivr2 r(:) = y(1:3) v(:) = y(4:6) j12(:,:) = reshape(y(7:12), (/3,2/)) t = trel + toffset do i=1,2 star12(:,i) = star_huntpol(t, i) end do m12 = star12(IMASS, :) s12 = star12(IRADI, :) k12 = star12(IAPSI, :) i12 = star12(IINER, :) q12 = star12(IQUAL, :) mr12(:) = m12(2:1:-1) mi12(:) = 1.d0 / m12 mrat12(:) = mr12(:) * mi12(:) qi12(:) = 1.d0 / q12(:) m = sum(m12(:)) mm1 = 1.d0 / m mu = product(m12) * mm1 do i=1,2 if (i12(i) >0.0) then w12(:,i) = j12(:,i) / i12(i) else w12(:,i) = 0.d0 endif end do wn212 = sum(w12**2, dim=1) rn = norm2(r) rnm1 = 1.d0 / rn rnm2 = rnm1**2 rnm3 = rnm2 * rnm1 rnm4 = rnm2**2 rnm8 = rnm4**2 re(:) = r(:) * rnm1 ! fail on merger if (rn < sum(s12(:))) then write(6, "(' [direct] object collision')") yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/1, 2/) return endif ! fail on escape if (has_cutoff.and.(rn > cutoff(1))) then write(6, "(' [direct] orbit escape')") yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/1, 0/) 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 (use_rxv) then nm1 = 1.d0 / norm2(w(:)) else eps = 0.5d0 * v2 - phi an = -0.5d0 * Gm / eps nm1 = sqrt(an**3 / Gm) end if s5k12(:) = s12(:)**5 * k12(:) vrn = dot_product(v(:), re(:)) eta = mu * mm1 f12(:,:) = 0.d0 torque(:,:) = 0.d0 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_td) then b = -2.d0 * h * wre12(i) fsd12(:,i) = & (h * (5.d0 * wre12(i)**2 - wn212(i))) * re(:) & + b * w12(:,i) f12(:,i) = f12(:,i) + fsd12(:,i) torque(:,i) = torque(:,i) + (b * mu) * wxr12(:,i) else ftd12(:,i) = 0.d0 endif end do else f12(:,:) = 0.d0 endif if (use_td .and. use_tf) then f = 6.d0 * Gm * nm1 * rnm8 vx(:) = (-2.d0 * vrn) * re(:) - v(:) do i=1,2 h = f * qi12(i) * s5k12(i) * mrat12(i) ftf12(:,i) = h * (vx(:) + wxr12(:,i)) f12(:,i) = f12(:,i) + ftf12(:,i) b = h * mu * rn torque(:,i) = torque(:,i) & + (b * wre12(i)) * r(:) & + (b * rn) * (w(:) - w12(:,i)) end do else ftf12(:,:) = 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 (local_store) then local_data( 1:12) = yp(1:12) local_data(13:18) = reshape(fsd12, (/6/)) local_data(19:24) = reshape(ftd12, (/6/)) local_data(25:30) = reshape(ftf12, (/6/)) local_data(31:33) = fpn(1:3) ! higher order GR terms local_data(34:36) = fso(1:3) local_data(37:39) = f2pn(1:3) local_data(40:42) = fss(1:3) local_data(43:45) = frr(1:3) local_data(46:51) = reshape(tgr12, (/6/)) endif end function direct ! ======================================================================= ! Secular code ! ======================================================================= function secular(trel, y) result(yp) use, intrinsic :: & ieee_arithmetic, only: & ieee_value, ieee_quiet_nan use types, only: & int32, real64 use physconst, only: & GRAV, CLIGHT use stardata, only: & IMASS, IRADI, IAPSI, IINER, IQUAL, & toffset, stars, & star_huntpol use vec, only: & cross_product use deriv_data, only: & local_store, local_data use parameters, only: & cutoff, has_cutoff use states, only: & status, status_stars, & STATUS_OK, STATUS_COLLIDE, STATUS_ESCAPE 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,q12, & mr12, mi12, qi12, & we12, wh12, wq12, & s5kmi12, x1, x2, 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 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, :) k12(:) = star12(IAPSI, :) i12(:) = star12(IINER, :) q12(:) = star12(IQUAL, :) mr12(:) = m12(2:1:-1) mi12(:) = 1.d0 / m12(:) qi12(:) = 1.d0 / q12(:) 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 merger: S1+S2 > rp = (1-e)*a = l/(1+e) if (1.d0 < sum(s12(:)) * c * (1.d0 + en)) then write(6, "(' [secular] object collision')") yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/1, 2/) return endif ! fail escape: cutoff < a = l/(1-e**2) if (has_cutoff.and.(1.d0 > cutoff(1) * c * f2)) then write(6, "(' [secular] orbit escape')") yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/1, 0/) 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 if (use_tf) 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 x2 = s5kmi12 * Qi12 * mr12 * (6.d0 * c5) do i=1,2 dTF12(:,i) = x2(i)*hn*( & ee(:)*(f_2*we12(i)) & +qe(:)*(g_2*wq12(i)) & +he(:)*(g_4*wh12(i)-g_3*hc2)) gTF12(:,i) = x2(i)*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) 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