1 files changed, 31 insertions, 38 deletions

diff --git a/src/brillq.F b/src/brillq.F
index bbd9220..df6a65d 100644
--- a/src/brillq.F
+++ b/src/brillq.F
@@ -50,93 +50,86 @@ c         q -> q [ 1 + d rho  cos  (n (phi)+phi0 ) / ( 1 + e rho ) ]
       implicit none
 
       DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
 
       logical firstcall
 
       integer qtype
-      integer b,c
 
       CCTK_REAL brillq,rho,z,phi
-      CCTK_REAL a,r0,sr,rho0,srho,sz
-      CCTK_REAL d,e,m,n,phi0
 
       data firstcall /.true./
 
       save firstcall,qtype
-      save a,b,c,r0,sr,rho0,srho,sz
-      save d,e,m,n,phi0
 
 c     Get parameters at first call.
 
       if (firstcall) then
 
-         qtype = brill_q
-
-         a = brill_a
-         b = brill_b
-         c = brill_c
-
-         r0   = brill_r0
-         sr   = brill_sr
-         rho0 = brill_rho0
-         srho = brill_srho
-         sz   = brill_sz
-
-         if (axisym.eq.0) then
-            d = brill_d
-            e = brill_e
-            m = brill_m
-            n = brill_n
-            phi0 = brill_phi0
+         if (CCTK_EQUALS(q_function,"exp")) then
+            qtype = 0
+         else if (CCTK_EQUALS(q_function,"eppley")) then
+            qtype = 1
+         else if (CCTK_EQUALS(q_function,"gundlach")) then
+            qtype = 2
+         else
+            call CCTK_WARN(0,"Brill wave data type not recognised")
          end if
 
          firstcall = .false.
 
       end if
 
+      if (rho < 0) then
+         rho = -rho
+      end if
+
 c     Calculate q(rho,z) from a choice of forms.
 
-c     brill_q = 0.
+c     "exp" brill data
 
       if (qtype.eq.0) then
 
-         brillq = a*dabs(rho)**(2 + b)
-     .          / dexp((rho - rho0)**2)/(rho**2 + z**2)
+         brillq = exp_a*dabs(rho)**(2 + exp_b)
+     .          / dexp((rho - exp_rho0)**2)/(rho**2 + z**2)
 
-         if (sz.ne.0.0D0) then
-            brillq = brillq/dexp((z/sz)**2)
+         if (exp_sigmaz.ne.0.0D0) then
+            brillq = brillq/dexp((z/exp_sigmaz)**2)
          end if
 
       else if (qtype.eq.1) then
 
-c     brill_q = 1.  This includes Eppleys choice of q. 
+c     "Eppley" brill data.  This includes Eppleys choice of q. 
 c     But note that q(Eppley) = 2q(Cactus).
 
-         brillq = a*(rho/srho)**b 
-     .          / ( 1.0D0 + ((rho**2 + z**2 - r0**2)/sr**2)**(c/2) )
+         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     brill_q = 2.   This includes my (Carstens) notion of what a
+c     "Gundlach" brill data   This includes my (Carstens) notion of what a
 c     smooth "pure quadrupole" q should be, plus the choice of
 c     Holz et al.
 
-         brillq = a*(rho/srho)**b
-     .          / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2))
+         brillq = gundlach_a*(rho/gundlach_sigmarho)**gundlach_b
+     &          / dexp(((rho**2 + z**2 - gundlach_r0**2)/
+     &            gundlach_sigmar**2)**(gundlach_c/2))
 
       else
 
 c     Unknown type for q function.
 
-         call CCTK_WARN(0,"Unknown type of Brill function q")
+         call CCTK_WARN(0,"brillq: Unknown type of Brill function q")
 
       end if
 
 c     If desired, multiply with a phi dependent factor.
 
-      if (axisym.eq.0) then
-         brillq = brillq*(1.0D0 + d*rho**m*cos(n*(phi-phi0))**2
-     .          / (1.0D0 + e*rho**m))
+      if (CCTK_EQUALS(initial_data,"brilldata")) 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