isol_hor.h

/*
 *  Definition of Lorene class Isol_Hor
 *
 */

/*
 *   Copyright (c) 2004-2005 Jose Luis Jaramillo
 *                      Francois Limousin
 *
 *   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 __Isol_hor_H_ 
#define __Isol_hor_H_ 

/*
 * $Id: isol_hor.h,v 1.54 2014/10/13 08:52:35 j_novak Exp $
 * $Log: isol_hor.h,v $
 * Revision 1.54  2014/10/13 08:52:35  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.53  2007/08/22 16:10:37  f_limousin
 * Correction of many errors in binhor_hh.C
 *
 * Revision 1.52  2007/04/13 15:29:44  f_limousin
 * Lots of improvements, generalisation to an arbitrary state of
 * rotation, implementation of the spatial metric given by Samaya.
 *
 * Revision 1.51  2006/08/01 14:36:25  f_limousin
 * New argument for the functions vv_bound_cart_bin( )
 *
 * Revision 1.50  2006/06/29 08:58:57  f_limousin
 * Boundary conditions in argument of write_global()
 *
 * Revision 1.49  2006/06/28 13:36:26  f_limousin
 * Convergence to a given irreductible mass
 *
 * Revision 1.48  2006/05/24 16:53:51  f_limousin
 * Funtion import(Bin_hor&)
 *
 * Revision 1.47  2006/02/22 16:55:42  f_limousin
 * boundary_nn_Neu_kk(int nn = 1)
 *
 * Revision 1.46  2006/02/22 16:20:54  jl_jaramillo
 * 2nc correction Valeur boundary_nn_Neu_kk(int nn)
 *
 * Revision 1.45  2006/02/20 16:48:14  jl_jaramillo
 * corrections numerical viscosity
 *
 * Revision 1.44  2006/02/20 15:01:37  jl_jaramillo
 * function for computing generalised Conformal THin Sandwich eqautions
 *
 * Revision 1.43  2006/01/18 12:29:18  jl_jaramillo
 * new method (chi,theta) for the spherically symmetric case
 *
 * Revision 1.42  2006/01/16 17:15:34  jl_jaramillo
 * function for solving the spherical case
 *
 * Revision 1.41  2005/11/02 16:10:38  jl_jaramillo
 * change in boundary_nn_Dir_lapl
 *
 * Revision 1.40  2005/10/24 16:45:16  jl_jaramillo
 * Cook boundary condition
 *
 * Revision 1.39  2005/10/23 12:32:00  f_limousin
 * Pure Dirichlet boundary condition for Psi.
 *
 * Revision 1.38  2005/09/13 18:33:17  f_limousin
 * New function vv_bound_cart_bin(double) for computing binaries with
 * berlin condition for the shift vector.
 * Suppress all the symy and asymy in the importations.
 *
 * Revision 1.37  2005/07/11 08:19:54  f_limousin
 * New function axi_break() to compute the departure to axisymmetry.
 *
 * Revision 1.36  2005/06/09 07:56:24  f_limousin
 * Implement a new function sol_elliptic_boundary() and
 * Vector::poisson_boundary(...) which solve the vectorial poisson
 * equation (method 6) with an inner boundary condition.
 *
 * Revision 1.35  2005/05/12 14:48:55  f_limousin
 * New boundary condition for the lapse : boundary_nn_lapl().
 *
 * Revision 1.34  2005/04/29 14:05:03  f_limousin
 * Important changes : manage the dependances between quantities (for
 * instance psi and psi4). New function write_global(ost).
 *
 * Revision 1.33  2005/04/15 11:20:39  jl_jaramillo
 * Function adapt_hor(double c_min, double c_max)  for adapting a given surface
 * to the excised surface
 *
 * Revision 1.32  2005/04/08 12:15:38  f_limousin
 * Function set_psi(). And dependance in phi.
 *
 * Revision 1.31  2005/04/03 19:48:38  f_limousin
 * Implementation of set_psi(psi_in).
 *
 * Revision 1.30  2005/04/02 15:50:08  f_limousin
 * New data member nz (number of zones). Delete ww.
 *
 * Revision 1.29  2005/03/31 09:48:04  f_limousin
 * New functions compute_ww(..) and aa_kerr_ww() and new data ww.
 *
 * Revision 1.28  2005/03/28 19:45:41  f_limousin
 * Implement Isol_hor::aa_kerr_perturb(...) and new member aa_quad_evol.
 *
 * Revision 1.27  2005/03/24 16:50:40  f_limousin
 * Add parameters solve_shift and solve_psi in par_isol.d and in function
 * init_dat(...). Implement Isolhor::kerr_perturb().
 *
 * Revision 1.26  2005/03/10 16:57:01  f_limousin
 * Improve the convergence of the code coal_bh.
 *
 * Revision 1.25  2005/03/10 10:19:42  f_limousin
 * Add the regularisation of the shift in the case of a single black hole
 * and lapse zero on the horizon.
 *
 * Revision 1.24  2005/03/09 10:28:37  f_limousin
 * Delete functions init_data_b_neumann(...) and init_data_berlin(...)
 * --> New parameter solve_lapse in the function init_data(...).
 * New function update_aa().
 *
 * Revision 1.23  2005/03/06 16:59:33  f_limousin
 * New function Isol_hor::aa() (the one belonging to the class
 * Time_slice_conf need to compute the time derivative of hh and thus
 * cannot work in the class Isol_hor).
 *
 * Revision 1.22  2005/03/04 09:39:31  f_limousin
 * Implement the constructor from a file, operator>>, operator<<
 * and function sauve in the class Bin_hor.
 *
 * Revision 1.21  2005/03/03 10:25:16  f_limousin
 * In the class Isol_hor :
 *    - Add the boost in x and z-direction (members boost_x and boost_z,
 *    and functions get_boost_x(), set_boost_x(double))
 *    - Add function area_hor()
 *    - Put the boundary conditions for the lapse, psi and beta in
 *    the parameter file.
 * In the class bin_hor :
 *    -  Introduce function to compute global quantities as ADM mass,
 *    Komar mass and angular momentum.
 *
 * Revision 1.20  2005/02/24 17:22:53  f_limousin
 * Suppression of the function beta_bound_cart().
 * The boundary conditions for psi, N and beta are now some parameters
 * in par_init.D and par_coal.d.
 *
 * Revision 1.19  2005/02/07 10:30:09  f_limousin
 * Add the regularisation in the case N=0 on the horizon.
 *
 * Revision 1.18  2004/12/31 15:33:37  f_limousin
 * Change the constructor from a file and the standard constructor.
 *
 * Revision 1.17  2004/12/29 16:30:00  f_limousin
 * Improve comments for doxygen
 *
 * Revision 1.16  2004/12/29 16:10:25  f_limousin
 * Add the new class Bin_hor.
 *
 * Revision 1.14  2004/11/24 19:32:05  jl_jaramillo
 * Method for initial data with Berlin boundary conditions
 *
 * Revision 1.13  2004/11/18 10:53:03  jl_jaramillo
 * Declarations for Berlin boundary conditions
 *
 * Revision 1.12  2004/11/05 11:01:13  f_limousin
 * Delete arguments ener_dens, mom_dens and trace_stress in all functions
 * source_nn, source_psi, source_beta, init_data. Delete also
 * argument partial_save in function save.
 *
 * Revision 1.11  2004/11/05 10:11:23  f_limousin
 * The member Metric met_gamt replace Sym_tensor gamt.
 *
 * Revision 1.10  2004/11/03 17:15:46  f_limousin
 * Change the standart constructor. Add 4 memebers : trK, trK_point,
 * gamt and gamt_point.
 * Add also a constructor from a file.
 *
 * Revision 1.9  2004/11/02 17:42:33  f_limousin
 * New method sauve(...) to save in a binary file.
 *
 * Revision 1.8  2004/11/02 16:15:12  f_limousin
 * Add new argument ang_vel in function init_dat(...).
 *
 * Revision 1.7  2004/10/29 15:46:14  jl_jaramillo
 * Remove 2 members, add ADM angular momentum and change name
 * of functions.
 *
 * Revision 1.6  2004/10/01 16:51:16  f_limousin
 * Pure Dirichlet boundary condition added
 *
 * Revision 1.5  2004/09/28 16:03:58  f_limousin
 * Add parameter niter in the parameter file par_hor.d. It appears in
 * argument of the function init_data_schwarz(...).
 *
 * Revision 1.4  2004/09/17 13:35:25  f_limousin
 * Introduction of relaxation in init_data_schwarz
 *
 * Revision 1.3  2004/09/16 08:36:57  f_limousin
 * New boundary conditions for lapse and psi.
 *
 * Revision 1.2  2004/09/09 17:04:27  jl_jaramillo
 * Elimination of _ih
 *
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/isol_hor.h,v 1.54 2014/10/13 08:52:35 j_novak Exp $
 *
 */

