star_QI.C

/*
 *  Methods of the class Star_QI
 *
 *    (see file compobj.h for documentation).
 *
 */
 
/*
 *   Copyright (c) 2012 Claire Some, Eric Gourgoulhon
 *
 *   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
 *
 */
 
char star_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $" ;
 
/*
 * $Id: star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
 * $Log: star_QI.C,v $
 * Revision 1.5  2014/10/13 08:52:49  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.4  2012/12/03 15:27:30  c_some
 * Small changes
 *
 * Revision 1.3  2012/11/22 16:04:51  c_some
 * Minor modifications
 *
 * Revision 1.2  2012/11/21 14:55:27  c_some
 * Added methods fait_shift and fait_nphi
 *
 * Revision 1.1  2012/11/20 16:29:50  c_some
 * New class Star_QI
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
 *
 */
 
 
// C headers
#include <cassert>
 
// Lorene headers
#include "compobj.h"
#include "unites.h"
 
                   //--------------//
                   // Constructors //
                   //--------------//
 
// Standard constructor
// --------------------
namespace Lorene {
Star_QI::Star_QI(Map& mpi) :
             Compobj_QI(mpi) ,
             logn(mpi), 
             tnphi(mpi), 
             nuf(mpi), 
             nuq(mpi), 
             dzeta(mpi), 
             tggg(mpi), 
             w_shift(mpi, CON, mp.get_bvect_cart()), 
             khi_shift(mpi), 
             ssjm1_nuf(mpi), 
             ssjm1_nuq(mpi), 
             ssjm1_dzeta(mpi), 
             ssjm1_tggg(mpi), 
             ssjm1_khi(mpi), 
             ssjm1_wshift(mpi, CON, mp.get_bvect_cart())    
{
    // Pointers of derived quantities initialized to zero : 
    set_der_0x0() ;
 
    // Initialization to a flat metric : 
    logn = 0 ;   
    tnphi = 0 ;   
    nuf = 0 ; 
    nuq = 0 ;
    dzeta = 0 ; 
    tggg = 0 ;   
 
    w_shift.set_etat_zero() ; 
    khi_shift =  0 ; 
 
    beta.set_etat_zero() ; 
 
    ssjm1_nuf = 0 ; 
    ssjm1_nuq = 0 ; 
    ssjm1_dzeta = 0 ; 
    ssjm1_tggg = 0 ; 
    ssjm1_khi = 0 ; 
 
    ssjm1_wshift.set_etat_zero() ; 
 
}
 
// Copy constructor
// --------------------
Star_QI::Star_QI(const Star_QI& st) :
             Compobj_QI(st), 
             logn(st.logn), 
             tnphi(st.tnphi), 
             nuf(st.nuf), 
             nuq(st.nuq),
             dzeta(st.dzeta), 
             tggg(st.tggg), 
             w_shift(st.w_shift), 
             khi_shift(st.khi_shift), 
             ssjm1_nuf(st.ssjm1_nuf), 
             ssjm1_nuq(st.ssjm1_nuq), 
             ssjm1_dzeta(st.ssjm1_dzeta), 
             ssjm1_tggg(st.ssjm1_tggg), 
             ssjm1_khi(st.ssjm1_khi), 
             ssjm1_wshift(st.ssjm1_wshift)           
{
    // Pointers of derived quantities initialized to zero : 
    set_der_0x0() ;
}
 
 
// Constructor from a file
// -----------------------
Star_QI::Star_QI(Map& mpi, FILE* fich) :
             Compobj_QI(mpi) , 
             logn(mpi), 
             tnphi(mpi), 
             nuf(mpi), 
             nuq(mpi), 
             dzeta(mpi), 
             tggg(mpi), 
             w_shift(mpi, CON, mp.get_bvect_cart()), 
             khi_shift(mpi), 
             ssjm1_nuf(mpi), 
             ssjm1_nuq(mpi), 
             ssjm1_dzeta(mpi), 
             ssjm1_tggg(mpi), 
             ssjm1_khi(mpi), 
             ssjm1_wshift(mpi, CON, mp.get_bvect_cart())             
{
    // Pointers of derived quantities initialized to zero : 
    set_der_0x0() ;
    
    // Read of the saved fields:
    // ------------------------
 
    Scalar nuf_file(mp, *(mp.get_mg()), fich) ; 
    nuf = nuf_file ; 
        
    Scalar nuq_file(mp, *(mp.get_mg()), fich) ; 
    nuq = nuq_file ; 
        
    logn = nuf + nuq ; //## to be checked !
    
    Scalar dzeta_file(mp, *(mp.get_mg()), fich) ; 
    dzeta = dzeta_file ; 
        
    Scalar tggg_file(mp, *(mp.get_mg()), fich) ; 
    tggg = tggg_file ; 
        
    Vector w_shift_file(mp, mp.get_bvect_cart(), fich) ; 
    w_shift = w_shift_file ;
    
    Scalar khi_shift_file(mp, *(mp.get_mg()), fich) ; 
    khi_shift = khi_shift_file ;
    
    fait_shift() ;      // constructs shift from w_shift and khi_shift
    fait_nphi() ;       // constructs N^phi from (N^x,N^y,N^z)
    
    Scalar ssjm1_nuf_file(mp, *(mp.get_mg()), fich) ; 
    ssjm1_nuf = ssjm1_nuf_file ; 
 
    Scalar ssjm1_nuq_file(mp, *(mp.get_mg()), fich) ; 
    ssjm1_nuq = ssjm1_nuq_file ; 
 
    Scalar ssjm1_dzeta_file(mp, *(mp.get_mg()), fich) ; 
    ssjm1_dzeta = ssjm1_dzeta_file ; 
 
    Scalar ssjm1_tggg_file(mp, *(mp.get_mg()), fich) ; 
    ssjm1_tggg = ssjm1_tggg_file ; 
 
    Scalar ssjm1_khi_file(mp, *(mp.get_mg()), fich) ; 
    ssjm1_khi = ssjm1_khi_file ; 
 
    Vector ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ; 
    ssjm1_wshift = ssjm1_wshift_file ; 
 
     
    // Initialization of N, A^2, B^2, gamma_ij, tkij and ak_car
    update_metric() ; 
    
}
 
