fixed of bug in schedule.ccl.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiOddBrillBH/trunk@51 b6f3ac56-194f-0410-8878-cdf6079d7f1b

5 files changed, 90 insertions, 39 deletions

diff --git a/interface.ccl b/interface.ccl
index d8c3ed9..2f29f81 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,7 +2,7 @@
 # $Header$
 
 implements: axiodd
-inherits: admbase staticconformal grid admanalysis driver
+inherits: ADMBase grid StaticConformal
 
 private:
 
diff --git a/param.ccl b/param.ccl
index 7db98ea..e321d87 100644
--- a/param.ccl
+++ b/param.ccl
@@ -20,17 +20,37 @@ EXTENDS KEYWORD initial_data
 
 private:
 
-REAL amp "Brill wave amplitude"
+REAL amp "Brill wave amplitude for extrinsic curveture"
 {
   *:* :: "No restriction"
 } 0.1
 
-REAL eta0 "Brill wave center (in eta coords)"
+REAL amp_me "Brill wave amplitude for metric"
 {
   *:* :: "No restriction"
 } 0.0
 
-REAL sigma "Brill wave width (in eta)"
+REAL xi "for recoil"
+{
+  *:* :: "No restriction"
+} 0.0
+
+REAL eta0 "Brill wave center (in eta coords) for extrisic curvature"
+{
+  *:* :: "No restriction"
+} 0.0
+
+REAL eta0_me "Brill wave center (in eta coords) for metric"
+{
+  *:* :: "No restriction"
+} 0.0
+
+REAL sigma "Brill wave width (in eta) for extrinsic curveture"
+{
+  *:* :: "No restriction"
+} 1.0
+
+REAL sigma_me "Brill wave width (in eta) for metric"
 {
   *:* :: "No restriction"
 } 1.0
@@ -46,7 +66,12 @@ REAL etamax "Eta value for outer edge of grid"
 } 5.0
 
 
-INT n "sin**n theta in brill wave"
+INT n "sin**n theta in brill wave for extrinsic curveture"
+{
+  *:* :: "No restriction"
+} 3
+
+INT n_me "sin**n theta in brill wave for metric"
 {
   *:* :: "No restriction"
 } 3
diff --git a/schedule.ccl b/schedule.ccl
index 79254f1..c9c43df 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -10,7 +10,7 @@ if (CCTK_Equals(initial_data,"axioddbh"))
     LANG: C
   } "Check Parameters"
 
-  schedule IDAxiOddBrillBH at CCTK_INITIAL
+  schedule IDAxiOddBrillBH IN ADMBase_InitialData
   {
     LANG: Fortran
   } "Construct IDAxiOddBrillBH"  
diff --git a/src/IDAxiOddBrillBH.F b/src/IDAxiOddBrillBH.F
index 838ac2f..781fd36 100644
--- a/src/IDAxiOddBrillBH.F
+++ b/src/IDAxiOddBrillBH.F
@@ -36,7 +36,8 @@ c@@*/
 
       real*8 deta,dq
       real*8, allocatable :: ac(:,:),ae(:,:),aw(:,:),an(:,:),as(:,:),
-     $     rhs(:,:),ksq(:,:), psi2dv(:,:),dpsi2dv(:,:),ddpsi2dv(:,:),
+     $     rhs(:,:),qfetaeta(:,:),qfqq(:,:),ksq(:,:), 
+     $     psi2dv(:,:),dpsi2dv(:,:),ddpsi2dv(:,:),
      $     detapsisph(:,:),dqpsisph(:,:),detaetapsisph(:,:),
      $     detaqpsisph(:,:),dqqpsisph(:,:),
      $     etagrd(:),qgrd(:)
@@ -85,11 +86,17 @@ c     --------------
          write(*,*) '"Solve Bradt-Seidel Initial Data Sets"'
       endif
 
-      write(*,*) 'Brill Wave: amplitude :',amp
-      write(*,*) 'Brill Wave: center :',eta0
-      write(*,*) 'Brill Wave: sigma :',sigma
+      write(*,*) 'Brill Wave: amplitude for K :',amp
+      write(*,*) 'Brill Wave: center for K:',eta0
+      write(*,*) 'Brill Wave: sigma for K:',sigma
       write(*,*) 'Bowen & York Momenta:  :',byJ
-      write(*,*) 'Brill Wave: sin^n theta :',n
+      write(*,*) 'Brill Wave: sin^n theta for K:',n
+
+      write(*,*) 'Brill Wave: amplitude for g :',amp_me
+      write(*,*) 'Brill Wave: center for g:',eta0_me
+      write(*,*) 'Brill Wave: sigma for g:',sigma_me
+      write(*,*) 'Brill Wave: sin^n theta for g:',n_me
+
 
 c Check if we should create and store conformal factor stuff
       
@@ -124,7 +131,8 @@ c     Add 2 zones for eta coordinate and 4 for theta
       nq = nq+2
 c     
       allocate(ac(ne,nq),ae(ne,nq),aw(ne,nq),an(ne,nq),as(ne,nq),
-     $     rhs(ne,nq),ksq(ne,nq),psi2dv(ne,nq),dpsi2dv(ne,nq),
+     $     rhs(ne,nq),qfetaeta(ne,nq),qfqq(ne,nq),ksq(ne,nq),
+     $     psi2dv(ne,nq),dpsi2dv(ne,nq),
      $     ddpsi2dv(ne,nq),detapsisph(ne,nq),dqpsisph(ne,nq),
      $     detaetapsisph(ne,nq),detaqpsisph(ne,nq),dqqpsisph(ne,nq),
      $     etagrd(ne),qgrd(nq))
