liblorene-dev 0.0.0~cvs20161116+dfsg-1ubuntu4 / usr / include / lorene / C++ / Include / hole_bhns.h

/*
 *  Definition of Lorene class Hole_bhns
 *
 */

/*
 *   Copyright (c) 2005-2007 Keisuke Taniguchi
 *
 *   This file is part of LORENE.
 *
 *   LORENE is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License version 2
 *   as published by the Free Software Foundation.
 *
 *   LORENE is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with LORENE; if not, write to the Free Software
 *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#ifndef __HOLE_BHNS_H_ 
#define __HOLE_BHNS_H_ 

/*
 * $Id: hole_bhns.h,v 1.3 2014/10/13 08:52:35 j_novak Exp $
 * $Log: hole_bhns.h,v $
 * Revision 1.3  2014/10/13 08:52:35  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.2  2008/05/15 18:55:04  k_taniguchi
 * Change of some parameters and introduction of some
 * computational routines.
 *
 * Revision 1.1  2007/06/22 01:04:14  k_taniguchi
 * *** empty log message ***
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/hole_bhns.h,v 1.3 2014/10/13 08:52:35 j_novak Exp $
 *
 */

// Headers Lorene
#include "blackhole.h"
#include "metric.h"

namespace Lorene {

// External classes which appear in the declaration of class Hole_bhns:
class Star_bhns ; 

                    //-------------------------//
                    //     Class Hole_bhns     //
                    //-------------------------//

/**
 * Class for black holes in black hole-neutron star binaries.
 * \ingroup(star)
 *
 */
class Hole_bhns : public Black_hole {

    // Data :
    // -----
    protected:
        /** \c true  for the Neumann type BC for the lapconf function,
	 *  \c false  for the Dirichlet type BH
	 */
        bool bc_lapconf_nd ;

	/** \c true  for the first type BC for the lapconf function,
	 *  \c false  for the second type BH
	 */
        bool bc_lapconf_fs ;

	/** \c true  for an irrotational black hole, \c false  for
	 *  a corotating one
	 */
	bool irrotational ;

	/// Spin angular velocity of the black hole
	double omega_spin ;

	/// Part of the lapconf function from the numerical computation
	Scalar lapconf_auto_rs ;

	/// Part of the lapconf function from the analytic background
	Scalar lapconf_auto_bh ;

        /// Lapconf function generated by the black hole
        Scalar lapconf_auto ;  // = lapconf_auto_rs + lapconf_auto_bh

	/// Lapconf function generated by the companion star
	Scalar lapconf_comp ;

	/// Total lapconf function
	Scalar lapconf_tot ;  // = lapconf_auto + lapconf_comp

	/// Lapse function of the "black hole" part
	Scalar lapse_auto ;  // = lapconf_auto / confo_tot

	/// Total lapse function
	Scalar lapse_tot ;  // = lapconf_tot / confo_tot

	/** Derivative of the part of the lapconf function
	 *  from the numerical computation
	 */
	Vector d_lapconf_auto_rs ;

	/** Derivative of the part of the lapconf function
	 *  from the analytic background
	 */
	Vector d_lapconf_auto_bh ;

	/// Derivative of the lapconf function generated by the black hole
	Vector d_lapconf_auto ;

	/// Derivative of the lapconf function generated by the companion star
	Vector d_lapconf_comp ;

	/// Part of the shift vector from the numerical computation
	Vector shift_auto_rs ;

	/// Part of the shift vector from the analytic background
	Vector shift_auto_bh ;

	/// Shift vector generated by the black hole
	Vector shift_auto ;

	/// Shift vector generated by the companion star
	Vector shift_comp ;

	/// Total shift vector ;
	Vector shift_tot ;

	/** Derivative of the part of the shift vector
	 *  from the numerical computation
	 */
	Tensor d_shift_auto_rs ;

	/** Derivative of the part of the shift vector
	 *  from the analytic background
	 */
	Tensor d_shift_auto_bh ;

	/// Derivative of the shift vector generated by the black hole
	Tensor d_shift_auto ;

	/// Derivative of the shift vector generated by the companion star
	Tensor d_shift_comp ;

	/// Part of the conformal factor from the numerical computation
	Scalar confo_auto_rs ;

	/// Part of the conformal factor from the analytic background
	Scalar confo_auto_bh ;

	/// Conformal factor generated by the black hole
	Scalar confo_auto ;

	/// Conformal factor generated by the companion star
	Scalar confo_comp ;

	/// Total conformal factor
	Scalar confo_tot ;

