star.h

/*
 *  Definition of Lorene classes Star, Star_bin and Star_bin_xcts
 *               
 */

/*
 *   Copyright (c) 2010 Michal Bejger
 *   Copyright (c) 2004 Francois Limousin
 *
 *   Copyright (c) 2000-2001 Eric Gourgoulhon (for preceding class Etoile)
 *   Copyright (c) 2000-2001 Keisuke Taniguchi (for preceding class Etoile)
 *
 *   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 as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   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_H_ 
#define __STAR_H_ 

/*
 * $Id: star.h,v 1.32 2014/10/13 08:52:36 j_novak Exp $
 * $Log: star.h,v $
 * Revision 1.32  2014/10/13 08:52:36  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.31  2010/12/09 10:36:42  m_bejger
 * Decouple method removed from declaration in Star_bin_xcts
 *
 * Revision 1.30  2010/10/26 18:47:36  m_bejger
 * Modification of Star_bin_xcts::equilibrium, added table fact_resize
 *
 * Revision 1.29  2010/10/18 19:11:53  m_bejger
 * Changes to Star::equilibrium_spher and Star_bin_xcts::equilibrium to allow for calculations with more than one domain in the star
 *
 * Revision 1.28  2010/06/15 08:17:43  m_bejger
 * Method Star_bin_xcts::set_chi_comp() declared
 *
 * Revision 1.27  2010/05/04 07:53:32  m_bejger
 * Class Star_bin_xcts added (initial version)
 *
 * Revision 1.26  2010/03/29 12:00:25  e_gourgoulhon
 * Removed lnq from documentation
 * (lnq has to be removed from base class Star; it is meaningfull only
 * for Star_bin).
 *
 * Revision 1.25  2010/01/24 16:07:45  e_gourgoulhon
 * New class Star_rot.
 *
 * Revision 1.24  2007/11/06 16:22:03  j_novak
 * The data member stress_euler is now a Sym_tensor instead of a Tensor.
 *
 * Revision 1.23  2007/06/21 19:48:25  k_taniguchi
 * Introduction of a method to compute ray_eq_3pis2.
 *
 * Revision 1.22  2006/05/31 09:25:47  f_limousin
 * Modif. of the size of the different domains
 *
 * Revision 1.21  2006/04/11 14:26:12  f_limousin
 * New version of the code : improvement of the computation of some
 * critical sources, estimation of the dirac gauge, helical symmetry...
 *
 * Revision 1.20  2005/09/15 15:56:28  e_gourgoulhon
 * Made the documentation compliant with Doxygen.
 *
 * Revision 1.19  2005/09/13 19:38:32  f_limousin
 * Reintroduction of the resolution of the equations in cartesian coordinates.
 *
 * Revision 1.18  2005/08/13 16:14:11  m_saijo
 * Corrected the document of the total shift vector
 *
 * Revision 1.17  2005/04/08 12:36:45  f_limousin
 * Just to avoid warnings...
 *
 * Revision 1.16  2005/02/24 16:09:29  f_limousin
 * Change the name of some variables (for instance dcov_logn --> dlogn).
 * Add also member dlnq but delete dlnpsi_auto and dlogn_auto.
 *
 * Revision 1.15  2005/02/17 17:28:18  f_limousin
 * Change the name of some quantities to be consistent with other classes
 * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
 *
 * Revision 1.14  2005/02/11 18:11:16  f_limousin
 * Introduction of a member Map_af in derived class Star_bin.
 *
 * Revision 1.13  2004/11/11 16:29:48  j_novak
 * set_der_0x0 is no longer virtual (to be coherent with Tensor/Scalar classes).
 *
 * Revision 1.12  2004/11/10 16:31:53  j_novak
 * Star is no longer an abstract class (mass_b and mass_g are no longer
 * pure virtual). Modified comments to be readable by doxygen.
 *
 * Revision 1.11  2004/07/21 11:48:30  f_limousin
 * Remove function sprod.
 *
 * Revision 1.10  2004/06/22 12:47:01  f_limousin
 * Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
 *
 * Revision 1.9  2004/05/25 14:48:57  f_limousin
 * Add a parameter for the function equilibrium.
 *
 * Revision 1.8  2004/03/23 09:53:50  f_limousin
 * Minor changes
 *
 * Revision 1.7  2004/02/27 09:41:52  f_limousin
 * Scalars ssjm1_logn, ssjm1_qq ... for all metric coefficients have been
 * in class Star_bin for the resolution of Poisson equations.
 * The class Star is now abstract : the computational routines mass_b()
 * and mass_g() = 0.
 *
 * Revision 1.6  2004/01/22 10:06:33  f_limousin
 * Add methods set_logn_comp() and set_shift_auto().
 *
 * Revision 1.5  2004/01/20 15:26:00  f_limousin
 * New class star and star_bin.
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/star.h,v 1.32 2014/10/13 08:52:36 j_novak Exp $
 *
 */

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