                //------------//
                // Destructor //
                //------------//
 
Star_QI::~Star_QI(){
 
    del_deriv() ; 
 
}
 
 
            //----------------------------------//
            // Management of derived quantities //
            //----------------------------------//
 
void Star_QI::del_deriv() const {
 
    Compobj_QI::del_deriv() ; 
 
    if (p_grv2 != 0x0) delete p_grv2 ; 
    if (p_grv3 != 0x0) delete p_grv3 ;
    if (p_mom_quad != 0x0) delete p_mom_quad ;
    if (p_mass_g != 0x0) delete p_mass_g ;
 
    Star_QI::set_der_0x0() ; 
}               
 
 
void Star_QI::set_der_0x0() const {
 
    p_grv2 = 0x0 ; 
    p_grv3 = 0x0 ;
    p_mom_quad = 0x0 ;
    p_mass_g = 0x0 ;
     
}               
 
                //--------------//
                //  Assignment  //
                //--------------//
 
// Assignment to another Star_QI
// --------------------------------
void Star_QI::operator=(const Star_QI& st) {
 
    // Assignment of quantities common to all the derived classes of Compobj_QI
    Compobj_QI::operator=(st) ;     
    
    logn = st.logn ; 
    tnphi = st.tnphi ; 
    nuf = st.nuf ; 
    nuq = st.nuq ; 
    dzeta = st.dzeta ; 
    tggg = st.tggg ; 
    w_shift = st.w_shift ;
    khi_shift = st.khi_shift ;
    ssjm1_nuf = st.ssjm1_nuf ;
    ssjm1_nuq = st.ssjm1_nuq ;
    ssjm1_dzeta = st.ssjm1_dzeta ; 
    ssjm1_tggg = st.ssjm1_tggg ;
    ssjm1_khi = st.ssjm1_khi ;
    ssjm1_wshift = st.ssjm1_wshift ; 
 
    del_deriv() ;  // Deletes all derived quantities
}   
 
                //--------------//
                //    Outputs   //
                //--------------//
 
// Save in a file
// --------------
void Star_QI::sauve(FILE* fich) const {
 
    nuf.sauve(fich) ; 
    nuq.sauve(fich) ; 
    dzeta.sauve(fich) ; 
    tggg.sauve(fich) ; 
    w_shift.sauve(fich) ; 
    khi_shift.sauve(fich) ; 
    
    ssjm1_nuf.sauve(fich) ; 
    ssjm1_nuq.sauve(fich) ; 
    ssjm1_dzeta.sauve(fich) ; 
    ssjm1_tggg.sauve(fich) ; 
    ssjm1_khi.sauve(fich) ; 
    ssjm1_wshift.sauve(fich) ; 
    
}
 
// Printing
// --------
 
ostream& Star_QI::operator>>(ostream& ost) const {
 
    using namespace Unites ;
    
    Compobj_QI::operator>>(ost) ; 
    
    ost << endl << "Axisymmetric stationary compact star in quasi-isotropic coordinates (class Star_QI) " << endl ; 
 
    ost << "Central values of various fields : " << endl ; 
    ost << "-------------------------------- " << endl ; 
    ost << "   ln(N) : " << logn.val_grid_point(0,0,0,0) << endl ; 
    ost << "   nuf : " << nuf.val_grid_point(0,0,0,0) << endl ; 
    ost << "   nuq : " << nuq.val_grid_point(0,0,0,0) << endl ; 
    ost << "   zeta = ln(AN): " << dzeta.val_grid_point(0,0,0,0) << endl << endl ;
    
    ost << "Error on the virial identity GRV2 : " <<  grv2() << endl ; 
    ost << "Error on the virial identity GRV3 : " <<  grv3(&ost) << endl ; 
 
    double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.)) 
                            / double(1.e38) ) ;
    ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
         << endl ; 
    ost << "c^4 Q / (G^2 M^3) :      " 
     << mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) ) 
     << endl ; 
    
    ost << "Total angular momentum J :      " 
     << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
     << endl ; 
    ost << "c J / (G M^2) :           " 
     << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ;           
        
    return ost ; 
      
}
 
                //-------------------------//
                //  Computational methods  //
                //-------------------------//
                