#include "time_slice.h"
#include "proto.h"
#include "headcpp.h"
#include "cmp.h"
#include "evolution.h"

namespace Lorene {
class Sym_tensor_trans ; 
class Sym_tensor ; 
class Vector ; 
class Scalar ; 
class Metric ; 
class Metric_flat ; 
class Base_vect ; 
class Map ; 
class Tbl ;
class Time_slice ;
class Time_slice_conf ;
                    
//----------------------------//
//       class Isol_Hor       //
//----------------------------//

class Isol_hor : public Time_slice_conf {

  // Data : 
  // -----
 protected: 
  Map_af& mp ;  

  int nz ;

  double radius ; 

  double omega ; 
  
  double boost_x ;
 
  double boost_z ;

  double regul ; 

  mutable Evolution_std<Scalar> n_auto_evol ; 

  mutable Evolution_std<Scalar> n_comp_evol ; 

  mutable Evolution_std<Scalar> psi_auto_evol ; 

  mutable Evolution_std<Scalar> psi_comp_evol ; 

  mutable Evolution_std<Vector> dn_evol ; 

  mutable Evolution_std<Vector> dpsi_evol ;   

  mutable Evolution_std<Vector> beta_auto_evol ; 

  mutable Evolution_std<Vector> beta_comp_evol ; 