namespace Lorene {
class Eos ;

            //---------------------------//
            //    base class Star        //
            //---------------------------//

class Star {

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

    int nzet ;
    
    const Eos& eos ;   
    
    // Fluid quantities with respect to the fluid frame
    // ------------------------------------------------

    Scalar ent ;       
    
    Scalar nbar ;      
    Scalar ener ;      
    Scalar press ;     

    // Fluid quantities with respect to the Eulerian frame
    // ---------------------------------------------------
    Scalar ener_euler ; 

    Scalar s_euler ;   

    Scalar gam_euler ; 
    
    Vector u_euler ; 
 
    Sym_tensor stress_euler ;
    
    // Metric potentials
    // -----------------
    
    Scalar logn ;

    Scalar nn ; 
    
    Vector beta ;
    
    // Scalar field \f$\ln Q = \ln(\psi^2  N)\f$ 
    //## to be removed from base class Star
    Scalar lnq ;
    
    Metric gamma ;


    // Derived data : 
    // ------------
    protected:
    mutable double* p_ray_eq ; 

    mutable double* p_ray_eq_pis2 ;

    mutable double* p_ray_eq_pi ;
    
    mutable double* p_ray_eq_3pis2 ;

    mutable double* p_ray_pole ;
    
    mutable Itbl* p_l_surf ; 
    
    mutable Tbl* p_xi_surf ; 
    
    mutable double* p_mass_b ;  
    mutable double* p_mass_g ;  


    // Constructors - Destructor
    // -------------------------
    public:
    
    Star(Map& mp_i, int nzet_i, const Eos& eos_i) ;         
    
    
    Star(const Star& ) ;        

    Star(Map& mp_i, const Eos& eos_i, FILE* fich) ;         

    virtual ~Star() ;           

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

    virtual void del_hydro_euler() ; 
    

    // Mutators / assignment
    // ---------------------
    public:
    void operator=(const Star&) ;   
    
    Map& set_mp() {return mp; } ; 

    void set_enthalpy(const Scalar& ) ; 
    
    void equation_of_state() ; 
    
     virtual void hydro_euler() ; 
    
     virtual void equilibrium_spher(double ent_c, double precis = 1.e-14, 
             const Tbl* pent_limit = 0x0 ) ; 

    // Accessors
    // ---------
    public:
    const Map& get_mp() const {return mp; } ; 

    int get_nzet() const {return nzet; } ; 

    const Eos& get_eos() const {return eos; } ; 

    const Scalar& get_ent() const {return ent;} ;

    const Scalar& get_nbar() const {return nbar;} ;

    const Scalar& get_ener() const {return ener;} ;

    const Scalar& get_press() const {return press;} ;

    const Scalar& get_ener_euler() const {return ener_euler;} ;

    const Scalar& get_s_euler() const {return s_euler;} ;

    const Scalar& get_gam_euler() const {return gam_euler;} ;

    const Vector& get_u_euler() const {return u_euler;} ;

    const Tensor& get_stress_euler() const {return stress_euler;} ;

    const Scalar& get_logn() const {return logn;} ;

    const Scalar& get_nn() const {return nn;} ;

    const Vector& get_beta() const {return beta;} ;

    // Returns the scalar field \f$\ln Q\f$.
    //## to be removed from base class Star
    const Scalar& get_lnq() const {return lnq;} ;

