module ternary use parameters, only: & verb_interact implicit none contains ! ======================================================================= ! Direct code, hierarchical ! ======================================================================= function direct3h(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, 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, & epslim, omegalim 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_pnf, use_sdh, use_tdh, use_tfh, use_pnh, & use_td3, use_rxv, & use_sdt, use_tdt, use_tft, use_pnt, & use_ori, use_orq, use_pqm, & use_pqh use orientation, only: & rotation_w2A, & rotation_f2A, & rotation_dodt, & rotation_dqdt implicit none real(kind=real64), intent(in) :: & trel real(kind=real64), intent(in), dimension(:) :: & y real(kind=real64), dimension(size(y, 1)) :: & yp integer(kind=int32), parameter :: & nstar = 3, & nquat = 4, & ndim = 3, & nbin = 2, & ! a binary contains 2 stars norb = nstar - 1 integer(kind=int32), parameter :: & nvorb = ndim * norb, & nvstar = ndim * nstar, & nqstar = nquat * nstar 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 real(kind=real64), dimension(nbin), parameter :: & sign_v3 = (/1.d0, -1.d0/) integer(kind=int32) :: & i real(kind=real64) :: & t, f, g, h, b, rnm8 real(kind=real64), dimension(ndim) :: & v, w, qr real(kind=real64), dimension(stars(1)%m, nstar) :: & star_v real(kind=real64), dimension(nstar) :: & m_v,s_v,k_v,t_v, & mi_v, wn2_v, s5k_v real(kind=real64), dimension(ndim, nstar) :: & i_v, j_v, w_v, torque_v, ori_v, rot_v, q_v real(kind=real64), dimension(nquat, nstar) :: & orq_v, roq_v logical, dimension(nstar) :: & lpqm_v real(kind=real64), dimension(ndim, ndim, nstar) :: & A_v real(kind=real64), dimension(ndim) :: & r_12, v_12, re_12, acc_12, & fpn_12, f_12 real(kind=real64) :: & m_12, mu_12, mi_12, mp_12, & rn_12, rnm1_12, rnm2_12, rnm3_12, rnm4_12, rnm8_12, & v2_12, phi_12, eps_12, an_12, nm1_12, Gm_12, & vrn_12, eta_12, ftgr_12, vr2_12 real(kind=real64), dimension(nbin) :: & mrat_v12, wre_v12, mr_v12, & mf_v3, wre_v3 real(kind=real64), dimension(ndim, nbin) :: & fsd_v12, ftd_v12, ftf_v12, & wxr_v12, ftd3_12, & fpq_v12, tpq_v12 real(kind=real64), dimension(ndim, nbin) :: & fsdt_v3, ftdt_v3, ftft_v3, wxr_v3 real(kind=real64) :: & m_123, mi_123, mui_123 real(kind=real64), dimension(ndim) :: & r_123, v_123, acc_123, f3_123 real(kind=real64) :: & rn_123, rnm1_123, rnm2_123, rnm3_123, rnm4_123, wre_123, mu_123, & v2_123, Gm_123, phi_123, eps_123, an_123, nm1_123, rnm8_123, vrn_123, & vr2_123, ftgr_123, eta_123 real(kind=real64), dimension(ndim) :: & f_123, re_123, fsd_123, ftd_123, & ftf_123, & fpn_123, wxr_123, w_123, ff, & fpq_123, tpq_123 real(kind=real64), dimension(ndim, nbin) :: & r_v3, v_v3, b_v3, & re_v3, fpn_v3 real(kind=real64), dimension(nbin) :: & rn_V3, rnm1_V3, rnm2_v3, rnm3_v3, vrn_v3, & eta_v3, m_v3, mp_v3, ftgr_v3, & rnm4_v3, mu_v3, v2_v3, vr2_v3, phi_v3 real(kind=real64) :: & nm1_i, phi_i, gm_i, eps_i, an_i ! [[star1, star2], star3] ! VARIABLE NAMES ! _v - values for 3 stars as last index ! _12 - scalar values for 2-star system ! _v12 - array values for 2-star system ! _123 - scalar values for 3rd-star system Jacobi coordinate ! _v3 - array values for 3rd-star system star 1 and 2 properties ! _i - local values 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, :) i_v(:,:) = star_v(IINER, :) k_v(:) = star_v(IAPSI, :) t_v(:) = star_v(ITAUQ, :) ! use t < 0 as flag for containing -1/Q, tau for q > 0 r_12(:) = y(1:3) r_123(:) = y(4:6) v_12(:) = y(7:9) v_123(:) = y(10:12) j_v(:,:) = reshape(y(13:21), (/3,3/)) if (use_ori) then if (use_orq) then orq_v(:,:) = reshape(y(1+2*nvorb+nvstar:2*nvorb+nvstar+nqstar), & (/nquat, nstar/)) roq_v(:,:) = 0.d0 else ori_v(:,:) = reshape(y(1+2*nvorb+nvstar:2*(nvorb+nvstar)), & (/ndim, nstar/)) rot_v(:,:) = 0.d0 endif end if if (use_ori) then do concurrent (i=1:nstar) ! alternate approach if use i_x == i_y == i_z ! if (abs(i_v(1,i)-i_v(2,i)) + abs(i_v(2,i)-i_v(3,i)) > 1.d-12 * i_v(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 (i_v(1,i) > 0.d0) then if (use_orq) then a_v(:,:,i) = rotation_f2A(orq_v(:,i)) else a_v(:,:,i) = rotation_w2A(ori_v(:,i)) endif w_v(:,i) = matmul(a_v(:,:,i), matmul(j_v(:,i), A_v(:,:,i)) / i_v(:,i)) if (i == 3) then lpqm_v(i) = use_pqh else lpqm_v(i) = use_pqm endif else if (i_v(3,i) > 0.d0) then w_v(:,i) = j_v(:,i) / i_v(3,i) lpqm_v(i) = .false. else w_v(:,i) = 0.d0 lpqm_v(i) = .false. endif if (i_v(3,i) > 0.d0) then if (use_orq) then roq_v(:,i) = rotation_dqdt(orq_v(:,i), w_v(:,i)) else rot_v(:,i) = rotation_dodt(ori_v(:,i), w_v(:,i)) endif endif ! TODO - modify rot(:,:) for GR/SR time delay if (lpqm_v(i)) then q_v(:, i) = ONETHIRD * sum(i_v(:,i)) - i_v(:,i) endif end do else do concurrent (i=1:nstar) if (i_v(3,i) > 0.d0) then w_v(:,i) = j_v(:,i) / i_v(3,i) else w_v(:,i) = 0.d0 endif lpqm_v(i) = .false. end do end if wn2_v(:) = sum(w_v(:,:)**2, dim=1) mi_v(:) = 1.d0 / m_v(:) s5k_v(:) = s_v(:)**5 * k_v(:) mr_v12(:) = m_v(2:1:-1) mrat_v12(:) = mr_v12(:) * mi_v(1:2) m_12 = sum(m_v(1:2)) mi_12 = 1.d0 / m_12 mp_12 = product(m_v(1:2)) mu_12 = mp_12 * mi_12 rn_12 = norm2(r_12) rnm1_12 = 1.d0 / rn_12 rnm2_12 = rnm1_12**2 rnm3_12 = rnm2_12 * rnm1_12 rnm4_12 = rnm2_12**2 rnm8_12 = rnm4_12**2 re_12(:) = r_12(:) * rnm1_12 ! fail on interaction if (has_interact.and.(rn_12 < interact(1,2))) then if (verb_interact) then write(6, "(' [direct3h] 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 binary merger if (rn_12 < sum(s_v(1:2))) then if (verb_interact) then write(6, "(' [direct3h] 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_12 > cutoff(1))) then if (verb_interact) then write(6, "(' [direct3h] orbit escaped ',i1)") 1 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 wre_v12(i) = dot_product(w_v(:,i), re_12(:)) enddo v2_12 = sum(v_12(:)**2) Gm_12 = GRAV * m_12 phi_12 = Gm_12 * rnm1_12 vrn_12 = dot_product(v_12(:), re_12(:)) if (use_sd .or. use_tf) then do i=1,2 wxr_v12(:,i) = cross_product(w_v(:, i), r_12(:)) enddo end if acc_12(:) = (- phi_12 * rnm2_12) * r_12(:) torque_v(:,:) = 0.d0 f_12(:) = 0.d0 if (use_pn) then eta_12 = mu_12 * mi_12 vr2_12 = vrn_12**2 fpn_12(:) = & phi_12 * rnm1_12 * CLIGHTm2 * ( & v_12(:) * ( & ( 4.d0 - 2.d0 * eta_12) * vrn_12) & + re_12(:) * ( & + ( 4.d0 + 2.d0 * eta_12) * phi_12 & + 1.5d0 * eta_12 * vr2_12 & + (-1.d0 - 3.d0 * eta_12) * v2_12) & ) f_12(:) = f_12(:) + fpn_12(:) ftgr_12 = mu_12 * (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_12 * (1.d0 - 3.d0 * eta_12) & + (3.d0 + eta_12) * phi_12)) else ftgr_12 = mu_12 fpn_12(:) = 0.d0 endif if (use_ori.and.use_pqm.and.any(lpqm_v(1:2))) then g = G3 * rnm2_12 * rnm3_12 / mu_12 do i=1,2 if (.not.lpqm_v(i)) then fpq_v12(:,i) = 0.d0 tpq_v12(:,i) = 0.d0 cycle endif f = g * mr_v12(i) qr(:) = matmul(a_v(:,:,i), q_v(:,i) * matmul(r_12(:), a_v(:,:,i))) fpq_v12(:,i) = f * ( & qr(:) - 2.5d0 * rnm2_12 * sum(r_12(:) * qr(:)) * r_12(:)) tpq_v12(:,i) = f * cross_product(qr(:), r_12(:)) * ftgr_12 f_12(:) = f_12(:) + fpq_v12(:,i) torque_v(:,i) = torque_v(:,i) + tpq_v12(:,i) end do else fpq_v12(:,:) = 0.d0 tpq_v12(:,:) = 0.d0 endif if (use_sd .or. use_td) then if (use_td) then f = -12.d0 * GRAV * rnm3_12 endif g = m_12 * rnm4_12 do i=1,2 h = s5k_v(i) * mi_v(i) * g if (use_td) then ftd_v12(:,i) = (h * f * mr_v12(i)) * re_12(:) f_12(:) = f_12(:) + ftd_v12(:,i) else ftd_v12(:,i) = 0.d0 endif if (use_sd) then b = -2.d0 * h * wre_v12(i) fsd_v12(:,i) = & (h * (5.d0 * wre_v12(i)**2 - wn2_v(i))) * re_12(:) & + b * w_v(:,i) f_12(:) = f_12(:) + fsd_v12(:,i) torque_v(:,i) = torque_v(:,i) & + (b * ftgr_12) * wxr_v12(:,i) else fsd_v12(:,i) = 0.d0 endif end do else fsd_v12(:,:) = 0.d0 ftd_v12(:,:) = 0.d0 endif if (use_tf) then w(:) = cross_product(re_12(:), v_12(:)) * rnm1_12 if (any(t_v(1:2) < 0.d0)) then if (use_rxv) then nm1_12 = 1.d0 / max(norm2(w(:)), & omegalim * sqrt(Gm_12 * rnm3_12)) else eps_12 = min(0.5d0 * v2_12 - phi_12, -epslim * phi_12) an_12 = -0.5d0 * Gm_12 / eps_12 nm1_12 = sqrt(an_12**3 / Gm_12) endif endif f = 12.d0 * Gm_12 * rnm8_12 v(:) = (-2.d0 * vrn_12) * re_12(:) - v_12(:) do i=1,2 h = t_v(i) if (h == 0) then ftf_v12(:,i) = 0.d0 cycle else if (h < 0.d0) then h = -0.5d0 * nm1_12 * h end if h = h * f * s5k_v(i) * mrat_v12(i) ftf_v12(:,i) = h * (v(:) + wxr_v12(:,i)) f_12(:) = f_12(:) + ftf_v12(:,i) b = h * ftgr_12 * rn_12 torque_v(:,i) = torque_v(:,i) & + (b * wre_v12(i)) * r_12(:) & + (b * rn_12) * (w(:) - w_v(:,i)) end do else ftf_v12(:,:) = 0.d0 endif ! 3rd body m_123 = sum(m_v(1:3)) mi_123 = 1.d0 / m_123 mu_123 = m_12 * m_v(3) * mi_123 do i=1,nbin mf_v3(i) = mr_v12(i) * mi_12 * sign_v3(i) r_v3(:,i) = r_123(:) + mf_v3(i) * r_12(:) rn_v3(i) = norm2(r_v3(:,i)) rnm1_v3(i) = 1.d0 / rn_v3(i) rnm2_v3(i) = rnm1_v3(i)**2 rnm3_v3(i) = rnm2_v3(i) * rnm1_v3(i) b_v3(:,i) = r_v3(:,i) * rnm3_v3(i) enddo if (use_sdh .or. use_tdh .or. use_tfh .or. use_pnh .or. use_pnt) then rn_123 = norm2(r_123) rnm1_123 = 1.d0 / rn_123 rnm2_123 = rnm1_123**2 rnm3_123 = rnm2_123*rnm1_123 rnm4_123 = rnm2_123**2 re_123(:) = r_123(:) * rnm1_123 wre_123 = dot_product(w_v(:,3), re_123(:)) vrn_123 = dot_product(v_123(:), re_123) Gm_123 = GRAV * m_123 phi_123 = Gm_123 * rnm1_123 endif ! fail on star interaction if (has_interact) then do i=1,nbin if (rn_v3(i) < interact(i,3)) then if (verb_interact) then write(6, "(' [direct3h] object interaction ',i1,' + ',i1)") i, 3 endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = (/i, 3/) terminate_integration = .true. return endif end do endif ! fail on star collisions do i=1,nbin if (rn_v3(i) < s_v(i) + s_v(3)) then if (verb_interact) then write(6, "(' [direct3h] object collision ',i1,' + ',i1)") i, 3 endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = (/i, 3/) terminate_integration = .true. return endif enddo ! fail on star escape if (has_cutoff.and.(rn_123 > cutoff(2))) then if (verb_interact) then write(6, "(' [direct3h] orbit escaped ',i1)") 2 endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/2, 0/) terminate_integration = .true. return endif f3_123(:) = GRAV * m_v(3) * (b_v3(:,2) - b_v3(:,1)) acc_12(:) = acc_12(:) + f3_123(:) acc_123(:) = (- GRAV * m_123 * mi_12) * ( & m_v(1) * b_v3(:,1) + m_v(2) * b_v3(:,2)) f_123(:) = 0.d0 if (use_pnf .or. use_tft) then do i=1,nbin v_v3(:,i) = v_123(:) + mf_v3(i) * v_12(:) vrn_v3(i) = dot_product(v_v3(:,i), r_v3(:,i)) * rnm1_v3(i) re_v3(:,i) = r_v3(:,i) * rnm1_v3(i) enddo endif ! hierarchical PN only if (use_sdh .or. use_tdh .or. use_tfh) then wxr_123(:) = cross_product(w_v(:, 3), r_123(:)) end if if (use_pnh .or. use_tfh) then v2_123 = sum(v_123(:)**2) endif if (use_pnh .or. use_pnt) then eta_123 = mu_123 * mi_123 ! use hierachical approximation, GR if (use_pnh) then vr2_123 = vrn_123**2 fpn_123(:) = & phi_123 * rnm1_123 * CLIGHTm2 * ( & v_123(:) * ( & ( 4.d0 - 2.d0 * eta_123) * vrn_123) & + re_123(:) * ( & + ( 4.d0 + 2.d0 * eta_123) * phi_123 & + 1.5d0 * eta_123 * vr2_123 & + (-1.d0 - 3.d0 * eta_123) * v2_123) & ) acc_123(:) = acc_123(:) + fpn_123(:) else fpn_123(:) = 0.d0 endif ftgr_123 = mu_123 * (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_123 * (1.d0 - 3.d0 * eta_123) & + (3.d0 + eta_123) * phi_123)) else ftgr_123 = mu_123 fpn_123(:) = 0.d0 endif if (use_sdt .or. use_tdt .or. use_tft .or. use_pqh .or. & (use_pnt .and. .not. use_pnf)) then mui_123 = 1.d0 / mu_123 endif ! use hierachical approximation, PQ if (use_ori.and.use_pqh.and.lpqm_v(3)) then g = G3 * rnm2_123 * rnm3_123 * mui_123 f = g * mr_v12(i) qr(:) = matmul(a_v(:,:,3), q_v(:,3) * matmul(r_123(:), a_v(:,:,3))) fpq_123(:) = f * ( & qr(:) - 2.5d0 * rnm2_123 * sum(r_123(:) * qr(:)) * r_123(:)) tpq_123(:) = f * cross_product(qr(:), r_123(:)) * ftgr_123 f_123(:) = f_123(:) + fpq_123(:) torque_v(:,i) = torque_v(:,i) + tpq_123(:) else fpq_123(:) = 0.d0 tpq_123(:) = 0.d0 endif ! use hierachical approximation, QD (TD + SD) if (use_sdh .or. use_tdh) then h = s5k_v(3) * mi_v(3) * m_123 * rnm4_123 if (use_tdh) then ftd_123(:) = (-12.d0 * h * GRAV * rnm3_123 * m_12) * re_123(:) else ftd_123(:) = 0.d0 endif if (use_sdh) then b = -2.d0 * h * wre_123 fsd_123(:) = & (h * ( 5.d0 * wre_123**2 - wn2_v(3))) * re_123(:) & + b * w_v(:,3) f_123(:) = f_123(:) + fsd_123(:) torque_v(:,3) = torque_v(:,3) + (b * ftgr_123) * wxr_123(:) else fsd_123(:) = 0.d0 endif else fsd_123(:) = 0.d0 ftd_123(:) = 0.d0 endif ! use hierachical approximation, TF if (use_tfh .and. t_v(3) /= 0.d0) then w_123(:) = cross_product(v_123(:), re_123(:)) * rnm1_123 rnm8_123 = rnm4_123**2 g = -2.d0 * vrn_123 v(:) = re_123(:) * g - v_123(:) f = 12.d0 * Gm_123 * rnm8_123 h = t_v(3) if (h < 0.d0) then if (use_rxv) then nm1_123 = 1.d0 / max(norm2(w_123(:)), & omegalim * sqrt(Gm_123 * rnm3_123)) else eps_123 = min(0.5d0 * v2_123 - phi_123, -epslim * phi_123) an_123 = -0.5d0 * Gm_123 / eps_123 nm1_123 = sqrt(an_123**3 / Gm_123) endif h = -0.5d0 * h * nm1_123 end if h = h * f * s5k_v(3) * m_12 * mi_v(3) ftf_123(:) = h * (v(:) + wxr_123(:)) f_123(:) = f_123(:) + ftf_123(:) b = h * ftgr_123 * rn_123 torque_v(:,3) = torque_v(:,3) + & + (b * wre_123) * r_123(:) & + (b * rn_123) * (w_123(:) - w_v(:,3)) else ftf_123(:) = 0.d0 endif ! impact of 3rd star on tides in inner binary if (use_td3) then b = m_12 * GRAV * rnm4_12 * rnm1_12 do i=1, nbin g = dot_product(r_12(:), r_v3(:,i)) f = b * s5k_v(i) * rnm3_v3(i) h = 6.d0 * f * g * rnm2_v3(i) ftd3_12(:,i) = & ((3.d0 - 15.d0 * g**2 * rnm2_v3(i) * rnm2_12) * f & ) * r_12(:) + & h * r_v3(:,i) f_12(:) = f_12(:) + ftd3_12(:,i) torque_v(:,i) = torque_v(:,i) & + (h * ftgr_12) * cross_product(r_v3(:,i), r_12(:)) enddo else ftd3_12(:,:) = 0.d0 endif ! full PN for gravity if (use_pnf) then do i=1,nbin mu_v3(i) = m_v(i) * m_v(3) / (m_v(i) + m_v(3)) m_v3(i) = m_v(i) + m_v(3) mp_v3(i) = m_v(i) * m_v(3) eta_v3(i) = mp_v3(i) / m_v3(i)**2 v2_v3(i) = sum(v_v3(:,i)**2) vr2_v3(i) = vrn_v3(i)**2 phi_v3(i) = GRAV * m_v3(i) * rnm1_v3(i) fpn_v3(:,i) = & phi_v3(i) * rnm1_v3(i) * CLIGHTm2 * ( & v_v3(:,i) * ( & ( 4.d0 - 2.d0 * eta_v3(i)) * vrn_v3(i)) & + re_v3(:,i) * ( & + ( 4.d0 + 2.d0 * eta_v3(i)) * phi_v3(i) & + 1.5d0 * eta_v3(i) * vr2_v3(i) & + (-1.d0 - 3.d0 * eta_v3(i)) * v2_v3(i)) & ) acc_123(:) = acc_123(:) + fpn_v3(:,i) * mu_v3(i) * & m_123 * mi_12 * mi_v(3) acc_12(:) = acc_12(:) + & fpn_v3(:,i) * mu_v3(i) * mi_v(i) * sign_v3(i) ! only use for direct tides tft, tdt, sdt if (use_tft .or. use_sdt) then ftgr_v3(i) = (1.d0 + CLIGHTm2 * ( & 0.5d0 * v2_v3(i) * (1.d0 - 3.d0 * eta_v3(i)) & + (3.d0 + eta_v3(i)) * phi_v3(i))) endif enddo else fpn_v3(:,:) = 0.d0 ftgr_v3(:) = 1.d0 endif ! approximate average PN corrections for 3rd star and binary CM if (use_pnt .and. .not. use_pnf) then ftgr_v3(:) = ftgr_123 * mui_123 endif ! impact of binary tides due to 3rd star on 3rd star orbit (SD, TD, TF, PN) if (use_sdt .or. use_tft) then do i=1,nbin wre_v3(i) = dot_product(w_v(:,3), r_v3(:,i)) * rnm1_v3(i) wxr_v3(:,i) = cross_product(w_v(:, i), r_v3(:,i)) enddo endif if (use_sdt .or. use_tft) then rnm4_v3(:) = rnm2_v3(:)**2 endif if (use_sdt .or. use_tdt) then if (use_tdt) then f = -12.d0 * GRAV * m_v(3) endif do i=1,2 h = s5k_v(i) * rnm4_v3(i) * m_v(3) if (use_tdt) then ff(:) = (h * f * rnm4_v3(i)) * r_v3(:,i) ftdt_v3(:,i) = mi_v(3) * ff(:) f_12(:) = f_12(:) + (mi_v(i) * sign_v3(i)) * ff(:) f_123(:) = f_123(:) + mui_123 * ff(:) else ftdt_v3(:,i) = 0.d0 endif if (use_sdt) then b = -2.d0 * h * wre_v3(i) ff(:) = & (h * rnm1_v3(i) * & (5.d0 * wre_v3(i)**2 - wn2_v(i))) * r_v3(:,i) & + b * w_v(:,i) fsdt_v3(:,i) = mi_v(3) * ff(:) f_12(:) = f_12(:) + (mi_v(i) * sign_v3(i)) * ff(:) f_123(:) = f_123(:) + mui_123 * ff(:) torque_v(:,i) = torque_v(:,i) & + (b * ftgr_v3(i)) * wxr_v3(:,i) else fsdt_v3(:,i) = 0.d0 endif end do else fsdt_v3(:,:) = 0.d0 ftdt_v3(:,:) = 0.d0 endif if (use_tft) then f = 12.d0 * GRAV * m_v(3)**2 do i=1, nbin w(:) = cross_product(v_v3(:,i), r_v3(:,i)) * rnm2_v3(i) rnm8 = rnm4_v3(i)**2 g = -2.d0 * vrn_v3(i) v(:) = re_v3(:,i) * g - v_v3(:,i) h = t_v(i) if (h == 0.d0) then ftft_v3(:,i) = 0.d0 cycle else if (h < 0.d0) then m_v3(i) = m_v(i) + m_v(3) ! also used in pnf if (use_rxv) then nm1_i = 1.d0 / max(norm2(w(:)), & omegalim * sqrt(GRAV * m_v3(i) * rnm3_v3(i))) else Gm_i = GRAV * m_v3(i) phi_i = -Gm_i * rnm1_v3(i) eps_i = min(0.5d0 * sum(v_v3(:,i)**2) + phi_i, epslim * phi_i) an_i = -0.5d0 * Gm_i / eps_i nm1_i = sqrt(an_i**3 / Gm_i) endif h = -0.5d0 * h * nm1_i endif h = h * f * s5k_v(i) * rnm8 ff(:) = h * (v(:) + wxr_v3(:,i)) ftft_v3(:,i)= mi_v(3) * ff(:) f_12(:) = f_12(:) + (mi_v(i) * sign_v3(i)) * ff(:) f_123(:) = f_123(:) + mui_123 * ff(:) b = h * ftgr_v3(i) * rn_v3(i) torque_v(:,i) = torque_v(:,i) & + (b * wre_v3(i)) * r_v3(:,i) & + (b * rn_v3(i)) * (w(:) - w_v(:,i)) end do else ftft_v3(:,:) = 0.d0 endif ! TODO - impact of in inner binary tides on 3rd star ! The effect may be small compared to quadrupole moment of ! binary orbit itself. ! TODO - impact of in inner binary permanent quadrupoles on 3rd star ! The effect may be small compared to quadrupole moment of ! binary orbit itself. acc_12(:) = acc_12(:) + f_12(:) acc_123(:) = acc_123(:) + f_123(:) yp( 1: 3) = v_12(1:3) yp( 4: 6) = v_123(1:3) yp( 7: 9) = acc_12(1:3) yp(10:12) = acc_123(1:3) yp(13:21) = reshape(torque_v, (/9/)) if (use_ori) then if (use_orq) then yp(22:33) = reshape(roq_v, (/nqstar/)) else yp(22:30) = reshape(rot_v, (/nvstar/)) endif endif if (local_store) then ! actual derivative data local_data( 1:size(yp,1)) = yp(:) ! basic forces and torques local_data( 34: 39) = reshape(fsd_v12, (/6/)) local_data( 40: 42) = fsd_123(1:3) local_data( 43: 48) = reshape(ftd_v12, (/6/)) local_data( 49: 51) = ftd_123(1:3) local_data( 52: 57) = reshape(ftf_v12, (/6/)) local_data( 58: 60) = ftf_123(1:3) local_data( 61: 63) = fpn_12(1:3) local_data( 64: 66) = fpn_123(1:3) local_data( 67: 69) = f3_123(1:3) local_data( 70: 75) = reshape(fpn_v3(1:3,1:2), (/6/)) local_data( 76: 81) = reshape(ftd3_12(1:3,1:2), (/6/)) local_data( 82: 87) = reshape(fsdt_v3(1:3,1:2), (/6/)) local_data( 88: 93) = reshape(ftdt_v3(1:3,1:2), (/6/)) local_data( 94: 99) = reshape(ftft_v3(1:3,1:2), (/6/)) ! permanent quadrupole momenent local_data(100:105) = reshape(fpq_v12, (/6/)) local_data(106:111) = reshape(tpq_v12, (/6/)) local_data(112:114) = fpq_123(1:3) local_data(115:117) = tpq_123(1:3) endif end function direct3h ! ======================================================================= ! Direct code, full ! ======================================================================= function direct3f(trel, y) result(yp) use, intrinsic :: & ieee_arithmetic, only: & ieee_value, ieee_quiet_nan, ieee_is_nan use typedef, only: & int32, real64 use physconst, only: & GRAV, CLIGHT use stardata, only: & IMASS, IRADI, IAPSI, IINER, ITAUQ, & toffset, stars, & star_huntpol use vector, only: & cross_product, & cross_product_tensor_contract use deriv_data, only: & local_store, local_data use integrator_flags, only: & terminate_integration use flags_data, only: & use_sd, use_td, use_tf, use_pn, & use_td3, & use_rxv, & use_ori, use_orq, use_ors, & use_pqm, use_pqq, use_cd use parameters, only: & cutoff, has_cutoff, & interact, has_interact, & epslim, omegalim use states, only: & status, status_stars, & STATUS_OK, STATUS_COLLIDE, STATUS_ESCAPE, STATUS_INTERACT use orientation, only: & rotation_w2A, & rotation_f2A, & rotation_dodt, & rotation_dqdt, & rotate_diag use collision, only: & collision_inspector ! , ifail implicit none integer(kind=int32), parameter :: & ndim = 3, & nquat = 4, & ncorn = 2, & nstar = 3, & norb = nstar - 1, & ncon = nstar - 1 integer(kind=int32), parameter :: & nstar1 = nstar - 1, & nstar2 = nstar - 2, & nedge = (nstar1 * nstar) / 2, & nvorb = ndim * norb, & nvstar = ndim * nstar, & nqstar = nquat * nstar ! stars belonging to each edge (ncorn=2) integer(kind=int32), dimension(ncorn, nedge), parameter :: & ies = reshape((/1,2, 1,3, 2,3/), (/ncorn, nedge/)) ! edge 1 2 3 ! edges belonging to each star (ncon=nstar-1) integer(kind=int32), dimension(ncon, nstar), parameter :: & ise = reshape((/1,2, 1,3, 2,3/), (/ncon, nstar/)) ! star 1 2 3 ! corners belonging to each star (ncon=nstar-1) integer(kind=int32), dimension(ncon, nstar), parameter :: & isc = reshape((/1,1, 2,1, 2,2/), (/ncon, nstar/)) ! star 1 2 3 ! edges belonging to each star pair combination (symmetric) integer(kind=int32), dimension(nstar, nstar), parameter :: & ipe = reshape((/0,1,2, & 1,0,3, & 2,3,0/), (/nstar, nstar/)) ! signs for corners real(kind=real64), dimension(ncorn), parameter :: & sc = (/1.d0, -1.d0/) ! reverse indices for corners integer(kind=int32), dimension(ncorn), parameter :: & irc = (/2, 1/) ! 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, 0.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, 3,3,2/), (/nstar, nstar1/)) integer(kind=int32), dimension(nstar, nstar1, nstar2), parameter :: & iijk = reshape((/3,3,2, 2,1,1/), (/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, & ONETHIRD = 1.d0 / 3.d0, & G3 = GRAV * 3.d0 ! local variables integer(kind=int32) :: & i, j, k, l, ir, ij, ik, j1, k1, ii real(kind=real64) :: & t real(kind=real64), dimension(stars(1)%m, nstar) :: & star_v real(kind=real64), dimension(norb) :: & rn_j real(kind=real64), dimension(ndim, norb) :: & r_j, v_j, acc_j real(kind=real64), dimension(ndim, nstar) :: & i_v, j_v, w_v, torque_v, ori_v, rot_v, q_v real(kind=real64), dimension(nquat, nstar) :: & orq_v, roq_v real(kind=real64), dimension(nstar) :: & m_v, s_v, k_v, t_v, mi_v, s5k_v, & wn2_v, m2_v logical, dimension(nstar) :: & lpqm_v real(kind=real64), dimension(ndim, ndim, nstar) :: & A_v, qq_v logical, dimension(nedge) :: & collide real(kind=real64), dimension(ncorn, nedge) :: & qrr_c real(kind=real64), dimension(ndim, ncorn, nedge) :: & qr_c ! coordinate setup real(kind=real64), dimension(ndim, nedge) :: & r_e, v_e real(kind=real64) :: & m_12, mi_12, m_123 real(kind=real64), dimension(ncorn) :: & mf_x ! 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, & Gm, eps, an, nm1 real(kind=real64), dimension(ndim) :: & v, w, rik, rij, ddd, tqqq real(kind=real64), dimension(ncorn, nedge) :: & wdr_c real(kind=real64), dimension(ndim, ncorn, nedge) :: & wxr_c ! torques real(kind=real64), dimension(ndim, ncorn) :: & qqqr_x ! for recording real(kind=real64), dimension(ndim, ncorn, nedge) :: & fsd_c, ftd_c, ftf_c, fpq_c, tpq_c, tqq_c real(kind=real64), dimension(ndim, nstar, nstar1, nstar2) :: & ftd3_d real(kind=real64), dimension(ndim, nedge) :: & fpn_e, fng_e, fqq_e ! [[star1, star2], star3] ! FUTURE - note on notations for hierarchies ! [[,],] +2,-1,-1 ! not supported (but redundant) ! [,[,]] +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 ! _c - (corner, edge) values r_j(:,:) = reshape(y(1 : nvorb), (/ndim, norb /)) if (ieee_is_nan(r_j(1,1))) then terminate_integration = .true. return endif v_j(:,:) = reshape(y(1+ nvorb:2*nvorb), (/ndim, norb /)) j_v(:,:) = reshape(y(1+2*nvorb:2*nvorb+nvstar), (/ndim, nstar/)) if (use_ori) then if (use_orq) then orq_v(:,:) = reshape(y(1+2*nvorb+nvstar:2*nvorb+nvstar+nqstar), & (/nquat, nstar/)) roq_v(:,:) = 0.d0 else ori_v(:,:) = reshape(y(1+2*nvorb+nvstar:2*(nvorb+nvstar)), & (/ndim, nstar/)) rot_v(:,:) = 0.d0 endif end if 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, :) i_v(:,:) = star_v(IINER, :) k_v(:) = star_v(IAPSI, :) t_v(:) = star_v(ITAUQ, :) ! use t < 0 as flag for containing -1/Q, tau for q > 0 mi_v(:) = 1.d0 / m_v(:) s5k_v(:) = s_v(:)**5 * k_v(:) m2_v(:) = m_v(:)**2 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 m_12 = m_e(1) mi_12 = 1.d0 / m_12 m_123 = m_12 + m_v(3) mf_x(1) = + m_v(2) * mi_12 mf_x(2) = - m_v(1) * mi_12 r_e(:,1) = r_j(:,1) r_e(:,2) = r_j(:,2) + mf_x(1) * r_j(:,1) r_e(:,3) = r_j(:,2) + mf_x(2) * r_j(:,1) v_e(:,1) = v_j(:,1) v_e(:,2) = v_j(:,2) + mf_x(1) * v_j(:,1) v_e(:,3) = v_j(:,2) + mf_x(2) * v_j(:,1) forall (l=1:nedge) rn_e(l) = norm2(r_e(:,l)) ! fail on star escape if (has_cutoff) then rn_j(1) = rn_e(1) rn_j(2) = norm2(r_j(:,2)) do j=1,norb if (rn_j(j) > cutoff(j)) then if (verb_interact) then write(6, "(' [direct3f] orbit escaped ',i1)") j endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_ESCAPE status_stars = (/j, 0/) terminate_integration = .true. return endif enddo endif ! fail on star interaction if (has_interact) then do l=1,nedge if (rn_e(l) < interact(ies(1,l),ies(2,l))) then if (verb_interact) then write(6, "(' [direct3f] object interaction ',i1,' + ',i1)") ies(:,l) endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_INTERACT status_stars = ies(:,l) terminate_integration = .true. return endif end do endif ! compute orientations and angular velocities if (use_ori) then do concurrent (i=1:nstar) ! do i=1, nstar ! alternate approach if use i_x == i_y == i_z ! if (abs(i_v(1,i)-i_v(2,i)) + abs(i_v(2,i)-i_v(3,i)) > 1.d-12 * i_v(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 (i_v(3,i) == 0.d0) then w_v(:,i) = 0.d0 lpqm_v(i) = .false. cycle endif if (i_v(1,i) > 0.d0) then if (use_orq) then a_v(:,:,i) = rotation_f2A(orq_v(:,i)) else a_v(:,:,i) = rotation_w2A(ori_v(:,i)) endif w_v(:,i) = matmul(a_v(:,:,i), matmul(j_v(:,i), a_v(:,:,i)) / i_v(:,i)) lpqm_v(i) = .true. else w_v(:,i) = j_v(:,i) / i_v(3,i) lpqm_v(i) = .false. endif if (.not.(lpqm_v(i) .or. use_ors)) & cycle if (use_orq) then roq_v(:,i) = rotation_dqdt(orq_v(:,i), w_v(:,i)) else rot_v(:,i) = rotation_dodt(ori_v(:,i), w_v(:,i)) endif ! TODO - modify rot(:,:) for GR/SR time delay if (.not.(lpqm_v(i) .and. (use_pqm .or. use_pqq))) & cycle q_v(:, i) = ONETHIRD * sum(i_v(:,i)) - i_v(:,i) do k=1, ncon l = ise(k,i) j = isc(k,i) qr_c(:,j,l) = matmul(a_v(:,:,i), q_v(:,i) * matmul(r_e(:,l), a_v(:,:,i))) qrr_c(j,l) = sum(r_e(:,l) * qr_c(:,j,l)) enddo end do else do concurrent (i=1:nstar) if (i_v(3,i) > 0.d0) then w_v(:,i) = j_v(:,i) / i_v(3,i) else w_v(:,i) = 0.d0 endif end do lpqm_v(:) = .false. end if wn2_v(:) = sum(w_v(:,:)**2, dim=1) ! fail on star collisions if (use_ori.and.use_cd) then collide(:) = collision_inspector(i_v, m_v, s_v, a_v, r_e, nstar, nedge) else collide(:) = rn_e(:) < s_v(ies(1,:)) + s_v(ies(2,:)) end if do l=1,nedge if (collide(l)) then if (verb_interact) then write(6, "(' [direct3f] object collision ',i1,' + ',i1)") ies(:,l) endif yp(:) = ieee_value(0.d0,ieee_quiet_nan) status = STATUS_COLLIDE status_stars = ies(:,l) terminate_integration = .true. return endif end do ! compute forces rn2_e(:) = rn_e(:)**2 rnm1_e(:) = 1.d0 / rn_e(:) rnm2_e(:) = rnm1_e(:)**2 rnm3_e(:) = rnm1_e(:) * rnm2_e(:) if (use_tf .or. use_td3 .or. use_sd .or. use_td) then rnm4_e(:) = rnm2_e(:)**2 rnm5_e(:) = rnm4_e(:) * rnm1_e(:) rnm7_e(:) = rnm4_e(:) * rnm3_e(:) rnm8_e(:) = rnm4_e(:)**2 endif ! used for qd or tf if (use_sd .or. use_tf) then do l=1,nedge do j=1,ncorn 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 ! permanent quadrupole-quadrupole interaction if (use_ori .and. use_pqq .and. any(lpqm_v)) then do i=1, nstar if (.not.lpqm_v(i)) & cycle qq_v(:,:,i) = rotate_diag(q_v(:,i),A_v(:,:,i)) enddo do l=1,nedge i = ies(1,l) ir = ies(2,l) if (.not.all(lpqm_v((/i,ir/)))) then fqq_e(:,l) = 0.d0 tqq_c(:,:,l) = 0.d0 cycle endif qqqr_x(:,1) = matmul(qq_v(:,:,i ), qr_c(:,2,l)) qqqr_x(:,2) = matmul(qq_v(:,:,ir), qr_c(:,1,l)) g = 7.5d0 * GRAV * rnm3_e(l) * rnm4_e(l) h = 7.d0 * rnm2_e(l) f = -31.5d0 * rnm4_e(l) * qrr_c(1,l) * qrr_c(2,l) & - sum(qq_v(:,:,i)*qq_v(:,:,ir)) & + 14.d0 * rnm2_e(l) * dot_product(qr_c(:,1,l), qr_c(:,2,l)) ddd(:) = -2.d0 * (qqqr_x(:,1) + qqqr_x(:,2)) & + (h * qrr_c(2,l)) * qr_c(:,1,l) & + (h * qrr_c(1,1)) * qr_c(:,2,l) fqq_e(:,l) = g * (ddd(:) + f * r_e(:,l)) f_e(:,l) = f_e(:,l) + fqq_e(:,l) tqqq(:) = 2.d0 * g * cross_product(qr_c(:,2,l), qr_c(:,1,l)) & + 3d0 * GRAV * rnm2_e(l) * rnm3_e(l) & * cross_product_tensor_contract(qq_v(:,:,i), qq_v(:,:,ir)) do j=1,ncorn ii = ies(j,l) tqq_c(:,j,l) = sc(j) * tqqq(:) & - 2.d0 * g * cross_product(qqqr_x(:,j),r_e(:,l)) & + g * h * qrr_c(irc(j),l) * cross_product(qr_c(:,j,l), r_e(:,l)) torque_v(:,ii) = torque_v(:,ii) + tqq_c(:,j,l) enddo enddo else fqq_e(:,:) = 0.d0 tqq_c(:,:,:) = 0.d0 endif ! permanent quadrupole-monopole interaction if (use_ori .and. use_pqm .and. any(lpqm_v)) then do l=1,nedge g = G3 * rnm2_e(l) * rnm3_e(l) do j=1,ncorn i = ies(j,l) if (.not.lpqm_v(i)) then fpq_c(:,j,l) = 0.d0 tpq_c(:,j,l) = 0.d0 cycle endif ir = ies(3-j,l) f = g * m_v(ir) fpq_c(:,j,l) = f * ( & qr_c(:,j,l) - 2.5d0 * rnm2_e(l) * qrr_c(j,l) * r_e(:,l)) tpq_c(:,j,l) = f * cross_product(qr_c(:,j,l), r_e(:,l)) * fgr_e(l) f_e(:,l) = f_e(:,l) + fpq_c(:,j,l) torque_v(:,i) = torque_v(:,i) + tpq_c(:,j,l) end do end do else fpq_c(:,:,:) = 0.d0 tpq_c(:,:,:) = 0.d0 endif 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,ncorn 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.d0 endif if (use_sd) then g = -2.d0 * rn2_e(l) * wdr_c(j,l) * f fsd_c(:,j,l) = & (f * (5.d0 * wdr_c(j,l)**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.d0 endif end do end do else ftd_c(:,:,:) = 0.d0 fsd_c(:,:,:) = 0.d0 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 = 6.d0 * f * g * rnm2_e(ij) v(:) = & ((3.d0 - 15.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 individual torques, cancel out overall due to symmetry. 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.d0 endif if (use_tf .and. any(t_v /= 0.d0)) then do l=1, nedge w(:) = cross_product(r_e(:,l), v_e(:,l)) * rnm2_e(l) g = 12.d0 * GRAV * rnm8_e(l) h = -2.d0 * dot_product(r_e(:,l), v_e(:,l)) * rnm2_e(l) if (any(t_v(ies(1:ncorn,l)) < 0.d0)) then if (use_rxv) then nm1 = 1.d0 / max(norm2(w), & omegalim * sqrt(GRAV * m_e(l) * rnm3_e(l))) else Gm = GRAV * m_e(l) phi = -Gm * rnm1_e(l) eps = min(0.5d0 * sum(v_e(:,l)**2) + phi, epslim * phi) an = -0.5d0 * Gm / eps nm1= sqrt(an**3 / Gm) endif endif do j=1,ncorn i = ies(j,l) f = t_v(i) if (f == 0.d0) then ftf_c(:,j,l) = 0.d0 continue else if (f < 0.d0) then f = -0.5d0 * f * nm1 endif ir = ies(3-j,l) f = f * g * s5k_v(i) * m2_v(ir) v(:) = f * (wxr_c(:,j,l) - v_e(:,l) + h * r_e(:,l)) ftf_c(:,j,l) = v(:) f_e(:,l) = f_e(:,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 point in opposite direction as \vec{r} as per ML02 acc_j(:,1) = & + ( mi_v(1) + mi_v(2)) * f_e(:,1) & + ( mi_v(1) ) * f_e(:,2) & + ( - mi_v(2)) * f_e(:,3) acc_j(:,2) = (m_123 * mi_12 * mi_v(3)) * sum(f_e(:,2:3), dim=2) yp( 1: 6) = reshape(v_j, (/ ndim*norb/)) yp( 7:12) = reshape(acc_j, (/ ndim*norb/)) yp(13:21) = reshape(torque_v, (/ ndim*nstar/)) if (use_ori) then if (use_orq) then yp(22:33) = reshape(roq_v, (/nqstar/)) else yp(22:30) = reshape(rot_v, (/nvstar/)) endif endif if (local_store) then ! actual derivaties local_data( 1: size(yp,1)) = yp(:) ! basic forces and torques local_data( 34: 51) = reshape(fsd_c, (/ndim*ncorn*nedge/)) !18 local_data( 52: 69) = reshape(ftd_c, (/ndim*ncorn*nedge/)) !18 local_data( 70: 87) = reshape(ftf_c, (/ndim*ncorn*nedge/)) !18 local_data( 88:105) = reshape(ftd3_d, (/ndim*ncorn*nedge/)) !18 local_data(106:114) = reshape(fpn_e, (/ndim*nedge/)) !9 local_data(115:123) = reshape(fng_e, (/ndim*nedge/)) !9 ! permanent quadrupole momenent local_data(124:141) = reshape(fpq_c, (/ndim*ncorn*nedge/)) !18 local_data(142:159) = reshape(tpq_c, (/ndim*ncorn*nedge/)) !18 local_data(160:168) = reshape(fqq_e, (/ndim*nedge/)) !9 local_data(169:186) = reshape(tqq_c, (/ndim*ncorn*nedge/)) !18 endif end function direct3f end module ternary