binary_global_xcts.C

/*
 * Methods of class Binary_xcts to compute global quantities
 * (see file binary_xcts.h for documentation)
 */
 
/*
 *   Copyright (c) 2010 Michal Bejger
 *
 *   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
 *
 */
 
 
 
/*
 * $Id: binary_global_xcts.C,v 1.12 2016/12/05 16:17:47 j_novak Exp $
 * $Log: binary_global_xcts.C,v $
 * Revision 1.12  2016/12/05 16:17:47  j_novak
 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
 *
 * Revision 1.11  2014/10/13 08:52:45  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.10  2010/12/21 11:15:46  m_bejger
 * Linear momentum properly defined
 *
 * Revision 1.9  2010/12/20 15:45:40  m_bejger
 * Spectral basis in lin_mom() was not defined
 *
 * Revision 1.8  2010/12/20 09:54:09  m_bejger
 * Angular momentum correction, stub for linear momentum added
 *
 * Revision 1.7  2010/12/09 10:39:41  m_bejger
 * Further corrections to integral quantities
 *
 * Revision 1.6  2010/10/26 19:16:26  m_bejger
 * Cleanup of some diagnostic messages
 *
 * Revision 1.5  2010/10/25 15:02:08  m_bejger
 * mass_kom_vol() corrected
 *
 * Revision 1.4  2010/10/24 21:45:24  m_bejger
 * mass_adm() corrected
 *
 * Revision 1.3  2010/06/17 14:48:14  m_bejger
 * Minor corrections
 *
 * Revision 1.2  2010/06/04 19:54:19  m_bejger
 * Minor corrections, mass volume integrals need to be checked out
 *
 * Revision 1.1  2010/05/04 07:35:54  m_bejger
 * Initial version
 *
 * $Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_global_xcts.C,v 1.12 2016/12/05 16:17:47 j_novak Exp $
 *
 */
 
// Headers C
#include "math.h"
 
// Headers Lorene
#include "nbr_spx.h"
#include "binary_xcts.h"
#include "unites.h"
#include "metric.h"
 
            //---------------------------------//
            //      ADM mass                   //
            //---------------------------------//
 
namespace Lorene {
double Binary_xcts::mass_adm() const {
    
  using namespace Unites ;
  
  if (p_mass_adm == 0x0) {      // a new computation is required
    
    p_mass_adm = new double ; 
    *p_mass_adm = 0 ; 
    
    const Map_af map0 (et[0]->get_mp()) ;
    const Metric& flat = (et[0]->get_flat()) ;
 
    Vector dpsi((et[0]->get_Psi()).derive_cov(flat)) ;
 
    dpsi.change_triad(map0.get_bvect_spher()) ;
 
    Scalar integrand ( dpsi(1) ) ;
 
    *p_mass_adm = - 2.* map0.integrale_surface_infini(integrand)/ qpig ;
    
    }
 
    return *p_mass_adm ; 
    
}
 
            //---------------------------------//
            //          ADM mass               //
            //          (volume integral)      // 
            //---------------------------------//
 
double Binary_xcts::mass_adm_vol() const {
 
  using namespace Unites ;
 
  double massadm = 0. ;
 
  for (int i=0; i<=1; i++) {        // loop on the stars
 
    // Declaration of all fields
 
      const Scalar& psi(et[i]->get_Psi()) ;
      Scalar psi5 = pow(psi, 5.) ;
      psi5.std_spectral_base() ;
      
      Scalar spsi7 = pow(psi, -7.) ; 
      spsi7.std_spectral_base() ;
 
      const Scalar& ener_euler = et[i]->get_ener_euler() ;
      const Scalar& hacar_auto = et[i]->get_hacar_auto() ;
      const Scalar& hacar_comp = et[i]->get_hacar_comp() ;
 
      Scalar source = psi5 % ener_euler 
                    + spsi7 % (hacar_auto + hacar_comp)/(4.*qpig) ;  
      
      source.std_spectral_base() ;
 
      massadm += source.integrale() ;
  }
  
  return massadm ;
}
 
            //---------------------------------//
            //      Komar mass                 //
            //---------------------------------//
 
double Binary_xcts::mass_kom() const {
    
  using namespace Unites ;
 
  if (p_mass_kom == 0x0) {      // a new computation is requireed
 
    p_mass_kom = new double ;    
    *p_mass_kom = 0 ; 
 
    const Scalar& logn = et[0]->get_logn() ;
    const Metric& flat = et[0]->get_flat() ; 
 
    Map_af map0 (et[0]->get_mp()) ; 
      
    Vector vect = logn.derive_con(flat) ; 
    vect.change_triad(map0.get_bvect_spher()) ;
 
    Scalar integrant (vect(1)) ;
    
    *p_mass_kom = map0.integrale_surface_infini (integrant) / qpig ;
       
  } // End of the case where a new computation was necessary
        
  return *p_mass_kom ; 
    
}
 
double Binary_xcts::mass_kom_vol() const {
    
  using namespace Unites ;
 
  double masskom = 0.;
 
  for (int i=0; i<=1; i++) {        // loop on the stars
 
      const Scalar& Psi  = et[i]->get_Psi() ;     
      const Scalar& psi4 = et[i]->get_psi4() ;
      const Scalar& chi  = et[i]->get_chi() ; 
      
      const Scalar& ener_euler = et[i]->get_ener_euler() ;
      const Scalar& s_euler = et[i]->get_s_euler() ;
      const Scalar& hacar_auto = et[i]->get_hacar_auto() ;
      const Scalar& hacar_comp = et[i]->get_hacar_comp() ;  
 
      Scalar psi4chi = psi4 % chi ; 
      psi4chi.std_spectral_base() ; 
 
      Scalar source = 0.5 * ener_euler * (psi4chi + psi4 % Psi) 
          + psi4chi * s_euler + pow(Psi, -7.) * (7.*chi/Psi + 1.) 
                    * (hacar_auto + hacar_comp) / (8.*qpig) ; 
      source.std_spectral_base() ;
      
      masskom += source.integrale() ;
      
  }
 
  return masskom ;
 
}
 
