aboutsummaryrefslogtreecommitdiff
path: root/src/qlm_multipoles.F90
diff options
context:
space:
mode:
Diffstat (limited to 'src/qlm_multipoles.F90')
-rw-r--r--src/qlm_multipoles.F90323
1 files changed, 323 insertions, 0 deletions
diff --git a/src/qlm_multipoles.F90 b/src/qlm_multipoles.F90
new file mode 100644
index 0000000..79b75d5
--- /dev/null
+++ b/src/qlm_multipoles.F90
@@ -0,0 +1,323 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+
+
+subroutine qlm_multipoles (CCTK_ARGUMENTS, hn)
+ use cctk
+ use constants
+ use qlm_boundary
+ use qlm_derivs
+ use qlm_variables
+ use ricci2
+ use tensor2
+ implicit none
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+ integer :: hn
+
+ CCTK_REAL, parameter :: one=1, two=2
+ CCTK_REAL, parameter :: half=one/2, fourth=one/4, eighth=one/8
+ CCTK_REAL, parameter :: sixteenth=one/16
+ CCTK_REAL, parameter :: o128=one/128
+
+ CCTK_REAL :: qq(2,2), dtq, rsc
+ CCTK_REAL :: ee(3,2)
+ CCTK_REAL :: kk(3,3)
+ CCTK_REAL :: ll(0:3), nn(0:3), ss(0:3)
+ CCTK_COMPLEX :: psi2
+ CCTK_REAL :: zz, dzz(2)
+ CCTK_REAL :: area, mass, spin
+
+ CCTK_REAL :: delta_space(2)
+
+ integer :: i, j
+!!$ integer :: a, b, c, d
+
+ if (veryverbose/=0) then
+ call CCTK_INFO ("Calculating multipole moments")
+ end if
+
+ delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /)
+
+ qlm_mp_m0(hn) = 0
+ qlm_mp_m1(hn) = 0
+ qlm_mp_m2(hn) = 0
+ qlm_mp_m3(hn) = 0
+ qlm_mp_m4(hn) = 0
+ qlm_mp_m5(hn) = 0
+ qlm_mp_m6(hn) = 0
+ qlm_mp_m7(hn) = 0
+ qlm_mp_m8(hn) = 0
+
+ qlm_mp_j0(hn) = 0
+ qlm_mp_j1(hn) = 0
+ qlm_mp_j2(hn) = 0
+ qlm_mp_j3(hn) = 0
+ qlm_mp_j4(hn) = 0
+ qlm_mp_j5(hn) = 0
+ qlm_mp_j6(hn) = 0
+ qlm_mp_j7(hn) = 0
+ qlm_mp_j8(hn) = 0
+
+ do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn)
+ do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(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)
+
+ 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)
+
+ 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)
+
+ 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)
+
+ ss = (ll - nn) / sqrt(two)
+
+ rsc = qlm_rsc(i,j,hn)
+ psi2 = qlm_psi2(i,j,hn)
+
+ zz = qlm_inv_z(i,j,hn)
+
+ dzz(:) = deriv (qlm_inv_z(:,:,hn), i, j, delta_space)
+
+ area = sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn)
+
+ mass = fourth * rsc
+
+ qlm_mp_m0(hn) = qlm_mp_m0(hn) + mass * p0(zz) * dzz(1) * area
+ qlm_mp_m1(hn) = qlm_mp_m1(hn) + mass * p1(zz) * dzz(1) * area
+ qlm_mp_m2(hn) = qlm_mp_m2(hn) + mass * p2(zz) * dzz(1) * area
+ qlm_mp_m3(hn) = qlm_mp_m3(hn) + mass * p3(zz) * dzz(1) * area
+ qlm_mp_m4(hn) = qlm_mp_m4(hn) + mass * p4(zz) * dzz(1) * area
+ qlm_mp_m5(hn) = qlm_mp_m5(hn) + mass * p5(zz) * dzz(1) * area
+ qlm_mp_m6(hn) = qlm_mp_m6(hn) + mass * p6(zz) * dzz(1) * area
+ qlm_mp_m7(hn) = qlm_mp_m7(hn) + mass * p7(zz) * dzz(1) * area
+ qlm_mp_m8(hn) = qlm_mp_m8(hn) + mass * p8(zz) * dzz(1) * area
+
+ spin = - aimag(psi2)
+
+ qlm_mp_j0(hn) = qlm_mp_j0(hn) + spin * p0(zz) * dzz(1) * area
+ qlm_mp_j1(hn) = qlm_mp_j1(hn) + spin * p1(zz) * dzz(1) * area
+ qlm_mp_j2(hn) = qlm_mp_j2(hn) + spin * p2(zz) * dzz(1) * area
+ qlm_mp_j3(hn) = qlm_mp_j3(hn) + spin * p3(zz) * dzz(1) * area
+ qlm_mp_j4(hn) = qlm_mp_j4(hn) + spin * p4(zz) * dzz(1) * area
+ qlm_mp_j5(hn) = qlm_mp_j5(hn) + spin * p5(zz) * dzz(1) * area
+ qlm_mp_j6(hn) = qlm_mp_j6(hn) + spin * p6(zz) * dzz(1) * area
+ qlm_mp_j7(hn) = qlm_mp_j7(hn) + spin * p7(zz) * dzz(1) * area
+ qlm_mp_j8(hn) = qlm_mp_j8(hn) + spin * p8(zz) * dzz(1) * area
+
+!!$ spin = 0
+!!$ do a=1,2
+!!$ do b=1,2
+!!$ do c=1,3
+!!$ do d=1,3
+!!$ spin = spin - half * area * epsilon2(a,b) * dzz(b) * ee(d,a) * kk(d,c) * ss(c)
+!!$ end do
+!!$ end do
+!!$ end do
+!!$ end do
+!!$
+!!$ qlm_mp_j0(hn) = qlm_mp_j0(hn) + spin * dp0(zz) * area
+!!$ qlm_mp_j1(hn) = qlm_mp_j1(hn) + spin * dp1(zz) * area
+!!$ qlm_mp_j2(hn) = qlm_mp_j2(hn) + spin * dp2(zz) * area
+!!$ qlm_mp_j3(hn) = qlm_mp_j3(hn) + spin * dp3(zz) * area
+!!$ qlm_mp_j4(hn) = qlm_mp_j4(hn) + spin * dp4(zz) * area
+!!$ qlm_mp_j5(hn) = qlm_mp_j5(hn) + spin * dp5(zz) * area
+!!$ qlm_mp_j6(hn) = qlm_mp_j6(hn) + spin * dp6(zz) * area
+!!$ qlm_mp_j7(hn) = qlm_mp_j7(hn) + spin * dp7(zz) * area
+!!$ qlm_mp_j8(hn) = qlm_mp_j8(hn) + spin * dp8(zz) * area
+
+ end do
+ end do
+
+ ! Normalise
+
+!!$ ! This is the normalisation for I_n and L_n
+!!$ qlm_mp_m0(hn) = qlm_mp_m0(hn) / sqrt(4*pi/ 1)
+!!$ qlm_mp_m1(hn) = qlm_mp_m1(hn) / sqrt(4*pi/ 3)
+!!$ qlm_mp_m2(hn) = qlm_mp_m2(hn) / sqrt(4*pi/ 5)
+!!$ qlm_mp_m3(hn) = qlm_mp_m3(hn) / sqrt(4*pi/ 7)
+!!$ qlm_mp_m4(hn) = qlm_mp_m4(hn) / sqrt(4*pi/ 9)
+!!$ qlm_mp_m5(hn) = qlm_mp_m5(hn) / sqrt(4*pi/11)
+!!$ qlm_mp_m6(hn) = qlm_mp_m6(hn) / sqrt(4*pi/13)
+!!$ qlm_mp_m7(hn) = qlm_mp_m7(hn) / sqrt(4*pi/15)
+!!$ qlm_mp_m8(hn) = qlm_mp_m8(hn) / sqrt(4*pi/17)
+!!$
+!!$ qlm_mp_j0(hn) = qlm_mp_j0(hn) / sqrt(4*pi/ 1)
+!!$ qlm_mp_j1(hn) = qlm_mp_j1(hn) / sqrt(4*pi/ 3)
+!!$ qlm_mp_j2(hn) = qlm_mp_j2(hn) / sqrt(4*pi/ 5)
+!!$ qlm_mp_j3(hn) = qlm_mp_j3(hn) / sqrt(4*pi/ 7)
+!!$ qlm_mp_j4(hn) = qlm_mp_j4(hn) / sqrt(4*pi/ 9)
+!!$ qlm_mp_j5(hn) = qlm_mp_j5(hn) / sqrt(4*pi/11)
+!!$ qlm_mp_j6(hn) = qlm_mp_j6(hn) / sqrt(4*pi/13)
+!!$ qlm_mp_j7(hn) = qlm_mp_j7(hn) / sqrt(4*pi/15)
+!!$ qlm_mp_j8(hn) = qlm_mp_j8(hn) / sqrt(4*pi/17)
+
+ ! This is the normalisation for M_n and J_n
+ qlm_mp_m0(hn) = qlm_mp_m0(hn) * qlm_mass(hn) * qlm_radius(hn)**0 / (2*pi)
+ qlm_mp_m1(hn) = qlm_mp_m1(hn) * qlm_mass(hn) * qlm_radius(hn)**1 / (2*pi)
+ qlm_mp_m2(hn) = qlm_mp_m2(hn) * qlm_mass(hn) * qlm_radius(hn)**2 / (2*pi)
+ qlm_mp_m3(hn) = qlm_mp_m3(hn) * qlm_mass(hn) * qlm_radius(hn)**3 / (2*pi)
+ qlm_mp_m4(hn) = qlm_mp_m4(hn) * qlm_mass(hn) * qlm_radius(hn)**4 / (2*pi)
+ qlm_mp_m5(hn) = qlm_mp_m5(hn) * qlm_mass(hn) * qlm_radius(hn)**5 / (2*pi)
+ qlm_mp_m6(hn) = qlm_mp_m6(hn) * qlm_mass(hn) * qlm_radius(hn)**6 / (2*pi)
+ qlm_mp_m7(hn) = qlm_mp_m7(hn) * qlm_mass(hn) * qlm_radius(hn)**7 / (2*pi)
+ qlm_mp_m8(hn) = qlm_mp_m8(hn) * qlm_mass(hn) * qlm_radius(hn)**8 / (2*pi)
+
+ qlm_mp_j0(hn) = qlm_mp_j0(hn) * qlm_radius(hn)**1 / (4*pi)
+ qlm_mp_j1(hn) = qlm_mp_j1(hn) * qlm_radius(hn)**2 / (4*pi)
+ qlm_mp_j2(hn) = qlm_mp_j2(hn) * qlm_radius(hn)**3 / (4*pi)
+ qlm_mp_j3(hn) = qlm_mp_j3(hn) * qlm_radius(hn)**4 / (4*pi)
+ qlm_mp_j4(hn) = qlm_mp_j4(hn) * qlm_radius(hn)**5 / (4*pi)
+ qlm_mp_j5(hn) = qlm_mp_j5(hn) * qlm_radius(hn)**6 / (4*pi)
+ qlm_mp_j6(hn) = qlm_mp_j6(hn) * qlm_radius(hn)**7 / (4*pi)
+ qlm_mp_j7(hn) = qlm_mp_j7(hn) * qlm_radius(hn)**8 / (4*pi)
+ qlm_mp_j8(hn) = qlm_mp_j8(hn) * qlm_radius(hn)**9 / (4*pi)
+
+contains
+
+ ! Legendre polynomials
+
+ function p0 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p0
+ p0 = 1
+ end function p0
+
+ function p1 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p1
+ p1 = z
+ end function p1
+
+ function p2 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p2
+ p2 = 3*half * z**2 - half
+ end function p2
+
+ function p3 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p3
+ p3 = 5*half * z**3 - 3*half * z
+ end function p3
+
+ function p4 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p4
+ p4 = 35*eighth * z**4 - 15*fourth * z**2 + 3*eighth
+ end function p4
+
+ function p5 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p5
+ p5 = 63*eighth * z**5 - 35*fourth * z**3 + 15*eighth * z
+ end function p5
+
+ function p6 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p6
+ p6 = 231*sixteenth * z**6 - 315*sixteenth * z**4 + 105*sixteenth * z**2 &
+ - 5*sixteenth
+ end function p6
+
+ function p7 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p7
+ p7 = 429*sixteenth * z**7 - 693*sixteenth * z**5 + 315*sixteenth * z**3 &
+ - 35*sixteenth * z
+ end function p7
+
+ function p8 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: p8
+ p8 = 6435*o128 * z**8 - 12012*o128 * z**6 + 6930*o128 * z**4 &
+ - 1260*o128 * z**2 + 35*o128
+ end function p8
+
+ ! Derivatives of the Legendre polynomials
+
+ function dp0 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp0
+ dp0 = 0
+ end function dp0
+
+ function dp1 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp1
+ dp1 = 1
+ end function dp1
+
+ function dp2 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp2
+ dp2 = 3 * z
+ end function dp2
+
+ function dp3 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp3
+ dp3 = 15*half * z**2 - 3*half
+ end function dp3
+
+ function dp4 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp4
+ dp4 = 35*half * z**3 - 15*half*z
+ end function dp4
+
+ function dp5 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp5
+ dp5 = 315*eighth * z**4 - 105*fourth * z**2 - 15*eighth
+ end function dp5
+
+ function dp6 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp6
+ dp6 = 693*eighth * z**5 - 315*fourth * z**3 + 105*eighth * z
+ end function dp6
+
+ function dp7 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp7
+ dp7 = 3003*sixteenth * z**6 - 3465*sixteenth * z**4 + 945*sixteenth * z**2 &
+ - 35*sixteenth
+ end function dp7
+
+ function dp8 (z)
+ CCTK_REAL, intent(in) :: z
+ CCTK_REAL :: dp8
+ dp8 = 51480*o128 * z**7 - 72072*o128 * z**5 + 27720*o128 * z**3 &
+ - 2520*o128 * z
+ end function dp8
+
+end subroutine qlm_multipoles
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 17:12:17 +0000