compobj_QI.C

/*
 *  Methods of the class Compobj_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 compobj_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $" ;
 
/*
 * $Id: compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
 * $Log: compobj_QI.C,v $
 * Revision 1.9  2014/10/13 08:52:49  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.8  2014/05/16 11:55:18  o_straub
 * fixed: GYOTO output from compobj & compobj_QI
 *
 * Revision 1.7  2013/07/25 19:44:11  o_straub
 * calculation of the marginally bound radius
 *
 * Revision 1.6  2013/04/04 15:32:32  e_gourgoulhon
 * Better computation of the ISCO
 *
 * Revision 1.5  2013/04/04 08:53:47  e_gourgoulhon
 * Minor improvements
 *
 * Revision 1.4  2013/04/03 12:10:13  e_gourgoulhon
 * Added member kk to Compobj; suppressed tkij
 *
 * Revision 1.3  2012/11/22 16:04:51  c_some
 * Minor modifications
 *
 * Revision 1.2  2012/11/20 16:28:48  c_some
 * -- tkij is created on the Cartesian triad.
 * -- implemented method extrinsic_curvature()
 *
 * Revision 1.1  2012/11/16 16:14:11  c_some
 * New class Compobj_QI
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
 *
 */
 
 
// C headers
#include <cassert>
#include <cmath>
#include <cstdio>
 
// Lorene headers
#include "compobj.h"
#include "nbr_spx.h"
#include "utilitaires.h"
 
 
 
                   //--------------//
                   // Constructors //
                   //--------------//
 
// Standard constructor
// --------------------
namespace Lorene {
Compobj_QI::Compobj_QI(Map& map_i) :
        Compobj(map_i) ,
        a_car(map_i) ,
        bbb(map_i) ,
        b_car(map_i) ,
        nphi(map_i) ,
        ak_car(map_i) 
{
    // Pointers of derived quantities initialized to zero : 
    set_der_0x0() ;
 
    // Initialization to a flat metric : 
    a_car = 1 ;
    a_car.std_spectral_base() ; 
    bbb = 1 ;
    bbb.std_spectral_base() ; 
    b_car = bbb*bbb ;
    nphi = 0 ;   
    ak_car = 0 ; 
 
}
 
// Copy constructor
// --------------------
Compobj_QI::Compobj_QI(const Compobj_QI& co) :
        Compobj(co), 
        a_car(co.a_car) , 
        bbb(co.bbb) , 
        b_car(co.b_car) ,
        nphi(co.nphi) ,
        ak_car(co.ak_car) 
{
    // Pointers of derived quantities initialized to zero : 
    set_der_0x0() ;
}
 
 
// Constructor from a file
// -----------------------
Compobj_QI::Compobj_QI(Map& map_i, FILE* fich) :
        Compobj(map_i) , 
        a_car(map_i, *(map_i.get_mg()), fich) , 
        bbb(map_i, *(map_i.get_mg()), fich) , 
        b_car(map_i) , 
        nphi(map_i, *(map_i.get_mg()), fich) , 
        ak_car(map_i) 
{
    // Pointers of derived quantities initialized to zero : 
    set_der_0x0() ;
    
    Scalar nn_file(mp, *(mp.get_mg()), fich) ; 
    nn = nn_file ;
    
    b_car = bbb*bbb ;
     
    // Initialization of gamma_ij:
    update_metric() ; 
    
    // Computation of K_ij and ak_car:
    extrinsic_curvature() ; 
}
 
                //------------//
                // Destructor //
                //------------//
 
Compobj_QI::~Compobj_QI(){
 
    del_deriv() ; 
 
}
 
 
            //----------------------------------//
            // Management of derived quantities //
            //----------------------------------//
 
void Compobj_QI::del_deriv() const {
 
    Compobj::del_deriv() ; 
 
    if (p_angu_mom != 0x0) delete p_angu_mom ; 
    if (p_r_mb != 0x0) delete p_r_mb ;
    if (p_r_isco != 0x0) delete p_r_isco ;
    if (p_f_isco != 0x0) delete p_f_isco ;
    if (p_lspec_isco != 0x0) delete p_lspec_isco ;
    if (p_espec_isco != 0x0) delete p_espec_isco ;
 
    Compobj_QI::set_der_0x0() ; 
}               
 
 
void Compobj_QI::set_der_0x0() const {
 
    p_angu_mom = 0x0 ; 
    p_r_mb = 0x0 ;
    p_r_isco = 0x0 ;
    p_f_isco = 0x0 ;
    p_lspec_isco = 0x0 ;
    p_espec_isco = 0x0 ;
     
}               
 
                //--------------//
                //  Assignment  //
                //--------------//
 
// Assignment to another Compobj_QI
// --------------------------------
void Compobj_QI::operator=(const Compobj_QI& co) {
 
    // Assignment of quantities common to all the derived classes of Compobj
    Compobj::operator=(co) ;        
    
    a_car = co.a_car ; 
    bbb = co.bbb ; 
    b_car = co.b_car ; 
    nphi = co.nphi ; 
    ak_car = co.ak_car ;
 
    del_deriv() ;  // Deletes all derived quantities
}   
 