                //-------------------------------------//
                //   Total angular momentum (z-axis)   //
                //-------------------------------------//
     
const Tbl& Binary_xcts::angu_mom() const {
 
  using namespace Unites ;
    
    if (p_angu_mom == 0x0) {        // a new computation is required
    
    p_angu_mom = new Tbl(3) ; 
    p_angu_mom->annule_hard() ; // fills the double array with zeros
 
    // Reference Cartesian vector basis of the Absolute frame
    Base_vect_cart bvect_ref(0.) ;  // 0. = parallel to the Absolute frame
    
    for (int i=0; i<=1; i++) {      // loop on the stars
 
        const Map& mp = et[i]->get_mp() ; 
        
        // Azimuthal vector d/dphi 
        Vector vphi(mp, CON, bvect_ref) ;       
        Scalar yya (mp) ; yya = mp.ya ;
        Scalar xxa (mp) ; xxa = mp.xa ; 
        vphi.set(1) = - yya ;   // phi^X
        vphi.set(2) = xxa ; 
        vphi.set(3) = 0 ;  
 
        vphi.std_spectral_base() ; 
        vphi.change_triad(mp.get_bvect_cart()) ; 
                
        // Matter part
        // -----------
        const Scalar& ee = et[i]->get_ener_euler() ;  
        const Scalar& pp = et[i]->get_press() ;
        
        Vector jmom = pow(et[i]->get_Psi(), 10) * (ee + pp) 
                    * (et[i]->get_u_euler()) ; 
        jmom.std_spectral_base() ;
        
        const Metric& flat = et[i]->get_flat() ;
        Vector vphi_cov = vphi.up_down(flat) ;
        
        Scalar integrand = contract(jmom, 0, vphi_cov, 0) ; 
              
        p_angu_mom->set(2) += integrand.integrale() ;
        
    }  // End of the loop on the stars
 
    }   // End of the case where a new computation was necessary
  
    return *p_angu_mom ; 
  
}
 
                //---------------------------------//
                //   Total linear momentum         //
                //---------------------------------//
    
const Tbl& Binary_xcts::lin_mom() const {
 
  using namespace Unites ;
    
    if (p_lin_mom == 0x0) {     // a new computation is required
    
    p_lin_mom = new Tbl(3) ; 
    p_lin_mom->annule_hard() ;
 
    // Reference Cartesian vector basis of the Absolute frame
    Base_vect_cart bvect_ref(0.) ;  // 0. = parallel to the Absolute frame
    
    for (int i=0; i<=1; i++) {          // loop on the stars
    
        const Scalar& ee = et[i]->get_ener_euler() ;  
        const Scalar& pp = et[i]->get_press() ;
        Vector lmom      = pow(et[i]->get_Psi(), 10) * (ee + pp) 
                         * ( et[i]->get_u_euler() ) ;  
    
        lmom.std_spectral_base() ;
        lmom.change_triad(bvect_ref) ; 
                
        // loop on the components               
        for (int j=1; j<=2; j++)        
            p_lin_mom->set(j-1) += lmom(j).integrale() ;
    
 
    } // End of the loop on the stars
 
    } // End of the case where a new computation was necessary
  
    return *p_lin_mom ; 
  
}
 
            //---------------------------------//
            //      Total energy               //
            //---------------------------------//
 
double Binary_xcts::total_ener() const {
   
    if (p_total_ener == 0x0) {      // a new computation is requireed
    
    p_total_ener = new double ; 
        
        *p_total_ener = mass_adm() - star1.mass_b() - star2.mass_b() ; 
        
    }   // End of the case where a new computation was necessary
    
    return *p_total_ener ; 
    
}
 
 
            //---------------------------------//
            //   Error on the virial theorem   //
            //---------------------------------//
 
double Binary_xcts::virial() const {
    
    if (p_virial == 0x0) {      // a new computation is requireed
    
    p_virial = new double ; 
        
        *p_virial = 1. - mass_kom() / mass_adm() ; 
        
    }
    
    return *p_virial ; 
    
}
 
            //---------------------------------//
            //   Error on the virial theorem   // 
            //   (volume version)              //
            //---------------------------------//
 
double Binary_xcts::virial_vol() const {
    
    if (p_virial_vol == 0x0) {      // a new computation is requireed
    
    p_virial_vol = new double ; 
        
        *p_virial_vol = 1. - mass_kom_vol() / mass_adm_vol() ; 
        
    }
    
    return *p_virial_vol ; 
    
}
}