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
29char 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
60namespace Lorene {
61double 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
74double 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
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
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
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
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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & dsdr() const
Returns of *this .
double radius_p(double phi)
Radius of the star to the direction of and .
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
double phi_local_min(double phi_ini)
Azimuthal angle when the indicator of the mass-shedding takes its local minimum.
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Scalar ent
Log-enthalpy.
Definition star.h:190
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
double ray_pole() const
Coordinate radius at [r_unit].
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Lorene prototypes.
Definition app_hor.h:64