  mutable Evolution_std<Sym_tensor> aa_auto_evol ; 

  mutable Evolution_std<Sym_tensor> aa_comp_evol ; 
  
  mutable Evolution_std<Sym_tensor> aa_nn ; 

  mutable Evolution_std<Scalar> aa_quad_evol ;

  Metric met_gamt ;
  
  Sym_tensor gamt_point ;
  
  Scalar trK ;
    
  Scalar trK_point ;

  Scalar decouple ;

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

  Isol_hor(Map_af& mpi, int depth_in = 3) ;

  Isol_hor(Map_af& mpi, const Scalar& lapse_in, const Scalar& psi_in,
       const Vector& shift_in, const Sym_tensor& aa_in, 
       const Metric& gamt, const Sym_tensor& gamt_point, 
       const Scalar& trK, const Scalar& trK_point, 
       const Metric_flat& ff_in, int depth_in = 3) ;    
  
  Isol_hor(const Isol_hor& ) ;  

  Isol_hor (Map_af& mp, FILE* fich, 
        bool partial_read, int depth_in = 3) ;  
  
  virtual ~Isol_hor() ;         
  

  // Mutators / assignment
  // ---------------------
 public:
  void operator=(const Isol_hor&) ; 
    

 public:
  const Map_af& get_mp() const {return mp;} ; 
  
  Map_af& set_mp() {return mp; } ;

  double get_radius() const {return radius;} ;
  
  void set_radius(double rad) {radius = rad ;} ;
    
  double get_omega() const {return omega ;} ;
  void set_omega(double ome) {omega = ome ;} ;

  double get_boost_x() const {return boost_x ;} ;
  void set_boost_x(double bo) {boost_x = bo ;} ;

  double get_boost_z() const {return boost_z ;} ;
  void set_boost_z(double bo) {boost_z = bo ;} ;


  
  // Accessors
  // ---------
 public:

  virtual const Scalar& n_auto() const ;

  virtual const Scalar& n_comp() const ;

  virtual const Scalar& psi_auto() const ;

  virtual const Scalar& psi_comp() const ;

  virtual const Vector& dnn() const ;

  virtual const Vector& dpsi() const ;

  virtual const Vector& beta_auto() const ;

  virtual const Vector& beta_comp() const ;

  virtual const Sym_tensor& aa_auto() const ; 

  virtual const Sym_tensor& aa_comp() const ; 

  virtual const Scalar& aa_quad() const ; 

  virtual const Metric& tgam() const {return met_gamt ;}

  const Scalar get_decouple() const {return decouple ;}


 public:
  void n_comp (const Isol_hor& comp) ;
  
  void psi_comp (const Isol_hor& comp) ;

  void beta_comp (const Isol_hor& comp) ;

  double viriel_seul () const ;

  void init_bhole () ;
  
  void init_met_trK() ;