    const Metric& get_gamma() const {return gamma;} ;

    // Outputs
    // -------
    public:
    virtual void sauve(FILE* ) const ;      
    
    friend ostream& operator<<(ostream& , const Star& ) ;   

    protected:
    virtual ostream& operator>>(ostream& ) const ;    

    // Global quantities
    // -----------------
    public:
    double ray_eq() const ; 
    
    double ray_eq_pis2() const ; 
    
    double ray_eq_pi() const ; 
    
    double ray_eq_3pis2() const ;

    double ray_pole() const ; 
    
    virtual const Itbl& l_surf() const ; 
    
    const Tbl& xi_surf() const ; 

        virtual double mass_b() const = 0  ;
    
        virtual double mass_g() const = 0  ;
    
};
ostream& operator<<(ostream& , const Star& ) ;  


            //---------------------------//
            //    class Star_bin         //
            //---------------------------//

class Star_bin : public Star {

    // Data : 
    // -----
    protected:
    bool irrotational ; 
    
    Scalar psi0 ; 

    Vector d_psi ; 
    
    Vector wit_w ; 
    
    Scalar loggam ; 

    Vector bsn ; 
    
    Scalar pot_centri ;     


    Scalar logn_auto ; 

    Scalar logn_comp ; 

    Vector dcov_logn ;

    Vector dcon_logn ;
    
    Scalar lnq_auto ;

    Scalar lnq_comp ;

    
    Scalar psi4 ;

    Vector dcov_phi ;
    Vector dcon_phi ;

    Metric_flat flat ;

    Metric gtilde ;

    Vector beta_auto ; 

    Vector beta_comp ; 


    Sym_tensor hij ; 


    Sym_tensor hij_auto ;

    Sym_tensor hij_comp ;

    Sym_tensor tkij_auto ;
    
    Sym_tensor tkij_comp ;
    
    Scalar kcar_auto ;
    
    Scalar kcar_comp ;

    Scalar ssjm1_logn ;

    Scalar ssjm1_lnq ;

    Scalar ssjm1_khi ;

    Vector ssjm1_wbeta ;
    
    Scalar ssjm1_h11 ;

    Scalar ssjm1_h21 ;

    Scalar ssjm1_h31 ;

    Scalar ssjm1_h22 ;

    Scalar ssjm1_h32 ;

    Scalar ssjm1_h33 ;

     Scalar decouple ;

     bool conf_flat ;
    
    // Derived data : 
    // ------------
    protected:
    mutable double* p_xa_barycenter ; 
     

    // Constructors - Destructor
    // -------------------------
    public:
    Star_bin(Map& mp_i, int nzet_i,  const Eos& eos_i,
           bool irrot, bool conf_flat) ;            
    
    
    Star_bin(const Star_bin& ) ;        

    Star_bin(Map& mp_i, const Eos& eos_i, FILE* fich) ;
            
    virtual ~Star_bin() ;           


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

    virtual void del_hydro_euler() ; 
    

    // Mutators / assignment
    // ---------------------
    public:
    void operator=(const Star_bin& ) ;  
    
    Scalar& set_pot_centri() ;
    
    Scalar& set_logn_comp() ;

    void set_logn_auto(const Scalar& logn_auto_new) {logn_auto = logn_auto_new ;
    return ;}

    void set_lnq_auto(const Scalar& lnq_auto_new) {lnq_auto = lnq_auto_new ;
    return ;}

    Vector& set_beta_auto() ;   

    Vector& set_beta() ;    
    
    void set_conf_flat(bool confflat) {conf_flat = confflat ; return ;}

    // Accessors
    // ---------
    public:
    bool is_irrotational() const {return irrotational; } ; 

    const Scalar& get_psi0() const {return psi0;} ;

    const Vector& get_d_psi() const {return d_psi;} ;

    const Vector& get_wit_w() const {return wit_w;} ;

    const Scalar& get_loggam() const {return loggam;} ;

    const Vector& get_bsn() const {return bsn;} ;

    const Scalar& get_pot_centri() const {return pot_centri;} ;

