






















































































































 


































































































































































 





































































































































































































 






















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                

















































































      function brillq(rho,z,phi)

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c

      implicit none

      REAL*8 brill3d_d
      integer, parameter :: cctki_use_brill3d_d = kind(brill3d_d)
      REAL*8 brill3d_e
      integer, parameter :: cctki_use_brill3d_e = kind(brill3d_e)
      REAL*8 brill3d_m
      integer, parameter :: cctki_use_brill3d_m = kind(brill3d_m)
      REAL*8 brill3d_n
      integer, parameter :: cctki_use_brill3d_n = kind(brill3d_n)
      REAL*8 brill3d_phi0
      integer, parameter :: cctki_use_brill3d_phi0 = kind(brill3d_phi0)
      REAL*8 const_v0
      integer, parameter :: cctki_use_const_v0 = kind(const_v0)
      REAL*8 eppley_a
      integer, parameter :: cctki_use_eppley_a = kind(eppley_a)
      REAL*8 eppley_r0
      integer, parameter :: cctki_use_eppley_r0 = kind(eppley_r0)
      REAL*8 eppley_sigmar
      integer, parameter :: cctki_use_eppley_sigmar = kind(eppley_sigmar
     &)
      REAL*8 eppley_sigmarho
      integer, parameter :: cctki_use_eppley_sigmarho = kind(eppley_sigm
     &arho)
      REAL*8 exp_a
      integer, parameter :: cctki_use_exp_a = kind(exp_a)
      REAL*8 exp_rho0
      integer, parameter :: cctki_use_exp_rho0 = kind(exp_rho0)
      REAL*8 exp_sigmaz
      integer, parameter :: cctki_use_exp_sigmaz = kind(exp_sigmaz)
      REAL*8 gundlach_a
      integer, parameter :: cctki_use_gundlach_a = kind(gundlach_a)
      REAL*8 gundlach_r0
      integer, parameter :: cctki_use_gundlach_r0 = kind(gundlach_r0)
      REAL*8 gundlach_sigmar
      integer, parameter :: cctki_use_gundlach_sigmar = kind(gundlach_si
     &gmar)
      REAL*8 gundlach_sigmarho
      integer, parameter :: cctki_use_gundlach_sigmarho = kind(gundlach_
     &sigmarho)
      REAL*8 rhofudge
      integer, parameter :: cctki_use_rhofudge = kind(rhofudge)
      REAL*8 robin_inf
      integer, parameter :: cctki_use_robin_inf = kind(robin_inf)
      REAL*8 thresh
      integer, parameter :: cctki_use_thresh = kind(thresh)
      integer*8 bound
      integer, parameter :: cctki_use_bound = kind(bound)
      integer*8 q_function
      integer, parameter :: cctki_use_q_function = kind(q_function)
      integer*8 solver
      integer, parameter :: cctki_use_solver = kind(solver)
      INTEGER*4 eppley_b
      integer, parameter :: cctki_use_eppley_b = kind(eppley_b)
      INTEGER*4 eppley_c
      integer, parameter :: cctki_use_eppley_c = kind(eppley_c)
      INTEGER*4 exp_b
      integer, parameter :: cctki_use_exp_b = kind(exp_b)
      INTEGER*4 gundlach_b
      integer, parameter :: cctki_use_gundlach_b = kind(gundlach_b)
      INTEGER*4 gundlach_c
      integer, parameter :: cctki_use_gundlach_c = kind(gundlach_c)
      INTEGER*4 output_coeffs
      integer, parameter :: cctki_use_output_coeffs = kind(output_coeffs
     &)
      INTEGER*4 robin_falloff
      integer, parameter :: cctki_use_robin_falloff = kind(robin_falloff
     &)
      INTEGER*4 sor_maxit
      integer, parameter :: cctki_use_sor_maxit = kind(sor_maxit)
      COMMON /IDBrillDatapriv/brill3d_d, brill3d_e, brill3d_m, brill3d_n
     &, brill3d_phi0, const_v0, eppley_a, eppley_r0, eppley_sigmar, eppl
     &ey_sigmarho, exp_a, exp_rho0, exp_sigmaz, gundlach_a, gundlach_r0,
     & gundlach_sigmar, gundlach_sigmarho, rhofudge, robin_inf, thresh, 
     &bound, q_function, solver, eppley_b, eppley_c, exp_b, gundlach_b, 
     &gundlach_c, output_coeffs, robin_falloff, sor_maxit
      integer*8 CCTKH0
      integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
      integer*8 CCTKH1
      integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
      integer*8 CCTKH2
      integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
      integer*8 CCTKH3
      integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
      integer*8 initial_data
      integer, parameter :: cctki_use_initial_data = kind(initial_data)
      integer*8 CCTKH4
      integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
      integer*8 CCTKH5
      integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
      integer*8 CCTKH6
      integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
      integer*8 CCTKH7
      integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
      integer*8 CCTKH8
      integer, parameter :: cctki_use_CCTKH8 = kind(CCTKH8)
      integer*8 CCTKH9
      integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
      integer*8 CCTKH11
      integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
      integer*8 metric_type
      integer, parameter :: cctki_use_metric_type = kind(metric_type)
      integer*8 CCTKH13
      integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
      integer*8 CCTKH14
      integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
      INTEGER*4 CCTKH10
      integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
      INTEGER*4 CCTKH12
      integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
      INTEGER*4 CCTKH15
      integer, parameter :: cctki_use_CCTKH15 = kind(CCTKH15)
      COMMON /ADMBASErest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, initial_data, 
     &CCTKH4, CCTKH5, CCTKH6, CCTKH7, CCTKH8, CCTKH9, CCTKH11, metric_ty
     &pe, CCTKH13, CCTKH14, CCTKH10, CCTKH12, CCTKH15
      integer*8 conformal_storage
      integer, parameter :: cctki_use_conformal_storage = kind(conformal
     &_storage)
      COMMON /STATICCONFORMALrest/conformal_storage
      interface
      integer*8 function CCTK_PointerTo (var)
      implicit none
      type(*),dimension(..),TARGET :: var
      end function CCTK_PointerTo
      end interface
      interface
      integer function CCTK_Equals (arg1, arg2)
      implicit none
      integer*8 arg1
      character(*) arg2
      end function CCTK_Equals
      integer function CCTK_MyProc (cctkGH)
      implicit none
      integer*8 cctkGH
      end function CCTK_MyProc
      integer function CCTK_nProcs (cctkGH)
      implicit none
      integer*8 cctkGH
      end function CCTK_nProcs
      integer function CCTK_IsThornActive (name)
      implicit none
      character(*) name
      end function CCTK_IsThornActive
      integer*8 function CCTK_NullPointer ()
      implicit none
      end function CCTK_NullPointer
      end interface

      logical firstcall

      integer qtype

      REAL*8 brillq,rho,z,phi

      data firstcall /.true./

      save firstcall,qtype

