Cleaning up.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@48 a678b1cf-93e1-4b43-a69d-d43939e66649

3 files changed, 99 insertions, 90 deletions

diff --git a/param.ccl b/param.ccl
index 5f15835..1e65c3b 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,4 +1,6 @@
 # Parameter definitions for thorn IDBrillData
+# $Header$
+
 
 shares:einstein
 
@@ -25,17 +27,17 @@ KEYWORD brill_bound "Which boundary condition to use"
   "robin" :: "Robin boundary: set brill_robin_falloff, brill_robin_inf"
 } "const"
 
-CCTK_REAL brill_const_v0 "Value of constant BC"
-{
-  : :: "anything goes"
-} 1.0
-
-CCTK_INT brill_robin_falloff "Fall-off of Robin BC"
+INT brill_robin_falloff "Fall-off of Robin BC"
 {
   0: :: "any positive integer value"
 } 1
 
-CCTK_REAL brill_robin_inf "Value at infinity of Robin BC"
+REAL brill_const_v0 "Value of constant BC"
+{
+  : :: "anything goes"
+} 1.0
+
+REAL brill_robin_inf "Value at infinity of Robin BC"
 {
   : :: "anything goes"
 } 1.0
@@ -55,22 +57,22 @@ BOOLEAN axisym "Axisymmetric initial data?"
 INT brill_q  "Form of function q [0,1,2]"
 {
   0:2 :: "Only cases 0,1,2 defined at the moment"
-} 1
+} 2
 
 REAL brill_a "Brill wave: Amplitude"
 {
   : :: "Anything"
 } 0.0
 
-REAL brill_b "Brill wave: rho^b"
+INT brill_b "Brill wave: rho^b"
 {
   : :: "Anything"
-} 2.0
+} 2
 
-REAL brill_c "Brill wave: (r^2 - r0^2)^(c/2)"
+INT brill_c "Brill wave: (r^2 - r0^2)^(c/2)"
 {  
   : :: "Anything"
-} 2.0
+} 2
 
 REAL brill_rho0 "Brill wave: radius of torus in rho"
 {
@@ -92,6 +94,11 @@ REAL brill_sr "Brill wave: sigma in r"
   : :: "Anything"
 } 1.0
 
+REAL brill_sz "Brill wave: sigma in z"
+{
+  : :: "Anything"
+} 1.0
+
 # 3D Brill wave parameters
 
 REAL brill_d "3D Brill wave:  d rho^m cos^2(n (phi + phi0))"
@@ -128,13 +135,14 @@ KEYWORD outputpsi "Output conformal factor?"
 
 # Additional parameters
 
-REAL brill_eps "epsilon for finite differencing"
-{
-  0: :: "Positive please"
-} 1.e-6
-
 REAL brill_rhofudge "delta rho for axis fudge"
 {
   0: :: "Positive please"
-} 1.e-6
+} 0.00001
+
+
+
+
+
+
 
diff --git a/src/brillq.F b/src/brillq.F
index 8d62064..ab299ec 100644
--- a/src/brillq.F
+++ b/src/brillq.F
@@ -10,44 +10,42 @@
 
 #include "cctk.h" 
 #include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
 
       function brillq(rho,z,phi)
 
-C     Calculates the function q that appear in the conformal
-C     metric for Brill waves:
-C
-C     ds^2  =  psi^4 ( e^(2q) (drho^2 + dz^2) + rho^2 dphi^2 )
-C
-C     There are three different choices for the form of q depending
-C     on the value of the parameter "brill_q":
-C
-C     brill_q = 0:
-C                                           2
-C                           2+b  - (rho-rho0)      2     2
-C                  q = a rho    e            (z/sz)  /  r
-C
-C
-C     brill_q = 1:  (includes Eppleys form as special case)
-C
-C
-C                                  b               2    2      2  c/2
-C                  q = a (rho/srho)  /  { 1 + [ ( r - r0 ) / sr  ]    }
-C
-C
-C     brill_q = 2:  (includes Holz et al form as special case)
-C
-C                                            2    2      2  c/2 
-C                                  b  - [ ( r - r0 ) / sr  ]
-C                  q = a (rho/srho)  e
-C
-C
-C     If we want 3D initial data, the function q is multiplied by an
-C     additional factor:
-C
-C                           m    2                              m
-C         q -> q [ 1 + d rho  cos  (n (phi)+phi0 ) / ( 1 + e rho ) ]
-
+c     Calculates the function q that appear in the conformal
+c     metric for Brill waves:
+c
+c     ds^2  =  psi^4 ( e^(2q) (drho^2 + dz^2) + rho^2 dphi^2 )
+c
+c     There are three different choices for the form of q depending
+c     on the value of the parameter "brill_q":
+c
+c     brill_q = 0:
+c                                           2
+c                           2+b  - (rho-rho0)      2     2
+c                  q = a rho    e            (z/sz)  /  r
+c
+c
+c     brill_q = 1:  (includes Eppleys form as special case)
+c
+c
+c                                  b               2    2      2  c/2
+c                  q = a (rho/srho)  /  { 1 + [ ( r - r0 ) / sr  ]    }
+c
+c
+c     brill_q = 2:  (includes Holz et al form as special case)
+c
+c                                            2    2      2  c/2 
+c                                  b  - [ ( r - r0 ) / sr  ]
+c                  q = a (rho/srho)  e
+c
+c
+c     If we want 3D initial data, the function q is multiplied by an
+c     additional factor:
+c
+c                           m    2                              m
+c         q -> q [ 1 + d rho  cos  (n (phi)+phi0 ) / ( 1 + e rho ) ]
 
       implicit none
 