    const Scalar& get_logn_auto() const {return logn_auto;} ;

    const Scalar& get_logn_comp() const {return logn_comp;} ;

    const Vector& get_beta_auto() const {return beta_auto;} ;

    const Vector& get_beta_comp() const {return beta_comp;} ;

    const Scalar& get_lnq_auto() const {return lnq_auto;} ;

    const Scalar& get_lnq_comp() const {return lnq_comp;} ;

    const Vector& get_dcov_logn() const {return dcov_logn;} ;

    const Vector& get_dcon_logn() const {return dcon_logn;} ;

    const Vector& get_dcov_phi() const {return dcov_phi;} ;

    const Vector& get_dcon_phi() const {return dcon_phi;} ;

    const Scalar& get_psi4() const {return psi4;} ;

    const Metric& get_flat() const {return flat;} ; 

    const Metric& get_gtilde() const {return gtilde;} ;
    
    const Sym_tensor& get_hij() const {return hij;} ;

    const Sym_tensor& get_hij_auto() const {return hij_auto;} ;

    const Sym_tensor& get_hij_comp() const {return hij_comp;} ;


    const Sym_tensor& get_tkij_auto() const {return tkij_auto;} ;

    const Sym_tensor& get_tkij_comp() const {return tkij_comp;} ;

    const Scalar& get_kcar_auto() const {return kcar_auto;} ;

    const Scalar& get_kcar_comp() const {return kcar_comp;} ;

    const Scalar get_decouple() const {return decouple ;}

    bool is_conf_flat() const {return conf_flat; } ; 

    // Outputs
    // -------
    public:
    virtual void sauve(FILE* ) const ;      
    
    protected:
    virtual ostream& operator>>(ostream& ) const ;    

    // Global quantities
    // -----------------
    public:
        virtual double mass_b() const ;
    
        virtual double mass_g() const ;
    
        virtual double xa_barycenter() const ;


    // Computational routines
    // ----------------------
    public: 

    virtual void hydro_euler() ; 
    
    void update_metric(const Star_bin& comp, double omega) ;
    
    void update_metric(const Star_bin& comp, const Star_bin& star_prev, 
               double relax, double omega) ; 
    
     void update_metric_der_comp(const Star_bin& comp, double omega) ;

    void kinematics(double omega, double x_axe) ; 
    
    void fait_d_psi() ; 
    
    void extrinsic_curvature(double omega) ; 
        
    
    void equilibrium(double ent_c, int mermax, int mermax_potvit, 
             int mermax_poisson, double relax_poisson, 
             double relax_potvit, double thres_adapt, Tbl& diff,
             double om) ;


     double velocity_potential(int mermax, double precis, double relax) ;

    void relaxation(const Star_bin& star_prev, double relax_ent, 
            double relax_met, int mer, int fmer_met) ;
    

    void test_K_Hi() const ;

    void helical(double omega) const ;

    friend class Binary ;

};

            //---------------------------//
            //    class Star_bin_xcts    //
            //---------------------------//

class Star_bin_xcts : public Star {

    // Data : 
    // -----
    protected:
    bool irrotational ; 
    
    Scalar psi0 ; 

    Vector d_psi ; 
    
    Vector wit_w ; 
    
    Scalar loggam ; 

    Vector bsn ; 
    
    Scalar pot_centri ;     

    Scalar chi_auto ;

    Scalar chi_comp ;

    Scalar Psi_auto ;

    Scalar Psi_comp ;

    Scalar Psi ;
    
    Scalar chi ;
    
    Scalar psi4 ;
        
    Vector w_beta ; 

    Scalar khi ;    

    Vector dcov_Psi ;

    Vector dcov_chi ;

    Metric_flat flat ;
    
    Vector beta_auto ; 
    
    Vector beta_comp ; 

    Sym_tensor haij_auto ;
    
    Sym_tensor haij_comp ;
    
    Scalar hacar_auto ;
    
    Scalar hacar_comp ;

    Scalar ssjm1_chi ;

    Scalar ssjm1_psi ;
        
    Scalar ssjm1_khi ;

