LORENE
star_bhns_chi.C
1 /*
2  * Method of class Star_bhns to compute a sensitve indicator of
3  * the mass-shedding and quantities related to the indicator
4  *
5  * (see file star_bhns.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2006 Keisuke Taniguchi
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 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 $" ;
30 
31 /*
32  * $Id: star_bhns_chi.C,v 1.3 2014/10/13 08:53:40 j_novak Exp $
33  * $Log: star_bhns_chi.C,v $
34  * Revision 1.3 2014/10/13 08:53:40 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.2 2014/10/06 15:13:16 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.1 2007/06/22 01:30:27 k_taniguchi
41  * *** empty log message ***
42  *
43  *
44  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_chi.C,v 1.3 2014/10/13 08:53:40 j_novak Exp $
45  *
46  */
47 
48 // C++ headers
49 //#include <>
50 
51 // C headers
52 #include <cmath>
53 
54 // Lorene headers
55 #include "star_bhns.h"
56 #include "param.h"
57 #include "tbl.h"
58 #include "utilitaires.h"
59 
60 namespace Lorene {
61 double Star_bhns::chi_rp(double radius, double phi) {
62 
63  const Scalar& dent = ent.dsdr() ;
64 
65  double dent_pole = dent.val_point(ray_pole(), 0., 0.) ;
66  double dent_eq = dent.val_point(radius, M_PI/2., phi) ;
67 
68  double chi = fabs( dent_eq / dent_pole ) ;
69 
70  return chi ;
71 
72 }
73 
74 double Star_bhns::radius_p(double phi) {
75 
76  double rad = mp.val_r(nzet-1, 1., M_PI/2., phi) ;
77  // We assume that the stellar surface is fitted to the domain (nzet-1)
78 
79  return rad ;
80 
81 }
82 
84 
85  const Mg3d* mg = mp.get_mg() ;
86  int np = mg->get_np(0) ;
87  int nps2 = np/2 ;
88 
89  Tbl phi(nps2+1) ;
90  phi.set_etat_qcq() ;
91 
92  for (int i=0; i<=nps2; i++) {
93  phi.set(i) = 2.*M_PI*i/np + 0.5*M_PI ;
94  }
95 
96  Tbl chi(nps2+1) ;
97  chi.set_etat_qcq() ;
98 
99  for (int i=0; i<=nps2; i++) {
100 
101  double phi_i = phi_local_min( phi(i) ) ;
102  double rad_i = radius_p( phi_i ) ;
103 
104  chi.set(i) = chi_rp(rad_i, phi_i) ;
105 
106  }
107 
108  for (int i=0; i<=nps2; i++) {
109 
110  cout.precision(16) ;
111  cout << "chi(" << i << ") = " << chi(i)
112  << " phi = " << phi_local_min( phi(i) ) / M_PI
113  << " [M_PI]" << endl ;
114  cout.precision(4) ;
115 
116  }
117 
118  double chi_ini = chi(0) ; // Initialization
119  double delta_chi ;
120  int jj = 0 ; // Initialization
121 
122  for (int i=0; i<nps2; i++) {
123 
124  if ( chi(i+1) < 1.e-12 )
125  chi.set(i+1) = 1. ;
126 
127  delta_chi = chi_ini - chi(i+1) ;
128 
129  if ( delta_chi > 0. ) {
130 
131  chi_ini = chi(i+1) ;
132  jj = i+1 ;
133 
134  }
135 
136  }
137 
138  double phi_glob_min = phi_local_min( phi(jj) ) ;
139 
140  return phi_glob_min ;
141 
142 }
143 
144 
145 double Star_bhns::phi_local_min(double phi_ini) {
146 
147  int mm ; // Number of steps to the phi direction
148  double ppp = phi_ini ; // Initial position of the phi coordinate
149  double diff ; // Difference between two succesive chi
150  double dp = M_PI/2. ; // Step interval to the phi direction
151  double ptmp ;
152 
153  double rad1, rad2 ;
154 
155  double init_check = chi_rp(radius_p(phi_ini), phi_ini)
156  - chi_rp(radius_p(phi_ini+1.e-10*dp), phi_ini+1.e-10*dp) ;
157 
158  if ( init_check >= 0. ) {
159 
160  while ( dp > 1.e-15 ) {
161 
162  diff = 1. ;
163  mm = 0 ;
164  dp = 0.1 * dp ;
165 
166  while ( diff > 0. && (ppp+mm*dp) < 2.*M_PI ) {
167 
168  mm++ ;
169  ptmp = ppp + mm * dp ;
170 
171  rad1 = radius_p(ptmp-dp) ;
172  rad2 = radius_p(ptmp) ;
173 
174  diff = chi_rp(rad1, ptmp-dp) - chi_rp(rad2, ptmp) ;
175 
176  }
177 
178  ppp += (mm - 2) * dp ;
179 
180  }
181 
182  if ( (ppp+2.*dp) >= 2.*M_PI ) {
183 
184  cout << "No minimum for phi > " << phi_ini / M_PI
185  << " [M_PI]" << endl ;
186 
187  }
188 
189  }
190  else {
191 
192  while ( dp > 1.e-15 ) {
193 
194  diff = 1. ;
195  mm = 0 ;
196  dp = 0.1 * dp ;
197 
198  while ( diff > 0. && (ppp-mm*dp) > 0. ) {
199 
200  mm++ ;
201  ptmp = ppp - mm * dp ;
202 
203  rad1 = radius_p(ptmp+dp) ;
204  rad2 = radius_p(ptmp) ;
205 
206  diff = chi_rp(rad1, ptmp+dp) - chi_rp(rad2, ptmp) ;
207 
208  }
209 
210  ppp -= (mm - 2) * dp ;
211 
212  }
213 
214  if ( (ppp-2.*dp) < 0. ) {
215 
216  cout << "No minimum for phi < " << phi_ini / M_PI
217  << " [M_PI]" << endl ;
218 
219  }
220 
221  }
222 
223  return ppp ;
224 
225 }
226 }
double phi_local_min(double phi_ini)
Azimuthal angle when the indicator of the mass-shedding takes its local minimum.
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Scalar ent
Log-enthalpy.
Definition: star.h:190
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
double radius_p(double phi)
Radius of the star to the direction of and .
Definition: star_bhns_chi.C:74
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
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Definition: star_bhns_chi.C:83
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
Definition: star_bhns_chi.C:61
Multi-domain grid.
Definition: grilles.h:273
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Basic array class.
Definition: tbl.h:161