@@ -56,17 +54,19 @@ C         q -> q [ 1 + d rho  cos  (n (phi)+phi0 ) / ( 1 + e rho ) ]
       logical firstcall
 
       integer qtype
+      integer b,c
 
-      real*8 brillq,rho,z,phi
-      real*8 a,b,c,r0,sr,rho0,srho,brill_sz
-      real*8 d,e,m,n,phi0
+      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,a,b,c,r0,sr,rho0,srho,brill_sz
+      save firstcall,qtype
+      save a,b,c,r0,sr,rho0,srho,sz
       save d,e,m,n,phi0
 
-C     Get parameters at first call.
+c     Get parameters at first call.
 
       if (firstcall) then
 
@@ -80,7 +80,7 @@ C     Get parameters at first call.
          sr   = brill_sr
          rho0 = brill_rho0
          srho = brill_srho
-         brill_sz = brill_sz
+         sz   = brill_sz
 
          if (axisym.eq.0) then
             d = brill_d
@@ -94,50 +94,51 @@ C     Get parameters at first call.
 
       end if
 
-      if (rho.lt.0) call CCTK_WARN(1,"Warning: negative rho in brillq:")
+c     Calculate q(rho,z) from a choice of forms.
 
-C     Calculate q(rho,z) from a choice of forms. Type 0, 1 and 2 are those
-C     of Bruegmann.
+c     brill_q = 0.
 
-C     brill_q = 0.
+      if (qtype.eq.0) then
 
-      if (qtype .eq. 0) then
+         brillq = a*dabs(rho)**(2 + b)
+     .          / dexp((rho - rho0)**2)/(rho**2 + z**2)
 
-         brillq = a * rho**(2.d0+b)
-     $        / dexp((rho - rho0)**2) / (rho**2 + z**2) 
-         if (brill_sz .ne. 0.d0) then
-            brillq = brillq / exp((z/brill_sz)**2)
+         if (sz.ne.0.0D0) then
+            brillq = brillq/dexp((z/sz)**2)
          end if
 
-      else if (qtype .eq. 1) then
+      else if (qtype.eq.1) then
+
+c     brill_q = 1.  This includes Eppleys choice of q. 
+c     But note that q(Eppley) = 2q(Cactus).
 
-C     brill_q = 1.  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 = a * (rho/srho)**b 
-     $        / ( 1.d0 + ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) )
+      else if (qtype.eq.2) then
 
-      else if (qtype .eq. 2) then
+c     brill_q = 2.   This includes my (Carstens) notion of what a
+c     smooth "pure quadrupole" q should be, plus the choice of
+c     Holz et al.
 
-C     brill_q = 2.   This includes my (Carstens) notion of what a
-C     smooth "pure quadrupole" q should be.
+         brillq = a*(rho/srho)**b
+     .          / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2))
 
-         brillq = a * (rho/srho)**b 
-     $        / dexp( ((rho**2 + z**2 - r0**2) / sr**2)**(c/2) )
+         print *,a,b,c,r0,sr,srho,z,rho,brillq
 
       else
 
-C     Unknown type for q function.
+c     Unknown type for q function.
 
-         call CCTK_WARN(0,"Unknown type of Brill function q.")
+         call CCTK_WARN(0,"Unknown type of Brill function q")
 
       end if
 
-C     If desired, multiply with phi dependent factor.
+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))
+     .          / (1.0D0 + e*rho**m))
       end if
 
       return
diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
index 36060e9..597c3ae 100644
--- a/src/finishbrilldata.F
+++ b/src/finishbrilldata.F
@@ -25,15 +25,15 @@
       integer nx,ny,nz
 
       CCTK_REAL x1,y1,z1,rho1,rho2
-      CCTK_REAL brillq,psi4,e2q,rhofudge
-      CCTK_REAL phi,phif
-      CCTK_REAL zero
+      CCTK_REAL phi,psi4,e2q,rhofudge
+      CCTK_REAL zero,one
 
-      external brillq
+      CCTK_REAL brillq,phif
 
 c     Numbers.
 
       zero = 0.0D0
+      one  = 1.0D0
 
 c     Set up grid size.
 
@@ -74,12 +74,12 @@ c              e^2q (dx^2 + dy^2 + dz^2) + (1-e^2q) (xdy-ydx)^2/rho^2
 c
 c              The individual coefficients can be read off as
 
-               if (rho1 .gt. rhofudge) then
+               if (rho1.gt.rhofudge) then
 
-                  gxx(i,j,k) = psi4*(e2q + (1.d0 - e2q)*y1**2/rho2)
-                  gyy(i,j,k) = psi4*(e2q + (1.d0 - e2q)*x1**2/rho2)
+                  gxx(i,j,k) = psi4*(e2q + (one - e2q)*y1**2/rho2)
+                  gyy(i,j,k) = psi4*(e2q + (one - e2q)*x1**2/rho2)
                   gzz(i,j,k) = psi4*e2q
-                  gxy(i,j,k) = - psi4*(1.d0 - e2q)*x1*y1/rho2
+                  gxy(i,j,k) = - psi4*(one - e2q)*x1*y1/rho2
 
                else