    Vector ssjm1_wbeta ;
     
    // Derived data : 
    // ------------
    protected:
    mutable double* p_xa_barycenter ; 

    // Constructors - Destructor
    // -------------------------
    public:
    Star_bin_xcts(Map& mp_i, int nzet_i, const Eos& eos_i, bool irrot) ;            
        
    Star_bin_xcts(const Star_bin_xcts& ) ;      

    Star_bin_xcts(Map& mp_i, const Eos& eos_i, FILE* fich) ;
            
    virtual ~Star_bin_xcts() ;          


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

    virtual void del_hydro_euler() ; 
    

    // Mutators / assignment
    // ---------------------
    public:
    void operator=(const Star_bin_xcts& ) ; 
    
    Scalar& set_pot_centri() ;
    
    Scalar& set_Psi_auto() ; 

    Scalar& set_Psi_comp() ; 

    Scalar& set_chi_auto() ; 

    Scalar& set_chi_comp() ; 
        
    Vector& set_beta_auto() ;   

    Vector& set_beta() ;    
    
    // Accessors
    // ---------
    public:
    bool is_irrotational() const { return irrotational; } ; 

    const Scalar& get_psi0() const {return psi0 ; } ;

    const Vector& get_d_psi() const {return d_psi ; } ;

    const Vector& get_wit_w() const {return wit_w ; } ;

    const Scalar& get_loggam() const {return loggam ; } ;

    const Vector& get_bsn() const {return bsn ; } ;

    const Scalar& get_pot_centri() const {return pot_centri ; } ;

    const Vector& get_beta_auto() const {return beta_auto ; } ;

    const Vector& get_beta_comp() const {return beta_comp ; } ;
    
    const Scalar& get_Psi_auto() const {return Psi_auto ; } ;

    const Scalar& get_Psi_comp() const {return Psi_comp ; } ;

    const Scalar& get_chi_auto() const {return chi_auto ; } ;

    const Scalar& get_chi_comp() const {return chi_comp ; } ;
        
    const Vector& get_dcov_chi() const {return dcov_chi ; } ;

    const Vector& get_dcov_Psi() const {return dcov_Psi ; } ;

    const Scalar& get_Psi() const {return Psi ; } ;

    const Scalar& get_chi() const {return chi ; } ;
    
    const Scalar& get_psi4() const {return psi4 ; } ;
    
    const Metric& get_flat() const {return flat ; } ;   

    const Sym_tensor& get_haij_auto() const {return haij_auto ; } ;

    const Sym_tensor& get_haij_comp() const {return haij_comp ; } ;

    const Scalar& get_hacar_auto() const {return hacar_auto ; } ;

    const Scalar& get_hacar_comp() const {return hacar_comp ; } ;
 
    // Outputs
    // -------
    public:
    virtual void sauve(FILE* ) const ;      
    
    protected:
    virtual ostream& operator>>(ostream& ) const ;    

    // Global quantities
    // -----------------
    public:
        virtual double mass_b() const ;
    
        virtual double mass_g() const ;
    
        virtual double xa_barycenter() const ;


    // Computational routines
    // ----------------------
    public: 

    virtual void hydro_euler() ; 
    
    void update_metric(const Star_bin_xcts& comp) ;
    
    void update_metric(const Star_bin_xcts& comp,
                       const Star_bin_xcts& star_prev, 
                       double relax) ; 
    
     void update_metric_der_comp(const Star_bin_xcts& comp) ;

    void kinematics(double omega, double x_axe) ; 
    
    void fait_d_psi() ; 
    
    void extrinsic_curvature() ; 
    
    void equilibrium(double ent_c, int mermax, int mermax_potvit, 
             int mermax_poisson, double relax_poisson, 
             double relax_potvit, double thres_adapt, 
             const Tbl& fact, const Tbl* pent_limit, 
             Tbl& diff) ;

     double velocity_potential(int mermax, double precis, double relax) ;

    void relaxation(const Star_bin_xcts& star_prev, double relax_ent, 
            double relax_met, int mer, int fmer_met) ;

    friend class Binary_xcts ;

}; 

}
#endif