star_bhns_chi.C
/*
* Method of class Star_bhns to compute a sensitve indicator of
* the mass-shedding and quantities related to the indicator
*
* (see file star_bhns.h for documentation).
*
*/
/*
* Copyright (c) 2006 Keisuke Taniguchi
*
* 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_bhns_chi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_chi.C,v 1.3 2014/10/13 08:53:40 j_novak Exp $" ;
/*
* $Id: star_bhns_chi.C,v 1.3 2014/10/13 08:53:40 j_novak Exp $
* $Log: star_bhns_chi.C,v $
* Revision 1.3 2014/10/13 08:53:40 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.2 2014/10/06 15:13:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.1 2007/06/22 01:30:27 k_taniguchi
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_chi.C,v 1.3 2014/10/13 08:53:40 j_novak Exp $
*
*/
// C++ headers
//#include <>
// C headers
#include <cmath>
// Lorene headers
#include "star_bhns.h"
#include "param.h"
#include "tbl.h"
#include "utilitaires.h"
namespace Lorene {
double Star_bhns::chi_rp(double radius, double phi) {
const Scalar& dent = ent.dsdr() ;
double dent_pole = dent.val_point(ray_pole(), 0., 0.) ;
double dent_eq = dent.val_point(radius, M_PI/2., phi) ;
double chi = fabs( dent_eq / dent_pole ) ;
return chi ;
}
double Star_bhns::radius_p(double phi) {
double rad = mp.val_r(nzet-1, 1., M_PI/2., phi) ;
// We assume that the stellar surface is fitted to the domain (nzet-1)
return rad ;
}
double Star_bhns::phi_min() {
const Mg3d* mg = mp.get_mg() ;
int np = mg->get_np(0) ;
int nps2 = np/2 ;
Tbl phi(nps2+1) ;
phi.set_etat_qcq() ;
for (int i=0; i<=nps2; i++) {
phi.set(i) = 2.*M_PI*i/np + 0.5*M_PI ;
}
Tbl chi(nps2+1) ;
chi.set_etat_qcq() ;
for (int i=0; i<=nps2; i++) {
double phi_i = phi_local_min( phi(i) ) ;
double rad_i = radius_p( phi_i ) ;
chi.set(i) = chi_rp(rad_i, phi_i) ;
}
for (int i=0; i<=nps2; i++) {
cout.precision(16) ;
cout << "chi(" << i << ") = " << chi(i)
<< " phi = " << phi_local_min( phi(i) ) / M_PI
<< " [M_PI]" << endl ;
cout.precision(4) ;
}
double chi_ini = chi(0) ; // Initialization
double delta_chi ;
int jj = 0 ; // Initialization
for (int i=0; i<nps2; i++) {
if ( chi(i+1) < 1.e-12 )
chi.set(i+1) = 1. ;
delta_chi = chi_ini - chi(i+1) ;
if ( delta_chi > 0. ) {
chi_ini = chi(i+1) ;
jj = i+1 ;
}
}
double phi_glob_min = phi_local_min( phi(jj) ) ;
return phi_glob_min ;
}
double Star_bhns::phi_local_min(double phi_ini) {
int mm ; // Number of steps to the phi direction
double ppp = phi_ini ; // Initial position of the phi coordinate
double diff ; // Difference between two succesive chi
double dp = M_PI/2. ; // Step interval to the phi direction
double ptmp ;
double rad1, rad2 ;
double init_check = chi_rp(radius_p(phi_ini), phi_ini)
- chi_rp(radius_p(phi_ini+1.e-10*dp), phi_ini+1.e-10*dp) ;
if ( init_check >= 0. ) {
while ( dp > 1.e-15 ) {
diff = 1. ;
mm = 0 ;
dp = 0.1 * dp ;
while ( diff > 0. && (ppp+mm*dp) < 2.*M_PI ) {
mm++ ;
ptmp = ppp + mm * dp ;
rad1 = radius_p(ptmp-dp) ;
rad2 = radius_p(ptmp) ;
diff = chi_rp(rad1, ptmp-dp) - chi_rp(rad2, ptmp) ;
}
ppp += (mm - 2) * dp ;
}
if ( (ppp+2.*dp) >= 2.*M_PI ) {
cout << "No minimum for phi > " << phi_ini / M_PI
<< " [M_PI]" << endl ;
}
}
else {
while ( dp > 1.e-15 ) {
diff = 1. ;
mm = 0 ;
dp = 0.1 * dp ;
while ( diff > 0. && (ppp-mm*dp) > 0. ) {
mm++ ;
ptmp = ppp - mm * dp ;
rad1 = radius_p(ptmp+dp) ;
rad2 = radius_p(ptmp) ;
diff = chi_rp(rad1, ptmp+dp) - chi_rp(rad2, ptmp) ;
}
ppp -= (mm - 2) * dp ;
}
if ( (ppp-2.*dp) < 0. ) {
cout << "No minimum for phi < " << phi_ini / M_PI
<< " [M_PI]" << endl ;
}
}
return ppp ;
}
}
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition scalar.C:890
Lorene::Star_bhns::radius_p
double radius_p(double phi)
Radius of the star to the direction of and .
Definition star_bhns_chi.C:74
Lorene::Star_bhns::chi_rp
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
Definition star_bhns_chi.C:61
Lorene::Star_bhns::phi_local_min
double phi_local_min(double phi_ini)
Azimuthal angle when the indicator of the mass-shedding takes its local minimum.
Definition star_bhns_chi.C:145
Lorene::Star_bhns::phi_min
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Definition star_bhns_chi.C:83
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition star.h:190
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition star.h:180
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Lorene::Star::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition star_global.C:278
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720