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

/*
 *  Definition of Lorene class Star_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 __STAR_BHNS_H_ 
#define __STAR_BHNS_H_ 

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

// Lorene headers
#include "star.h"

namespace Lorene {

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

/**
 * Class for stars in black hole-neutron star binaries.
 * \ingroup(star)
 * 
 */
class Star_bhns : public Star {

    // Data : 
    // -----
    protected:
        /** Affine mapping for solving poisson's equations of
	 *  metric quantities
	 */
        Map_af mp_aff ;

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

	/** Scalar potential \f$\Psi_0\f$ of the non-translational part of
	 *  fluid 4-velocity (in the irrotational case)
	 */
	Scalar psi0 ;

	/** Gradient of \f$\Psi\f$ (in the irrotational case)
	 *  (Spherical components with respect to the mapping of the star)
	 */
	Vector d_psi ;

	/** Spatial projection of the fluid 3-velocity with respect to
	 *  the co-orbiting observer.
	 *  (Spherical components with respect to the mapping of the star)
	 */
	Vector wit_w ;

	/** Logarithm of the Lorentz factor between the fluid and
	 *  the co-orbiting observer.
	 */
	Scalar loggam ;

	/** 3-vector shift, divided by \e N , of the rotating coordinates,
	 *  \f$\beta^i/N\f$.
	 *  (Spherical components with respect to the mapping of the star)
	 */
	Vector bsn ;

	/// Lorentz factor between the fluid and the co-orbiting observer
	Scalar gam ;

	/** Lorentz factor between the co-orbiting observer
	 *  and the Eulerian one
	 */
	Scalar gam0 ;

	/// Centrifugal potential
	Scalar pot_centri ;

	/// Lapconf function generated by the star
        Scalar lapconf_auto ;

	/// Lapconf function generated by the companion black hole
	Scalar lapconf_comp ;

	/// Total lapconf function
	Scalar lapconf_tot ;

	/// Lapse function generated by the "star"
        Scalar lapse_auto ;  // = lapconf_auto / confo_tot

	/// Total lapse function
	Scalar lapse_tot ;

	/** Derivative of the lapconf function generated by the star
	 *  \f$ \partial_j \alpha \f$
	 */
	Vector d_lapconf_auto ;

	/** Derivative of the lapconf function generated by the companion
	 *  black hole
	 */
	Vector d_lapconf_comp ;

	/// Shift vector generated by the star
	Vector shift_auto ;

	/// Shift vector generated by the companion black hole
	Vector shift_comp ;

	/// Total shift vector
	Vector shift_tot ;

	/** Derivative of the shift vector generated by the star
	 *  \f$ \eta^{ik} \partial_k \beta^j \f$
	 */
	Tensor d_shift_auto ;

	/** Derivative of the shift vector generated by the companion
	 *  black hole
	 */
	Tensor d_shift_comp ;

	/// Conformal factor generated by the star
	Scalar confo_auto ;

	/// Conformal factor generated by the companion black hole
	Scalar confo_comp ;

	/// Total conformal factor
	Scalar confo_tot ;

	/** Derivative of the conformal factor generated by the star
	 *  \f$ \partial_j \psi \f$
	 */
	Vector d_confo_auto ;

	/** Derivative of the conformal factor generated by the companion
	 *  black hole
	 */
	Vector d_confo_comp ;

	/// Fourth power of the total conformal factor
	Scalar psi4 ;  // psi4 = pow(confo_tot, 4.)

	/** Part of the extrinsic curvature tensor \f$ A^{ij}\f$
	 *  generated by \c shift_auto , \c lapse_auto , and
	 *  \c confo_auto .
	 */
	Sym_tensor taij_auto ;

	/** Part of the scalar \f$\eta_{ik} \eta_{jl} A^{ij} A^{kl}\f$
	 *  generated by \f$A_{ij}^{\rm auto}\f$
	 */
	Scalar taij_quad_auto ;

	/** Flat metric defined on the mapping (Spherical components
	 *  with respect to the mapping of the star ).
	 */
	Metric_flat flat ;