  void init_bhole_seul () ;

  void set_psi(const Scalar& psi_in) ; 

  void set_nn(const Scalar& nn_in) ; 
  
  void set_gamt(const Metric& gam_tilde) ;

  // Physical parameters
  //--------------------
 public:
 
  
  const Vector radial_vect_hor() const ;

  const Vector tradial_vect_hor() const ;

  const Scalar b_tilde() const ;

  const Scalar darea_hor() const ;

  double area_hor() const ;

  double radius_hor() const ;
  
  double ang_mom_hor() const ;
  
  double mass_hor() const ;
  
  double kappa_hor() const ;
  
  double omega_hor() const ;
  
  double ang_mom_adm() const ;

  Scalar expansion() const ;


  //Computational methods
  //---------------------
 public:
  
  /* function to compute initial data for a single black hole
   *  @param bound_nn boundary condition for the lapse
   *  @param lim_nn value of the boundary condition for the lapse 
   *  @param bound_psi boundary condition for \f$ \Psi \f$
   *  @param bound_beta boundary condition for the shift
   *  @param solve_lapse do we solve the equation for the lapse ?
   *  @param precis precision for the convergence
   *  @param relax relaxation
   *  @param niter number of iterations
   */
  void init_data(int bound_nn, double lim_nn, int bound_psi, int bound_beta,
         int solve_lapse, int solve_psi, int solve_shift, 
         double precis = 1.e-12,
         double relax_nn = 0.5, double relax_psi = 0.5,  
         double relax_beta = 0.5, int niter = 100) ; 

  void init_data_loop(int bound_nn, double lim_nn, int bound_psi, 
              int bound_beta, int solve_lapse, int solve_psi,
              int solve_shift, double precis= 1.e-12, 
              double precis_loop= 1.e-12,
              double relax_nn = 1., double relax_psi= 1.,  
              double relax_beta = 1., double relax_loop = 1., 
              int niter = 100) ;



  void init_data_spher(int bound_nn, double lim_nn, int bound_psi, 
               int bound_beta, int solve_lapse, int solve_psi, 
               int solve_shift, double precis = 1.e-12, 
               double relax = 1., int niter = 100) ; 

  void init_data_alt(int bound_nn, double lim_nn, int bound_psi, 
               int bound_beta, int solve_lapse, int solve_psi, 
               int solve_shift, double precis = 1.e-12, 
               double relax = 1., int niter = 100) ; 

  void init_data_CTS_gen(int bound_nn, double lim_nn, int bound_psi, int bound_beta,
         int solve_lapse, int solve_psi, int solve_shift, 
         double precis = 1.e-12, double relax_nn = 1.,
         double relax_psi = 1., double relax_beta = 1., 
         int niter = 100, double a = 1., double zeta = 4.) ; 




  //Sources
  //-------
  
  const Scalar source_psi() const ;

  const Scalar source_nn() const ;

  const Vector source_beta() const ;

  const Scalar source_b_tilde() const ;

  const Vector source_vector_b() const ;
  
  
  // BOUNDARY CONDITIONS
  //--------------------
  
  const Valeur boundary_psi_Dir_evol() const ;

  const Valeur boundary_psi_Neu_evol() const ;

  const Valeur boundary_psi_Dir_spat() const ;

  const Valeur boundary_psi_Neu_spat() const ;

  const Valeur boundary_psi_app_hor() const ;

  const Valeur boundary_psi_Dir() const ;

  const Valeur boundary_nn_Dir_kk() const ;

  const Valeur boundary_nn_Neu_kk(int nn = 1) const ;   

  const Valeur boundary_nn_Neu_Cook() const ;   

  const Valeur boundary_nn_Dir_eff(double aa) const ;

  const Valeur boundary_nn_Dir_lapl(int mer = 1) const ;

  const Valeur boundary_nn_Neu_eff(double aa) const ;   

  const Valeur boundary_nn_Dir(double aa) const ;

  const Valeur boundary_beta_r() const ;
  
  const Valeur boundary_beta_theta() const ;
  
  const Valeur boundary_beta_phi(double om) const ;
  
  const Valeur boundary_beta_x(double om) const ;
  
  const Valeur boundary_beta_y(double om) const ;
  
  const Valeur boundary_beta_z() const ;

  const Valeur beta_boost_x() const ;  

  const Valeur beta_boost_z() const ;  

  const Vector vv_bound_cart(double om) const ;

