1 files changed, 116 insertions, 0 deletions
diff --git a/src/Bin_BH.cc b/src/Bin_BH.cc
new file mode 100644
index 0000000..693cdbc
--- /dev/null
+++ b/src/Bin_BH.cc
@@ -0,0 +1,116 @@
+#include <cassert>
+#include <cstdio>
+#include <vector>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <bin_bh.h>
+
+using namespace std;
+
+
+
+extern "C"
+void ID_Bin_BH_initialise (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INFO ("Setting up LORENE Bin_BH initial data");
+  
+  
+  
+  CCTK_INFO ("Setting up coordinates");
+  
+  int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+  
+  vector<double> xx(npoints), yy(npoints), zz(npoints);
+  
+#pragma omp parallel for
+  for (int i=0; i<npoints; ++i) {
+    xx[i] = x[i];
+    yy[i] = y[i];
+    zz[i] = z[i];
+  }
+  
+  
+  
+  CCTK_VInfo (CCTK_THORNSTRING, "Reading from file \"%s\"", filename);
+  
+  Bin_BH bin_bh (npoints, &xx[0], &yy[0], &zz[0], 1, filename);
+  
+  CCTK_VInfo (CCTK_THORNSTRING, "Omega [1/a]: %g", bin_bh.omega);
+  CCTK_VInfo (CCTK_THORNSTRING, "dist [a]:    %g", bin_bh.dist);
+  CCTK_VInfo (CCTK_THORNSTRING, "radius1 [a]: %g", 1.0);
+  CCTK_VInfo (CCTK_THORNSTRING, "radius2 [a]: %g", bin_bh.radius2);
+  assert (bin_bh.np == npoints);
+  
+  
+  
+  CCTK_INFO ("Filling in Cactus grid points");
+  
+#pragma omp parallel for
+  for (int i=0; i<npoints; ++i) {
+    
+    alp[i] = bin_bh.nnn[i];
+    
+    betax[i] = bin_bh.beta_x[i];
+    betay[i] = bin_bh.beta_y[i];
+    betaz[i] = bin_bh.beta_z[i];
+    
+    CCTK_REAL g[3][3];
+    g[0][0] = bin_bh.g_xx[i];
+    g[0][1] = bin_bh.g_xy[i];
+    g[0][2] = bin_bh.g_xz[i];
+    g[1][1] = bin_bh.g_yy[i];
+    g[1][2] = bin_bh.g_yz[i];
+    g[2][2] = bin_bh.g_zz[i];
+    g[1][0] = g[0][1];
+    g[2][0] = g[0][2];
+    g[2][1] = g[1][2];
+    
+    CCTK_REAL ku[3][3];
+    ku[0][0] = bin_bh.k_xx[i];
+    ku[0][1] = bin_bh.k_xy[i];
+    ku[0][2] = bin_bh.k_xz[i];
+    ku[1][1] = bin_bh.k_yy[i];
+    ku[1][2] = bin_bh.k_yz[i];
+    ku[2][2] = bin_bh.k_zz[i];
+    ku[1][0] = ku[0][1];
+    ku[2][0] = ku[0][2];
+    ku[2][1] = ku[1][2];
+    
+    CCTK_REAL k[3][3];
+    for (int a=0; a<3; ++a) {
+      for (int b=0; b<3; ++b) {
+        k[a][b] = 0.0;
+        for (int c=0; c<3; ++c) {
+          for (int d=0; d<3; ++d) {
+            k[a][b] += g[a][c] * g[b][d] * ku[c][d];
+          }
+        }
+      }
+    }
+    
+    gxx[i] = g[0][0];
+    gxy[i] = g[0][1];
+    gxz[i] = g[0][2];
+    gyy[i] = g[1][1];
+    gyz[i] = g[1][2];
+    gzz[i] = g[2][2];
+    
+    kxx[i] = k[0][0];
+    kxy[i] = k[0][1];
+    kxz[i] = k[0][2];
+    kyy[i] = k[1][1];
+    kyz[i] = k[1][2];
+    kzz[i] = k[2][2];
+    
+  } // for i
+  
+  
+  
+  CCTK_INFO ("Done.");
+}