	/** Effective source at the previous step for the resolution of
	 *  the Poisson equation for \c lapconf_auto .
	 */
	Scalar ssjm1_lapconf ;

	/** Effective source at the previous step for the resolution of
	 *  the Poisson equation for \c confo_auto .
	 */
	Scalar ssjm1_confo ;

	/** Effective source at the previous step for the resolution of
	 *  the Poisson equation for the scalar \f$\chi\f$ by means of
	 *  \c Map_et::poisson .
	 *  \f$\chi\f$ is an intermediate quantity for the resolution of the
	 *  elliptic equation for the shift vector \f$N^i\f$
	 */
	Scalar ssjm1_khi ;

	/** Effective source at the previous step for the resolution of
	 *  the vector Poisson equation for \f$W^i\f$ by means of
	 *  \c Map_et::poisson .
	 *  \f$W^i\f$ is an intermediate quantity for the resolution of the
	 *  elliptic equation for the shift vector \f$N^i\f$
	 *  (Components with respect to the Cartesian triad associated with
	 *  the mapping \c mp )
	 */
	Vector ssjm1_wshift ;

    // Derived data
    // ------------
    protected:
	mutable double* p_mass_b_bhns ;  ///< Baryon mass
	mutable double* p_mass_g_bhns ;  ///< Gravitational mass

    // Constructors - Destructor
    // -------------------------
    public:
	/** Standard constructor
	 *
	 *  @param mp_i Mapping on which the star will be defined
	 *  @param nzet_i Number of domains occupied by the star
	 *  @param eos_i Equation of state of the stellar matter
	 *  @param irrot_i should be \c true  for an irrotational star,
	 *                 \c false  for a corotating one
	 */
	Star_bhns(Map& mp_i, int nzet_i, const Eos& eos_i, bool irrot_i) ;

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

	/** Constructor from a file (see \c sauve(FILE*) )
	 *   @param mp_i  Mapping on which the star will be defined
	 *   @param eos_i Equation of state of the stellar matter
	 *   @param fich  input file (must have been created by the function
	 *                \c sauve )
	 */
	Star_bhns(Map& mp_i, const Eos& eos_i, FILE* fich) ;

	virtual ~Star_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 Star_bhns
	void operator=(const Star_bhns&) ;

	/// Read/write the centrifugal potential
	Scalar& set_pot_centri() ;

	/// Read/write of the lapconf function generated by the neutron star
	Scalar& set_lapconf_auto() ;

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

	/// Read/write of the shift vector generated by the neutron star
	Vector& set_shift_auto() ;

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

	/// Read/write of the conformal factor generated by the neutron star
	Scalar& set_confo_auto() ;

	/** Read/write of the conformal factor generated by the companion
	 *  black hole
	 */
	Scalar& set_confo_comp() ;
	
    // Accessors
    // ---------
    public:
	/** Returns \c true  for an irrotational motion, \c false  for
	 *  a corotating one.
	 */
	bool is_irrotational() const {return irrotational; } ;

	/// Returns the non-translational part of the velocity potential
	const Scalar& get_psi0() const {return psi0; } ;

	/** Returns the covariant derivative of the velocity potential 
	 *  (Spherical components with respect to the mapping of the star)
	 */
	const Vector& get_d_psi() const {return d_psi; } ;

	/** Returns the spatial projection of the fluid 3-velocity with
	 *  respect to the co-orbiting observer.
	 *  (Spherical components with respect to the mapping of the star)
	 */
	const Vector& get_wit_w() const {return wit_w; } ;

	/** Returns the logarithm of the Lorentz factor between the fluid and
	 *  the co-orbiting observer.
	 */
	const Scalar& get_loggam() const {return loggam; } ;

	/** Returns the shift vector, divided by \e N , of the rotating
	 *  coordinates, \f$\beta^i/N\f$.
	 *  (Spherical components with respect to the mapping of the star)
	 */
	const Vector& get_bsn() const {return bsn; } ;

	/// Returns the Lorentz factor gam
	const Scalar& get_gam() const {return gam; } ;