  const Vector vv_bound_cart_bin(double om, int hole = 0) const ;

  const Valeur boundary_vv_x(double om) const ;

  const Valeur boundary_vv_y(double om) const ;
  
  const Valeur boundary_vv_z(double om) const ;

  const Valeur boundary_vv_x_bin(double om, int hole = 0) const ;

  const Valeur boundary_vv_y_bin(double om, int hole = 0) const ;

  const Valeur boundary_vv_z_bin(double om, int hole = 0) const ;

  const Valeur boundary_b_tilde_Neu() const ;

  const Valeur boundary_b_tilde_Dir() const ;
 
  void update_aa() ; 
  
  double regularisation (const Vector& shift_auto, const Vector& shift_comp, 
             double ang_vel) ;

  double regularise_one() ;

  void met_kerr_perturb() ;
  
  /* Perturbation of Kerr using  \f$ A^{ij}A_{ij} \f$ from 
   * equation (14) of Dain (2002).
   * @param mm mass of the Kerr black hole metric.
   * @param aa rotation parameter of the Kerr black hole metric.
   */
  void aa_kerr_ww(double mm, double aa) ;
  
    
  double axi_break() const ;

  /*  Calculation of the outermost trapped surface and adaptation
   *  of all necessary quantities
   */

  void adapt_hor(double c_min, double c_max) ;
  

  // Outputs
  // -------
 protected:
  virtual ostream& operator>>(ostream& ) const ;    
  

  public :
    virtual void sauve(FILE* fich, bool partial_save) const ; 
    
  friend class Bin_hor ; 

};

class Single_hor {

  // Data : 
  // -----
 protected: 
  Map_af& mp ;  

  int nz ;

  double radius ; 

  double omega ; 
  
  double regul ; 

  Scalar n_auto ; 

  Scalar n_comp ; 

  Scalar nn ; 

  Scalar psi_auto ; 

  Scalar psi_comp ; 

  Scalar psi ; 

  mutable Scalar* p_psi4 ; 

  Vector dn ; 

  Vector dpsi ;   

  Vector beta_auto ; 

  Vector beta_comp ; 

  Vector beta ; 

  mutable Metric* p_gam ;

  Sym_tensor aa_auto ; 

  Sym_tensor aa_comp ;
 
  Sym_tensor aa ; 

  mutable Sym_tensor* p_k_dd ; 

  Metric tgam ;
  
  Metric_flat ff ;
  
  Sym_tensor hh ;
  
  Sym_tensor gamt_point ;
  
  Scalar trK ;
    
  Scalar trK_point ;

  Scalar decouple ;

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

  Single_hor(Map_af& mpi) ;
 
  Single_hor(const Single_hor& ) ;  

  Single_hor (Map_af& mp, FILE* fich) ;  
  
  virtual ~Single_hor() ;           
  

  // Mutators / assignment
  // ---------------------
 public:
  void operator=(const Single_hor&) ;   
    

 public:
  const Map_af& get_mp() const {return mp;} ; 
  
  Map_af& set_mp() {return mp; } ;

  double get_radius() const {return radius;} ;
  
  void set_radius(double rad) {radius = rad ;} ;
    
  double get_omega() const {return omega ;} ;
  void set_omega(double ome) {omega = ome ;} ;

  // Memory management
  // -----------------
 protected:
  
  void del_deriv() const ; 
  
  void set_der_0x0() const ; 
  

  
  // Accessors
  // ---------
 public:

  const Scalar& get_n_auto() const ;

  const Scalar& get_n_comp() const ;

  const Scalar& get_nn() const ;

  const Scalar& get_psi_auto() const ;

  const Scalar& get_psi_comp() const ;

  const Scalar& get_psi() const ;

  const Scalar& get_psi4() const ;

  const Vector& get_dn() const ;

  const Vector& get_dpsi() const ;

  const Vector& get_beta_auto() const ;

  const Vector& get_beta_comp() const ;

  const Vector& get_beta() const ;

  const Sym_tensor& get_aa_auto() const ; 

  const Sym_tensor& get_aa_comp() const ; 

  const Sym_tensor& get_aa() const ; 

  const Metric& get_tgam() const {return tgam ;}

  const Metric& get_gam() const ;

  const Sym_tensor& get_k_dd() const ;

  const Scalar get_decouple() const {return decouple ;}


 public:
  void n_comp_import (const Single_hor& comp) ;
  