@@ -141,7 +149,16 @@ c
       do i=1,ne
          etagrd(i) = (i-2)*deta
       enddo
-      
+c     
+c Initialize q-function and its derivatives: should be generalized
+c     
+       do j = 1,nq
+          do i = 1,ne
+#include "qfunc_even.x"
+          enddo
+       enddo
+
+
 c     
 c     Specify the initial guess for the conformal factor:
 c     
@@ -190,7 +207,9 @@ c
      $              (dpsi2dv(i+1,j)-2.*dpsi2dv(i,j)+dpsi2dv(i-1,j))/
      $              deta**2+(dpsi2dv(i,j+1)-2.*dpsi2dv(i,j)+
      $              dpsi2dv(i,j-1))/dq**2+0.5*(dpsi2dv(i,j+1)-
-     $              dpsi2dv(i,j-1))/(dq*tan(qgrd(j)))-0.25*dpsi2dv(i,j)+
+     $              dpsi2dv(i,j-1))/(dq*tan(qgrd(j)))
+     $              +0.25*dpsi2dv(i,j)*
+     $              (qfetaeta(i,j)+qfqq(i,j)-1.0d0)+
      $              0.125*ksq(i,j)/(psi2dv(i,j)+dpsi2dv(i,j))**7
             enddo
          enddo
@@ -220,8 +239,9 @@ c     Compute stencil coefficients:
 c     
          do j = 2,nq-1
             do i = 2,ne-1
-               ac(i,j) = -2./deta**2-2./dq**2-sin(qgrd(j))**2-0.25-
-     $              0.875*ksq(i,j)/(psi2dv(i,j)+ 
+               ac(i,j) = -2./deta**2-2./dq**2-sin(qgrd(j))**2+0.25
+     $              *(qfetaeta(i,j)+qfqq(i,j) - 1.0d0)
+     $              -0.875*ksq(i,j)/(psi2dv(i,j)+ 
      $              dpsi2dv(i,j))**8
 c     
                ae(i,j) = 1./deta**2
@@ -549,7 +569,8 @@ c      print *,'Angular momentum: ',Jmom
          alp = (2.*r - adm)/(2.*r+adm)
       endif
       
-      deallocate(ac,ae,aw,an,as,rhs,ksq,psi2dv,dpsi2dv,ddpsi2dv,
+      deallocate(ac,ae,aw,an,as,rhs,qfetaeta,qfqq,
+     $     ksq,psi2dv,dpsi2dv,ddpsi2dv,
      $     detapsisph,dqpsisph,detaetapsisph,detaqpsisph,dqqpsisph,
      $     etagrd,qgrd)
       
diff --git a/src/gij.x b/src/gij.x
index 693b57d..1f57c56 100644
--- a/src/gij.x
+++ b/src/gij.x
@@ -1,24 +1,29 @@
-        gxx(i,j,k) = 1.00000000000000d0
-c        dxgxx(i,j,k) = 0
-c        dygxx(i,j,k) = 0
-c        dzgxx(i,j,k) = 0
-        gxy(i,j,k) = 0
-c        dxgxy(i,j,k) = 0
-c        dygxy(i,j,k) = 0
-c        dzgxy(i,j,k) = 0
+        o1 = 2.d0*phi(i,j,k)
+        o2 = cos(o1)
+        o3 = -xi
+        o4 = cos(q(i,j,k))
+        o5 = o4*xi
+        o6 = 1.d0 + o3 + o5
+        o7 = -eta0_me
+        o8 = eta(i,j,k) + o7
+        o9 = o8**2
+        o10 = sigma_me**2
+        o11 = 1/o10
+        o12 = -(o11*o9)
+        o13 = exp(o12)
+        o14 = eta(i,j,k) + eta0_me
+        o15 = o14**2
+        o16 = -(o11*o15)
+        o17 = exp(o16)
+        o18 = o13 + o17
+        o19 = sin(q(i,j,k))
+        o20 = o19**n_me
+        o21 = 2.d0*amp_me*o18*o20*o6
+        o22 = exp(o21)
+        o23 = -1.d0 + o22
+        gxx(i,j,k) = 5.d-1*(1.d0 + o22 + o2*o23)
+        gxy(i,j,k) = o23*cos(phi(i,j,k))*sin(phi(i,j,k))
         gxz(i,j,k) = 0
-c        dxgxz(i,j,k) = 0
-c        dygxz(i,j,k) = 0
-c        dzgxz(i,j,k) = 0
-        gyy(i,j,k) = 1.00000000000000d0
-c        dxgyy(i,j,k) = 0
-c        dygyy(i,j,k) = 0
-c        dzgyy(i,j,k) = 0
+        gyy(i,j,k) = 5.d-1*(1.d0 + o22 - o2*o23)
         gyz(i,j,k) = 0
-c        dxgyz(i,j,k) = 0
-c        dygyz(i,j,k) = 0
-c        dzgyz(i,j,k) = 0
-        gzz(i,j,k) = 1.00000000000000d0
-c        dxgzz(i,j,k) = 0
-c        dygzz(i,j,k) = 0
-c        dzgzz(i,j,k) = 0
+        gzz(i,j,k) = o22