Globals.h

//---------------------------------------------------------------------------
//
// $Id: Globals.h,v 1.6 2018/08/02 23:05:02 sahughes Exp $
//
//---------------------------------------------------------------------------
#ifndef _Globals_H
#define _Globals_H

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <complex>

using std::cout;
using std::cerr;
using std::endl;
using std::ios;

// For convenience, can change the floating-point type with a
// single preprocessor symbol
#define USE_FLOAT       0
#define USE_DOUBLE      1
#define USE_LONG_DOUBLE 2

#if !defined(REAL_TYPE) || REAL_TYPE == USE_DOUBLE
typedef double Real;
#elif REAL_TYPE == USE_FLOAT
typedef float Real;
#elif REAL_TYPE == USE_LONG_DOUBLE
typedef long double Real;
#else
Illegal value!
#endif

typedef std::complex<Real> Complex;
const Complex II = Complex(0.,1.);

// Reasonable maximum
template<class NUMBER>
inline NUMBER Max(const NUMBER a, const NUMBER b) {
  return (a > b ? a : b);
};

// Reasonable minimum
template<class NUMBER>
inline NUMBER Min(const NUMBER a, const NUMBER b) {
  return (a < b ? a : b);
};


class Kerr {
 public:
  static Real rplus(const Real a)
    {
      return(1. + sqrt((1. - a)*(1. + a)));
    }; 

  static Real rminus(const Real a)
    {
      return(1. - sqrt((1. - a)*(1. + a)));
    };

  static Real rstar(const Real r, const Real a)
    {
      const Real rm = rminus(a);
      const Real rp = rplus(a);
      
      return(r + ((2.*rp)/(rp - rm))*log((r - rp)/2.) -
       ((2.*rm)/(rp - rm))*log((r - rm)/2.));
    };

  static Real Delta(const Real r, const Real a)
    {
      return(r*r - 2.*r + a*a);
    }

  static Real dr_Delta(const Real r)
    {
      return(2.*(r - 1.));
    }

  static Real ddr_Delta()
    {
      return(2.);
    }

  static Real Sigma(const Real r, const Real a, const Real z)
    {
      return(r*r + a*a*z);
    };

  static Real dr_Sigma(const Real r)
    {
      return(2.*r);
    };

  static Real ddr_Sigma()
    {
      return(2.);
    };

  static Real Eeqpro(const Real r, const Real a)
    {
      const Real v = 1./sqrt(r);

      const Real numer = 1. - v*v*(2. - a*v);
      const Real denom = 1. - v*v*(3. - 2.*a*v);

      return(numer/sqrt(denom));
    };

  static Real Eeqret(const Real r, const Real a)
    {
      const Real v = 1./sqrt(r);
      
      const Real numer = 1. - v*v*(2. + a*v);
      const Real denom = 1. - v*v*(3. + 2.*a*v);
      
      return(numer/sqrt(denom));
    };
  
  static Real Lzeqpro(const Real r, const Real a)
    {
      const Real v = 1./sqrt(r);
      
      const Real numer = 1. - a*v*v*v*(2. - a*v);
      const Real denom = 1. - v*v*(3. - 2.*a*v);
      
      return(r*v*numer/sqrt(denom));
    };

  static Real Lzeqret(const Real r, const Real a)
    {
      const Real v = 1./sqrt(r);
      
      const Real numer = 1. + a*v*v*v*(2. + a*v);
      const Real denom = 1. - v*v*(3. + 2.*a*v);
      
      return(-r*v*numer/sqrt(denom));
    };

  static Real Omega_phi_eqpro(const Real r, const Real a)
    {
      return(1./(sqrt(r*r*r) + a));
    };

  static Real Omega_phi_eqret(const Real r, const Real a)
    {
      return(-1./(sqrt(r*r*r) - a));
    };

  static Real isco_pro(const Real a)
    {
      const Real Z1 = 1. + (pow(1. + a, 1./3.) + pow(1. - a, 1./3.))*
  pow((1. + a)*(1. - a), 1./3.);
      const Real Z2 = sqrt(3.*a*a + Z1*Z1);

      return(3. + Z2 - sqrt((3. - Z1)*(3. + Z1 + 2.*Z2)));
    };

  static Real isco_ret(const Real a)
    {
      const Real Z1 = 1. + (pow(1. + a, 1./3.) + pow(1. - a, 1./3.))*
  pow((1. + a)*(1. - a), 1./3.);
      const Real Z2 = sqrt(3.*a*a + Z1*Z1);

      return(3. + Z2 + sqrt((3. - Z1)*(3. + Z1 + 2.*Z2)));
    };
};


struct DataHolder {
  int l, m, k;
  Real r, a, Lz, E, Q;
  Real w, p, lambda;
  Complex ZI, ZH;
  Real EdotInf, LzdotInf, QdotInf, rdotInf;
  Real EdotH, LzdotH, QdotH, rdotH;
};
//
// Index range
//
#define L_MAX 100
#define K_MAX 100

#endif