	/** Derivative of the part of the conformal factor
	 *  from the numerical computation
	 */
	Vector d_confo_auto_rs ;

	/** Derivative of the part of the conformal factor
	 *  from the analytic background
	 */
	Vector d_confo_auto_bh ;

	/// Derivative of the conformal factor generated by the black hole
	Vector d_confo_auto ;

	/// Derivative of the conformal factor generated by the companion star
	Vector d_confo_comp ;

	/** Part of the extrinsic curvature tensor
	 *  from the numerical computation
	 */
	Sym_tensor taij_tot_rs ;

	/** Part of the extrinsic curvature tensor
	 *  from the rotation shift vector
	 */
	Sym_tensor taij_tot_rot ;

	/** Part of the extrinsic curvature tensor
	 *  from the analytic background
	 */
	Sym_tensor taij_tot_bh ;

	/** Total extrinsic curvature tensor \f$\ A^{ij}\f$
	 *  generated by \c shift_tot , \c lapse_tot , and
	 *  \c confo_tot .
	 */
	Sym_tensor taij_tot ;

	/** Part of the extrinsic curvature tensor numericalty computed
	 *  for the black hole.
	 */
	Sym_tensor taij_auto_rs ;

	/** Part of the extrinsic curvature tensor generated by
	 *  the black hole.
	 */
	Sym_tensor taij_auto ;

	/** Part of the extrinsic curvature tensor generated by
	 *  the companion star.
	 */
	Sym_tensor taij_comp ;

	/// Part of the scalar from the numerical computation
	Scalar taij_quad_tot_rs ;

	/// Part of the scalar from the rotation shift vector
	Scalar taij_quad_tot_rot ;

	/// Part of the scalar from the analytic background
	Scalar taij_quad_tot_bh ;

	/** Total scalar \f$\eta_{ik} \eta_{jl} A^{ij} A^{kl}\f$
	 *  generated by \f$A_{ij}^{\rm tot}\f$
	 */
	Scalar taij_quad_tot ;

	/// Part of the scalar from the black hole.
	Scalar taij_quad_auto ;

	/// Part of the scalar from the companion star.
	Scalar taij_quad_comp ;

    // Derived data : 
    // ------------
    protected:
	mutable double* p_mass_irr_bhns ; /// Irreducible mass of BH
	/// Quasilocal spin angular momentum of BH
	mutable double* p_spin_am_bhns ;

    // Constructors - Destructor
    // -------------------------
    public:

	/** Standard constructor.
	 *
	 *  @param mp_i Mapping on which the black hole will be defined
	 */
	Hole_bhns(Map& mp_i, bool kerrschild_i, bool bc_lapconf_nd,
		  bool bc_lapconf_fs, bool irrotational, double massbh) ;

	Hole_bhns(const Hole_bhns& ) ;       ///< Copy constructor

	/** Constructor from a file (see \c sauve(FILE*) )
	 *
	 *  @param mp_i Mapping on which the black hole will be defined
	 *  @param fich input file (must have been created by the function
	 *         \c sauve )
	 */
	Hole_bhns(Map& mp_i, FILE* fich) ;    		

	virtual ~Hole_bhns() ;			///< Destructor

    // Memory management
    // -----------------
    protected:
	/// Deletes all the derived quantities
	virtual void del_deriv() const ;

	/// Sets to \c 0x0 all the pointers on derived quantities
	void set_der_0x0() const ;


    // Mutators / assignment
    // ---------------------
    public:
	/// Assignment to another \c Hole_bhns
	void operator=(const Hole_bhns&) ;

	/// Sets the spin angular velocity of the black hole [{\tt f\_unit}]
	double& set_omega_spin() { return omega_spin ; } ;

	/// Read/write of the lapconf function generated by the black hole
	Scalar& set_lapconf_auto_rs() ;

	/// Read/write of the lapconf function generated by the black hole
	Scalar& set_lapconf_auto_bh() ;

	/// Read/write of the lapconf function generated by the black hole
	Scalar& set_lapconf_auto() ;

	/// Read/write of the lapconf function generated by the companion star
	Scalar& set_lapconf_comp() ;

	/// Read/write of the total lapconf function
	Scalar& set_lapconf_tot() ;

	/// Read/write of the lapse function generated by the black hole
	Scalar& set_lapse_auto() ;

	/// Read/write of the total lapse function
	Scalar& set_lapse_tot() ;

	/// Read/write of the shift vector generated by the black hole
	Vector& set_shift_auto_rs() ;

	/// Read/write of the shift vector generated by the black hole
	Vector& set_shift_auto_bh() ;

	/// Read/write of the shift vector generated by the black hole
	Vector& set_shift_auto() ;