	/// Returns the Lorentz factor gam0
	const Scalar& get_gam0() const {return gam0; } ;

	/// Returns the centrifugal potential
	const Scalar& get_pot_centri() const {return pot_centri; } ;

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

	/** Returns the part of the lapconf function generated by the
	 *  companion black hole
	 */
	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 star
	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 lapse function generated by the star
	const Vector& get_d_lapconf_auto() const {return d_lapconf_auto; } ;

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

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

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

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

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

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

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

	/** Returns the part of the conformal factor generated by the
	 *  companion black hole
	 */
	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 star
	 */
	const Vector& get_d_confo_auto() const {return d_confo_auto; } ;

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

	/// Returns the fourth power of the total conformal factor
	const Scalar& get_psi4() const {return psi4; } ;

	/** Returns the part of the extrinsic curvature tensor
	 *  \f$\tilde A^{ij}\f$ generated by the neutron star part.
	 */
	const Sym_tensor& get_taij_auto() const {return taij_auto; } ;

	/** Returns the part of the scalar
	 *  \f$\eta_{ik} \eta_{jl} A^{ij} A^{kl}\f$
	 *  generated by \f$A_{ij}^{\rm auto}\f$
	 */
	const Scalar& get_taij_quad_auto() const {return taij_quad_auto; } ;

    // 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:
	/// Baryon mass
	virtual double mass_b() const ;

	virtual double mass_b_bhns(bool kerrschild, const double& mass_bh,
				   const double& sepa) const ;

	/// Gravitational mass
	virtual double mass_g() const ;

	virtual double mass_g_bhns() const ;

    // Computational routines
    // ----------------------
    public:
	/** Computes the hydrodynamical quantities relative to the Eulerian
	 *  observer from those in the fluid frame, as well as
	 *  \c wit_w  and \c loggam .
	 *
	 *  The calculation is performed starting from the quantities
	 *  \c ent , \c ener , \c press , \c a_car  and \c bsn ,
	 *  which are supposed to be up to date.
	 *  From these,  the following fields are updated:
	 *  \c gam_euler , \c u_euler , \c ener_euler ,
	 *  \c s_euler , \c stress_euler ,
	 *  \c wit_w  and \c loggam .
	 *
	 *  @param kerrschild should be \c true  for a Kerr-Schild background,
	 *                    \c false  for a Conformally flat one
	 *  @param mass_bh BH mass in the background metric
	 *  @param sepa Separation between NS and BH
	 *
	 */
	void hydro_euler_bhns(bool kerrschild, const double& mass_bh,
			      const double& sepa) ;

	/** 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 hole companion black hole
	 *  @param star_prev previous value of the star
	 *  @param relax relaxation parameter
	 *
	 */
	void update_metric_bhns(const Hole_bhns& hole,
				const Star_bhns& star_prev,
				double relax) ;

	/** Computes derivative of metric quantities from
	 *  the companion black hole
	 *
	 *  @param hole companion black hole
	 *
	 */
	void update_met_der_comp_bhns(const Hole_bhns& hole) ;

	/** Computes the quantities \c bsn  and \c pot_centri .
	 *
	 *  The calculation is performed starting from the quantities
	 *  \c lapse_tot , \c shift_tot ,
	 *  which are supposed to be up to date.
	 *
	 *  @param kerrschild should be \c true  for a Kerr-Schild background,
	 *                    \c false  for a Conformally flat one
	 *  @param mass_bh BH mass in the background metric
	 *  @param sepa Separation between NS and BH
	 *  @param omega angular velocity with respect to an asymptotically
	 *               inertial observer
	 *  @param x_rot absolute X coordinate of the rotation axis
	 *  @param y_rot absolute Y coordinate of the rotation axis
	 *
	 */
	void kinema_bhns(bool kerrschild, const double& mass_bh,
			 const double& sepa, double omega,
			 double x_rot, double y_rot) ;

	/// Computes the gradient of the total velocity potential \f$\psi\f$.
	void fait_d_psi_bhns() ;

