Equations.cpp

Go to the documentation of this file.

/* TwoPunctures:  File  "Equations.c"*/
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include "libtwopunctures/TwoPunctures.h"
 
namespace TP {
using namespace Utilities;
 
/* U.d0[ivar]   = U[ivar];  (ivar = 0..nvar-1) */
/* U.d1[ivar]   = U[ivar]_x;  */
/* U.d2[ivar]   = U[ivar]_y;  */
/* U.d3[ivar]   = U[ivar]_z;  */
/* U.d11[ivar]  = U[ivar]_xx; */
/* U.d12[ivar]  = U[ivar]_xy; */
/* U.d13[ivar]  = U[ivar]_xz;*/
/* U.d22[ivar]  = U[ivar]_yy;*/
/* U.d23[ivar]  = U[ivar]_yz;*/
/* U.d33[ivar]  = U[ivar]_zz;*/
 
double
TwoPunctures::BY_KKofxyz (double x, double y, double z)
{
  
  int i, j;
  double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
    Aij, AijAij, n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
 
  r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
  r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
  r_plus = sqrt (r2_plus);
  r_minus = sqrt (r2_minus);
  r3_plus = r_plus * r2_plus;
  r3_minus = r_minus * r2_minus;
 
  n_plus[0] = (x - par_b) / r_plus;
  n_minus[0] = (x + par_b) / r_minus;
  n_plus[1] = y / r_plus;
  n_minus[1] = y / r_minus;
  n_plus[2] = z / r_plus;
  n_minus[2] = z / r_minus;
 
  /* dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) */
  np_Pp = 0;
  nm_Pm = 0;
  for (i = 0; i < 3; i++)
  {
    np_Pp += n_plus[i] * par_P_plus[i];
    nm_Pm += n_minus[i] * par_P_minus[i];
  }
  /* cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i*/
  np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
  np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
  np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
  nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
  nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
  nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
  AijAij = 0;
  for (i = 0; i < 3; i++)
  {
    for (j = 0; j < 3; j++)
    {               /* Bowen-York-Curvature :*/
      Aij =
    + 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
                 + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
    + 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
         + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
    - 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
    - 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
      if (i == j)
    Aij -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
      AijAij += Aij * Aij;
    }
  }
 
  return AijAij;
}
 
void
TwoPunctures::BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
{
  
  int i, j;
  double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
    n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
 
  const double TP4 = TP_epsilon*TP_epsilon*TP_epsilon*TP_epsilon;
  r2_plus  = (x - par_b) * (x - par_b) + y * y + z * z;
  r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
  r2_plus  = sqrt (r2_plus*r2_plus   + TP4);
  r2_minus = sqrt (r2_minus*r2_minus + TP4);
 
  const double TPT2 = TP_Tiny*TP_Tiny;
  if (r2_plus < TPT2)
    r2_plus = TPT2;
  if (r2_minus < TPT2)
    r2_minus = TPT2;
  r_plus   = sqrt (r2_plus);
  r_minus  = sqrt (r2_minus);
  r3_plus  = r_plus * r2_plus;
  r3_minus = r_minus * r2_minus;
 
  n_plus[0]  = (x - par_b) / r_plus;
  n_minus[0] = (x + par_b) / r_minus;
  n_plus[1]  = y / r_plus;
  n_minus[1] = y / r_minus;
  n_plus[2]  = z / r_plus;
  n_minus[2] = z / r_minus;
 
  /* dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) */
  np_Pp = 0;
  nm_Pm = 0;
  for (i = 0; i < 3; i++)
  {
    np_Pp += n_plus[i] * par_P_plus[i];
    nm_Pm += n_minus[i] * par_P_minus[i];
  }
  /* cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i*/
  np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
  np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
  np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
  nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
  nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
  nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
  for (i = 0; i < 3; i++)
  {
    for (j = 0; j < 3; j++)
    {               /* Bowen-York-Curvature :*/
      Aij[i][j] =
        + 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
         + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
    + 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
         + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
    - 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
    - 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
    }
  }
 
  Aij[0][0] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
  Aij[1][1] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
  Aij[2][2] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
 
 
}
 
/*-----------------------------------------------------------*/
/********           Nonlinear Equations                ***********/
/*-----------------------------------------------------------*/
void
TwoPunctures::NonLinEquations (double rho_adm,
     double A, double B, double X, double R,
         double x, double r, double phi,
         double y, double z, derivs U, double *values)
{
  
  double r_plus, r_minus, psi, psi2, psi4, psi7;
  double mu;
 
  r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
  r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
 
  psi =
    1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
  psi2 = psi * psi;
  psi4 = psi2 * psi2;
  psi7 = psi * psi2 * psi4;
 
  values[0] =
    U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7 +
    2.0 * Pi / psi2/psi * rho_adm;
 
}
 
/*-----------------------------------------------------------*/
/********               Linear Equations                ***********/
/*-----------------------------------------------------------*/
void
TwoPunctures::LinEquations (double A, double B, double X, double R,
          double x, double r, double phi,
          double y, double z, derivs dU, derivs U, double *values)
{
  
  double r_plus, r_minus, psi, psi2, psi4, psi8;
 
  r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
  r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
 
  psi =
    1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
  psi2 = psi * psi;
  psi4 = psi2 * psi2;
  psi8 = psi4 * psi4;
 
  values[0] = dU.d11[0] + dU.d22[0] + dU.d33[0]
    - 0.875 * BY_KKofxyz (x, y, z) / psi8 * dU.d0[0];
}
 
/*-----------------------------------------------------------*/
 
} // namespace TP