	/// Read/write of the shift vector generated by the companion star
	Vector& set_shift_comp() ;

	/// Read/write of the total shift vector
	Vector& set_shift_tot() ;

	/// Read/write of the conformal factor generated by the black hole
	Scalar& set_confo_auto_rs() ;

	/// Read/write of the conformal factor generated by the black hole
	Scalar& set_confo_auto_bh() ;

	/// Read/write of the conformal factor generated by the black hole
	Scalar& set_confo_auto() ;

	/// Read/write of the conformal factor generated by the companion star
	Scalar& set_confo_comp() ;

	/// Read/write of the total conformal factor
	Scalar& set_confo_tot() ;

    // Accessors
    // ---------
    public:
	/** Returns \c true  for the Neumann type BC for the lapconf function,
	 *  \c false  for the Dirichlet type BH
	 */
	bool has_bc_lapconf_nd() const {return bc_lapconf_nd; } ;

	/** Returns \c true  for the first type BC for the lapconf function,
	 *  \c false  for the second type BH
	 */
	bool has_bc_lapconf_fs() const {return bc_lapconf_fs; } ;

	/** Returns \c true  for an irrotational black hole, \c false  for
	 *  a corotating one
	 */
	bool is_irrotational() const {return irrotational; } ;

	/// Returns the spin angular velocity of the black hole [{\tt f\_unit}]
	double get_omega_spin() const {return omega_spin; } ;

	/// Returns the part of the lapconf function from numerical computation
	const Scalar& get_lapconf_auto_rs() const {return lapconf_auto_rs; } ;

	/** Returns the part of the lapconf function generated by
	 *  the black hole
	 */
	const Scalar& get_lapconf_auto() const {return lapconf_auto; } ;

	/** Returns the part of the lapconf function generated by the
	 *  companion star
	 */
	const Scalar& get_lapconf_comp() const {return lapconf_comp; } ;

	/// Returns the total lapconf function
	const Scalar& get_lapconf_tot() const {return lapconf_tot; } ;

	/** Returns the part of the lapse function generated by
	 *  the black hole
	 */
	const Scalar& get_lapse_auto() const {return lapse_auto; } ;

	/// Returns the total lapse function
	const Scalar& get_lapse_tot() const {return lapse_tot; } ;

	/** Returns the derivative of the lapconf function generated by
	 *  the black hole
	 */
	const Vector& get_d_lapconf_auto_rs() const
	  {return d_lapconf_auto_rs; } ;

	/** Returns the derivative of the lapconf function generated by
	 *  the black hole
	 */
	const Vector& get_d_lapconf_auto() const {return d_lapconf_auto; } ;

	/** Returns the derivative of the lapconf function generated by the
	 *  companion star
	 */
	const Vector& get_d_lapconf_comp() const {return d_lapconf_comp; } ;

	/// Returns the part of the shift vector from numerical computation
	const Vector& get_shift_auto_rs() const {return shift_auto_rs; } ;

	/// Returns the part of the shift vector generated by the black hole
	const Vector& get_shift_auto() const {return shift_auto; } ;

	/** Returns the part of the shift vector generated by the
	 *  companion star
	 */
	const Vector& get_shift_comp() const {return shift_comp; } ;

	/// Returns the total shift vector
	const Vector& get_shift_tot() const {return shift_tot; } ;

	/** Returns the derivative of the shift vector generated by the
	 *  black hole
	 */
	const Tensor& get_d_shift_auto_rs() const {return d_shift_auto_rs; } ;

	/** Returns the derivative of the shift vector generated by the
	 *  black hole
	 */
	const Tensor& get_d_shift_auto() const {return d_shift_auto; } ;

	/** Returns the derivative of the shift vector generated by the
	 *  companion star
	 */
	const Tensor& get_d_shift_comp() const {return d_shift_comp; } ;

	/// Returns the part of the conformal factor from numerical computation
	const Scalar& get_confo_auto_rs() const {return confo_auto_rs; } ;

	/** Returns the part of the conformal factor generated by the
	 *  black hole
	 */
	const Scalar& get_confo_auto() const {return confo_auto; } ;

	/** Returns the part of the conformal factor generated by the
	 *  companion star
	 */
	const Scalar& get_confo_comp() const {return confo_comp; } ;

	/// Returns the total conformal factor
	const Scalar& get_confo_tot() const {return confo_tot; } ;

	/** Returns the derivative of the conformal factor generated by
	 *  the black hole
	 */
	const Vector& get_d_confo_auto_rs() const {return d_confo_auto_rs; } ;

