module quaternary ! based on binary use typedef, only: & real64, int32 implicit none integer(kind=int32), parameter :: & VECLEN = 30 contains ! ======================================================================= ! Direct code, pair of pairs ! ======================================================================= function direct4p(trel, y) result(yp) use, intrinsic :: & ieee_arithmetic, only: & ieee_value, ieee_quiet_nan, ieee_is_nan, ieee_signaling_nan use typedef, only: & int32, real64 use physconst, only: & GRAV, CLIGHT use stardata, only: & IMASS, IRADI, IAPSI, IINER, IQUAL, & 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_OK, 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_pnq, use_rxv, use_td3, & use_pnf, & use_pnt, use_sdt, use_tdt, use_tft, use_t3t, use_t3h, use_t3s 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 integer(kind=int32), parameter :: & nstar = 4, & norb = 3, & ndim = 3, & nbin = 2, & ! a binary contains 2 stars ndoub = 2 ! number of binaries integer(kind=int32), dimension(ndoub, nbin) :: & idstar = reshape((/1,3,2,4/), (/2,2/)) real(kind=real64), dimension(nbin,nbin, ndoub) :: & sgn_f = reshape((/1.d0,1.d0,-1.d0,-1.d0,-1.d0,1.d0,-1.d0,1.d0/), & (/nbin,nbin, ndoub/)) real(kind=real64), dimension(2) :: & sgni = (/-1.d0, 1.d0/) ! local variables real(kind=real64), dimension(stars(1)%m, nstar) :: & star_v integer(kind=int32) :: & i, j, k, i1, j1, k1 real(kind=real64) :: & t, f, g, h, b real(kind=real64), dimension(ndim) :: & v, ff, ft real(kind=real64), dimension(ndim, norb) :: & r_r, v_r real(kind=real64), dimension(ndim, nbin, ndoub) :: & j_p, w_p, & mc_p, fsd_p, ftd_p, ftf_p, ftd3_p, torque_p, wxr_p, & mv_p real(kind=real64), dimension(ndim, nbin, nbin) :: & r_c, b_c, v_c, fpn_c, re_c real(kind=real64), dimension(nbin, nbin) :: & m_c, phi_c, rnm1_c, rnm2_c, rnm3_c, v2_c, eta_c, vr2_c, vrn_c, & rnm4_c, ftgr_c real(kind=real64), dimension(ndim, ndoub) :: & b_d, fng_d, acc_d, fpn_d, & re_d, f_d, w_d real(kind=real64), dimension(nbin, ndoub) :: & m_p, s_p, k_p, i_p, q_p, & mi_p, wn2_p, qi_p, & mr_p, mrat_p, mf_p, s5k_p, wre_p real(kind=real64), dimension(nbin, nbin) :: & mp_c, rn_c real(kind=real64), dimension(ndoub) :: & m_d, mm1_d, mu_d, & rn_d, rnm1_d, rnm2_d, rnm3_d, rnm4_d, rnm7_d, rnm8_d, & vrn_d, v2_d, Gm_d, an_d, phi_d, eps_d, nm1_d, & eta_d, vr2_d, ftgr_d real(kind=real64), dimension(ndim) :: & re_q, fng_q, acc_q, fpn_q, & rij real(kind=real64) :: & m_q, mi_q, eta_q, mp_q, phi_q, & rn2_q, rnm1_q, rnm2_q, rn_q, v2_q, vrn_q, & vr2_q, mpi_q, mui_q, & rnm1_ij, rnm2_ij, rnm3_ij real(kind=real64), dimension(ndim, nbin, nbin, ndoub) :: & ftdt_f, fsdt_f, ftft_f, wxr_f real(kind=real64), dimension(nbin, nbin, ndoub) :: & wre_f ! [[star1, star2], [star3, star4]] ! VARIABLE NAMES ! _v - values for 4 stars as last index ( ..., 4) ! _r - values for 3 relative coordinates ! last coordinate is for outer double ( ..., 3) ! _p - pair of 2-star values (..., 2, 2) ! _c - pair of 2-star coordinate (..., 2, 2) ! _d - pair of bianry values (..., 2) ! _q - outer double values (..., 1) ! _f - full combinations (...., 2, 2, 2) ! if (.not.all(ieee_is_finite(y))) then ! yp(:) = ieee_value(0.d0,ieee_quiet_nan) ! return ! endif r_r(:,:) = reshape(y( 1: 9), (/3, 3/)) v_r(:,:) = reshape(y(10:18), (/3, 3/)) j_p(:,:,:) = reshape(y(19:30), (/3, 2, 2/)) t = trel + toffset do i=1,4 star_v(:,i) = star_huntpol(t, i) end do m_p = reshape(star_v(IMASS, :), (/2, 2/)) s_p = reshape(star_v(IRADI, :), (/2, 2/)) k_p = reshape(star_v(IAPSI, :), (/2, 2/)) i_p = reshape(star_v(IINER, :), (/2, 2/)) q_p = reshape(star_v(IQUAL, :), (/2, 2/)) mi_p(:,:) = 1.d0 / m_p(:,:) qi_p(:,:) = 1.d0 / q_p(:,:) s5k_p(:,:) = s_p(:,:)**5 * k_p(:,:) mr_p(:,:) = m_p(2:1:-1,:) mrat_p(:,:) = mr_p(:,:) * mi_p(:,:) m_d(:) = sum(m_p, dim=1) mm1_d(:) = 1.d0 / m_d(:) mu_d (:)= product(m_p(:,:), dim=1) * mm1_d(:) do i=1, nbin do j=1, ndoub if (i_p(i,j) >0.d0) then w_p(:,i,j) = j_p(:,i,j) / i_p(i,j) else w_p(:,i,j) = 0.d0 end if end do end do wn2_p(:,:) = sum(w_p(:,:,:)**2, dim=1) do j=1,ndoub rn_d(j) = norm2(r_r(:,j)) enddo rnm1_d(:) = 1.d0 / rn_d(:) rnm2_d(:) = rnm1_d(:)**2 rnm3_d(:) = rnm2_d(:) * rnm1_d(:) rnm4_d(:) = rnm2_d(:)**2 rnm7_d(:) = rnm4_d(:) * rnm3_d(:) rnm8_d(:) = rnm4_d(:)**2 do j=1, ndoub re_d(:,j) = r_r(:,j) * rnm1_d(j) enddo ! fail on star interaction if (has_interact) then do j=1,ndoub if (rn_d(j) < interact(idstar(j,1), idstar(j,2))) then write(6, "(' [direct4p] object interaction ',i1,' + ',i1)") idstar(j,:) yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = idstar(j,:) terminate_integration = .true. return endif end do endif do j=1,ndoub ! fail on binary merger if (rn_d(j) < sum(s_p(:,j))) then write(6, "(' [direct4p] object collision ',i1,' + ',i1)") idstar(j,:) yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = idstar(j,:) terminate_integration = .true. return endif enddo ! fail on binary escape if (has_cutoff) then do j=1,norb if (rn_d(j) > cutoff(j)) then write(6, "(' [direct4p] orbit escape ',i1)") j yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/j, 0/) terminate_integration = .true. return endif enddo endif if (use_sd .or. use_tf) then do j=1,ndoub do i=1,nbin wre_p(i,j) = dot_product(w_p(:,i,j), re_d(:,j)) wxr_p(:,i,j) = cross_product(w_p(:,i,j), r_r(:,j)) enddo enddo end if do j=1,ndoub v2_d(j) = sum(v_r(:,j)**2) enddo Gm_d(:) = GRAV * m_d(:) phi_d(:) = Gm_d(:) * rnm1_d(:) do j=1,ndoub vrn_d(j) = dot_product(v_r(:,j), re_d(:,j)) enddo ! gravity m_q = sum(m_d) mi_q = 1.d0 / m_q mp_q = product(m_d) mpi_q = 1.d0 / mp_q mui_q = m_q * mpi_q do i=1,nbin do j=1,ndoub mp_c(i,j) = m_p(i,1) * m_p(j,2) mf_p(i,j) = mr_p(i,j) * mm1_d(j) enddo enddo mf_p(2,1) = - mf_p(2,1) mf_p(1,2) = - mf_p(1,2) do i=1,nbin do j=1,ndoub mc_p(:,i,j) = mf_p(i,j) * r_r(:,j) enddo enddo do i=1,nbin do j=1,nbin r_c(:,i,j) = r_r(:,3) + mc_p(:,i,1) + mc_p(:,j,2) rn_c(i,j) = norm2(r_c(:,i,j)) rnm1_c(i,j) = 1.d0 / rn_c(i,j) rnm2_c(i,j) = rnm1_c(i,j)**2 rnm3_c(i,j) = rnm1_c(i,j) * rnm2_c(i,j) enddo enddo ! fail on star interaction if (has_interact) then do i=1,nbin do j=1,ndoub if (rn_c(i,j) < interact(i,j+2)) then write(6, "(' [direct4p] object interaction ',i1,' + ',i1)") i, j+2 yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = (/i, j+2/) terminate_integration = .true. return endif end do end do endif ! fail on star collisions do i=1,nbin do j=1,ndoub if (rn_c(i,j) < s_p(i,1) + s_p(j,2)) then write(6, "(' [direct4p] object collision ',i1,' + ',i1)") i, j+2 yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/i, j+2/) terminate_integration = .true. return endif enddo enddo rn_q = norm2(r_r(:,3)) do i=1,nbin do j=1,nbin b_c(:,i,j) = r_c(:,i,j) * rnm3_c(i,j) enddo enddo do j=1,ndoub b_d(:,j) = r_r(:,j) * rnm3_d(j) enddo fng_d(:,1) = GRAV * ( & m_p(1,2) * (b_c(:,2,1) - b_c(:,1,1)) & + m_p(2,2) * (b_c(:,2,2) - b_c(:,1,2)) & - m_d(1) * b_d(:,1) & ) fng_d(:,2) = GRAV * ( & m_p(1,1) * (b_c(:,1,1) - b_c(:,1,2)) & + m_p(2,1) * (b_c(:,2,1) - b_c(:,2,2)) & - m_d(2) * b_d(:,2) & ) fng_q(:) = (- GRAV * mui_q) * ( & mp_c(1,1) * b_c(:,1,1) & + mp_c(1,2) * b_c(:,1,2) & + mp_c(2,1) * b_c(:,2,1) & + mp_c(2,2) * b_c(:,2,2) & ) acc_d(:,:) = fng_d(:,:) acc_q(:) = fng_q(:) f_d(:,:) = 0.d0 torque_p(:,:,:) = 0.d0 ! first order post-Newtonian GR corrections if (use_pn) then eta_d(:) = mu_d(:) * mm1_d(:) vr2_d(:) = vrn_d(:)**2 do j=1,ndoub fpn_d(:,j) = & phi_d(j) * rnm1_d(j) * CLIGHTm2 * ( & v_r(:,j) * ( & ( 4.d0 - 2.d0 * eta_d(j)) * vrn_d(j)) & + re_d(:,j) * ( & + ( 4.d0 + 2.d0 * eta_d(j)) * phi_d(j) & + 1.5d0 * eta_d(j) * vr2_d(j) & + (-1.d0 - 3.d0 * eta_d(j)) * v2_d(j)) & ) f_d(:,j) = f_d(:,j) + fpn_d(:,j) ftgr_d(j) = mu_d(j) * (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_d(j) * (1.d0 - 3.d0 * eta_d(j)) & + (3.d0 + eta_d(j)) * phi_d(j))) end do else ftgr_d(:) = mu_d(:) fpn_d(:,:) = 0.d0 endif if (use_pnq) then if (use_pnf) stop ' [direct4p] cannot use use_pnf and use_pnq' rn2_q = rn_q**2 rnm1_q = 1.d0 / rn_q rnm2_q = rnm1_q**2 re_q = r_r(:,3) * rnm1_q eta_q = mp_q * mi_q**2 vrn_q = dot_product(v_r(:,3), re_q(:)) vr2_q = vrn_q**2 v2_q = sum(v_r(:,3)**2) phi_q = GRAV * m_q * rnm1_q fpn_q(:) = (Gcm2 * m_q * rnm2_q) * ( & v_r(:,3) * ( & ( 4.d0 - 2.d0 * eta_q) * vrn_q) & + re_q(:) * ( & + ( 4.d0 + 2.d0 * eta_q) * phi_q & + 1.5d0 * eta_q * vr2_q & + (-1.d0 - 3.d0 * eta_q) * v2_q) & ) acc_q(:) = acc_q(:) + fpn_q(:) else fpn_q(:) = 0.d0 endif ! indidividual binaries SD, TD, and TF if (use_sd .or. use_td) then do j=1,ndoub if (use_td) then f = -12.d0 * GRAV * rnm3_d(j) endif g = m_d(j) * rnm4_d(j) do i=1,nbin h = s5k_p(i,j) * mi_p(i,j) * g if (use_td) then ftd_p(:,i,j) = (h * f * mr_p(i,j)) * re_d(:,j) f_d(:,j) = f_d(:,j) + ftd_p(:,i,j) else ftd_p(:,i,j) = 0.d0 end if if (use_sd) then b = -2.d0 * h * wre_p(i,j) fsd_p(:,i,j) = & (h * (5.d0 * wre_p(i,j)**2 - wn2_p(i,j))) * re_d(:,j) & + b * w_p(:,i,j) f_d(:,j) = f_d(:,j) + fsd_p(:,i,j) torque_p(:,i,j) = torque_p(:,i,j) & + (b * ftgr_d(j)) * wxr_p(:,i,j) else fsd_p(:,i,j) = 0.d0 end if end do end do else fsd_p(:,:,:) = 0.d0 ftd_p(:,:,:) = 0.d0 endif if (use_tf) then do j=1,ndoub w_d(:,j) = cross_product(re_d(:,j), v_r(:,j)) * rnm1_d(j) enddo if (use_rxv) then do j=1,2 nm1_d(j) = 1.d0 / norm2(w_d(:,j)) enddo else eps_d(:) = 0.5d0 * v2_d(:) - phi_d(:) an_d(:) = max(0.d0, -0.5d0 * Gm_d(:) / eps_d(:)) nm1_d(:) = sqrt(an_d(:)**3 / Gm_d(:)) end if do j=1,ndoub f = 6.d0 * Gm_d(j) * nm1_d(j) * rnm8_d(j) v(:) = (-2.d0 * vrn_d(j)) * re_d(:,j) - v_r(:,j) do i=1,nbin h = f * qi_p(i,j) * s5k_p(i,j) * mrat_p(i,j) ftf_p(:,i,j) = h * (v(:) + wxr_p(:,i,j)) f_d(:,j) = f_d(:,j) + ftf_p(:,i,j) b = h * ftgr_d(j) * rn_d(j) torque_p(:,i,j) = torque_p(:,i,j) & + (b * wre_p(i,j)) * r_r(:,j) & + (b * rn_d(j)) * (w_d(:,J) - w_p(:,i,j)) end do end do else ftf_p(:,:,:) = 0.d0 endif ! impact of other binary CM gravity on tides in a binary if (use_td3) then do j=1,ndoub b = m_d(j) * GRAV * rnm4_d(j) * rnm1_d(j) do i=1, nbin if (j == 1) then rij(:) = r_r(:,3) + mc_p(:,i,1) else rij(:) = -r_r(:,3) - mc_p(:,2,i) endif rnm1_ij = 1.d0 / norm2(rij) rnm2_ij = rnm1_ij**2 rnm3_ij = rnm2_ij * rnm1_ij g = dot_product(r_r(:,j), rij(:)) f = b * s5k_p(i,j) * rnm3_ij h = 2.d0 * f * g * rnm2_ij ftd3_p(:,i,j) = & ((3.d0 + 11.d0 * g**2 * rnm2_ij * rnm2_d(j)) * f) * r_r(:,j) + & h * rij(:) f_d(:,j) = f_d(:,j) + ftd3_p(:,i,j) torque_p(:,i,j) = torque_p(:,i,j) & + (h * ftgr_d(j)) * cross_product(rij(:), r_r(:,j)) enddo enddo else ftd3_p(:,:,:) = 0.0d0 endif if (use_pnf .or. use_sdt .or. use_tdt .or. use_tft) then do i=1,nbin do j=1,nbin re_c(:,i,j) = r_c(:,i,j) * rnm1_c(i,j) enddo enddo endif if (use_pnf) then if (.not. use_pn) stop ' [direct4p] require use_pn' do i=1,nbin do j=1,ndoub mv_p(:,i,j) = mf_p(i,j) * v_r(:,j) enddo enddo do i=1,nbin do j=1,nbin m_c(i,j) = m_p(i,1) + m_p(j,2) eta_c(i,j) = mp_c(i,j) / m_c(j,j)**2 v_c(:,i,j) = v_r(:,3) + mv_p(:,i,1) + mv_p(:,j,2) vrn_c(i,j) = dot_product(v_c(:,i,j), re_c(:,i,j)) vr2_c(i,j) = vrn_c(i,j)**2 v2_c(i,j) = sum(v_c(:,i,j)**2) phi_c(i,j) = GRAV * m_c(i,j) * rnm1_c(i,j) fpn_c(:,i,j) = GRAV * mp_c(i,j) * rnm2_c(i,j) * CLIGHTm2 * ( & v_c(:,i,j) * ( & ( 4.d0 - 2.d0 * eta_c(i,j)) * vrn_c(i,j)) & + re_c(:,i,j) * ( & + ( 4.d0 + 2.d0 * eta_c(i,j)) * phi_c(i,j) & + 1.5d0 * eta_c(i,j) * vr2_c(i,j) & + (-1.d0 - 3.d0 * eta_c(i,j)) * v2_c(i,j))& ) acc_q(:) = acc_q(:) + mui_q * fpn_c(:,i,j) acc_d(:,1) = acc_d(:,1) + mi_p(i,1) * fpn_c(:,i,j) * sgni(3-i) acc_d(:,2) = acc_d(:,2) + mi_p(i,2) * fpn_c(:,i,j) * sgni(i) ! for recoding if (local_store) then fpn_c(:,i,j) = fpn_c(:,i,j) * (m_c(i,j) / mp_c(i,j)) endif if ((use_sdt .or. use_tdt .or. use_tft) .and. use_pnt) then ftgr_c(i,j) = (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_c(i,j) * (1.d0 - 3.d0 * eta_c(i,j)) & + (3.d0 + eta_c(i,j)) * phi_c(i,j))) else ftgr_c(i,j) = 1.d0 end if enddo enddo else fpn_c(:,:,:) = 0.d0 ftgr_c(:,:) = 1.d0 endif ! tide interaction on resolved star-basis (may be overkill) if (use_sdt .or. use_tdt .or. use_tft) then ! TODO - add star-wise tide-orbit feedback? ! what to do about triple star interaction in this case? stop ' [direct4p] not implemented: use_sdt, use_tdt, use_tft' do i=1,nbin do j=1, nbin rnm4_c(i,j) = rnm3_c(i,j) * rnm1_c(i,j) enddo enddo if (use_sd) then do k=1, ndoub do i=1, nbin do j=1,nbin wre_f(i,j,k) = dot_product(w_p(:,i,k), re_c(:,i,j)) wxr_f(:,i,j,k) = cross_product(w_p(:,i,k), r_c(:,i,j)) enddo enddo enddo endif if (use_tft) then if (.not. use_rxv) & error stop ' [direct4p] use_tft requires use_rxv' ! compute nm1_c endif if (use_pnt .and. (.not. use_pnf)) then ! need to fill in hierarchical calculation ftgr_c(:,:) = 1.d0 end if do k=1, ndoub ! k to correspond to i, k1 to j k1 = 3 - k do i=1, nbin do j=1, nbin if (k == 1) then i1 = i j1 = j else i1 = j j1 = i endif f = -12.d0 * GRAV * rnm3_c(i1,j1) * m_p(j,k1) h = s5k_p(i,k) * m_p(j,k1) * rnm4_c(i1,j1) if (use_tdt) then ff(:) = (h * f) * re_c(:,i1,j1) ftdt_f(:,i1,j1,k) = ff(:) * mi_p(j,k1) ft(:) = ff(:) else ftdt_f(:,i1,j1,k) = 0.d0 ft(:) = 0.d0 end if if (use_sd) then b = -2.d0 * h * wre_f(i1,j1,k) ff(:) = & (h * (5.d0 * wre_f(i1,j1,k)**2 - wn2_p(i,k))) * re_c(:,i1,j1) & + b * w_p(:,i,k) fsdt_f(:,i1,j1,k) = ff(:) * mi_p(j,k1) ft(:) = ft(:) + ff(:) torque_p(:,i,k) = torque_p(:,i,k) & + (b * ftgr_c(i1,j1)) * wxr_f(:,i1,j1,k) else fsdt_f(:,i1,j1,k) = 0.d0 end if acc_d(:,k ) = acc_d(:,k ) + ft(:) * (mi_p(i,k ) * sgn_f(i1,j1,k )) acc_d(:,k1) = acc_d(:,k1) + ft(:) * (mi_p(i,k1) * sgn_f(i1,j1,k1)) acc_q(:) = acc_q(:) + ft(:) * mui_q enddo enddo enddo else ftdt_f(:,:,:,:) = 0.d0 fsdt_f(:,:,:,:) = 0.d0 ftft_f(:,:,:,:) = 0.d0 endif if (use_t3t) then ! 3-tide of individual star in other binary stop ' [direct4p] not implemented: use_t3t' endif if (use_t3h) then ! 3-tide of binay on cms of other binary stop ' [direct4p] not implemented: use_t3h' endif if (use_t3s) then ! 3-tide of binay on individual star in other binary stop ' [direct4p] not implemented: use_t3s' endif acc_d(:,:) = acc_d(:,:) + f_d(:,:) yp( 1: 9) = reshape(v_r(1:3,1:3), (/9/)) yp(10:15) = reshape(acc_d(1:3,1:2), (/6/)) yp(16:18) = acc_q(1:3) yp(19:30) = reshape(torque_p(1:3,1:2,1:2), (/12/)) if (local_store) then local_data(1:30) = yp(1:30) local_data(31:42) = reshape(fsd_p, (/12/)) local_data(43:54) = reshape(ftd_p, (/12/)) local_data(55:66) = reshape(ftf_p, (/12/)) local_data(67:72) = reshape(fpn_d(1:3,1:2), (/6/)) local_data(73:84) = reshape(fpn_c, (/12/)) local_data(85:87) = fpn_q(1:3) local_data(88:99) = reshape(ftd3_p, (/12/)) endif end function direct4p ! ======================================================================= ! Direct code, hierarchical ! ======================================================================= function direct4h(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, IQUAL, & 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_OK, 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_rxv, & use_pnf, use_pnh, use_tfh, use_sdh, use_tdh, use_td3, & use_sdt, use_tdt, use_tft, use_pnt 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 integer(kind=int32), parameter :: & nstar = 4, & norb = 3, & ndim = 3, & nbin = 2 real(kind=real64), dimension(nbin), parameter :: & sign_b1 = (/1.d0, -1.d0/) integer(kind=int32) :: & i, j, & imax, jmax real(kind=real64) :: & t, b, h, f, g, f1, g1, & rnm2_ij, rnm3_ij, rnm4_ij, rnm8_ij, nm1_ij integer(kind=int32), dimension(nstar), parameter :: & v2b_v = (/(max(1, j-1),j=1,nstar)/) real(kind=real64), dimension(stars(1)%m, nstar) :: & star_v real(kind=real64), dimension(nstar) :: & m_v,s_v,k_v,i_v,q_v, & mi_v, wn2_v, s5k_v, qi_v, & wre_v, mr_v, mrat_v real(kind=real64), dimension(ndim, nstar) :: & j_v, w_v, torque_v, & fsd_v, ftd_v, ftf_v, wxr_v real(kind=real64), dimension(ndim, norb) :: & r_b, v_b, re_b, & acc_b, fpn_b, f_b, w_b real(kind=real64), dimension(norb) :: & rn_b, rnm1_b, rnm2_b, rnm3_b, rnm4_b, rnm8_b, & m_b, mi_b, mu_b, v2_b, Gm_b, phi_b, & nm1_b, eps_b, an_b, & vrn_b, eta_b, vr2_b, ftgr_b real(kind=real64), dimension(nbin, norb) :: & m_p real(kind=real64), dimension(ndim, nstar) :: & rx_x, vx_x real(kind=real64), dimension(nstar) :: & mf_x real(kind=real64), dimension(ndim, nstar-1, 2:nstar) :: & r_c, b_c, & v_c, re_c, fpn_c real(kind=real64), dimension(nstar-1, 2:nstar) :: & rn_c, rnm1_c, rnm2_c, rnm3_c, & m_c, phi_c, eta_c, vrn_c, vr2_c, v2_c real(kind=real64), dimension(ndim) :: & v, rij, ff, ft, & w_ij real(kind=real64), dimension(ndim,nbin,nbin+1:nstar) :: & ftd3_c real(kind=real64), dimension(ndim,nstar-1,nbin+1:nstar) :: & ftdt_c, fsdt_c, ftft_c, & wxr_c real(kind=real64), dimension(nstar-1,nbin+1:nstar) :: & wre_c, ftgr_c, mu_c real(kind=real64), dimension(2:norb) :: & mui_b ! [[[star1, star2], star3], star4] ! VARIABLE NAMES ! _v - values for stars as last index (nstar) ! _b - values for binary systems (norb) ! _p - pairs of values for binary systems (2,norb) ! _c - coordinate pairs ! _x - other data r_b(:,:) = reshape(y(1 : ndim*norb), (/ndim, norb/)) v_b(:,:) = reshape(y(1+ ndim*norb:2*ndim*norb), (/ndim, norb/)) j_v(:,:) = reshape(y(1+2*ndim*norb:ndim*(2*norb+nstar)), (/ndim, nstar/)) t = trel + toffset do i=1,nstar star_v(:,i) = star_huntpol(t, i) end do m_v(:) = star_v(IMASS, :) s_v(:) = star_v(IRADI, :) k_v(:) = star_v(IAPSI, :) i_v(:) = star_v(IINER, :) q_v(:) = star_v(IQUAL, :) mi_v(:) = 1.d0 / m_v(:) qi_v(:) = 1.d0 / q_v(:) s5k_v(:) = s_v(:)**5 * k_v(:) m_b(1) = sum(m_v(1:2)) m_p(1:2,1) = m_v(1:2) m_p(2,2:norb) = m_v(3:nstar) do j=2, norb m_b(j) = m_b(j-1) + m_v(j+1) m_p(1,j) = m_b(j-1) enddo mi_b(:) = 1.d0 / m_b(:) mu_b(:) = product(m_p(1:2,:), dim=1) * mi_b(:) mr_v(1:2) = m_v(2:1:-1) mr_v(3) = sum(mr_v(1:2)) do i=4,nstar mr_v(i) = mr_v(i-1) + m_v(i-1) enddo mrat_v(:) = mr_v(:) * mi_v(:) do i=1,nstar if (i_v(i) > 0.d0) then w_v(:,i) = j_v(:,i) / i_v(i) else w_v(:,i) = 0.d0 endif enddo wn2_v = sum(w_v**2, dim=1) do j=1,norb rn_b(j) = norm2(r_b(:,j)) enddo ! fail on binary escape if (has_cutoff) then do j=1,norb if (rn_b(j) > cutoff(j)) then write(6, "(' [direct4h] orbit escape ',i1)") j yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/j, 0/) terminate_integration = .true. return endif enddo end if rnm1_b(:) = 1.d0 / rn_b(:) rnm2_b(:) = rnm1_b(:)**2 rnm3_b(:) = rnm2_b(:) * rnm1_b(:) rnm4_b(:) = rnm2_b(:)**2 rnm8_b(:) = rnm4_b(:)**2 do j=1,norb re_b(:,j) = r_b(:,j) * rnm1_b(j) enddo do i=1,nstar j = v2b_v(i) wre_v( i) = dot_product(w_v(:,i), re_b(:,j)) wxr_v(:,i) = cross_product(w_v(:,i), r_b(:,j)) end do v2_b(:) = sum(v_b(:,:)**2, dim=1) Gm_b(:) = GRAV * m_b(:) phi_b(:) = Gm_b(:) * rnm1_b(:) do j=1,norb vrn_b(j) = dot_product(v_b(:,j), re_b(:,j)) enddo f_b(:,:) = 0.d0 torque_v(:,:) = 0.d0 ! gravity mf_x(1) = m_v(2) * mi_b(1) mf_x(2) = -m_v(1) * mi_b(1) mf_x(3) = m_v(3) * mi_b(2) mf_x(4) = -m_b(1) * mi_b(2) rx_x(:,1) = mf_x(1) * r_b(:,1) rx_x(:,2) = mf_x(2) * r_b(:,1) rx_x(:,3) = mf_x(3) * r_b(:,2) + r_b(:,3) r_c(:,1,2) = r_b(:,1) r_c(:,1,3) = rx_x(:,1) + r_b(:,2) r_c(:,2,3) = rx_x(:,2) + r_b(:,2) r_c(:,1,4) = rx_x(:,1) + rx_x(:,3) r_c(:,2,4) = rx_x(:,2) + rx_x(:,3) r_c(:,3,4) = mf_x(4) * r_b(:,2) + r_b(:,3) do i=1,nstar-1 do j=i+1,nstar rn_c(i,j) = norm2(r_c(:,i,j)) ! fail on star interaction if (has_interact.and.(rn_c(i,j) < interact(i,j))) then write(6, "(' [direct4h] object interaction ',i1,' + ',i1)") i,j yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = (/i, j/) terminate_integration = .true. return endif ! fail on binary merger if (rn_c(i,j) < s_v(i) + s_v(j)) then write(6, "(' [direct4h] object collision ',i1,' + ',i1)") i, j yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/i, j/) terminate_integration = .true. return endif rnm1_c(i,j) = 1.d0 / rn_c(i,j) rnm2_c(i,j) = rnm1_c(i,j) **2 rnm3_c(i,j) = rnm2_c(i,j) * rnm1_c(i,j) b_c(:,i,j) = r_c(:,i,j) * rnm3_c(i,j) enddo enddo rx_x(:,4) = & + m_v(1) * b_c(:,1,4) & + m_v(2) * b_c(:,2,4) & + m_v(3) * b_c(:,3,4) acc_b(:,1) = (- GRAV) * ( & m_b(1) * b_c(:,1,2) & + m_v(3) * (b_c(:,1,3) - b_c(:,2,3)) & + m_v(4) * (b_c(:,1,4) - b_c(:,2,4))) acc_b(:,2) = (GRAV * mi_b(1)) * ( & m_b(2) * ( & - m_v(1) * b_c(:,1,3) & - m_v(2) * b_c(:,2,3)) & + m_v(4) * ( & - m_v(1) * b_c(:,1,4) & - m_v(2) * b_c(:,2,4))) & + (GRAV * m_v(4)) * b_c(:,3,4) acc_b(:,3) = (- GRAV) * m_b(3) * mi_b(2) * rx_x(:,4) if ((use_pn .and. use_pnf) .or. use_tft) then vx_x(:,1) = mf_x(1) * v_b(:,1) vx_x(:,2) = mf_x(2) * v_b(:,1) vx_x(:,3) = mf_x(3) * v_b(:,2) + v_b(:,3) v_c(:,1,2) = v_b(:,1) v_c(:,1,3) = rx_x(:,1) + v_b(:,2) v_c(:,2,3) = rx_x(:,2) + v_b(:,2) v_c(:,1,4) = rx_x(:,1) + vx_x(:,3) v_c(:,2,4) = rx_x(:,2) + vx_x(:,3) v_c(:,3,4) = mf_x(4) * v_b(:,2) + v_b(:,3) do i=1,nstar-1 do j=i+1, nstar re_c(:,i,j) = r_c(:,i,j) * rnm1_c(i,j) vrn_c(i,j) = dot_product(re_c(:,i,j), v_c(:,i,j)) enddo enddo endif ! GR if (use_pn) then if (use_pnf) then do i=1,nstar-1 do j=i+1, nstar m_c(i,j) = m_v(i) + m_v(j) phi_c(i,j) = GRAV * m_c(i,j) eta_c(i,j) = m_v(i) * m_v(j) / m_c(i,j)**2 vr2_c(i,j) = vrn_c(i,j)**2 v2_c(i,j) = sum(v_c(:,i,j)**2) fpn_c(:,i,j) = & phi_c(i,j) * rnm1_c(i,j) * CLIGHTm2 * ( & v_c(:,i,j) * ( & ( 4.d0 - 2.d0 * eta_c(i,j)) * vrn_c(i,j)) & + re_c(:,i,j) * ( & + ( 4.d0 + 2.d0 * eta_c(i,j)) * phi_c(i,j) & + 1.5d0 * eta_c(i,j) * vr2_c(i,j) & + (-1.d0 - 3.d0 * eta_c(i,j)) * v2_c(i,j) ) & ) if (use_sdt .or. use_tdt .or. use_tft) then if (j < nbin+1) cycle mu_c(i,j) = eta_c(i,j) * m_c(i,j) if (use_pnt) then ftgr_c(i,j) = mu_c(i,j) * (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_c(i,j) * (1.d0 - 3.d0 * eta_c(i,j)) & + (3.d0 + eta_c(i,j)) * phi_c(i,j))) else ftgr_c(i,j) = mu_c(i,j) endif endif end do end do else if (use_pnh .or. use_pnt) then jmax = norb else jmax = 1 endif eta_b(1:jmax) = mu_b(1:jmax) * mi_b(1:jmax) vr2_b(1:jmax) = vrn_b(1:jmax)**2 do j=1,jmax fpn_b(:,j) = & phi_b(j) * rnm1_b(j) * CLIGHTm2 * ( & v_b(:,j) * ( & ( 4.d0 - 2.d0 * eta_b(j)) * vrn_b(j)) & + re_b(:,j) * ( & + ( 4.d0 + 2.d0 * eta_b(j)) * phi_b(j) & + 1.5d0 * eta_b(j) * vr2_b(j) & + (-1.d0 - 3.d0 * eta_b(j)) * v2_b(j) ) & ) f_b(:,j) = f_b(:,j) + fpn_b(:,j) ftgr_b(j) = mu_b(j) * (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_b(j) * (1.d0 - 3.d0 * eta_b(j)) & + (3.d0 + eta_b(j)) * phi_b(j))) end do do j = jmax+1,norb fpn_b(:,j) = 0.d0 ftgr_b(j) = mu_b(j) enddo end if else fpn_b(:,:) = 0.d0 ftgr_b(:) = mu_b(:) endif if (use_sd .or. use_td) then if (use_sdh .or. use_tdh) then imax = nstar else imax = nbin endif do i=1,imax j = v2b_v(i) if (i /= 2) then g = m_b(j) * rnm4_b(j) endif h = s5k_v(i) * mi_v(i) * g if ((use_td) .and. ((i <= nbin) .or. use_tdh)) then if (i /= 2) then f = (-12.d0 * GRAV) * rnm3_b(j) endif ftd_v(:,i) = (f * h * mr_v(i)) * re_b(:,j) f_b(:,j) = f_b(:,j) + ftd_v(:,i) else ftd_v(:,j) = 0.0d0 end if if ((use_sd) .and. ((i <= nbin) .or. use_sdh)) then b = -2.d0 * h * wre_v(i) fsd_v(:,i) = & (h * (5.d0 * wre_v(i)**2 - wn2_v(i))) * re_b(:,j) & + b * w_v(:,i) f_b(:,j) = f_b(:,j) + ftd_v(:,i) torque_v(:,i) = torque_v(:,i) + & (b * ftgr_b(j)) * wxr_v(:,i) else fsd_v(:,i) = 0.0d0 end if end do fsd_v(:,imax+1:nstar) = 0.d0 ftd_v(:,imax+1:nstar) = 0.d0 else fsd_v(:,:) = 0.d0 ftd_v(:,:) = 0.d0 endif if (use_tf) then if (use_tfh) then imax = nstar jmax = norb else imax = nbin jmax = 1 endif do j=1,jmax w_b(:,j) = cross_product(v_b(:,j), re_b(:,j)) * rnm1_b(j) enddo if (use_rxv) then do j=1, jmax nm1_b(j) = 1.d0 / norm2(w_b(:,j)) enddo else do j=1, jmax eps_b(j) = 0.5d0 * v2_b(j) - phi_b(j) an_b(j) = max(0.d0, -0.5d0 * Gm_b(j) / eps_b(j)) nm1_b(j) = sqrt(an_b(j)**3 / Gm_b(j)) end do endif do i=1,imax j = v2b_v(i) if (i /= 2) then f = 6.d0 * Gm_b(j) * nm1_b(j) * rnm8_b(j) v(:) = (-2.d0 * vrn_b(j)) * re_b(:,j) - v_b(:,j) end if h = f * qi_v(i) * s5k_v(i) * mrat_v(i) ftf_v(:,i) = h * (v(:) + wxr_v(:,i)) f_b(:,j) = f_b(:,j) + ftf_v(:,i) b = h * ftgr_b(j) * rn_b(j) torque_v(:,i) = torque_v(:,i) & + (b * wre_v(i)) * r_b(:,j) & + (b * rn_b(j)) * (w_b(:,j) - w_v(:,i)) end do ftf_v(:,imax+1:nstar) = 0.d0 else ftf_v(:,:) = 0.d0 endif ! impact of outer binaries on tides in inner binary ! TODO - add 3-star impact 4-->3-->(1,2) ? if (use_td3) then b = m_b(1) * GRAV * rnm4_b(1) * rnm1_b(1) do j=3, nstar do i=1, 2 rij(:) = r_c(:,i,j) rnm2_ij = rnm2_c(i,j) rnm3_ij = rnm3_c(i,j) g = dot_product(r_b(:,1), rij(:)) f = b * s5k_v(i) * rnm3_ij h = 2.d0 * f * g * rnm2_ij ftd3_c(:,i,j) = & ((3.d0 + 11.d0 * g**2 * rnm2_ij * rnm2_b(1)) * f) * r_b(:,1) + & h * rij(:) f_b(:,1) = f_b(:,1) + ftd3_c(:,i,j) torque_v(:,i) = torque_v(:,i) & + (h * ftgr_b(1)) * cross_product(rij(:), r_b(:,1)) enddo enddo else ftd3_c(:,:,:) = 0.d0 endif ! impact of tides in inner binaries due to 3rd and 4th star on their ! star orbits (SD, TD, TF) and corresponding torques on 3rd/4th object if (use_sdt .or. use_tdt .or. use_tft) then ! The direct action of individual stars, may be ! excessive relative to full model. ! But maybe indspensible simple extension of direct3h if (use_tdh .and. use_tdt) STOP 'conflict use_tdh and use_tdt' if (use_sdh .and. use_sdt) STOP 'conflict use_sdh and use_sdt' if (use_tfh .and. use_tft) STOP 'conflict use_tfh and use_tft' mui_b(2:norb) = 1.d0 / mu_b(2:norb) if (use_sdt .or. use_tft) then do j=nbin+1, nstar do i=1,j-1 wre_c(i,j) = dot_product(w_v(:,j), r_c(:,i,j)) * rnm1_c(i,j) wxr_c(:,i,j) = cross_product(w_v(:,i), r_c(:,i,j)) enddo enddo endif if (use_pnt .and. .not. use_pnf) then do j=nbin+1, nstar ftgr_c(1:nstar-1,j) = ftgr_b(j-1) / mu_b(j-1) enddo endif if ((use_tft) .and. (.not. use_rxv)) stop '[tft] require use_rxv -- True.' do j=nbin+1, nstar if (use_tdt) then f = -12.d0 * GRAV * m_v(j) endif if (use_tft) then f1 = 6.d0 * GRAV * m_v(j)**2 endif do i=1,j-1 rnm4_ij = rnm2_c(i,j)**2 h = s5k_v(i) * rnm4_ij * m_v(j) if (use_tdt) then ft(:) = (h * f * rnm4_ij) * r_c(:,i,j) ftdt_c(:,i,j) = mi_v(j) * ft(:) else ftdt_c(:,i,j) = 0.d0 ft = 0.d0 endif if (use_sdt) then b = -2.d0 * h * wre_c(i,j) ff(:) = & (h * rnm1_c(i,j) * & (5.d0 * wre_c(i,j)**2 - wn2_v(i))) * r_c(:,i,j) & + b * w_v(:,i) fsdt_c(:,i,j) = mi_v(3) * ff(:) ft(:) = ft(:) + ff(:) torque_v(:,i) = torque_v(:,i) & + (b * ftgr_c(i,j)) * wxr_c(:,i,j) else fsdt_c(:,i,j) = 0.d0 endif if (use_tft) then w_ij(:) = cross_product(v_c(:,i,j), r_c(:,i,j)) * rnm2_c(i,j) nm1_ij = 1.d0 / norm2(w_ij(:)) rnm8_ij = rnm4_ij**2 g1 = -2.d0 * vrn_c(i,j) v(:) = re_c(:,i,j) * g1 - v_c(:,i,j) h = f1 * qi_v(i) * s5k_v(i) * nm1_ij * rnm8_ij ff(:) = h * (v(:) + wxr_c(:,i,j)) ftft_c(:,i,j)= mi_v(j) * ff(:) ft(:) = ft(:) + ff(:) b = h * ftgr_c(i,j) * rn_c(i,j) torque_v(:,i) = torque_v(:,i) & + (b * wre_c(i,j)) * r_c(:,i,j) & + (b * rn_c(i,j)) * (w_ij(:) - w_v(:,i)) else ftft_c(:,i,j) = 0.d0 endif ! impact on orbits is tricky ... f_b(:,1) = f_b(:,1) + (mi_v(i) * sign_b1(i)) * ff(:) if (j == 3) then f_b(:,2) = f_b(:,2) + mui_b(2) * ff(:) else if (i < 3) then f_b(:,2) = f_b(:,2) + mi_b(1) * ff(:) f_b(:,3) = f_b(:,3) + mui_b(3) * ff(:) else f_b(:,2) = f_b(:,2) - mi_v(3) * ff(:) f_b(:,3) = f_b(:,3) + mui_b(3) * ff(:) endif enddo enddo ! STOP ' [direct4h] not implemented: use_sdt, use_tdt, use_tft' else ftdt_c(:,:,:) = 0.d0 fsdt_c(:,:,:) = 0.d0 ftft_c(:,:,:) = 0.d0 endif acc_b(:,:) = acc_b(:,:) + f_b(:,:) yp( 1: 9) = reshape(v_b(:,:), (/ndim*norb/)) yp(10:18) = reshape(acc_b(:,:), (/ndim*norb/)) yp(19:30) = reshape(torque_v(:,:), (/ndim*nstar/)) if (local_store) then local_data( 1: 30) = yp(:) local_data( 31: 42) = reshape(fsd_v (:,:), (/ndim*nstar/)) local_data( 43: 54) = reshape(ftd_v (:,:), (/ndim*nstar/)) local_data( 55: 66) = reshape(ftf_v (:,:), (/ndim*nstar/)) local_data( 67: 75) = reshape(fpn_b (:,:), (/ndim*norb/)) local_data( 76: 93) = reshape(fpn_c (:,:,:), (/ndim*(nstar*(nstar-1)/2)/)) local_data( 94:105) = reshape(ftd3_c(:,:,:), (/ndim*nbin*(nstar-nbin)/)) local_data(121:126) = reshape(fsdt_c(:,1:2,3), (/ndim*nbin/)) local_data(127:135) = reshape(fsdt_c(:,1:3,4), (/ndim*(nbin+1)/)) local_data(106:111) = reshape(ftdt_c(:,1:2,3), (/ndim*nbin/)) local_data(112:120) = reshape(ftdt_c(:,1:3,4), (/ndim*(nbin+1)/)) local_data(136:141) = reshape(ftft_c(:,1:2,3), (/ndim*nbin/)) local_data(142:150) = reshape(ftft_c(:,1:3,4), (/ndim*(nbin+1)/)) endif end function direct4h ! ======================================================================= ! Direct code, full ! ======================================================================= function direct4f(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, IQUAL, & 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_OK, 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_td3, & use_rxv, use_qp implicit none integer(kind=int32), parameter :: & ndim = 3, & ncorn = 2, & nstar = 4, & norb = nstar - 1 integer(kind=int32), parameter :: & nstar1 = nstar - 1, & nstar2 = nstar - 2, & nedge = (nstar1 * nstar) / 2, & nvorb = ndim * norb, & nvstar = ndim * nstar ! stars belonging to each edge (ncorn=2) integer(kind=int32), dimension(ncorn, nedge), parameter :: & ies = reshape((/1,2, 1,3, 1,4, 2,3, 2,4, 3,4/), (/ncorn, nedge/)) ! edge 1 2 3 4 5 6 ! edges belongig to each star pair combination (symmetric) integer(kind=int32), dimension(nstar, nstar), parameter :: & ipe = reshape((/0, 1, 2, 3, & 1, 0, 4, 5, & 2, 4, 0, 6, & 3, 5, 6, 0 /), (/nstar, nstar/)) ! DEBUG - not currently used ! signs of edges belongig to each star pair combination (symmetric) ! real(kind=real64), dimension(nstar, nstar), parameter :: & ! spe = reshape((/ 0.d0, 1.d0, 1.d0, 1.d0, & ! -1.d0, 0.d0, 1.d0, 1.d0, & ! -1.d0,-1.d0, 0.d0, 1.d0, & ! -1.d0,-1.d0,-1.d0, 0.d0 /), (/nstar,nstar/)) ! seems to be basically sign(j-i) ! indices for td3 combinations integer(kind=int32), dimension(nstar, nstar1), parameter :: & iij = reshape((/2,1,1,1, 3,3,2,2, 4,4,4,3/), (/nstar, nstar1/)) integer(kind=int32), dimension(nstar, nstar1, nstar2), parameter :: & iijk = reshape((/3,3,2,2, 2,1,1,1, 2,1,1,1, & 4,4,4,3, 4,4,4,3, 3,3,2,2 /), (/nstar, nstar1, nstar2/)) 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 ! local variables integer(kind=int32) :: & i, j, k, l, ir, ij, ik, j1, k1 real(kind=real64) :: & t real(kind=real64), dimension(stars(1)%m, nstar) :: & star_v real(kind=real64), dimension(ndim, norb) :: & r_j, v_j, acc_j real(kind=real64), dimension(ndim, nstar) :: & j_v, w_v, torque_v real(kind=real64), dimension(nstar) :: & m_v, s_v, k_v, i_v, q_v, mi_v, qi_v, s5k_v, & wn2_v, m2_v, & s5kqi_v ! coordinate setup real(kind=real64), dimension(ndim, nedge) :: & r_e, v_e real(kind=real64) :: & m_12, mi_12, m_123, mi_123, m_1234 real(kind=real64), dimension(2) :: & mi_d real(kind=real64), dimension(2,2) :: & mc_y real(kind=real64), dimension(ndim, 2,2) :: & rc_y, vc_y ! forces real(kind=real64), dimension(nstar1, nstar) :: & gjk_x, mjk_x real(kind=real64), dimension(ndim, nedge) :: & f_e real(kind=real64), dimension(nedge) :: & rn_e, rn2_e, rnm1_e, rnm2_e, rnm3_e, rnm4_e, rnm5_e, rnm7_e, rnm8_e, & m_e, mp_e, mm2_e, fgr_e, ggr_e real(kind=real64) :: & f, g, h, & eta, phi, vn2, vr real(kind=real64), dimension(ndim) :: & v, w, rik, rij real(kind=real64), dimension(ncorn, nedge) :: & wdr_c real(kind=real64), dimension(ndim, ncorn, nedge) :: & wxr_c ! for recording real(kind=real64), dimension(ndim, ncorn, nedge) :: & fsd_c, ftd_c, ftf_c real(kind=real64), dimension(ndim, nstar, nstar1, nstar2) :: & ftd3_d real(kind=real64), dimension(ndim, nedge) :: & fpn_e, fng_e ! [[[star1, star2], star3], star4] or [[star1, star2], [star3, star4]] ! FUTURE - note on notations for hierarchies ! [[[,],],] +3,-1,-1,-1 [[[,],],] ! [[,],[,]] +2,-1,+1,-2 [[,],[,]] ! [,[,]] +1,+1,-2 ! VARIABLE NAMES ! _v - values for stars as last index (nstar) ! _j - values for relative star coordnates in Jacobi coordinates !? _x - helper arrays of arbitrary dimensions ! _e - array of edges ! _y - helper arrays for coordinate construction if (.not.use_rxv) ERROR STOP 'require rxv' r_j(:,:) = reshape(y(1 : ndim*norb), (/ndim, norb /)) v_j(:,:) = reshape(y(1+ ndim*norb:2*ndim*norb), (/ndim, norb /)) j_v(:,:) = reshape(y(1+2*ndim*norb:ndim*(2*norb+nstar)), (/ndim, nstar/)) t = trel + toffset do i=1,nstar star_v(:,i) = star_huntpol(t, i) end do m_v(:) = star_v(IMASS, :) s_v(:) = star_v(IRADI, :) k_v(:) = star_v(IAPSI, :) i_v(:) = star_v(IINER, :) q_v(:) = star_v(IQUAL, :) mi_v(:) = 1.d0 / m_v(:) qi_v(:) = 1.d0 / q_v(:) s5k_v(:) = s_v(:)**5 * k_v(:) m2_v(:) = m_v(:)**2 do concurrent (i=1:nstar) if (i_v(i) > 0.d0) then w_v(:,i) = j_v(:,i) / i_v(i) else w_v(:,i) = 0.d0 end if wn2_v(i) = sum(w_v(:,i)**2) end do mp_e(:) = m_v(ies(1,:)) * m_v(ies(2,:)) m_e(:) = m_v(ies(1,:)) + m_v(ies(2,:)) mm2_e(:) = 1.d0 / m_e(:)**2 ! convert from Jacobi coordinates to real relative vectors if (use_qp) then mi_d(1) = 1.d0 / m_e(1) mi_d(2) = 1.d0 / m_e(6) mc_y(1,1) = m_v(2) * mi_d(1) mc_y(2,1) = -m_v(1) * mi_d(1) mc_y(1,2) = -m_v(4) * mi_d(2) mc_y(2,2) = m_v(3) * mi_d(2) ! do i=1,2 ! do j=1,2 ! k = 1 - i + 2 * j ! use map ! mc_y(i,j) = float((3 - 2 * j) * (3 - 2 * i), kind=real64) * mi_d(j) ! enddo ! enddo do i=1,2 do j=1,2 rc_y(:,i,j) = mc_y(i,j) * r_j(:,j) vc_y(:,i,j) = mc_y(i,j) * v_j(:,j) enddo enddo r_e(:,1) = r_j(:,1) r_e(:,6) = r_j(:,2) r_e(:,2) = r_j(:,3) + rc_y(:,1,1) + rc_y(:,1,2) r_e(:,3) = r_j(:,3) + rc_y(:,1,1) + rc_y(:,2,2) r_e(:,4) = r_j(:,3) + rc_y(:,2,1) + rc_y(:,1,2) r_e(:,5) = r_j(:,3) + rc_y(:,2,1) + rc_y(:,2,2) v_e(:,1) = v_j(:,1) v_e(:,6) = v_j(:,2) v_e(:,2) = v_j(:,3) + vc_y(:,1,1) + vc_y(:,1,2) v_e(:,3) = v_j(:,3) + vc_y(:,1,1) + vc_y(:,2,2) v_e(:,4) = v_j(:,3) + vc_y(:,2,1) + vc_y(:,1,2) v_e(:,5) = v_j(:,3) + vc_y(:,2,1) + vc_y(:,2,2) ! do i=1,2 ! do j=1,2 ! k = 2 * i + j - 1 ! use map ! r_e(:,k) = r_j(:,3) + rc_y(:,i,1) + rc_y(:,j,2) ! v_e(:,k) = v_j(:,3) + vc_y(:,i,1) + vc_y(:,j,2) ! enddo ! enddo else m_12 = m_e(1) mi_12 = 1.d0 / m_12 m_123 = m_12 + m_v(3) mi_123 = 1.d0 / m_123 m_1234 = m_123 + m_v(4) mc_y(1,1) = + m_v(2) * mi_12 mc_y(1,2) = - m_v(1) * mi_12 mc_y(2,1) = + m_v(3) * mi_123 mc_y(2,2) = - m_12 * mi_123 do i=1,2 rc_y(:,1,i) = mc_y(1,i) * r_j(:,1) rc_y(:,2,i) = mc_y(2,i) * r_j(:,2) + r_j(:,3) vc_y(:,1,i) = mc_y(1,i) * v_j(:,1) vc_y(:,2,i) = mc_y(2,i) * v_j(:,2) + v_j(:,3) enddo r_e(:,1) = r_j(:,1) r_e(:,2) = r_j(:,2) + rc_y(:,1,1) r_e(:,4) = r_j(:,2) + rc_y(:,1,2) r_e(:,3) = rc_y(:,2,1) + rc_y(:,1,1) r_e(:,5) = rc_y(:,2,1) + rc_y(:,1,2) r_e(:,6) = rc_y(:,2,2) v_e(:,1) = v_j(:,1) v_e(:,2) = v_j(:,2) + vc_y(:,1,1) v_e(:,4) = v_j(:,2) + vc_y(:,1,2) v_e(:,3) = vc_y(:,2,1) + vc_y(:,1,1) v_e(:,5) = vc_y(:,2,1) + vc_y(:,1,2) v_e(:,6) = vc_y(:,2,2) endif do l=1,nedge rn_e(l) = norm2(r_e(:,l)) enddo ! fail on star interaction if (has_interact) then do l=1,nedge if (rn_e(l) < interact(ies(1,l),ies(2,l))) then write(6, "(' [direct4f] object interaction ',i1,' + ',i1)") ies(:,l) yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = ies(:,l) terminate_integration = .true. return endif end do endif ! fail on star collisions do l=1,nedge if (rn_e(l) < s_v(ies(1,l)) + s_v(ies(2,l))) then write(6, "(' [direct4f] object collision ',i1,' + ',i1)") ies(:,l) yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = ies(:,l) terminate_integration = .true. return endif end do ! fail on star escape if (has_cutoff) then do j=1, norb if (norm2(r_j(:,j)) > cutoff(j)) then write(6, "(' [direct4f] orbit escape ',i1)") j yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/j, 0/) terminate_integration = .true. return endif enddo endif ! compute forces rn2_e(:) = rn_e(:)**2 rnm1_e(:) = 1.d0 / rn_e(:) rnm2_e(:) = rnm1_e(:)**2 rnm3_e(:) = rnm1_e(:) * rnm2_e(:) rnm4_e(:) = rnm2_e(:)**2 rnm5_e(:) = rnm4_e(:) * rnm1_e(:) rnm7_e(:) = rnm4_e(:) * rnm3_e(:) rnm8_e(:) = rnm4_e(:)**2 ! used for tf if (use_tf) then s5kqi_v(:) = s5k_v(:) * qi_v(:) endif ! used for qd or tf if (use_sd .or. use_tf) then do l=1,nedge do j=1,2 i = ies(j,l) wdr_c(j,l) = dot_product(w_v(:,i), r_e(:,l)) wxr_c(:,j,l) = cross_product(w_v(:,i), r_e(:,l)) end do end do end if ! direct gravity do l=1,nedge v(:) = ((- GRAV) * rnm3_e(l) * mp_e(l)) * r_e(:,l) fng_e(:,l) = v(:) f_e(:,l) = v(:) enddo if (use_pn) then do l=1, nedge eta = mp_e(l) * mm2_e(l) phi = GRAV * rnm1_e(l) * m_e(l) vn2 = sum(v_e(:,l)**2) vr = dot_product(v_e(:,l), r_e(:,l)) g = Gcm2 * mp_e(l) * rnm3_e(l) v(:) = & ((( 4.d0 + 2.d0 * eta) * phi + & (-1.d0 - 3.d0 * eta) * vn2 + & ( 1.5d0 * eta) * vr**2 * rnm2_e(l)) * g) * r_e(:,l) + & (( 4.d0 - 2.d0 * eta) * vr * g) * v_e(:,l) f_e(:,l) = f_e(:,l) + v(:) fpn_e(:,l) = v(:) ggr_e(l) = CLIGHTm2 * ( & (1.d0 - 3.d0 * eta) * 0.5d0 * vn2 & + (3.d0 + eta) * phi) fgr_e(l) = 1.d0 + ggr_e(l) enddo else ggr_e(:) = 0.d0 fgr_e(:) = 1.d0 fpn_e(:,:) = 0.d0 endif torque_v(:,:) = 0.d0 if (use_sd .or. use_td) then do l=1,nedge if (use_td) then h = (-12.d0 * GRAV) * rnm1_e(l) endif do j=1,2 i = ies(j,l) ir = ies(3-j,l) f = s5k_v(i) * m_v(ir) * rnm7_e(l) if (use_td) then ftd_c(:,j,l) = (f * h * m_v(ir)) * r_e(:,l) f_e(:,l) = f_e(:,l) + ftd_c(:,j,l) else ftd_c(:,j,l) = 0.0d0 endif if (use_sd) then g = -2.d0 * rn_e(l) * wdr_c(j,l) * f fsd_c(:,j,l) = & (f * (5.d0 * wdr_c(i,j)**2 - rn2_e(l) * wn2_v(i))) * r_e(:,l) & + g * w_v(:,i) f_e(:,l) = f_e(:,l) + fsd_c(:,j,l) torque_v(:,i) = torque_v(:,i) + (g * fgr_e(l)) * wxr_c(:,j,l) else fsd_c(:,j,l) = 0.0d0 endif end do end do else ftd_c(:,:,:) = 0.0d0 fsd_c(:,:,:) = 0.0d0 endif if (use_td3) then do i=1, nstar do j1=1, nstar1 j = iij(i,j1) ij = ipe(i,j) rij(:) = r_e(:,ij) do k1=1, nstar2 k = iijk(i,j1,k1) ! DEBUG if ((i==j).or.(j==k).or.(k==i)) STOP 'index bug' ik = ipe(i,k) rik(:) = r_e(:,ik) ! all signs of r cancel out for f_e if (j < k) then ! no need to redo symmetric (dot) products g = dot_product(rij(:), rik(:)) gjk_x(j,k) = g f = m_v(j) * m_v(k) mjk_x(j,k) = f else g = gjk_x(k,j) f = mjk_x(k,j) endif f = f * s5k_v(i) * GRAV * rnm3_e(ij) * rnm5_e(ik) h = 2.d0 * f * g * rnm2_e(ij) v(:) = & ((3.d0 + 11.d0 * g**2 * rnm2_e(ij) * rnm2_e(ik)) * f & ) * rik(:) + & h * rij(:) f_e(:, ik) = f_e(:, ik) + v(:) ftd3_d(:,i,j1,k1) = v(:) ! DEVEL - only record net torque, which is due to PN if (.not. use_pn) cycle if (j < k) then ! no need to redo antisymmetric cross product v(:) = (h * (ggr_e(ik) - ggr_e(ij))) * & cross_product(rij(:), rik(:)) torque_v(:,i) = torque_v(:,i) + v(:) endif enddo enddo enddo else ftd3_d(:,:,:,:) = 0.0d0 endif if (use_tf) then do l=1, nedge w(:) = cross_product(r_e(:,l), v_e(:,l)) * rnm2_e(l) g = 6.d0 * GRAV * rnm8_e(l) / norm2(w) h = -2.d0 * dot_product(r_e(:,l), v_e(:,l)) * rnm2_e(l) do j=1, 2 i = ies(j,l) ir = ies(3-j,l) f = g * s5kqi_v(i) * m2_v(ir) v(:) = f * (wxr_c(:,j,l) - v_e(:,l)) & + (f * h) * r_e(:,l) f_e(:,l) = f_e(:,l) + v(:) ftf_c(:,j,l) = v(:) f = f * fgr_e(l) v(:) = (f * wdr_c(j,l)) * r_e(:,l) & + (f * rn2_e(l)) * (w(:) - w_v(:,i)) torque_v(:,i) = torque_v(:,i) + v(:) end do end do else ftf_c(:,:,:) = 0.d0 endif ! Note that f_e points in opposite direction as r_e as per ML02 acc_j(:,1) = & + ( mi_v(1) + mi_v(2)) * f_e(:,1) & + mi_v(1) * (f_e(:,2) + f_e(:,3)) & + (-mi_v(2)) * (f_e(:,4) + f_e(:,5)) if (use_qp) then acc_j(:,2) = & + ( mi_v(3) + mi_v(4)) * f_e(:,6) & + (-mi_v(3)) * (f_e(:,2) + f_e(:,4)) & + mi_v(4) * (f_e(:,3) + f_e(:,5)) acc_j(:,3) = (mi_d(1) + mi_d(2)) * sum(f_e(:,2:5), dim=2) else acc_j(:,2) = (m_123 * mi_12 * mi_v(3)) * (f_e(:,2) + f_e(:,4)) & + mi_12 * (f_e(:,3) + f_e(:,5)) + (-mi_v(3)) * f_e(:,6) acc_j(:,3) = (m_1234 * mi_123 * mi_v(4)) * (f_e(:,3) + f_e(:,5) + f_e(:,6)) endif yp( 1: 9) = reshape(v_j, (/ndim*norb/)) yp(10:18) = reshape(acc_j, (/ndim*norb/)) yp(19:30) = reshape(torque_v, (/ndim*nstar/)) if (local_store) then local_data( 1: 30) = yp(1:30) local_data( 31: 66) = reshape(fsd_c, (/ndim*ncorn*nedge/)) !36 local_data( 67:102) = reshape(ftd_c, (/ndim*ncorn*nedge/)) !36 local_data(103:138) = reshape(ftf_c, (/ndim*ncorn*nedge/)) !36 local_data(139:210) = reshape(ftd3_d, (/ndim*nstar*nstar1*nstar2/)) !72 local_data(211:228) = reshape(fpn_e, (/ndim*nedge/)) !18 local_data(229:246) = reshape(fng_e, (/ndim*nedge/)) !18 endif end function direct4f end module quaternary