diff options
Diffstat (limited to 'src/spin_integrand3D.F')
| -rw-r--r-- | src/spin_integrand3D.F | 193 |
1 files changed, 193 insertions, 0 deletions
diff --git a/src/spin_integrand3D.F b/src/spin_integrand3D.F new file mode 100644 index 0000000..3f1afc3 --- /dev/null +++ b/src/spin_integrand3D.F @@ -0,0 +1,193 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz, + & gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & spin_int,Psi,Psi_power) + +c ------------------------------------------------------------------------ +c +c Estimates the spin at a given radius using Equation (12) from +c <Time-asymmetric initial data for black holes and black-hole +c collisions>, Bowen and York, Phys. Rev. D, 21 2047(1980) +c +c ----------------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & spin_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip,count + + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz, + & Pi,idet,p, + & term1,term2 + + data count / 1 / + save count +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + + +c Integrands +c ---------- + + if (count.eq.1) then + term1 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)) + & +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxy*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uyz*(ux*hxz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (uy*term1-uz*term2)*rad**3 + + else if (count.eq.2) then + term1 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uxy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uxz*(uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)) + & +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (uz*term1-ux*term2)*rad**3 + + else if (count.eq.3) then + term1 = uxy*(ux*hxx(i,j,k)+uz*hxz(i,j,k)) + & +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uyz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uxy*(uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uxz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (ux*term1-uy*term2)*rad**3 + + end if + + + ELSE + + spin_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + spin_int(1,j,k) = spin_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + spin_int(i,1,k) = spin_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + spin_int(i,j,1) = spin_int(i,j,2) + ENDDO + ENDDO + + +c setting counter +c --------------- + + if (count.eq.3) then + count = 1 + else + count = count + 1 + end if + + +c end +c --- + + END SUBROUTINE spin_integrand3D |