	/** Returns the derivative of the conformal factor generated by the
	 *  black hole
	 */
	const Vector& get_d_confo_auto() const {return d_confo_auto; } ;

	/** Returns the derivative of the conformal factor generated by the
	 *  companion star
	 */
	const Vector& get_d_confo_comp() const {return d_confo_comp; } ;

	/// Returns the part of rs of the extrinsic curvature tensor
	const Sym_tensor& get_taij_tot_rs() const {return taij_tot_rs; } ;

	/// Returns the total extrinsic curvature tensor
	const Sym_tensor& get_taij_tot() const {return taij_tot; } ;

	/** Returns the part of rs of the extrinsic curvature tensor
	 *  generated by the black hole
	 */
	const Sym_tensor& get_taij_auto_rs() const {return taij_auto_rs; } ;

	/** Returns the part of the extrinsic curvature tensor
	 *  generated by the black hole
	 */
	const Sym_tensor& get_taij_auto() const {return taij_auto; } ;

	/** Returns the part of the extrinsic curvature tensor
	 *  generated by the companion star
	 */
	const Sym_tensor& get_taij_comp() const {return taij_comp; } ;

	/// Returns the part of rs
	const Scalar& get_taij_quad_tot_rs() const {return taij_quad_tot_rs; } ;
	/// Returns the part of rot
	const Scalar& get_taij_quad_tot_rot() const {return taij_quad_tot_rot; } ;

	/// Returns the total scalar \f$\eta_{ik} \eta_{jl} A^{ij} A^{kl}\f$
	const Scalar& get_taij_quad_tot() const {return taij_quad_tot; } ;

	/// Returns the part of rs generated by the black hole
	const Scalar& get_taij_quad_auto() const {return taij_quad_auto; } ;

	/// Returns the part of rs generated by the companion star
	const Scalar& get_taij_quad_comp() const {return taij_quad_comp; } ;

    // Outputs
    // -------
    public:
	virtual void sauve(FILE *) const ;	    ///< Save in a file
    
    protected:
	/// Operator >> (virtual function called by the operator <<)
	virtual ostream& operator>>(ostream& ) const ;

    // Global quantities
    // -----------------
    public:
	/// Irreducible mass of the black hole
	virtual double mass_irr_bhns() const ;

	/// Spin angular momentum of the black hole
	double spin_am_bhns(const Tbl& xi_i, const double& phi_i,
			    const double& theta_i, const int& nrk_phi,
			    const int& nrk_theta) const ;

    // Computational routines
    // ----------------------
    public:
	/** Boundary condition on the apparent horizon of the black hole
	 *  for the lapconf function: 2-D \c Valeur
	 */
	const Valeur bc_lapconf() const ;

	/** Boundary condition on the apparent horizon of the black hole
	 *  for the shift vector of the \fx\f direction: 2-D \c Valeur
	 */
	const Valeur bc_shift_x(double ome_orb, double y_rot) const ;

	/** Boundary condition on the apparent horizon of the black hole
	 *  for the shift vector of the \fy\f direction: 2-D \c Valeur
	 */
	const Valeur bc_shift_y(double ome_orb, double x_rot) const ;

	/** Boundary condition on the apparent horizon of the black hole
	 *  for the shift vector of the \fz\f direction: 2-D \c Valeur
	 */
	const Valeur bc_shift_z() const ;

	/** Boundary condition on the apparent horizon of the black hole
	 *  for the conformal factor: 2-D \c Valeur
	 *
	 *  @param ome_orb orbital angular velocity
	 *  @param x_rot position of the rotation axis
	 */
	const Valeur bc_confo(double ome_orb, double x_rot,
			      double y_rot) const ;

	/** Computes metric coefficients from known potentials
	 *  with relaxation when the companion is a black hole.
	 *
	 *  The calculation is performed starting from the quantities
	 *  \c lapse_auto , \c shift_auto , \c confo_auto ,
	 *  \c comp.lapse_auto , \c comp.confo_auto
	 *  which are supposed to be up to date.
	 *  From these, the following fields are updated:
	 *  \c lapse_comp , \c lapse_tot, \c confo_comp,
	 *  \c confo_tot , \c psi4 ,
	 *
	 *  @param star companion neutron star
	 *  @param hole_prev previous value of the hole
	 *  @param relax relaxation parameter
	 *
	 */
	void update_metric_bhns(const Star_bhns& star,
				const Hole_bhns& hole_prev,
				double relax) ;

	/** Computes derivative of metric quantities from
	 *  the companion neutron star
	 *
	 *  @param star companion neutron star
	 *
	 */
	void update_met_der_comp_bhns(const Star_bhns& star) ;

