#include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" #include "cctk_Parameters.h" subroutine qlm_analyse (CCTK_ARGUMENTS, hn) use adm_metric_simple use cctk use constants use qlm_derivs use qlm_variables use tensor use tensor2 implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS DECLARE_CCTK_PARAMETERS integer :: hn CCTK_REAL, parameter :: zero=0, two=2 integer, parameter :: rk = kind(zero) CCTK_REAL :: theta, phi CCTK_REAL :: ll(0:3), nn(0:3) CCTK_COMPLEX :: mm(0:3) CCTK_REAL :: ss(0:3), tt(0:3) CCTK_COMPLEX :: npalpha, npbeta CCTK_REAL :: gg(3,3), dgg(3,3,3), kk(3,3), gg_dot(3,3) CCTK_REAL :: g4(0:3,0:3), dg4(0:3,0:3,0:3) CCTK_REAL :: h4(0:3,0:3), dh4(0:3,0:3,0:3) CCTK_REAL :: dtg, gu(3,3), trk, dgu(3,3,3) CCTK_REAL :: xx(3), ee(3,2) CCTK_REAL :: xi(2), xi1(3) CCTK_REAL :: qq(2,2), dtq CCTK_REAL :: adm_energy, adm_mom(3), adm_amom(3) CCTK_REAL :: ll_energy, ll_mom(3), ll_amom(3,3) CCTK_COMPLEX :: ev CCTK_REAL :: spin CCTK_REAL :: npspin CCTK_REAL :: wsspin CCTK_REAL :: tx, ty, tz CCTK_REAL :: xi1_x(3), xi1_y(3), xi1_z(3) CCTK_REAL :: coordspinx, coordspiny, coordspinz CCTK_REAL :: delta_space(2) CCTK_REAL, allocatable :: weights(:) CCTK_REAL :: sum1 integer :: ntheta_inner, nphi_inner integer :: i_eq, j_p0, j_p2 integer :: i, j, l, m integer :: a, b, c character :: msg*1000 if (veryverbose/=0) then call CCTK_INFO ("Calculating spin") end if delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) ! Equatorial circumference i_eq = (qlm_ntheta(hn) + 1) / 2 ! Polar circumference at phi=0 j_p0 = 1+qlm_nghostsphi(hn) ! Polar circumference at phi=pi/2 j_p2 = 1+qlm_nghostsphi(hn) + (qlm_nphi(hn) - 2*qlm_nghostsphi(hn) - 1) / 4 ! Initial values qlm_equatorial_circumference(hn) = 0 qlm_polar_circumference_0(hn) = 0 qlm_polar_circumference_pi_2(hn) = 0 qlm_area(hn) = 0 qlm_spin(hn) = 0 qlm_cvspin(hn) = 0 qlm_npspin(hn) = 0 qlm_wsspin(hn) = 0 qlm_coordspinx(hn) = 0 qlm_coordspiny(hn) = 0 qlm_coordspinz(hn) = 0 qlm_adm_energy(hn) = 0 qlm_adm_momentum_x(hn) = 0 qlm_adm_momentum_y(hn) = 0 qlm_adm_momentum_z(hn) = 0 qlm_adm_angular_momentum_x(hn) = 0 qlm_adm_angular_momentum_y(hn) = 0 qlm_adm_angular_momentum_z(hn) = 0 qlm_ll_energy(hn) = 0 qlm_ll_momentum_x(hn) = 0 qlm_ll_momentum_y(hn) = 0 qlm_ll_momentum_z(hn) = 0 qlm_ll_angular_momentum_x(hn) = 0 qlm_ll_angular_momentum_y(hn) = 0 qlm_ll_angular_momentum_z(hn) = 0 ! Compute weights for spherical integration (see Driscoll and Healy). ! These are the correct weights in a Gauss-Legendre-Senc. ! Also compare qlm_area with AHFinderDirect's area, they are in much ! better agreement than with the previous method using ! delta_theta*delta_phi. ntheta_inner = qlm_ntheta(hn) - 2*qlm_nghoststheta(hn) nphi_inner = qlm_nphi(hn) - 2*qlm_nghostsphi(hn) allocate (weights(1+qlm_nghoststheta(hn) : qlm_ntheta(hn)-qlm_nghoststheta(hn))) do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn) sum1 = 0 do l = 0, (ntheta_inner-1)/2 sum1 = sum1 + sin((2*l+1)*theta)/(2*l+1) end do weights(i) = 8*pi * sum1 / (nphi_inner * ntheta_inner) end do do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn) phi = qlm_origin_phi(hn) + (j-1)*qlm_delta_phi(hn) ! 2-metric on the horizon qq(1,1) = qlm_qtt(i,j,hn) qq(1,2) = qlm_qtp(i,j,hn) qq(2,2) = qlm_qpp(i,j,hn) qq(2,1) = qq(1,2) call calc_2det (qq, dtq) ll(0) = qlm_l0(i,j,hn) ll(1) = qlm_l1(i,j,hn) ll(2) = qlm_l2(i,j,hn) ll(3) = qlm_l3(i,j,hn) nn(0) = qlm_n0(i,j,hn) nn(1) = qlm_n1(i,j,hn) nn(2) = qlm_n2(i,j,hn) nn(3) = qlm_n3(i,j,hn) mm(0) = qlm_m0(i,j,hn) mm(1) = qlm_m1(i,j,hn) mm(2) = qlm_m2(i,j,hn) mm(3) = qlm_m3(i,j,hn) tt = (ll + nn) / sqrt(two) ss = (ll - nn) / sqrt(two) npalpha = qlm_npalpha(i,j,hn) npbeta = qlm_npbeta(i,j,hn) gg(1,1) = qlm_gxx(i,j) gg(1,2) = qlm_gxy(i,j) gg(1,3) = qlm_gxz(i,j) gg(2,2) = qlm_gyy(i,j) gg(2,3) = qlm_gyz(i,j) gg(3,3) = qlm_gzz(i,j) gg(2,1) = gg(1,2) gg(3,1) = gg(1,3) gg(3,2) = gg(2,3) dgg(1,1,1) = qlm_dgxxx(i,j) dgg(1,2,1) = qlm_dgxyx(i,j) dgg(1,3,1) = qlm_dgxzx(i,j) dgg(2,2,1) = qlm_dgyyx(i,j) dgg(2,3,1) = qlm_dgyzx(i,j) dgg(3,3,1) = qlm_dgzzx(i,j) dgg(2,1,1) = dgg(1,2,1) dgg(3,1,1) = dgg(1,3,1) dgg(3,2,1) = dgg(2,3,1) dgg(1,1,2) = qlm_dgxxy(i,j) dgg(1,2,2) = qlm_dgxyy(i,j) dgg(1,3,2) = qlm_dgxzy(i,j) dgg(2,2,2) = qlm_dgyyy(i,j) dgg(2,3,2) = qlm_dgyzy(i,j) dgg(3,3,2) = qlm_dgzzy(i,j) dgg(2,1,2) = dgg(1,2,2) dgg(3,1,2) = dgg(1,3,2) dgg(3,2,2) = dgg(2,3,2) dgg(1,1,3) = qlm_dgxxz(i,j) dgg(1,2,3) = qlm_dgxyz(i,j) dgg(1,3,3) = qlm_dgxzz(i,j) dgg(2,2,3) = qlm_dgyyz(i,j) dgg(2,3,3) = qlm_dgyzz(i,j) dgg(3,3,3) = qlm_dgzzz(i,j) dgg(2,1,3) = dgg(1,2,3) dgg(3,1,3) = dgg(1,3,3) dgg(3,2,3) = dgg(2,3,3) kk(1,1) = qlm_kxx(i,j) kk(1,2) = qlm_kxy(i,j) kk(1,3) = qlm_kxz(i,j) kk(2,2) = qlm_kyy(i,j) kk(2,3) = qlm_kyz(i,j) kk(3,3) = qlm_kzz(i,j) kk(2,1) = kk(1,2) kk(3,1) = kk(1,3) kk(3,2) = kk(2,3) call calc_det (gg, dtg) call calc_inv (gg, dtg, gu) call calc_trace (gu, kk, trk) ! Calculate time derivatives call calc_3metricdot_simple (kk, gg_dot) ! Calculate 4-metric call calc_4metricderivs_simple (gg, dgg, gg_dot, g4,dg4) ! Calculate derivative of 4-metric inverse call calc_invderiv (gu, dgg, dgu) ! "Perturbative" 4-metric h4 = g4 - eta4 dh4 = dg4 xx(1) = qlm_x(i,j,hn) xx(2) = qlm_y(i,j,hn) xx(3) = qlm_z(i,j,hn) ee(1,1:2) = deriv (qlm_x(:,:,hn), i, j, delta_space) ee(2,1:2) = deriv (qlm_y(:,:,hn), i, j, delta_space) ee(3,1:2) = deriv (qlm_z(:,:,hn), i, j, delta_space) xi(1) = qlm_xi_t(i,j,hn) xi(2) = qlm_xi_p(i,j,hn) if (i == i_eq) then qlm_equatorial_circumference(hn) = qlm_equatorial_circumference(hn) & + sqrt(qq(2,2)) * qlm_delta_phi(hn) end if if (j == j_p0) then qlm_polar_circumference_0(hn) = qlm_polar_circumference_0(hn) & + sqrt(qq(1,1)) * qlm_delta_theta(hn) end if if (j == j_p2) then qlm_polar_circumference_pi_2(hn) = qlm_polar_circumference_pi_2(hn) & + sqrt(qq(1,1)) * qlm_delta_theta(hn) end if qlm_area(hn) = qlm_area(hn) & + sqrt(dtq) * weights(i) ! s^i: outward spacelike normal ! K_ij: extrinsic curvature ! phi^i: rotational Killing vector do a=1,3 xi1(a) = 0 do b=1,2 xi1(a) = xi1(a) + ee(a,b) * xi(b) end do end do ! phi^i omega_i = - phi^i s^j K_ij spin = 0 do a=1,3 do b=1,3 spin = spin + xi1(a) * ss(b) * kk(a,b) end do end do qlm_spin_density(i,j) = spin ! phi^i omega_i = (alpha + ~beta) phi^i m_i + complex conjugate npspin = 0 do a=1,3 do b=1,3 npspin = npspin + 2 * real((npalpha + conjg(npbeta)) * xi1(a) * gg(a,b) * mm(b)) end do end do ! phi^i omega_i = 1/2 f Im Psi_2 ! (or is it phi^i omega_i = - f Im Psi_2 ?) wsspin = - qlm_inv_z(i,j,hn) * aimag(qlm_psi2(i,j,hn)) ! x = sin theta cos phi ! y = sin theta sin phi ! z = cos theta ! xi_x = (0,-z,y) ! xi_y = (z,0,-x) ! xi_z = (-y,x,0) tx = qlm_x(i,j,hn) - qlm_origin_x(hn) ty = qlm_y(i,j,hn) - qlm_origin_y(hn) tz = qlm_z(i,j,hn) - qlm_origin_z(hn) xi1_x(:) = (/ zero, -tz, ty /) xi1_y(:) = (/ tz, zero, -tx /) xi1_z(:) = (/ -ty, tx, zero /) coordspinx = 0 coordspiny = 0 coordspinz = 0 do a=1,3 do b=1,3 coordspinx = coordspinx + xi1_x(a) * ss(b) * kk(a,b) coordspiny = coordspiny + xi1_y(a) * ss(b) * kk(a,b) coordspinz = coordspinz + xi1_z(a) * ss(b) * kk(a,b) end do end do xi1_x(:) = (/ 1, 0, 0 /) xi1_y(:) = (/ 0, 1, 0 /) xi1_z(:) = (/ 0, 0, 1 /) qlm_spin(hn) = qlm_spin(hn) & + spin * sqrt(dtq) * weights(i) qlm_npspin(hn) = qlm_npspin(hn) & + npspin * sqrt(dtq) * weights(i) qlm_wsspin(hn) = qlm_wsspin(hn) & + wsspin * sqrt(dtq) * weights(i) qlm_coordspinx(hn) = qlm_coordspinx(hn) & + coordspinx * sqrt(dtq) * weights(i) qlm_coordspiny(hn) = qlm_coordspiny(hn) & + coordspiny * sqrt(dtq) * weights(i) qlm_coordspinz(hn) = qlm_coordspinz(hn) & + coordspinz * sqrt(dtq) * weights(i) ! ADM quantities ! Alcubierre, Appendix A, p. 402 ff, eqns. (A.5) - (A.7): ! ADM energy ! E_adm = (1/16 pi) int_S(r) [delta^jk h_ik,j - (trace(h)),i] n^i r^2 dOmega adm_energy = 0 do a=1,3 do b=1,3 do c=1,3 adm_energy = adm_energy + (delta3(a,b) * dh4(a,c,b) & & - gu(a,b) * dh4(a,b,c) - h4(a,b) * dgu(a,b,c) ) * ss(c) end do end do end do qlm_adm_energy(hn) = qlm_adm_energy(hn) & & + adm_energy / (16*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) ! ADM momentum ! P_adm^i = (1/8 pi) int_S(r) [K^i_j - K delta^i_j] n^j r^2 dOmega do a=1,3 adm_mom(a) = 0 do b=1,3 do c=1,3 adm_mom(a) = adm_mom(a) & & + gu(a,c) * kk(c,b) * ss(b) end do adm_mom(a) = adm_mom(a) & + delta3(a,b) * trk * ss(b) end do end do qlm_adm_momentum_x(hn) = qlm_adm_momentum_x(hn) & & + adm_mom(1) / (8*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) qlm_adm_momentum_y(hn) = qlm_adm_momentum_y(hn) & & + adm_mom(2) / (8*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) qlm_adm_momentum_z(hn) = qlm_adm_momentum_z(hn) & & + adm_mom(3) / (8*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) ! ADM angular momentum ! J_adm^i = (1/8 pi) int_S(r) ! epsilon^ijk x_j (K_kl - delta_kl K) * n^l dS do a=1,3 adm_amom(a) = 0 do b=1,3 do c=1,3 do l=1,3 do m=1,3 adm_amom(a) = adm_amom(a) & & + epsilon3(a,b,c) * xx(b) * gu(c,m) * kk(m,l) * ss(l) end do adm_amom(a) = adm_amom(a) & & - epsilon3(a,b,c) * xx(b) * delta3(c,l) * trk * ss(l) end do end do end do end do qlm_adm_angular_momentum_x(hn) = qlm_adm_angular_momentum_x(hn) & & + adm_amom(1) / (8*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) qlm_adm_angular_momentum_y(hn) = qlm_adm_angular_momentum_y(hn) & & + adm_amom(2) / (8*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) qlm_adm_angular_momentum_z(hn) = qlm_adm_angular_momentum_z(hn) & & + adm_amom(3) / (8*pi) & & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) ! Landau-Lifshitz quantities ! Weinberg, chapter 7.6, pp. 165 ff, eqns. (7.6.22) - (7.6.24): ! Landau-Lifshitz energy ! E_ll = (-1/16 pi) int_S(r) [h_jj,i - hij,j] n_i r^2 dOmega ll_energy = 0 do a=1,3 do b=1,3 ll_energy = ll_energy + (dh4(b,b,a) - dh4(a,b,b)) * ss(a) end do end do qlm_ll_energy(hn) = qlm_ll_energy(hn) & & + ll_energy / (-16*pi) & & * sqrt(dtq) * weights(i) ! Landau-Lifshitz momentum ! P_ll^j = (-1/16 pi) int_S(r) ! [- h_kk,0 delta_ij + h_k0,k delta_ij ! - h_j0,i + h_ij,0 ] n_i r^2 dOmega ! ! P_ll^j = (1/8 pi) int_S(r) [K_ij - K delta_ij] n_i r^2 dOmega do a=1,3 ll_mom(a) = 0 do b=1,3 ll_mom(a) = ll_mom(a) & + (- dh4(b,b,0) + dh4(b,0,b)) * ss(a) & + (- dh4(a,0,b) + dh4(b,a,0)) * ss(b) end do end do qlm_ll_momentum_x(hn) = qlm_ll_momentum_x(hn) & & + ll_mom(1) / (-16*pi) & & * sqrt(dtq) * weights(i) qlm_ll_momentum_y(hn) = qlm_ll_momentum_y(hn) & & + ll_mom(2) / (-16*pi) & & * sqrt(dtq) * weights(i) qlm_ll_momentum_z(hn) = qlm_ll_momentum_z(hn) & & + ll_mom(3) / (-16*pi) & & * sqrt(dtq) * weights(i) ! Landau-Lifshitz angular momentum ! J_ll^jk = (-1/16 pi) int_S(r) ! [- x_j h_0k,i + x_k h_0j,i ! + x_j h_ki,0 - x_k h_ji,0 ! + h_0k delta_ij - h_0j delta_ik] n_i r^2 dOmega do a=1,3 do b=1,3 ll_amom(a,b) = 0 do c=1,3 ll_amom(a,b) = ll_amom(a,b) + ss(c) * ( & + (- xx(a) * dh4(0,b,c) + xx(b) * dh4(0,a,c)) & + (+ xx(a) * dh4(b,c,0) - xx(b) * dh4(a,c,0)) & + (h4(0,b) * delta4(c,a) - h4(0,a) * delta4(c,b))) end do end do end do qlm_ll_angular_momentum_x(hn) = qlm_ll_angular_momentum_x(hn) & & + ll_amom(2,3) / (-16*pi) & & * sqrt(dtq) * weights(i) qlm_ll_angular_momentum_y(hn) = qlm_ll_angular_momentum_y(hn) & & + ll_amom(3,1) / (-16*pi) & & * sqrt(dtq) * weights(i) qlm_ll_angular_momentum_z(hn) = qlm_ll_angular_momentum_z(hn) & & + ll_amom(1,2) / (-16*pi) & & * sqrt(dtq) * weights(i) end do end do deallocate(weights) qlm_polar_circumference_0(hn) = qlm_polar_circumference_0(hn) * 2 qlm_polar_circumference_pi_2(hn) = qlm_polar_circumference_pi_2(hn) * 2 ! A = 4 pi R^2 ! R = 2 M qlm_radius(hn) = sqrt(qlm_area(hn) / (4*pi)) qlm_irreducible_mass(hn) = qlm_radius(hn) / 2 qlm_spin(hn) = qlm_spin(hn) / (8*pi) qlm_mass(hn) = 1/(2*qlm_radius(hn)) * sqrt(qlm_radius(hn)**4 + 4*qlm_spin(hn)**2) qlm_cvspin(hn) = qlm_cvspin(hn) / (8*pi) qlm_npspin(hn) = qlm_npspin(hn) / (-8*pi) qlm_wsspin(hn) = qlm_wsspin(hn) / (-4*pi) qlm_coordspinx(hn) = qlm_coordspinx(hn) / (8*pi) qlm_coordspiny(hn) = qlm_coordspiny(hn) / (8*pi) qlm_coordspinz(hn) = qlm_coordspinz(hn) / (8*pi) ! The event horizon is at r = M + sqrt (M^2 - a^2) ! with x^2 + y^2 + z^2 = rho^2 = r^2 + a^2 (1 - z^2 / r^2) call guess_mass_spin & (qlm_area(hn), qlm_equatorial_circumference(hn), qlm_mass_guess(hn), qlm_spin_guess(hn)) if (verbose/=0) then write (msg, '("Geometric quantities for surface ",i4,":")') hn-1 call CCTK_INFO (msg) write (msg, '(" Area A: ",g16.6)') qlm_area(hn) call CCTK_INFO (msg) write (msg, '(" Irreducible mass M = R/2: ",g16.6)') qlm_irreducible_mass(hn) call CCTK_INFO (msg) write (msg, '(" Areal radius R = sqrt(A/4pi): ",g16.6)') qlm_radius(hn) call CCTK_INFO (msg) write (msg, '("Coordinate-dependent quantities for surface ",i4,":")') hn-1 call CCTK_INFO (msg) write (msg, '(" Equatorial circumference: ",g16.6)') & qlm_equatorial_circumference(hn) call CCTK_INFO (msg) write (msg, '(" Polar circumference at phi=0: ",g16.6)') & qlm_polar_circumference_0(hn) call CCTK_INFO (msg) write (msg, '(" Polar circumference at phi=pi/2: ",g16.6)') & qlm_polar_circumference_pi_2(hn) call CCTK_INFO (msg) if (qlm_spin_guess(hn) >= 0) then write (msg, '(" Spin guess J from distortion: ",g16.6)') & qlm_spin_guess(hn) call CCTK_INFO (msg) else call CCTK_INFO (" No valid spin guess from distortion.") call CCTK_INFO (" (Spin guess is imaginary.)") write (msg, '(" Magnitude of invalid spin guess: ",g16.6)') & abs(qlm_spin_guess(hn)) call CCTK_INFO (msg) end if write (msg, '(" Mass guess M from distortion: ",g16.6)') & qlm_mass_guess(hn) call CCTK_INFO (msg) write (msg, '("Isolated Horizon quantities for surface ",i4,":")') hn-1 call CCTK_INFO (msg) ev = cmplx(qlm_killing_eigenvalue_re(hn), qlm_killing_eigenvalue_im(hn),rk) write (msg, '(" Killing vector field eigenvalue norm: ",g14.6)') abs(ev) call CCTK_INFO (msg) write (msg, '(" Spin J: ",g14.6)') qlm_spin(hn) call CCTK_INFO (msg) write (msg, '(" Kerr spin parameter a = J/M: ",g14.6)') qlm_spin(hn) / qlm_mass(hn) call CCTK_INFO (msg) write (msg, '(" Dimensionless spin parameter a = J/M^2: ",g14.6)') qlm_spin(hn) / qlm_mass(hn)**2 call CCTK_INFO (msg) write (msg, '(" Spin J from NP: ",g14.6)') qlm_npspin(hn) call CCTK_INFO (msg) write (msg, '(" Spin J from phi-coordinate-vector: ",g14.6)') qlm_cvspin(hn) call CCTK_INFO (msg) write (msg, '(" Mass M: ",g14.6)') qlm_mass(hn) call CCTK_INFO (msg) write (msg, '("Global quantities for surface ",i4,":")') hn-1 call CCTK_INFO (msg) write (msg, '(" ADM energy: ",g16.6)') qlm_adm_energy(hn) call CCTK_INFO (msg) write (msg, '(" ADM momentum x: ",g16.6)') qlm_adm_momentum_x(hn) call CCTK_INFO (msg) write (msg, '(" ADM momentum y: ",g16.6)') qlm_adm_momentum_y(hn) call CCTK_INFO (msg) write (msg, '(" ADM momentum z: ",g16.6)') qlm_adm_momentum_z(hn) call CCTK_INFO (msg) write (msg, '(" ADM angular momentum x: ",g16.6)') qlm_adm_angular_momentum_x(hn) call CCTK_INFO (msg) write (msg, '(" ADM angular momentum y: ",g16.6)') qlm_adm_angular_momentum_y(hn) call CCTK_INFO (msg) write (msg, '(" ADM angular momentum z: ",g16.6)') qlm_adm_angular_momentum_z(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz energy: ",g16.6)') qlm_ll_energy(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz momentum x: ",g16.6)') qlm_ll_momentum_x(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz momentum y: ",g16.6)') qlm_ll_momentum_y(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz momentum z: ",g16.6)') qlm_ll_momentum_z(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz angular momentum x: ",g16.6)') qlm_ll_angular_momentum_x(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz angular momentum y: ",g16.6)') qlm_ll_angular_momentum_y(hn) call CCTK_INFO (msg) write (msg, '(" Landau-Lifshitz angular momentum z: ",g16.6)') qlm_ll_angular_momentum_z(hn) call CCTK_INFO (msg) end if contains subroutine guess_mass_spin (area, circumference, mass, spin) CCTK_REAL, intent(in) :: area, circumference CCTK_REAL, intent(out) :: mass, spin CCTK_REAL :: radius, amom2, amom ! equatorial circumference L, area A ! L = 2 pi (r^2 + a^2) / r ! A = 4 pi (r^2 + a^2) ! r = M + sqrt (M^2 - a^2) ! r = A / (2 L) ! a^2 = A / (4 pi) - r^2 ("spin" a = J/M = specific angular momentum) ! M = (r^2 + a^2) / (2 r) ! J = a M (angular momentum) radius = area / (2*circumference) amom2 = area / (4*pi) - radius**2 amom = sign(sqrt(abs(amom2)), amom2) mass = (radius**2 + amom2) / (2*radius) spin = amom * mass ! equatorial circumference L, area A ! mass M, areal radius R, spin a = J/M^2 ! A = 4 pi R^2 ! L = 4 pi M ! a^2 = (R/M)^2 - 1/4 (R/M)^4 ! J = a M^2 end subroutine guess_mass_spin end subroutine qlm_analyse