  void psi_comp_import (const Single_hor& comp) ;

  void beta_comp_import (const Single_hor& comp) ;

  double viriel_seul () const ;

  void init_bhole () ;
  
  void init_met_trK() ;

  void init_bhole_seul () ;

  void set_psi_auto(const Scalar& psi_in) ; 

  void set_n_auto(const Scalar& nn_in) ; 
  
  void set_beta_auto(const Scalar& shift_in) ; 

  void set_aa_auto(const Scalar& aa_auto_in) ; 

  void set_aa_comp(const Scalar& aa_comp_in) ; 

  void set_aa(const Scalar& aa_in) ; 


  // Physical parameters
  //--------------------
 public:
 

  const Scalar b_tilde() const ;

  const Scalar darea_hor() const ;

  double area_hor() const ;

  double radius_hor() const ;
  
  double ang_mom_hor() const ;
  
  double mass_hor() const ;
  
  double kappa_hor() const ;
  
  double omega_hor() const ;
  
  double ang_mom_adm() const ;

  Scalar expansion() const ;

  
  // BOUNDARY CONDITIONS
  //--------------------

  const Valeur boundary_psi_app_hor() const ;

  const Valeur boundary_nn_Dir(double aa) const ;

  const Valeur boundary_nn_Neu(double aa) const ;   
  
  const Valeur boundary_beta_x(double om_orb, double om_loc) const ;
  
  const Valeur boundary_beta_y(double om_orb, double om_loc) const ;
  
  const Valeur boundary_beta_z() const ;

  
  double regularisation (const Vector& shift_auto, const Vector& shift_comp, 
             double ang_vel) ;

  double regularise_one() ;  

  public :
    virtual void sauve(FILE* fich) const ; 
    
  friend class Bin_hor ; 

};

class Bin_hor {
    
    // data :
    private:
    // lThe two black holes
    Single_hor hole1 ;  
    Single_hor hole2 ;  
    
    Single_hor* holes[2] ; 
    
    double omega ;  
    
    public:

    Bin_hor(Map_af& mp1, Map_af& mp2) ; 

    Bin_hor(const Bin_hor& ) ;  

    Bin_hor (Map_af& mp1, Map_af& mp2, FILE* fich) ;  
    
    virtual ~Bin_hor() ;  
    
     public :
      void sauve(FILE* fich) const ; 
      
    void write_global(ostream&, double lim_nn, int bound_nn, 
              int bound_psi, int bound_beta, double alpha) const  ;

       public:
    
    void operator=(const Bin_hor&) ; 
    
    Single_hor& set(int i) 
        { assert( (i==1) || (i==2) ); 
          return *holes[i-1] ;} ; 
    void set_omega(double ome) {omega = ome ; 
                 hole1.set_omega (ome) ;
                 hole2.set_omega (ome) ;} ;
    
    public: const Single_hor& operator()(int i) const 
        { assert( (i==1) || (i==2) ); 
          return *holes[i-1] ;} ;
          
    double get_omega() const {return omega; } ; 
    
    void init_bin_hor() ;
    
    double viriel() const ;

    void extrinsic_curvature () ;
    
    void decouple () ;
    
    public:
    void set_statiques (double precis, double relax, int bound_nn,
                double lim_nn, int bound_psi) ;
     
    double coal (double ang_vel, double relax, int nb_om,
             int nb_it, int bound_nn, double lim_nn, 
             int bound_psi, int bound_beta, double omega_eff, 
             double alpha,
             ostream& fich_iteration, ostream& fich_correction,
             ostream& fich_viriel, ostream& fich_kss, 
             int step, int search_mass, double mass_irr, 
             const int sortie = 0) ;

      
      void solve_lapse (double precis, double relax, int bound_nn,
                double lim_nn) ;
      
      void solve_psi (double precis, double relax, int bound_psi) ;
      
      void solve_shift (double precis, double relax, int bound_beta,
                double omega_eff) ;
    
      void import_bh (const Bin_hor& bin) ;

      double adm_mass() const ;
      
      double komar_mass() const ;
      
      double ang_mom_hor() const ;

      double ang_mom_adm() const ;

      double proper_distance(const int nr = 65) const ;

      Sym_tensor hh_Samaya_hole1() ;
      
      Sym_tensor hh_Samaya_hole2() ;

      void set_hh_Samaya() ;



} ;


}
#endif