Import initial version of Marcus Ansorg's code, converted into a

Cactus thorn.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@2 b2a53a04-0f4f-0410-87ed-f9f25ced00cf

1 files changed, 744 insertions, 0 deletions

diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
new file mode 100644
index 0000000..568ac13
--- /dev/null
+++ b/src/FuncAndJacobian.c
@@ -0,0 +1,744 @@
+// TwoPunctures:  File  "FuncAndJacobian.c"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <time.h>
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+#define FAC sin(al)*sin(be)*sin(al)*sin(be)*sin(al)*sin(be)
+//#define FAC sin(al)*sin(be)*sin(al)*sin(be)
+//#define FAC 1
+
+// -------------------------------------------------------------------------------
+int
+Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
+{
+  int i1 = i, j1 = j, k1 = k;
+
+  if (i1 < 0)
+    i1 = -(i1 + 1);
+  if (i1 >= n1)
+    i1 = 2 * n1 - (i1 + 1);
+
+  if (j1 < 0)
+    j1 = -(j1 + 1);
+  if (j1 >= n2)
+    j1 = 2 * n2 - (j1 + 1);
+
+  if (k1 < 0)
+    k1 = k1 + n3;
+  if (k1 >= n3)
+    k1 = k1 - n3;
+
+  return ivar + nvar * (i1 + n1 * (j1 + n2 * k1));
+}
+
+// -------------------------------------------------------------------------------
+void
+allocate_derivs (derivs * v, int n)
+{
+  int m = n - 1;
+  (*v).d0 = dvector (0, m);
+  (*v).d1 = dvector (0, m);
+  (*v).d2 = dvector (0, m);
+  (*v).d3 = dvector (0, m);
+  (*v).d11 = dvector (0, m);
+  (*v).d12 = dvector (0, m);
+  (*v).d13 = dvector (0, m);
+  (*v).d22 = dvector (0, m);
+  (*v).d23 = dvector (0, m);
+  (*v).d33 = dvector (0, m);
+}
+
+// -------------------------------------------------------------------------------
+void
+free_derivs (derivs * v, int n)
+{
+  int m = n - 1;
+  free_dvector ((*v).d0, 0, m);
+  free_dvector ((*v).d1, 0, m);
+  free_dvector ((*v).d2, 0, m);
+  free_dvector ((*v).d3, 0, m);
+  free_dvector ((*v).d11, 0, m);
+  free_dvector ((*v).d12, 0, m);
+  free_dvector ((*v).d13, 0, m);
+  free_dvector ((*v).d22, 0, m);
+  free_dvector ((*v).d23, 0, m);
+  free_dvector ((*v).d33, 0, m);
+}
+
+// -------------------------------------------------------------------------------
+void
+Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v)
+{
+  int i, j, k, ivar, N, *indx;
+  double *p, *dp, *d2p, *q, *dq, *r, *dr;
+
+  N = maximum3 (n1, n2, n3);
+  p = dvector (0, N);
+  dp = dvector (0, N);
+  d2p = dvector (0, N);
+  q = dvector (0, N);
+  dq = dvector (0, N);
+  r = dvector (0, N);
+  dr = dvector (0, N);
+  indx = ivector (0, N);
+
+  for (ivar = 0; ivar < nvar; ivar++)
+  {
+    for (k = 0; k < n3; k++)
+    {				// Calculation of Derivatives w.r.t. A-Dir. 
+      for (j = 0; j < n2; j++)
+      {				// (Chebyshev_Zeros)
+	for (i = 0; i < n1; i++)
+	{
+	  indx[i] = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  p[i] = v.d0[indx[i]];
+	}
+	chebft_Zeros (p, n1, 0);
+	chder (p, dp, n1);
+	chder (dp, d2p, n1);
+	chebft_Zeros (dp, n1, 1);
+	chebft_Zeros (d2p, n1, 1);
+	for (i = 0; i < n1; i++)
+	{
+	  v.d1[indx[i]] = dp[i];
+	  v.d11[indx[i]] = d2p[i];
+	}
+      }
+    }
+    for (k = 0; k < n3; k++)
+    {				// Calculation of Derivatives w.r.t. A-Dir. 
+      for (i = 0; i < n1; i++)
+      {				// (Chebyshev_Zeros)
+	for (j = 0; j < n2; j++)
+	{
+	  indx[j] = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  p[j] = v.d0[indx[j]];
+	  q[j] = v.d1[indx[j]];
+	}
+	chebft_Zeros (p, n2, 0);
+	chebft_Zeros (q, n2, 0);
+	chder (p, dp, n2);
+	chder (dp, d2p, n2);
+	chder (q, dq, n2);
+	chebft_Zeros (dp, n2, 1);
+	chebft_Zeros (d2p, n2, 1);
+	chebft_Zeros (dq, n2, 1);
+	for (j = 0; j < n2; j++)
+	{
+	  v.d2[indx[j]] = dp[j];
+	  v.d22[indx[j]] = d2p[j];
+	  v.d12[indx[j]] = dq[j];
+	}
+      }
+    }
+    for (i = 0; i < n1; i++)
+    {				// Calculation of Derivatives w.r.t. phi-Dir. (Fourier)
+      for (j = 0; j < n2; j++)
+      {
+	for (k = 0; k < n3; k++)
+	{
+	  indx[k] = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  p[k] = v.d0[indx[k]];
+	  q[k] = v.d1[indx[k]];
+	  r[k] = v.d2[indx[k]];
+	}
+	fourft (p, n3, 0);
+	fourder (p, dp, n3);
+	fourder2 (p, d2p, n3);
+	fourft (dp, n3, 1);
+	fourft (d2p, n3, 1);
+	fourft (q, n3, 0);
+	fourder (q, dq, n3);
+	fourft (dq, n3, 1);
+	fourft (r, n3, 0);
+	fourder (r, dr, n3);
+	fourft (dr, n3, 1);
+	for (k = 0; k < n3; k++)
+	{
+	  v.d3[indx[k]] = dp[k];
+	  v.d33[indx[k]] = d2p[k];
+	  v.d13[indx[k]] = dq[k];
+	  v.d23[indx[k]] = dr[k];
+	}
+      }
+    }
+  }
+  free_dvector (p, 0, N);
+  free_dvector (dp, 0, N);
+  free_dvector (d2p, 0, N);
+  free_dvector (q, 0, N);
+  free_dvector (dq, 0, N);
+  free_dvector (r, 0, N);
+  free_dvector (dr, 0, N);
+  free_ivector (indx, 0, N);
+}
+
+// -------------------------------------------------------------------------------
+void
+F_of_v (int nvar, int n1, int n2, int n3, derivs v, double *F,
+	derivs u)
+{				//      Calculates the left hand sides of the non-linear equations F_m(v_n)=0
+  //      and the function u (u.d0[]) as well as its derivatives
+  //      (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])
+  //      at interior points and at the boundaries "+/-"
+  int i, j, k, ivar, indx;
+  double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
+  derivs U;
+
+  values = dvector (0, nvar - 1);
+  allocate_derivs (&U, nvar);
+
+  Derivatives_AB3 (nvar, n1, n2, n3, v);
+
+  for (i = 0; i < n1; i++)
+  {
+    for (j = 0; j < n2; j++)
+    {
+      for (k = 0; k < n3; k++)
+      {
+
+	al = Pih * (2 * i + 1) / n1;
+	A = -cos (al);
+	be = Pih * (2 * j + 1) / n2;
+	B = -cos (be);
+	phi = 2. * Pi * k / n3;
+
+	Am1 = A - 1;
+	for (ivar = 0; ivar < nvar; ivar++)
+	{
+	  indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  U.d0[ivar] = Am1 * v.d0[indx];	// U
+	  U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx];	// U_A
+	  U.d2[ivar] = Am1 * v.d2[indx];	// U_B
+	  U.d3[ivar] = Am1 * v.d3[indx];	// U_3
+	  U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx];	// U_AA
+	  U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx];	// U_AB
+	  U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx];	// U_AB
+	  U.d22[ivar] = Am1 * v.d22[indx];	// U_BB
+	  U.d23[ivar] = Am1 * v.d23[indx];	// U_B3
+	  U.d33[ivar] = Am1 * v.d33[indx];	// U_33
+	}
+	// Calculation of (X,R) and
+	// (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)
+	AB_To_XR (nvar, A, B, &X, &R, U);
+	// Calculation of (x,r) and
+	// (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)
+	C_To_c (nvar, X, R, &x, &r, U);
+	// Calculation of (y,z) and
+	// (U, U_x, U_y, U_z, U_xx, U_xy, U_xz, U_yy, U_yz, U_zz)
+	rx3_To_xyz (nvar, x, r, phi, &y, &z, U);
+	NonLinEquations (A, B, X, R, x, r, phi, y, z, U, values);
+	for (ivar = 0; ivar < nvar; ivar++)
+	{
+	  indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  F[indx] = values[ivar] * FAC;
+	  u.d0[indx] = U.d0[ivar];	//  U
+	  u.d1[indx] = U.d1[ivar];	//      U_x
+	  u.d2[indx] = U.d2[ivar];	//      U_y
+	  u.d3[indx] = U.d3[ivar];	//      U_z
+	  u.d11[indx] = U.d11[ivar];	//      U_xx
+	  u.d12[indx] = U.d12[ivar];	//      U_xy
+	  u.d13[indx] = U.d13[ivar];	//      U_xz
+	  u.d22[indx] = U.d22[ivar];	//      U_yy
+	  u.d23[indx] = U.d23[ivar];	//      U_yz
+	  u.d33[indx] = U.d33[ivar];	//      U_zz
+	}
+      }
+    }
+  }
+  free_dvector (values, 0, nvar - 1);
+  free_derivs (&U, nvar);
+}
+
+// -------------------------------------------------------------------------------
+void
+J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
+	    double *Jdv, derivs u)
+{				//      Calculates the left hand sides of the non-linear equations F_m(v_n)=0
+  //      and the function u (u.d0[]) as well as its derivatives
+  //      (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])
+  //      at interior points and at the boundaries "+/-"
+  DECLARE_CCTK_PARAMETERS;
+  int i, j, k, ivar, indx;
+  double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
+  derivs dU, U;
+
+
+  values = dvector (0, nvar - 1);
+  allocate_derivs (&dU, nvar);
+  allocate_derivs (&U, nvar);
+
+  Derivatives_AB3 (nvar, n1, n2, n3, dv);
+
+  for (i = 0; i < n1; i++)
+  {
+    for (j = 0; j < n2; j++)
+    {
+      for (k = 0; k < n3; k++)
+      {
+
+	al = Pih * (2 * i + 1) / n1;
+	A = -cos (al);
+	be = Pih * (2 * j + 1) / n2;
+	B = -cos (be);
+	phi = 2. * Pi * k / n3;
+
+	Am1 = A - 1;
+	for (ivar = 0; ivar < nvar; ivar++)
+	{
+	  indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  dU.d0[ivar] = Am1 * dv.d0[indx];	// dU
+	  dU.d1[ivar] = dv.d0[indx] + Am1 * dv.d1[indx];	// dU_A
+	  dU.d2[ivar] = Am1 * dv.d2[indx];	// dU_B
+	  dU.d3[ivar] = Am1 * dv.d3[indx];	// dU_3
+	  dU.d11[ivar] = 2 * dv.d1[indx] + Am1 * dv.d11[indx];	// dU_AA
+	  dU.d12[ivar] = dv.d2[indx] + Am1 * dv.d12[indx];	// dU_AB
+	  dU.d13[ivar] = dv.d3[indx] + Am1 * dv.d13[indx];	// dU_AB
+	  dU.d22[ivar] = Am1 * dv.d22[indx];	// dU_BB
+	  dU.d23[ivar] = Am1 * dv.d23[indx];	// dU_B3
+	  dU.d33[ivar] = Am1 * dv.d33[indx];	// dU_33
+	  U.d0[ivar] = u.d0[indx];	// U   
+	  U.d1[ivar] = u.d1[indx];	// U_x
+	  U.d2[ivar] = u.d2[indx];	// U_y
+	  U.d3[ivar] = u.d3[indx];	// U_z
+	  U.d11[ivar] = u.d11[indx];	// U_xx
+	  U.d12[ivar] = u.d12[indx];	// U_xy
+	  U.d13[ivar] = u.d13[indx];	// U_xz
+	  U.d22[ivar] = u.d22[indx];	// U_yy
+	  U.d23[ivar] = u.d23[indx];	// U_yz
+	  U.d33[ivar] = u.d33[indx];	// U_zz
+	}
+	// Calculation of (X,R) and
+	// (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)
+	AB_To_XR (nvar, A, B, &X, &R, dU);
+	// Calculation of (x,r) and
+	// (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)
+	C_To_c (nvar, X, R, &x, &r, dU);
+	// Calculation of (y,z) and
+	// (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)
+	rx3_To_xyz (nvar, x, r, phi, &y, &z, dU);
+	LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values);
+	for (ivar = 0; ivar < nvar; ivar++)
+	{
+	  indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  Jdv[indx] = values[ivar] * FAC;
+	}
+      }
+    }
+  }
+  free_dvector (values, 0, nvar - 1);
+  free_derivs (&dU, nvar);
+  free_derivs (&U, nvar);
+}
+
+// -------------------------------------------------------------------------------
+void
+JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2,
+	      int n3, derivs dv, derivs u, double *values)
+{				// Calculates rows of the vector 'J(FD)*dv'.
+  // First row to be calculated: row = Index(0,      i, j, k; nvar, n1, n2, n3)
+  // Last  row to be calculated: row = Index(nvar-1, i, j, k; nvar, n1, n2, n3)
+  // These rows are stored in the vector JFDdv[0] ... JFDdv[nvar-1].
+  DECLARE_CCTK_PARAMETERS;
+  int ivar, indx;
+  double al, be, A, B, X, R, x, r, phi, y, z, Am1;
+  double sin_al, sin_al_i1, sin_al_i2, sin_al_i3, cos_al;
+  double sin_be, sin_be_i1, sin_be_i2, sin_be_i3, cos_be;
+  double dV0, dV1, dV2, dV3, dV11, dV12, dV13, dV22, dV23, dV33,
+    ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
+  derivs dU, U;
+
+  allocate_derivs (&dU, nvar);
+  allocate_derivs (&U, nvar);
+
+  if (k < 0)
+    k = k + n3;
+  if (k >= n3)
+    k = k - n3;
+
+  ha = Pi / n1;			// ha: Stepsize with respect to (al)
+  al = ha * (i + 0.5);
+  A = -cos (al);
+  ga = 1 / ha;
+  ga2 = ga * ga;
+
+  hb = Pi / n2;			// hb: Stepsize with respect to (be)
+  be = hb * (j + 0.5);
+  B = -cos (be);
+  gb = 1 / hb;
+  gb2 = gb * gb;
+  gagb = ga * gb;
+
+  hp = 2 * Pi / n3;		// hp: Stepsize with respect to (phi)
+  phi = hp * j;
+  gp = 1 / hp;
+  gp2 = gp * gp;
+  gagp = ga * gp;
+  gbgp = gb * gp;
+
+
+  sin_al = sin (al);
+  sin_be = sin (be);
+  sin_al_i1 = 1 / sin_al;
+  sin_be_i1 = 1 / sin_be;
+  sin_al_i2 = sin_al_i1 * sin_al_i1;
+  sin_be_i2 = sin_be_i1 * sin_be_i1;
+  sin_al_i3 = sin_al_i1 * sin_al_i2;
+  sin_be_i3 = sin_be_i1 * sin_be_i2;
+  cos_al = -A;
+  cos_be = -B;
+
+  Am1 = A - 1;
+  for (ivar = 0; ivar < nvar; ivar++)
+  {
+    int iccc = Index (ivar, i, j, k, nvar, n1, n2, n3),
+      ipcc = Index (ivar, i + 1, j, k, nvar, n1, n2, n3),
+      imcc = Index (ivar, i - 1, j, k, nvar, n1, n2, n3),
+      icpc = Index (ivar, i, j + 1, k, nvar, n1, n2, n3),
+      icmc = Index (ivar, i, j - 1, k, nvar, n1, n2, n3),
+      iccp = Index (ivar, i, j, k + 1, nvar, n1, n2, n3),
+      iccm = Index (ivar, i, j, k - 1, nvar, n1, n2, n3),
+      icpp = Index (ivar, i, j + 1, k + 1, nvar, n1, n2, n3),
+      icmp = Index (ivar, i, j - 1, k + 1, nvar, n1, n2, n3),
+      icpm = Index (ivar, i, j + 1, k - 1, nvar, n1, n2, n3),
+      icmm = Index (ivar, i, j - 1, k - 1, nvar, n1, n2, n3),
+      ipcp = Index (ivar, i + 1, j, k + 1, nvar, n1, n2, n3),
+      imcp = Index (ivar, i - 1, j, k + 1, nvar, n1, n2, n3),
+      ipcm = Index (ivar, i + 1, j, k - 1, nvar, n1, n2, n3),
+      imcm = Index (ivar, i - 1, j, k - 1, nvar, n1, n2, n3),
+      ippc = Index (ivar, i + 1, j + 1, k, nvar, n1, n2, n3),
+      impc = Index (ivar, i - 1, j + 1, k, nvar, n1, n2, n3),
+      ipmc = Index (ivar, i + 1, j - 1, k, nvar, n1, n2, n3),
+      immc = Index (ivar, i - 1, j - 1, k, nvar, n1, n2, n3);
+    // Derivatives of (dv) w.r.t. (al,be,phi):
+    dV0 = dv.d0[iccc];
+    dV1 = 0.5 * ga * (dv.d0[ipcc] - dv.d0[imcc]);
+    dV2 = 0.5 * gb * (dv.d0[icpc] - dv.d0[icmc]);
+    dV3 = 0.5 * gp * (dv.d0[iccp] - dv.d0[iccm]);
+    dV11 = ga2 * (dv.d0[ipcc] + dv.d0[imcc] - 2 * dv.d0[iccc]);
+    dV22 = gb2 * (dv.d0[icpc] + dv.d0[icmc] - 2 * dv.d0[iccc]);
+    dV33 = gp2 * (dv.d0[iccp] + dv.d0[iccm] - 2 * dv.d0[iccc]);
+    dV12 =
+      0.25 * gagb * (dv.d0[ippc] - dv.d0[ipmc] + dv.d0[immc] - dv.d0[impc]);
+    dV13 =
+      0.25 * gagp * (dv.d0[ipcp] - dv.d0[imcp] + dv.d0[imcm] - dv.d0[ipcm]);
+    dV23 =
+      0.25 * gbgp * (dv.d0[icpp] - dv.d0[icpm] + dv.d0[icmm] - dv.d0[icmp]);
+    // Derivatives of (dv) w.r.t. (A,B,phi):
+    dV11 = sin_al_i3 * (sin_al * dV11 - cos_al * dV1);
+    dV12 = sin_al_i1 * sin_be_i1 * dV12;
+    dV13 = sin_al_i1 * dV13;
+    dV22 = sin_be_i3 * (sin_be * dV22 - cos_be * dV2);
+    dV23 = sin_be_i1 * dV23;
+    dV1 = sin_al_i1 * dV1;
+    dV2 = sin_be_i1 * dV2;
+    // Derivatives of (dU) w.r.t. (A,B,phi):
+    dU.d0[ivar] = Am1 * dV0;
+    dU.d1[ivar] = dV0 + Am1 * dV1;
+    dU.d2[ivar] = Am1 * dV2;
+    dU.d3[ivar] = Am1 * dV3;
+    dU.d11[ivar] = 2 * dV1 + Am1 * dV11;
+    dU.d12[ivar] = dV2 + Am1 * dV12;
+    dU.d13[ivar] = dV3 + Am1 * dV13;
+    dU.d22[ivar] = Am1 * dV22;
+    dU.d23[ivar] = Am1 * dV23;
+    dU.d33[ivar] = Am1 * dV33;
+
+    indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+    U.d0[ivar] = u.d0[indx];	// U   
+    U.d1[ivar] = u.d1[indx];	// U_x
+    U.d2[ivar] = u.d2[indx];	// U_y
+    U.d3[ivar] = u.d3[indx];	// U_z
+    U.d11[ivar] = u.d11[indx];	// U_xx
+    U.d12[ivar] = u.d12[indx];	// U_xy
+    U.d13[ivar] = u.d13[indx];	// U_xz
+    U.d22[ivar] = u.d22[indx];	// U_yy
+    U.d23[ivar] = u.d23[indx];	// U_yz
+    U.d33[ivar] = u.d33[indx];	// U_zz
+  }
+  // Calculation of (X,R) and
+  // (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)
+  AB_To_XR (nvar, A, B, &X, &R, dU);
+  // Calculation of (x,r) and
+  // (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)
+  C_To_c (nvar, X, R, &x, &r, dU);
+  // Calculation of (y,z) and
+  // (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)
+  rx3_To_xyz (nvar, x, r, phi, &y, &z, dU);
+  LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values);
+  for (ivar = 0; ivar < nvar; ivar++)
+    values[ivar] *= FAC;
+
+  free_derivs (&dU, nvar);
+  free_derivs (&U, nvar);
+}
+
+// -------------------------------------------------------------------------------
+void
+SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u,
+	       int *ncols, int **cols, double **Matrix)
+{
+  DECLARE_CCTK_PARAMETERS;
+  int column, row, mcol;
+  int i, i1, i_0, i_1, j, j1, j_0, j_1, k, k1, k_0, k_1, N1, N2, N3,
+    ivar, ivar1, ntotal = nvar * n1 * n2 * n3;
+  double *values;
+  derivs dv;
+
+  values = dvector (0, nvar - 1);
+  allocate_derivs (&dv, ntotal);
+
+  N1 = n1 - 1;
+  N2 = n2 - 1;
+  N3 = n3 - 1;
+
+  for (i = 0; i < n1; i++)
+  {
+    for (j = 0; j < n2; j++)
+    {
+      for (k = 0; k < n3; k++)
+      {
+	for (ivar = 0; ivar < nvar; ivar++)
+	{
+	  row = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  ncols[row] = 0;
+	  dv.d0[row] = 0;
+	}
+      }
+    }
+  }
+  for (i = 0; i < n1; i++)
+  {
+    for (j = 0; j < n2; j++)
+    {
+      for (k = 0; k < n3; k++)
+      {
+	for (ivar = 0; ivar < nvar; ivar++)
+	{
+	  column = Index (ivar, i, j, k, nvar, n1, n2, n3);
+	  dv.d0[column] = 1;
+
+	  i_0 = maximum2 (0, i - 1);
+	  i_1 = minimum2 (N1, i + 1);
+	  j_0 = maximum2 (0, j - 1);
+	  j_1 = minimum2 (N2, j + 1);
+	  k_0 = k - 1;
+	  k_1 = k + 1;
+/*					i_0 = 0;
+					i_1 = N1;
+					j_0 = 0;
+					j_1 = N2;
+					k_0 = 0;
+					k_1 = N3;*/
+
+	  for (i1 = i_0; i1 <= i_1; i1++)
+	  {
+	    for (j1 = j_0; j1 <= j_1; j1++)
+	    {
+	      for (k1 = k_0; k1 <= k_1; k1++)
+	      {
+		JFD_times_dv (i1, j1, k1, nvar, n1, n2, n3,
+			      dv, u, values);
+		for (ivar1 = 0; ivar1 < nvar; ivar1++)
+		{
+		  if (values[ivar1] != 0)
+		  {
+		    row = Index (ivar1, i1, j1, k1, nvar, n1, n2, n3);
+		    mcol = ncols[row];
+		    cols[row][mcol] = column;
+		    Matrix[row][mcol] = values[ivar1];
+		    ncols[row] += 1;
+		  }
+		}
+	      }
+	    }
+	  }
+
+	  dv.d0[column] = 0;
+	}
+      }
+    }
+  }
+  free_derivs (&dv, ntotal);
+  free_dvector (values, 0, nvar - 1);
+}
+
+// -------------------------------------------------------------------------------
+// Calculates the value of v at an arbitrary position (A,B,phi)
+double
+PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
+			  int n1, int n2, int n3)
+{
+  int i, j, k, N;
+  double *p, *values1, **values2, result;
+
+  N = maximum3 (n1, n2, n3);
+  p = dvector (0, N);
+  values1 = dvector (0, N);
+  values2 = dmatrix (0, N, 0, N);
+
+  for (k = 0; k < n3; k++)
+  {
+    for (j = 0; j < n2; j++)
+    {
+      for (i = 0; i < n1; i++)
+	p[i] = v[i + n1 * (j + n2 * k)];
+      chebft_Zeros (p, n1, 0);
+      values2[j][k] = chebev (-1, 1, p, n1, A);
+    }
+  }
+
+  for (k = 0; k < n3; k++)
+  {
+    for (j = 0; j < n2; j++)
+      p[j] = values2[j][k];
+    chebft_Zeros (p, n2, 0);
+    values1[k] = chebev (-1, 1, p, n2, B);
+  }
+
+  fourft (values1, n3, 0);
+  result = fourev (values1, n3, phi);
+
+  free_dvector (p, 0, N);
+  free_dvector (values1, 0, N);
+  free_dmatrix (values2, 0, N, 0, N);
+
+  return result;
+}
+
+// -------------------------------------------------------------------------------
+void
+calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2,
+		  int n3, derivs v, derivs vv)
+{
+  double al = Pih * (2 * i + 1) / n1, be = Pih * (2 * j + 1) / n2,
+    sin_al = sin (al), sin2_al = sin_al * sin_al, cos_al = cos (al),
+    sin_be = sin (be), sin2_be = sin_be * sin_be, cos_be = cos (be);
+
+  vv.d0[0] = v.d0[Index (ivar, i, j, k, nvar, n1, n2, n3)];
+  vv.d1[0] = v.d1[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al;
+  vv.d2[0] = v.d2[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_be;
+  vv.d3[0] = v.d3[Index (ivar, i, j, k, nvar, n1, n2, n3)];
+  vv.d11[0] = v.d11[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin2_al
+    + v.d1[Index (ivar, i, j, k, nvar, n1, n2, n3)] * cos_al;
+  vv.d12[0] =
+    v.d12[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al * sin_be;
+  vv.d13[0] = v.d13[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al;
+  vv.d22[0] = v.d22[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin2_be
+    + v.d2[Index (ivar, i, j, k, nvar, n1, n2, n3)] * cos_be;
+  vv.d23[0] = v.d23[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_be;
+  vv.d33[0] = v.d33[Index (ivar, i, j, k, nvar, n1, n2, n3)];
+}
+
+// -------------------------------------------------------------------------------
+double
+interpol (double a, double b, double c, derivs v)
+{
+  return v.d0[0]
+    + a * v.d1[0] + b * v.d2[0] + c * v.d3[0]
+    + 0.5 * a * a * v.d11[0] + a * b * v.d12[0] + a * c * v.d13[0]
+    + 0.5 * b * b * v.d22[0] + b * c * v.d23[0] + 0.5 * c * c * v.d33[0];
+}
+
+// -------------------------------------------------------------------------------
+// Calculates the value of v at an arbitrary position (A,B,phi)
+double
+PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
+			    int n2, int n3, derivs v, double x, double y,
+			    double z)
+{
+  DECLARE_CCTK_PARAMETERS;
+  double xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
+    result, Us, Ui;
+  int i, j, k;
+  derivs vv;
+  allocate_derivs (&vv, 1);
+
+  xs = x / par_b;
+  ys = y / par_b;
+  zs = z / par_b;
+  rs2 = ys * ys + zs * zs;
+  phi = atan2 (z, y);
+  if (phi < 0)
+    phi += 2 * Pi;
+
+  aux1 = 0.5 * (xs * xs + rs2 - 1);
+  aux2 = sqrt (aux1 * aux1 + rs2);
+  X = asinh (sqrt (aux1 + aux2));
+  R = asin (sqrt (-aux1 + aux2));
+  if (x < 0)
+    R = Pi - R;
+
+  A = 2 * tanh (0.5 * X) - 1;
+  B = tan (0.5 * R - Piq);
+  al = Pi - acos (A);
+  be = Pi - acos (B);
+
+  i = rint (al * n1 / Pi - 0.5);
+  j = rint (be * n2 / Pi - 0.5);
+  k = rint (0.5 * phi * n3 / Pi);
+
+  a = al - Pi * (i + 0.5) / n1;
+  b = be - Pi * (j + 0.5) / n2;
+  c = phi - 2 * Pi * k / n3;
+
+  calculate_derivs (i, j, k, ivar, nvar, n1, n2, n3, v, vv);
+  result = interpol (a, b, c, vv);
+  free_derivs (&vv, 1);
+
+//      Us = (A-1)*PunctEvalAtArbitPosition(v.d0, A, B, phi, n1, n2, n3);
+  Ui = (A - 1) * result;
+//      printf("%e %e %e Us=%e Ui=%e 1-Ui/Us=%e\n",x,y,z,Us,Ui,1-Ui/Us);
+
+  return Ui;
+
+/*	calculate_derivs(i+1, j,   k,   ivar, nvar, n1, n2, n3, v, vv, 2);
+	calculate_derivs(i,   j+1, k,   ivar, nvar, n1, n2, n3, v, vv, 3);
+	calculate_derivs(i+1, j+1, k,   ivar, nvar, n1, n2, n3, v, vv, 4);
+	calculate_derivs(i,   j,   k+1, ivar, nvar, n1, n2, n3, v, vv, 5);
+	calculate_derivs(i+1, j,   k+1, ivar, nvar, n1, n2, n3, v, vv, 6);
+	calculate_derivs(i,   j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 7);
+	calculate_derivs(i+1, j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 8);
+	v_intpol[2]=interpol(ha*(a-1),hb*b,    hp*c,    vv,2);
+	v_intpol[3]=interpol(ha*a,    hb*(b-1),hp*c,    vv,3);
+	v_intpol[4]=interpol(ha*(a-1),hb*(b-1),hp*c,    vv,4);
+	v_intpol[5]=interpol(ha*a,    hb*b,    hp*(c-1),vv,5);
+	v_intpol[6]=interpol(ha*(a-1),hb*b,    hp*(c-1),vv,6);
+	v_intpol[7]=interpol(ha*a,    hb*(b-1),hp*(c-1),vv,7);
+	v_intpol[8]=interpol(ha*(a-1),hb*(b-1),hp*(c-1),vv,8);
+	for(j=1;j<=8;j++)
+		printf("Value U[V] at origin:%16.15f\tj=%d\n",-2*v_intpol[j],j);*/
+/*	result = 0;
+	for(mi=i-npol; mi<=i+npol; mi++){
+		al_m  = Pih*(2*mi+1)/n1;
+		for(mj=j-npol; mj<=j+npol; mj++){
+			be_m  = Pih*(2*mj+1)/n2;  
+			for(mk=k-npol; mk<=k+npol; mk++){
+				phi_m = 2.*Pi*mk/n3;
+				g_m   = 1;
+				for(ni=i-npol; ni<=i+npol; ni++){
+					al_n  = Pih*(2*ni+1)/n1;
+					if(ni != mi) g_m *= (al - al_n)/(al_m - al_n);
+				}
+				for(nj=j-npol; nj<=j+npol; nj++){
+					be_n  = Pih*(2*nj+1)/n2;  
+					if(nj != mj) g_m *= (be - be_n)/(be_m - be_n);
+				}
+				for(nk=k-npol; nk<=k+npol; nk++){
+					phi_n = 2.*Pi*nk/n3;
+					if(nk != mk) g_m *= (phi-phi_n)/(phi_m - phi_n);
+				}
+				result += g_m*v.d0[Index(ivar,mi,mj,mk,nvar,n1,n2,n3)];
+			}
+		}
+	}	*/
+}
+
+// -------------------------------------------------------------------------------