c

      if (firstcall) then

         if ((CCTK_Equals(q_function,"exp") .ne. 0)) then
            qtype = 0
         else if ((CCTK_Equals(q_function,"eppley") .ne. 0)) then
            qtype = 1
         else if ((CCTK_Equals(q_function,"gundlach") .ne. 0)) then
            qtype = 2
         else
            call CCTK_Warn(0,77,"brillq.F","IDBrillData","Brill wave dat
     &a type not recognised")
         end if

         firstcall = .false.

      end if

      if (rho < 0) then
         rho = -rho
      end if

c

c

      if (qtype.eq.0) then

         brillq = exp_a*abs(rho)**(2 + exp_b)
     .          / exp((rho - exp_rho0)**2)/(rho**2 + z**2)

         if (exp_sigmaz.ne.0.0D0) then
            brillq = brillq/exp((z/exp_sigmaz)**2)
         end if

      else if (qtype.eq.1) then

c
c

         brillq = eppley_a*(rho/eppley_sigmarho)**eppley_b 
     &          / ( 1.0D0 + ((rho**2 + z**2 - eppley_r0**2)/
     &          eppley_sigmar**2)**(eppley_c/2) )

      else if (qtype.eq.2) then

c
c
c

         brillq = gundlach_a*(rho/gundlach_sigmarho)**gundlach_b
     &          / exp(((rho**2 + z**2 - gundlach_r0**2)/
     &            gundlach_sigmar**2)**(gundlach_c/2))

      else

c

         call CCTK_Warn(0,124,"brillq.F","IDBrillData","brillq: Unknown 
     &type of Brill function q")

      end if

c

      if ((CCTK_Equals(initial_data,"brilldata") .ne. 0)) then
         brillq = brillq*(1.0D0 + brill3d_d*rho**brill3d_m*
     &        cos(brill3D_n*(phi-brill3d_phi0))**2
     &        / (1.0D0 + brill3d_e*rho**brill3d_m))
      end if

      return
      end