	/** Computes \c taij_tot , \c taij_quad_tot  from
	 *  \c shift_tot , \c lapse_tot , \c confo_tot .
	 *
	 *  @param omega_orb [input] Orbital angular velocity
	 *  @param x_rot [input] Position of the rotation axis
	 */
	void extr_curv_bhns(double omega_orb, double x_rot, double y_rot) ;

	/** Computes a black-hole part in a black hole-neutron star binary
	 *  by giving boundary conditions on the apparent horizon.
	 *
	 *  @param mer [input] Number of iteration
	 *  @param mermax_bh [input] Maximum number of iteration steps
	 *  @param filter_r [input] No. points to be deleted for (r): lap,conf
	 *  @param filter_r_s [input] No. points to be deleted for (r): shift
	 *  @param filter_p_s [input] No. points to be deleted for (phi): shift
	 *  @param x_rot [input] Position of the rotation axis
	 *  @param y_rot [input] Position of the rotation axis
	 *  @param precis [input] threshold in the relative difference of
	 *                a metric quantity between two succesive steps
	 *                to stop the iterative procedure
	 *                (default value: 1.e-14)
	 *  @param omega_orb [input] Orbital angular velocity
	 *  @param resize_bh [input] Resize factor for the first shell
	 *  @param fact_resize [input] 1-D \c Tbl  for the input of some
	 *                             factors : \\
	 *              \c fact(0)  : A resizing factor for the first shell
	 *  @param diff [output] 1-D \c Tbl  for the storage of some
	 *                       error indicators :
	 */
	void equilibrium_bhns(int mer, int mermax_bh,
			      int filter_r, int filter_r_s, int filter_p_s,
			      double x_rot, double y_rot, double precis,
			      double omega_orb, double resize_bh,
			      const Tbl& fact_resize, Tbl& diff) ;

	/** Performs a relaxation on \c lapse_auto ,
	 *  \c shift_auto , \c confo_auto .
	 *
	 *  @param hole_prev [input] black hole at the previous step
	 *  @param relax_met [input] Relaxation factor for \c lapse_auto ,
	 *                           \c shift_auto , \c confo_auto ,
	 *                           only if \c (mer%fmer_met=0) .
	 *  @param mer       [input] Step number
	 *  @param fmer_met  [input] Step interval between metric updates
	 */
	void relax_bhns(const Hole_bhns& hole_prev,
			double relax_met, int mer, int fmer_met) ;

	/** Compute a forth-order Runge-Kutta integration to the phi
	 *  direction for the solution of the Killing vectors on the
	 *  equator
	 *
	 *  @param xi_i [input] initial set of the Killing vectors
	 *              at the starting point on the equator
	 *  @param phi_i [input] initial phi coordinate at which the
	 *               integration starts
	 *  @param nrk [input] number of the Runge-Kutta integration
	 *             between two successive collocation points
	 */
	Tbl runge_kutta_phi(const Tbl& xi_i, const double& phi_i,
			    const int& nrk) const ;

	/** Compute a forth-order Runge-Kutta integration to the theta
	 *  direction for the solution of the Killing vectors
	 *
	 *  @param xi_i [input] initial set of the Killing vectors
	 *              at the starting point on the equator
	 *  @param theta_i [input] initial theta coordinate at which the
	 *                 integration starts
	 *  @param phi [input] fixed phi coordinate during integration to
	 *             the theta direction
	 *  @param nrk [input] number of the Runge-Kutta integration
	 *             between two successive collocation points
	 */

	Tbl runge_kutta_theta(const Tbl& xi_i, const double& theta_i,
			      const double& phi, const int& nrk) const ;

	/** Compute the Killing vector of a black hole normalized so that
	 *  its affine length is 2 M_PI
	 *
	 *  @param xi_i [input] initial set of the Killing vectors
	 *              at the starting point on the equator
	 *  @param phi_i [input] initial phi coordinate at which the
	 *               integration starts
	 *  @param theta_i [input] initial theta coordinate at which the
	 *                 integration starts
	 *  @param nrk_phi [input] number of the Runge-Kutta integration
	 *                 between two successive collocation points
	 *                 for the phi direction
	 *  @param nrk_theta [input] number of the Runge-Kutta integration
	 *                   between two successive collocation points
	 *                   for the theta direction
	 */
	Vector killing_vect(const Tbl& xi_i, const double& phi_i,
			    const double& theta_i, const int& nrk_phi,
			    const int& nrk_theta) const ;

	friend class Bin_bhns ;

};

}
#endif