module quaternary use parameters, only: & verb_interact implicit none 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, 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_pnq, use_rxv, use_td3, & use_pnf, use_pnh, use_tfh, use_sdh, use_tdh, & use_pnt, use_sdt, use_tdt, use_tft, use_t3t, use_t3h, use_t3s, & 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 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 = 4, & nquat = 4, & norb = 3, & ndim = 3, & nbin = 2, & ! a binary contains 2 stars ndoub = 2 ! number of binaries integer(kind=int32), parameter :: & nvstar = ndim*nstar, & nqstar = nquat*nstar 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, qr real(kind=real64), dimension(ndim, norb) :: & r_r, v_r real(kind=real64), dimension(ndim, nbin, ndoub) :: & i_p, j_p, w_p, & mc_p, fsd_p, ftd_p, ftf_p, ftd3_p, torque_p, wxr_p, & mv_p, & ori_p, rot_p, q_p, & fpq_p, tpq_p real(kind=real64), dimension(nquat, nbin, ndoub) :: & orq_p, roq_p real(kind=real64), dimension(ndim, ndim, nbin, ndoub) :: & A_p logical, dimension(nbin,ndoub) :: & lpqm_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, t_p, & mi_p, wn2_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, phi_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) :: & an_x, eps_x 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) [which pair, which star] ! _c - pair of 2-star coordinate (..., 2, 2) [star in bin 1, star in bin 2] ! _d - pair of bianry values (..., 2) ! _q - outer double values (..., 1) ! _f - full combinations (...., 2, 2, 2) ! _x - local sacala variable ! if (.not.all(ieee_is_finite(y))) then ! yp(:) = ieee_value(0.d0,ieee_quiet_nan) ! return ! endif 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/)) i_p(:,:,:) = reshape(star_v(IINER, :), (/3, 2, 2/)) k_p(:,:) = reshape(star_v(IAPSI, :), (/2, 2/)) t_p(:,:) = reshape(star_v(ITAUQ, :), (/2, 2/)) ! use t < 0 as flag for containing -1/Q, tau for q > 0 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/)) if (use_ori) then if (use_orq) then orq_p(:,:,:) = reshape(y(31:46), (/4, 2, 2/)) roq_p(:,:,:) = 0.d0 else ori_p(:,:,:) = reshape(y(31:42), (/3, 2, 2/)) rot_p(:,:,:) = 0.d0 endif end if if (use_ori) then do concurrent (i=1:nbin) do concurrent (j=1:ndoub) ! alternate approach if use i_x == i_y == i_z ! if (abs(i_p(1,i,j)-i_p(2,i,j)) + abs(i_p(2,i,j)-i_p(3,i,j)) > 1.d-12 * i_p(3,i,j)) 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_p(1,i,j) > 0.d0) then if (use_orq) then a_p(:,:,i,j) = rotation_f2A(orq_p(:,i,j)) else a_p(:,:,i,j) = rotation_w2A(ori_p(:,i,j)) endif w_p(:,i,j) = matmul(a_p(:,:,i,j), matmul(j_p(:,i,j), A_p(:,:,i,j)) / i_p(:,i,j)) lpqm_p(i,j) = use_pqm .or. use_pqh else if (i_p(3,i,j) > 0.d0) then w_p(:,i,j) = j_p(:,i,j) / i_p(3,i,j) lpqm_p(i,j) = .false. else w_p(:,i,j) = 0.d0 lpqm_p(i,j) = .false. endif if (i_p(3,i,j) > 0.d0) then if (use_orq) then roq_p(:,i,j) = rotation_dqdt(orq_p(:,i,j), w_p(:,i,j)) else rot_p(:,i,j) = rotation_dodt(ori_p(:,i,j), w_p(:,i,j)) endif endif ! TODO - modify rot(:,:,:) for GR/SR time delay if (lpqm_p(i,j)) then q_p(:, i,j) = ONETHIRD * sum(i_p(:,i,j)) - i_p(:,i,j) endif end do enddo else do concurrent (i=1:nbin) do concurrent (j=1:ndoub) if (i_p(3,i,j) > 0.d0) then w_p(:,i,j) = j_p(:,i,j) / i_p(3,i,j) else w_p(:,i,j) = 0.d0 endif lpqm_p(i,j) = .false. end do end do end if wn2_p(:,:) = sum(w_p(:,:,:)**2, dim=1) mi_p(:,:) = 1.d0 / m_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 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 if (verb_interact) then write(6, "(' [direct4p] object interaction ',i1,' + ',i1)") idstar(j,:) endif 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 if (verb_interact) then write(6, "(' [direct4p] object collision ',i1,' + ',i1)") idstar(j,:) endif 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 if (verb_interact) then write(6, "(' [direct4p] orbit escape ',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 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 if (verb_interact) then write(6, "(' [direct4p] object interaction ',i1,' + ',i1)") i, j+2 endif 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 if (verb_interact) then write(6, "(' [direct4p] object collision ',i1,' + ',i1)") i, j+2 endif 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 in each binary 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 ! first order newtonian corrections for for force between binaries 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 PQ, SD, TD, TF if (use_ori .and. use_pqm .and. any(lpqm_p)) then do j=1,ndoub if (.not. any(lpqm_p(:,j))) then fpq_p(:,:,j) = 0.d0 tpq_p(:,:,j) = 0.d0 cycle endif g = G3 * rnm2_d(j) * rnm3_d(j) / mu_d(j) do i=1,nbin if (.not.lpqm_p(i,j)) then fpq_p(:,i,j) = 0.d0 tpq_p(:,i,j) = 0.d0 cycle endif f = g * mr_p(i,j) qr(:) = matmul(a_p(:,:,i,j),q_p(:,i,j)*matmul(r_r(:,j),a_p(:,:,i,j))) fpq_p(:,i,j) = f * ( & qr(:) - 2.5d0 * rnm2_d(j) * sum(r_r(:,j) * qr(:)) * r_r(:,j)) tpq_p(:,i,j) = f * cross_product(qr(:), r_r(:,j)) * ftgr_d(j) f_d(:,j) = f_d(:,j) + fpq_p(:,i,j) torque_p(:,i,j) = torque_p(:,i,j) + tpq_p(:,i,j) end do end do else fpq_p(:,:,:) = 0.d0 tpq_p(:,:,:) = 0.d0 endif 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 .and. any(t_p /= 0.d0)) then do j=1,ndoub w_d(:,j) = cross_product(re_d(:,j), v_r(:,j)) * rnm1_d(j) if (any(t_p(:,j) < 0.d0)) then if (use_rxv) then nm1_d(j) = 1.d0 / max(norm2(w_d(:,j)), & omegalim * sqrt(Gm_d(j) * rnm3_d(j))) else eps_x = min(0.5d0 * v2_d(j) - phi_d(j), -epslim * phi_d(j)) an_x = -0.5d0 * Gm_d(j) / eps_x nm1_d(j) = sqrt(an_x**3 / Gm_d(j)) end if endif enddo do j=1,ndoub f = 12.d0 * Gm_d(j) * rnm8_d(j) v(:) = (-2.d0 * vrn_d(j)) * re_d(:,j) - v_r(:,j) do i=1,nbin h = t_p(i,j) if (h == 0.d0) then ftf_p(:,i,j) = 0.d0 continue else if (h < 0.d0) then h = -0.5d0 * nm1_d(j) endif h = h * f * 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 tide in a binary (3-star tide) ! - does only require torque GR correction within binary if (use_td3) then do j=1,ndoub b = m_d(j) * GRAV * rnm4_d(j) * rnm1_d(j) do i=1, nbin ! check if sign is really as desired. if (j == 1) then rij(:) = r_r(:,3) + mc_p(:,i,1) else rij(:) = -r_r(:,3) - mc_p(:,i,2) endif rnm2_ij = 1.d0 / sum(rij**2) rnm1_ij = sqrt(rnm1_ij) rnm3_ij = rnm2_ij * rnm1_ij g = dot_product(r_r(:,j), rij(:)) f = b * s5k_p(i,j) * rnm3_ij h = 6.d0 * f * g * rnm2_ij ftd3_p(:,i,j) = & ((3.d0 - 15.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] use_pnf requires 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 ! TODO - individual star QP, SD, TD, TF with CM of other binary ! (the hierarchical cases this routine is about) ! PNH is only used for torque corrections if (use_pnh .and. (use_sdh .or. use_tdh .or. use_tfh .or. use_pqh)) then ! TODO - compute torque factors stop ' [direct4p] not implemented: use_pnh' else if (use_sdh .or. use_tdh .or. use_tfh .or. use_pqh) then ! TODO - set unit ftgrh factors endif if (use_ori .and. use_pqh .and. any(lpqm_p)) then stop ' [direct4p] not implemented: use_sdh, tdh' endif if (use_tfh .and. any(t_p /= 0.d0)) then stop ' [direct4p] not implemented: use_tfh' endif if (use_sdh .or. use_tdh) then stop ' [direct4p] not implemented: use_sdh, tdh' endif ! tide interaction on resolved star-basis (may be overkill) ! TODO - this seems incomplete if (use_sdt .or. use_tdt .or. use_tft) then ! TODO - add star-wise tide-orbit feedback? stop ' [direct4p] not implemented: use_sdt, use_tdt, use_tft' ! TODO - ensure we are not using hierarchical approach concurrently 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 ! TODO - add tidal friction? 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 ! say 1-->3-->4 stop ' [direct4p] not implemented: use_t3t' endif if (use_t3h) then ! 3-tide of binay on cms of other binary ! say 1-->2-->(3+4) - NO - we do not have QP on orbit stop ' [direct4p] not implemented: use_t3h' endif if (use_t3s) then ! 3-tide of binay on individual star in other binary ! say 1-->2-->3 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 (use_ori) then if (use_orq) then yp(31:46) = reshape(roq_p, (/nqstar/)) else yp(31:42) = reshape(rot_p, (/nvstar/)) endif endif if (local_store) then ! actual derivaties local_data( 1: size(yp,1)) = yp(:) ! basic forces and torques local_data( 47: 58) = reshape(fsd_p, (/12/)) local_data( 59: 70) = reshape(ftd_p, (/12/)) local_data( 71: 82) = reshape(ftf_p, (/12/)) local_data( 83: 88) = reshape(fpn_d(1:3,1:2), (/6/)) local_data( 89:100) = reshape(fpn_c, (/12/)) local_data(101:103) = fpn_q(1:3) local_data(104:115) = reshape(ftd3_p, (/12/)) ! permanent quadrupole momenent ! TODO 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, 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_rxv, & use_pnf, use_pnh, use_tfh, use_sdh, use_tdh, use_td3, & 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 = 4, & nquat = 4, & norb = 3, & ndim = 3, & nbin = 2 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_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, an_ij, phi_ij, eps_ij, Gm_ij integer(kind=int32), dimension(nstar), parameter :: & v2b_v = (/(max(1, j-1),j=1,nstar)/) integer(kind=int32), dimension(nbin,norb), parameter :: & b2v_b = reshape((/1, 2, 1, 3, 1, 4/), (/2, 3/)) 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, & wre_v, mr_v, mrat_v real(kind=real64), dimension(ndim, nstar) :: & i_v, j_v, w_v, torque_v, ori_v, rot_v, q_v, & fsd_v, ftd_v, ftf_v, wxr_v, & fpq_v, tpq_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, 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, qr 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 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_b(:,:) = reshape(y(1 : nvorb ), (/ndim, norb/)) v_b(:,:) = 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 do concurrent (i=1:nstar) if (use_ori) then ! 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 else 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 if end do wn2_v = sum(w_v**2, dim=1) mi_v(:) = 1.d0 / m_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 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 if (verb_interact) then write(6, "(' [direct4h] orbit escape ',i1)") j endif 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 if (verb_interact) then write(6, "(' [direct4h] object interaction ',i1,' + ',i1)") i,j endif 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 if (verb_interact) then write(6, "(' [direct4h] object collision ',i1,' + ',i1)") i, j endif 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_ori .and. & ((use_pqm .and. any(lpqm_v(1:2))) .or. & (use_pqh .and. any(lpqm_v(3:4))))) then if (use_pqh) then imax = nstar else imax = nbin endif do i=1,imax if (.not.lpqm_v(i)) then fpq_v(:,i) = 0.d0 tpq_v(:,i) = 0.d0 cycle endif j = v2b_v(i) if (i /= 2) then g = G3 * rnm2_b(j) * rnm3_b(j) / mu_b(j) endif f = g * mr_v(i) qr(:) = matmul(a_v(:,:,i), q_v(:,i) * matmul(r_b(:,j), a_v(:,:,i))) fpq_v(:,i) = f * ( & qr(:) - 2.5d0 * rnm2_b(j) * sum(r_b(:,j) * qr(:)) * r_b(:,j)) tpq_v(:,i) = f * cross_product(qr(:), r_b(:,j)) * ftgr_b(j) f_b(:,j) = f_b(:,j) + fpq_v(:,i) torque_v(:,i) = torque_v(:,i) + tpq_v(:,i) end do fpq_v(:,imax+1:nstar) = 0.d0 fpq_v(:,imax+1:nstar) = 0.d0 else fpq_v(:,:) = 0.d0 tpq_v(:,:) = 0.d0 endif ! TODO - impact of PQ in inner binary stars on outer stars ! TODO - impact of PQ in star 3 on star 4 and vise versa 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 .and. any(t_v /= 0.d0)) then if (use_tfh) then imax = nstar jmax = norb else imax = nbin jmax = 1 endif do j=1,jmax w_b(:,j) = cross_product(re_b(:,j), v_b(:,j)) * rnm1_b(j) enddo if (use_rxv) then do j=1, jmax if (any(t_v(b2v_b(:,j)) < 0.d0)) then nm1_b(j) = 1.d0 / max(norm2(w_b(:,j)), & omegalim * sqrt(Gm_b(j) * rnm3_b(j))) endif enddo else do j=1, jmax if (any(t_v(b2v_b(:,j)) < 0.d0)) then eps_b(j) = min(0.5d0 * v2_b(j) - phi_b(j), -epslim * phi_b(j)) an_b(j) = -0.5d0 * Gm_b(j) / eps_b(j) nm1_b(j) = sqrt(an_b(j)**3 / Gm_b(j)) endif end do endif do i=1,imax h = t_v(i) if (h == 0.d0) then ftf_v(:,i) = 0.d0 continue else if (h < 0.d0) then h = -0.5d0 * h * nm1_b(j) endif j = v2b_v(i) if (i /= 2) then ! check signs f = 12.d0 * Gm_b(j) * rnm8_b(j) v(:) = (-2.d0 * vrn_b(j)) * re_b(:,j) - v_b(:,j) end if h = h * f * 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 stop '[td3] needs fixing' 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 = 6.d0 * f * g * rnm2_ij ftd3_c(:,i,j) = & ((3.d0 - 15.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) * mui_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 = 12.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 h = t_v(i) if (use_tft .and. h /= 0.d0) then w_ij(:) = cross_product(v_c(:,i,j), r_c(:,i,j)) * rnm2_c(i,j) if (h < 0.d0) then Gm_ij = GRAV * (m_v(i) + m_v(j)) if (use_rxv) then nm1_ij = 1.d0 / max(norm2(w_ij(:)), & omegalim * sqrt(Gm_ij * rnm3_c(i,j))) else phi_ij = -Gm_ij * rnm1_c(i,j) eps_ij = min(0.5d0 * sum(v_c(:,i,j)**2) + phi_ij, & epslim * phi_ij) an_ij = -0.5d0 * Gm_ij / eps_ij nm1_ij = sqrt(an_ij**3 / Gm_ij) endif h = -0.5d0 * nm1_ij * h endif rnm8_ij = rnm4_ij**2 g1 = -2.d0 * vrn_c(i,j) v(:) = re_c(:,i,j) * g1 - v_c(:,i,j) h = h * f1 * s5k_v(i) * 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 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(:,:), (/nvorb/)) yp(10:18) = reshape(acc_b(:,:), (/nvorb/)) yp(19:30) = reshape(torque_v(:,:), (/nvstar/)) if (use_ori) then if (use_orq) then yp(31:46) = reshape(roq_v, (/nqstar/)) else yp(31:42) = 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( 47: 58) = reshape(fsd_v (:,:), (/ndim*nstar/)) local_data( 59: 70) = reshape(ftd_v (:,:), (/ndim*nstar/)) local_data( 71: 82) = reshape(ftf_v (:,:), (/ndim*nstar/)) local_data( 83: 91) = reshape(fpn_b (:,:), (/ndim*norb/)) local_data( 92:109) = reshape(fpn_c (:,:,:), (/ndim*(nstar*(nstar-1)/2)/)) local_data(110:121) = reshape(ftd3_c(:,:,:), (/ndim*nbin*(nstar-nbin)/)) local_data(137:142) = reshape(fsdt_c(:,1:2,3), (/ndim*nbin/)) local_data(143:151) = reshape(fsdt_c(:,1:3,4), (/ndim*(nbin+1)/)) local_data(122:127) = reshape(ftdt_c(:,1:2,3), (/ndim*nbin/)) local_data(128:136) = reshape(ftdt_c(:,1:3,4), (/ndim*(nbin+1)/)) local_data(152:157) = reshape(ftft_c(:,1:2,3), (/ndim*nbin/)) local_data(158:166) = reshape(ftft_c(:,1:3,4), (/ndim*(nbin+1)/)) ! permanent quadrupole momenent local_data(167:178) = reshape(fpq_v, (/nvstar/)) local_data(179:190) = reshape(tpq_v, (/nvstar/)) 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, 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_td3, & use_rxv, use_qp, & 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 implicit none integer(kind=int32), parameter :: & ndim = 3, & nquat = 4, & ncorn = 2, & nstar = 4, & 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, 1,4, 2,3, 2,4, 3,4/), (/ncorn, nedge/)) ! edge 1 2 3 4 5 6 ! edges belonging to each star (ncon=nstar-1) integer(kind=int32), dimension(ncon, nstar), parameter :: & ise = reshape((/1,2,3, 1,4,5, 2,4,6, 3,5,6/), (/ncon, nstar/)) ! star 1 2 3 4 ! corners belonging to each star (ncon=nstar-1) integer(kind=int32), dimension(ncon, nstar), parameter :: & isc = reshape((/1,1,1, 2,1,1, 2,2,1, 2,2,2/), (/ncon, nstar/)) ! star 1 2 3 4 ! edges belonging 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, & ONETHIRD = 1.d0 / 3.d0, & G3 = GRAV * 3.d0 ! 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(nquat, nstar) :: & orq_v, roq_v real(kind=real64), dimension(ndim, ndim, nstar) :: & a_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(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(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, 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, & Gm, eps, an, nm1 real(kind=real64), dimension(ndim) :: & v, w, rik, rij, qr 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, fpq_c, tpq_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 ! not supported (but redundant) ! [[,[,]],] +2,+1,-2,-1 ! [,[,[,]]] +1,+1,+1,-3 ! 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' 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 q < 0 as flag for containing -1/Q, tau for q > 0 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/)) 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 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 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 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(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 enddo 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) ! 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, "(' [direct4f] 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 ! fail on star collisions do l=1,nedge if (rn_e(l) < s_v(ies(1,l)) + s_v(ies(2,l))) then if (verb_interact) then write(6, "(' [direct4f] 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 ! fail on star escape if (has_cutoff) then do j=1, norb if (norm2(r_j(:,j)) > cutoff(j)) then if (verb_interact) then write(6, "(' [direct4f] orbit escape ',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 ! 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 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 ! permanent quadrupole-quadrupole interaction if (use_ori .and. use_pqq .and. any(lpqm_v)) then error stop '[direct4f] use_pqq not implemented' 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,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 * 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.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 = 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 net torque, cancels out 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.0d0 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 end if do j=1, ncorn i = ies(j,l) f = t_v(i) if (f == 0.0d0) then ftf_c(:,j,l) = 0.d0 continue else if (f < 0.0d0) 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 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 (use_ori) then if (use_orq) then yp(31:46) = reshape(roq_v, (/nqstar/)) else yp(31:42) = 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( 47: 82) = reshape(fsd_c, (/ndim*ncorn*nedge/)) !36 local_data( 83:118) = reshape(ftd_c, (/ndim*ncorn*nedge/)) !36 local_data(119:154) = reshape(ftf_c, (/ndim*ncorn*nedge/)) !36 local_data(155:226) = reshape(ftd3_d, (/ndim*nstar*nstar1*nstar2/)) !72 local_data(227:244) = reshape(fpn_e, (/ndim*nedge/)) !18 local_data(245:262) = reshape(fng_e, (/ndim*nedge/)) !18 ! permanent quadrupole momenent local_data(263:298) = reshape(fpq_c, (/ndim*ncorn*nedge/)) !36 local_data(299:334) = reshape(tpq_c, (/ndim*ncorn*nedge/)) !36 endif end function direct4f end module quaternary