	/** Computes \c taij_auto , \c taij_quad_auto  from
	 *  \c shift_auto , \c lapse_auto , \c confo_auto .
	 */
	void extr_curv_bhns() ;

	/** Computes an equilibrium configuration
	 *
	 *  @param ent_c  [input] Central enthalpy
	 *  @param mass_bh [input] BH mass in the background metric
	 *  @param sepa [input] Separation between NS and BH
	 *  @param kerrschild should be \c true  for a Kerr-Schild background,
	 *                    \c false  for a Conformally flat one
	 *  @param mer [input] Number of iteration
	 *  @param mermax_ns [input] Maximum number of iteration steps
	 *  @param mermax_potvit [input] Maximum number of steps in
	 *                               Map_radial::poisson_compact
	 *  @param mermax_poisson [input] Maximum number of steps in
	 *                                poisson scalar
	 *  @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 relax_poisson [input] Relaxation factor in poisson scalar
	 *  @param relax_potvit [input] Relaxation factor in
	 *                              Map_radial::poisson_compact
	 *  @param thres_adapt [input] Threshold on dH/dr for the adaptation
	 *                             of the mapping
	 *  @param resize_ns [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(double ent_c, const double& mass_bh,
			      const double& sepa, bool kerrschild,
			      int mer, int mermax_ns, int mermax_potvit,
			      int mermax_poisson, int filter_r, int filter_r_s,
			      int filter_p_s, double relax_poisson,
			      double relax_potvit, double thres_adapt,
			      double resize_ns,
			      const Tbl& fact_resize, Tbl& diff) ;

	/** Computes the non-translational part of the velocity scalar
	 *  potential \f$\psi0\f$ by solving the continuity equation.
	 *
	 *  @param mass_bh [input] BH mass in the background metric
	 *  @param sepa [input] Separation between NS and BH
	 *  @param kerrschild should be \c true  for a Kerr-Schild background,
	 *                    \c false  for a Conformally flat one
	 *  @param mermax [input] Maximum number of steps in the iteration
	 *  @param precis [input] Required precision: the iteration will
	 *                        be stopped when the relative difference
	 *                        on \f$\psi0\f$ between two successive steps
	 *                        is lower than \c precis .
	 *  @param relax [input] Relaxation factor.
	 *
	 *  @return Relative error of the resolution obtained by comparing
	 *          the operator acting on the solution with the source.
	 */
	double velo_pot_bhns(const double& mass_bh, const double& sepa,
			     bool kerrschild,
			     int mermax, double precis, double relax) ;

	/** Sensitive indicator of the mass-shedding to the direction of
	 *  \f$r\f$, \f$\theta=\pi/2\f$, \f$\phi\f$.
	 *
	 *  @param radius [input] Radial coordinate
	 *  @param phi [input] Azimuthal angle
	 */
	double chi_rp(double radius, double phi) ;

	/** Radius of the star to the direction of
	 *  \f$\theta=\pi/2\f$ and \f$\phi\f$.
	 *
	 *  @param phi [input] Azimuthal angle
	 */
	double radius_p(double phi) ;

	/** Azimuthal angle when the indicator of the mass-shedding
	 *  takes its minimum chi_min
	 */
	double phi_min() ;

	/** Azimuthal angle when the indicator of the mass-shedding
	 *  takes its local minimum
	 *
	 *  @param phi_ini [input] Initial azumuthal angle to search minimum
	 */
	double phi_local_min(double phi_ini) ;

	/** Performs a relaxation on \c ent , \c lapse_auto ,
	 *  \c shift_auto , \c confo_auto .
	 *
	 *  @param star_prev [input] star at the previous step
	 *  @param relax_ent [input] Relaxation factor for \c ent
	 *  @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 Star_bhns& star_prev, double relax_ent,
			double relax_met, int mer, int fmer_met) ;

	/** Computes a spherical configuration
	 *
	 *  @param ent_c  [input] Central enthalpy
	 *  @param precis [input] precision
	 */
	void equil_spher_bhns(double ent_c, double precis) ;

	friend class Bin_bhns ;

};

}
#endif