void Star_QI::fait_shift() {
 
    Vector d_khi = khi_shift.derive_con( mp.flat_met_cart() ) ; 
    
    d_khi.dec_dzpuis(2) ;   // divide by r^2 in the external compactified domain  
 
    // x_k dW^k/dx_i
      Scalar xx(mp) ; 
      Scalar yy(mp) ; 
      Scalar zz(mp) ; 
      Scalar sintcosp(mp) ;
      Scalar sintsinp(mp) ;
      Scalar cost(mp) ;
      xx = mp.x ; 
      yy = mp.y ; 
      zz = mp.z ; 
      sintcosp = mp.sint * mp.cosp ;
      sintsinp = mp.sint * mp.sinp ;
      cost = mp.cost ;
 
      int nz = mp.get_mg()->get_nzone() ;
      Vector xk(mp, COV, mp.get_bvect_cart()) ;
      xk.set(1) = xx ;
      xk.set(1).set_domain(nz-1) = sintcosp.domain(nz-1) ;
      xk.set(1).set_dzpuis(-1) ;
      xk.set(2) = yy ;
      xk.set(2).set_domain(nz-1) = sintsinp.domain(nz-1) ;
      xk.set(2).set_dzpuis(-1) ;
      xk.set(3) = zz ;
      xk.set(3).set_domain(nz-1) = cost.domain(nz-1) ;
      xk.set(3).set_dzpuis(-1) ;
      xk.std_spectral_base() ;
      
      Tensor d_w = w_shift.derive_con( mp.flat_met_cart() ) ;
      
      Vector x_d_w = contract(xk, 0, d_w, 0) ;
      x_d_w.dec_dzpuis() ;
 
      double lambda = double(1) / double(3) ; 
 
      beta = - (lambda+2)/2./(lambda+1) * w_shift
        + (lambda/2./(lambda+1))  * (d_khi + x_d_w) ;      
    
} 
 
 
void Star_QI::fait_nphi() {
 
    if ( (beta(1).get_etat() == ETATZERO) && (beta(2).get_etat() == ETATZERO) ) {
        tnphi = 0 ; 
        nphi = 0 ;
    return ;  
    }
 
    // Computation of tnphi
    // --------------------
    
    mp.comp_p_from_cartesian( -beta(1), -beta(2), tnphi ) ; 
 
    // Computation of nphi
    // -------------------
    
    nphi = tnphi ; 
    nphi.div_rsint() ; 
    
}
 
 
// Updates the 3-metric and the shift
 
void Star_QI::update_metric() {
 
    // Lapse function N
    // ----------------
    
    nn = exp( logn ) ; 
 
    nn.std_spectral_base() ;   // set the bases for spectral expansions
    
    
    // Metric factor A^2
    // -----------------
    
    a_car = exp( 2*( dzeta - logn ) ) ; 
 
    a_car.std_spectral_base() ;   // set the bases for spectral expansions
 
    // Metric factor B
    // ---------------
    
    Scalar tmp = tggg ; 
    tmp.div_rsint() ;           //... Division of tG by r sin(theta)
 
    bbb = (1 + tmp) / nn ; 
 
    bbb.std_spectral_base() ;   // set the bases for spectral expansions
        
    b_car = bbb * bbb ; 
    
 
    Compobj_QI::update_metric() ;  // updates gamma_{ij} 
    
    
    // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
    // -------------------------------------------------
    
    //## extrinsic_curvature() ; // should be done by Compobj_QI::update_metric()
    
  
    // The derived quantities are no longer up to date : 
    // -----------------------------------------------
 
    del_deriv() ;  
    
}
 
 
 
}