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 
30 
31 /*
32  * $Id: star_bhns_chi.C,v 1.4 2016/12/05 16:18:16 j_novak Exp $
33  * $Log: star_bhns_chi.C,v $
34  * Revision 1.4 2016/12/05 16:18:16 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.3 2014/10/13 08:53:40 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.2 2014/10/06 15:13:16 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.1 2007/06/22 01:30:27 k_taniguchi
44  * *** empty log message ***
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_chi.C,v 1.4 2016/12/05 16:18:16 j_novak Exp $
48  *
49  */
50 
51 // C++ headers
52 //#include <>
53 
54 // C headers
55 #include <cmath>
56 
57 // Lorene headers
58 #include "star_bhns.h"
59 #include "param.h"
60 #include "tbl.h"
61 #include "utilitaires.h"
62 
63 namespace Lorene {
64 double Star_bhns::chi_rp(double radius, double phi) {
65 
66  const Scalar& dent = ent.dsdr() ;
67 
68  double dent_pole = dent.val_point(ray_pole(), 0., 0.) ;
69  double dent_eq = dent.val_point(radius, M_PI/2., phi) ;
70 
71  double chi = fabs( dent_eq / dent_pole ) ;
72 
73  return chi ;
74 
75 }
76 
77 double Star_bhns::radius_p(double phi) {
78 
79  double rad = mp.val_r(nzet-1, 1., M_PI/2., phi) ;
80  // We assume that the stellar surface is fitted to the domain (nzet-1)
81 
82  return rad ;
83 
84 }
85 
87 
88  const Mg3d* mg = mp.get_mg() ;
89  int np = mg->get_np(0) ;
90  int nps2 = np/2 ;
91 
92  Tbl phi(nps2+1) ;
93  phi.set_etat_qcq() ;
94 
95  for (int i=0; i<=nps2; i++) {
96  phi.set(i) = 2.*M_PI*i/np + 0.5*M_PI ;
97  }
98 
99  Tbl chi(nps2+1) ;
100  chi.set_etat_qcq() ;
101 
102  for (int i=0; i<=nps2; i++) {
103 
104  double phi_i = phi_local_min( phi(i) ) ;
105  double rad_i = radius_p( phi_i ) ;
106 
107  chi.set(i) = chi_rp(rad_i, phi_i) ;
108 
109  }
110 
111  for (int i=0; i<=nps2; i++) {
112 
113  cout.precision(16) ;
114  cout << "chi(" << i << ") = " << chi(i)
115  << " phi = " << phi_local_min( phi(i) ) / M_PI
116  << " [M_PI]" << endl ;
117  cout.precision(4) ;
118 
119  }
120 
121  double chi_ini = chi(0) ; // Initialization
122  double delta_chi ;
123  int jj = 0 ; // Initialization
124 
125  for (int i=0; i<nps2; i++) {
126 
127  if ( chi(i+1) < 1.e-12 )
128  chi.set(i+1) = 1. ;
129 
130  delta_chi = chi_ini - chi(i+1) ;
131 
132  if ( delta_chi > 0. ) {
133 
134  chi_ini = chi(i+1) ;
135  jj = i+1 ;
136 
137  }
138 
139  }
140 
141  double phi_glob_min = phi_local_min( phi(jj) ) ;
142 
143  return phi_glob_min ;
144 
145 }
146 
147 
148 double Star_bhns::phi_local_min(double phi_ini) {
149 
150  int mm ; // Number of steps to the phi direction
151  double ppp = phi_ini ; // Initial position of the phi coordinate
152  double diff ; // Difference between two succesive chi
153  double dp = M_PI/2. ; // Step interval to the phi direction
154  double ptmp ;
155 
156  double rad1, rad2 ;
157 
158  double init_check = chi_rp(radius_p(phi_ini), phi_ini)
159  - chi_rp(radius_p(phi_ini+1.e-10*dp), phi_ini+1.e-10*dp) ;
160 
161  if ( init_check >= 0. ) {
162 
163  while ( dp > 1.e-15 ) {
164 
165  diff = 1. ;
166  mm = 0 ;
167  dp = 0.1 * dp ;
168 
169  while ( diff > 0. && (ppp+mm*dp) < 2.*M_PI ) {
170 
171  mm++ ;
172  ptmp = ppp + mm * dp ;
173 
174  rad1 = radius_p(ptmp-dp) ;
175  rad2 = radius_p(ptmp) ;
176 
177  diff = chi_rp(rad1, ptmp-dp) - chi_rp(rad2, ptmp) ;
178 
179  }
180 
181  ppp += (mm - 2) * dp ;
182 
183  }
184 
185  if ( (ppp+2.*dp) >= 2.*M_PI ) {
186 
187  cout << "No minimum for phi > " << phi_ini / M_PI
188  << " [M_PI]" << endl ;
189 
190  }
191 
192  }
193  else {
194 
195  while ( dp > 1.e-15 ) {
196 
197  diff = 1. ;
198  mm = 0 ;
199  dp = 0.1 * dp ;
200 
201  while ( diff > 0. && (ppp-mm*dp) > 0. ) {
202 
203  mm++ ;
204  ptmp = ppp - mm * dp ;
205 
206  rad1 = radius_p(ptmp+dp) ;
207  rad2 = radius_p(ptmp) ;
208 
209  diff = chi_rp(rad1, ptmp+dp) - chi_rp(rad2, ptmp) ;
210 
211  }
212 
213  ppp -= (mm - 2) * dp ;
214 
215  }
216 
217  if ( (ppp-2.*dp) < 0. ) {
218 
219  cout << "No minimum for phi < " << phi_ini / M_PI
220  << " [M_PI]" << endl ;
221 
222  }
223 
224  }
225 
226  return ppp ;
227 
228 }
229 }
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:479
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar ent
Log-enthalpy.
Definition: star.h:190
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
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:77
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:896
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Definition: star_bhns_chi.C:86
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:281
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
Definition: star_bhns_chi.C:64
Multi-domain grid.
Definition: grilles.h:279
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Basic array class.
Definition: tbl.h:164