                //--------------//
                //    Outputs   //
                //--------------//
 
 
// Save in a file
// --------------
void Compobj_QI::sauve(FILE* fich) const {
 
    a_car.sauve(fich) ; 
    bbb.sauve(fich) ;
    nphi.sauve(fich) ;
 
    nn.sauve(fich) ;     
 
}
 
            
// Save in a file for GYOTO input  
// -------------------------------
 
// Redefinition of Compobj::gyoto_data
 
void Compobj_QI::gyoto_data(const char* file_name) const {
    
    FILE* file_out = fopen(file_name, "w") ;
    double total_time = 0. ; // for compatibility
    double RISCO=r_isco(0) ;
    double RMB=r_mb(0);
 
    fwrite_be(&total_time, sizeof(double), 1, file_out) ;
    mp.get_mg()->sauve(file_out) ;
    mp.sauve(file_out) ;
    nn.sauve(file_out) ;
    beta.sauve(file_out) ;
    gamma.cov().sauve(file_out) ;
    gamma.con().sauve(file_out) ;
    kk.sauve(file_out) ;
    fwrite_be(&RISCO, sizeof(double), 1, file_out) ;
    fwrite_be(&RMB, sizeof(double), 1, file_out) ;
    fclose(file_out) ;    
 
    cout << "WRITING RISCO TO GYOTO FILE : " << RISCO << endl ; 
    cout << "WRITING RMB TO GYOTO FILE   : " << RMB << endl ; 
    cout << "WRITING TO GYOTO FILE - end of part " << endl ; 
}
 
 
 
 
 
 
// Printing
// --------
 
ostream& Compobj_QI::operator>>(ostream& ost) const {
   
    Compobj::operator>>(ost) ; 
    
    ost << endl << "Axisymmetric stationary compact object in quasi-isotropic coordinates (class Compobj_QI) " << endl ; 
 
    ost << "Central values of various fields : " << endl ; 
    ost << "-------------------------------- " << endl ; 
    ost << "   metric coefficient A^2   : " << a_car.val_grid_point(0,0,0,0) << endl ; 
    ost << "   metric coefficient B^2   : " << b_car.val_grid_point(0,0,0,0) << endl ; 
    ost << "   metric coefficient N^phi : " << nphi.val_grid_point(0,0,0,0) << endl ; 
    ost << "   A^2 K_{ij} K^{ij}        = " << ak_car.val_grid_point(0,0,0,0) << endl << endl ; 
 
    
    double RISCO = r_isco(0, &ost) ; 
    ost << "Coordinate r at the innermost stable circular orbit (ISCO) : " << 
        RISCO  << endl ;  
     // ost << "Circumferential radius of the innermost stable circular orbit (ISCO) : " << 
     // bbb.val_point(risco, M_PI/2, 0)*risco << endl ;  
     // ost << "Orbital frequency at the ISCO : " << f_isco(0) << endl ; 
     // ost << "Specific energy of a particle on the ISCO : " << espec_isco(0) << endl ;    
     // ost << "Specific angular momentum of a particle on the ISCO : " << lspec_isco(0) << endl ;  
    
 
    double RMB = r_mb(0, &ost) ; 
    ost << "Coordinate r at the marginally bound circular orbit (R_mb) : " << RMB  << endl ; 
 
 
  //  ost << "A^2 : " << a_car << endl ; 
  //  ost << "B^2 : " << b_car << endl ; 
  //  ost << "nphi : " << nphi << endl ; 
    
    return ost ; 
      
}
 
// Updates the 3-metric and the shift
 
void Compobj_QI::update_metric() {
 
    Sym_tensor gam(mp, COV, mp.get_bvect_spher()) ; 
    gam.set(1,1) = a_car ; 
    gam.set(1,2) = 0 ; 
    gam.set(1,3) = 0 ; 
    gam.set(2,2) = a_car ; 
    gam.set(2,3) = 0 ; 
    gam.set(3,3) = b_car ;
    
    gamma = gam ;
 
    assert(*(beta.get_triad()) == mp.get_bvect_spher()) ; 
    
    beta.set(1) = 0 ;
    beta.set(2) = 0 ;
    Scalar nphi_ortho(nphi) ; 
    nphi_ortho.mult_rsint() ;
    beta.set(3) = - nphi_ortho ; 
    
    // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
    // -------------------------------------------------
    
    extrinsic_curvature() ; 
    
  
    // The derived quantities are no longer up to date : 
    // -----------------------------------------------
 
    del_deriv() ;  
    
}
 
 
// Updates the extrinsic curvature
// -------------------------------
 
void Compobj_QI::extrinsic_curvature() {
 
    // Special treatment for axisymmetric case:
    
    if ( (mp.get_mg())->get_np(0) == 1) {
    
        Scalar dnpdr = nphi.dsdr() ;        // d/dr (N^phi)
        Scalar dnpdt = nphi.dsdt() ;        // d/dtheta (N^phi)
        
        // What follows is valid only for a mapping of class Map_radial :   
        assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
        
        dnpdr.mult_rsint() ;    // multiplication by r sin(theta)
        kk.set(1,3) = - b_car * dnpdr / (2*nn) ; 
        
        dnpdt.mult_sint() ; // multiplication by sin(theta)
        kk.set(2,3) = - b_car * dnpdt / (2*nn) ; 
        kk.set(2,3).inc_dzpuis(2) ;  // to have the same dzpuis as kk(1,3)
        
        kk.set(1,1) = 0 ; 
        kk.set(1,2) = 0 ; 
        kk.set(2,2) = 0 ; 
        kk.set(3,3) = 0 ; 
    }
    else {
 
    // General case:
 
        Compobj::extrinsic_curvature() ; 
   }
    
    // Computation of A^2 K_{ij} K^{ij}
    // --------------------------------
        
    ak_car = 2 * ( kk(1,3)*kk(1,3) +  kk(2,3)*kk(2,3) ) / b_car ;
    
    del_deriv() ; 
 
}
}