LORENE
compobj_QI.C
1/*
2 * Methods of the class Compobj_QI
3 *
4 * (see file compobj.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2012 Claire Some, Eric Gourgoulhon
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char compobj_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $" ;
29
30/*
31 * $Id: compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
32 * $Log: compobj_QI.C,v $
33 * Revision 1.9 2014/10/13 08:52:49 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.8 2014/05/16 11:55:18 o_straub
37 * fixed: GYOTO output from compobj & compobj_QI
38 *
39 * Revision 1.7 2013/07/25 19:44:11 o_straub
40 * calculation of the marginally bound radius
41 *
42 * Revision 1.6 2013/04/04 15:32:32 e_gourgoulhon
43 * Better computation of the ISCO
44 *
45 * Revision 1.5 2013/04/04 08:53:47 e_gourgoulhon
46 * Minor improvements
47 *
48 * Revision 1.4 2013/04/03 12:10:13 e_gourgoulhon
49 * Added member kk to Compobj; suppressed tkij
50 *
51 * Revision 1.3 2012/11/22 16:04:51 c_some
52 * Minor modifications
53 *
54 * Revision 1.2 2012/11/20 16:28:48 c_some
55 * -- tkij is created on the Cartesian triad.
56 * -- implemented method extrinsic_curvature()
57 *
58 * Revision 1.1 2012/11/16 16:14:11 c_some
59 * New class Compobj_QI
60 *
61 *
62 * $Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
63 *
64 */
65
66
67// C headers
68#include <cassert>
69#include <cmath>
70#include <cstdio>
71
72// Lorene headers
73#include "compobj.h"
74#include "nbr_spx.h"
75#include "utilitaires.h"
76
77
78
79 //--------------//
80 // Constructors //
81 //--------------//
82
83// Standard constructor
84// --------------------
85namespace Lorene {
87 Compobj(map_i) ,
88 a_car(map_i) ,
89 bbb(map_i) ,
90 b_car(map_i) ,
91 nphi(map_i) ,
92 ak_car(map_i)
93{
94 // Pointers of derived quantities initialized to zero :
95 set_der_0x0() ;
96
97 // Initialization to a flat metric :
98 a_car = 1 ;
100 bbb = 1 ;
102 b_car = bbb*bbb ;
103 nphi = 0 ;
104 ak_car = 0 ;
105
106}
107
108// Copy constructor
109// --------------------
111 Compobj(co),
112 a_car(co.a_car) ,
113 bbb(co.bbb) ,
114 b_car(co.b_car) ,
115 nphi(co.nphi) ,
116 ak_car(co.ak_car)
117{
118 // Pointers of derived quantities initialized to zero :
119 set_der_0x0() ;
120}
121
122
123// Constructor from a file
124// -----------------------
125Compobj_QI::Compobj_QI(Map& map_i, FILE* fich) :
126 Compobj(map_i) ,
127 a_car(map_i, *(map_i.get_mg()), fich) ,
128 bbb(map_i, *(map_i.get_mg()), fich) ,
129 b_car(map_i) ,
130 nphi(map_i, *(map_i.get_mg()), fich) ,
131 ak_car(map_i)
132{
133 // Pointers of derived quantities initialized to zero :
134 set_der_0x0() ;
135
136 Scalar nn_file(mp, *(mp.get_mg()), fich) ;
137 nn = nn_file ;
138
139 b_car = bbb*bbb ;
140
141 // Initialization of gamma_ij:
142 update_metric() ;
143
144 // Computation of K_ij and ak_car:
146}
147
148 //------------//
149 // Destructor //
150 //------------//
151
153
154 del_deriv() ;
155
156}
157
158
159 //----------------------------------//
160 // Management of derived quantities //
161 //----------------------------------//
162
164
166
167 if (p_angu_mom != 0x0) delete p_angu_mom ;
168 if (p_r_mb != 0x0) delete p_r_mb ;
169 if (p_r_isco != 0x0) delete p_r_isco ;
170 if (p_f_isco != 0x0) delete p_f_isco ;
171 if (p_lspec_isco != 0x0) delete p_lspec_isco ;
172 if (p_espec_isco != 0x0) delete p_espec_isco ;
173
175}
176
177
179
180 p_angu_mom = 0x0 ;
181 p_r_mb = 0x0 ;
182 p_r_isco = 0x0 ;
183 p_f_isco = 0x0 ;
184 p_lspec_isco = 0x0 ;
185 p_espec_isco = 0x0 ;
186
187}
188
189 //--------------//
190 // Assignment //
191 //--------------//
192
193// Assignment to another Compobj_QI
194// --------------------------------
196
197 // Assignment of quantities common to all the derived classes of Compobj
198 Compobj::operator=(co) ;
199
200 a_car = co.a_car ;
201 bbb = co.bbb ;
202 b_car = co.b_car ;
203 nphi = co.nphi ;
204 ak_car = co.ak_car ;
205
206 del_deriv() ; // Deletes all derived quantities
207}
208
209 //--------------//
210 // Outputs //
211 //--------------//
212
213
214// Save in a file
215// --------------
216void Compobj_QI::sauve(FILE* fich) const {
217
218 a_car.sauve(fich) ;
219 bbb.sauve(fich) ;
220 nphi.sauve(fich) ;
221
222 nn.sauve(fich) ;
223
224}
225
226
227// Save in a file for GYOTO input
228// -------------------------------
229
230// Redefinition of Compobj::gyoto_data
231
232void Compobj_QI::gyoto_data(const char* file_name) const {
233
234 FILE* file_out = fopen(file_name, "w") ;
235 double total_time = 0. ; // for compatibility
236 double RISCO=r_isco(0) ;
237 double RMB=r_mb(0);
238
239 fwrite_be(&total_time, sizeof(double), 1, file_out) ;
240 mp.get_mg()->sauve(file_out) ;
241 mp.sauve(file_out) ;
242 nn.sauve(file_out) ;
243 beta.sauve(file_out) ;
244 gamma.cov().sauve(file_out) ;
245 gamma.con().sauve(file_out) ;
246 kk.sauve(file_out) ;
247 fwrite_be(&RISCO, sizeof(double), 1, file_out) ;
248 fwrite_be(&RMB, sizeof(double), 1, file_out) ;
249 fclose(file_out) ;
250
251 cout << "WRITING RISCO TO GYOTO FILE : " << RISCO << endl ;
252 cout << "WRITING RMB TO GYOTO FILE : " << RMB << endl ;
253 cout << "WRITING TO GYOTO FILE - end of part " << endl ;
254}
255
256
257
258
259
260
261// Printing
262// --------
263
264ostream& Compobj_QI::operator>>(ostream& ost) const {
265
266 Compobj::operator>>(ost) ;
267
268 ost << endl << "Axisymmetric stationary compact object in quasi-isotropic coordinates (class Compobj_QI) " << endl ;
269
270 ost << "Central values of various fields : " << endl ;
271 ost << "-------------------------------- " << endl ;
272 ost << " metric coefficient A^2 : " << a_car.val_grid_point(0,0,0,0) << endl ;
273 ost << " metric coefficient B^2 : " << b_car.val_grid_point(0,0,0,0) << endl ;
274 ost << " metric coefficient N^phi : " << nphi.val_grid_point(0,0,0,0) << endl ;
275 ost << " A^2 K_{ij} K^{ij} = " << ak_car.val_grid_point(0,0,0,0) << endl << endl ;
276
277
278 double RISCO = r_isco(0, &ost) ;
279 ost << "Coordinate r at the innermost stable circular orbit (ISCO) : " <<
280 RISCO << endl ;
281 // ost << "Circumferential radius of the innermost stable circular orbit (ISCO) : " <<
282 // bbb.val_point(risco, M_PI/2, 0)*risco << endl ;
283 // ost << "Orbital frequency at the ISCO : " << f_isco(0) << endl ;
284 // ost << "Specific energy of a particle on the ISCO : " << espec_isco(0) << endl ;
285 // ost << "Specific angular momentum of a particle on the ISCO : " << lspec_isco(0) << endl ;
286
287
288 double RMB = r_mb(0, &ost) ;
289 ost << "Coordinate r at the marginally bound circular orbit (R_mb) : " << RMB << endl ;
290
291
292 // ost << "A^2 : " << a_car << endl ;
293 // ost << "B^2 : " << b_car << endl ;
294 // ost << "nphi : " << nphi << endl ;
295
296 return ost ;
297
298}
299
300// Updates the 3-metric and the shift
301
303
304 Sym_tensor gam(mp, COV, mp.get_bvect_spher()) ;
305 gam.set(1,1) = a_car ;
306 gam.set(1,2) = 0 ;
307 gam.set(1,3) = 0 ;
308 gam.set(2,2) = a_car ;
309 gam.set(2,3) = 0 ;
310 gam.set(3,3) = b_car ;
311
312 gamma = gam ;
313
314 assert(*(beta.get_triad()) == mp.get_bvect_spher()) ;
315
316 beta.set(1) = 0 ;
317 beta.set(2) = 0 ;
318 Scalar nphi_ortho(nphi) ;
319 nphi_ortho.mult_rsint() ;
320 beta.set(3) = - nphi_ortho ;
321
322 // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
323 // -------------------------------------------------
324
326
327
328 // The derived quantities are no longer up to date :
329 // -----------------------------------------------
330
331 del_deriv() ;
332
333}
334
335
336// Updates the extrinsic curvature
337// -------------------------------
338
340
341 // Special treatment for axisymmetric case:
342
343 if ( (mp.get_mg())->get_np(0) == 1) {
344
345 Scalar dnpdr = nphi.dsdr() ; // d/dr (N^phi)
346 Scalar dnpdt = nphi.dsdt() ; // d/dtheta (N^phi)
347
348 // What follows is valid only for a mapping of class Map_radial :
349 assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
350
351 dnpdr.mult_rsint() ; // multiplication by r sin(theta)
352 kk.set(1,3) = - b_car * dnpdr / (2*nn) ;
353
354 dnpdt.mult_sint() ; // multiplication by sin(theta)
355 kk.set(2,3) = - b_car * dnpdt / (2*nn) ;
356 kk.set(2,3).inc_dzpuis(2) ; // to have the same dzpuis as kk(1,3)
357
358 kk.set(1,1) = 0 ;
359 kk.set(1,2) = 0 ;
360 kk.set(2,2) = 0 ;
361 kk.set(3,3) = 0 ;
362 }
363 else {
364
365 // General case:
366
368 }
369
370 // Computation of A^2 K_{ij} K^{ij}
371 // --------------------------------
372
373 ak_car = 2 * ( kk(1,3)*kk(1,3) + kk(2,3)*kk(2,3) ) / b_car ;
374
375 del_deriv() ;
376
377}
378}
Base class for axisymmetric stationary compact objects in Quasi-Isotropic coordinates (under developm...
Definition compobj.h:280
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition compobj_QI.C:178
double * p_lspec_isco
Specific angular momentum of a particle at the ISCO.
Definition compobj.h:327
virtual void extrinsic_curvature()
Computes the extrinsic curvature and ak_car from nphi , nn and b_car .
Definition compobj_QI.C:339
double * p_r_mb
Coordinate r of the marginally bound orbit.
Definition compobj.h:328
void operator=(const Compobj_QI &)
Assignment to another Compobj_QI.
Definition compobj_QI.C:195
double * p_r_isco
Coordinate r of the ISCO.
Definition compobj.h:322
virtual void sauve(FILE *) const
Save in a file.
Definition compobj_QI.C:216
void gyoto_data(const char *file_name) const
Save in a file for GYOTO.
Definition compobj_QI.C:232
virtual double r_isco(int lmin, ostream *ost=0x0) const
Coordinate r of the innermost stable circular orbit (ISCO).
Scalar nphi
Metric coefficient .
Definition compobj.h:296
virtual void del_deriv() const
Deletes all the derived quantities.
Definition compobj_QI.C:163
Compobj_QI(Map &map_i)
Standard constructor.
Definition compobj_QI.C:86
Scalar ak_car
Scalar .
Definition compobj.h:315
virtual ~Compobj_QI()
Destructor.
Definition compobj_QI.C:152
virtual void update_metric()
Updates the 3-metric from A and B and the shift vector from .
Definition compobj_QI.C:302
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition compobj_QI.C:264
Scalar b_car
Square of the metric factor B.
Definition compobj.h:293
Scalar bbb
Metric factor B.
Definition compobj.h:290
double * p_espec_isco
Specific energy of a particle at the ISCO.
Definition compobj.h:325
Scalar a_car
Square of the metric factor A.
Definition compobj.h:287
double * p_angu_mom
Angular momentum.
Definition compobj.h:321
double * p_f_isco
Orbital frequency of the ISCO.
Definition compobj.h:323
virtual double r_mb(int lmin, ostream *ost=0x0) const
Coordinate r of the marginally bound circular orbit (R_mb).
Base class for stationary compact objects (under development).
Definition compobj.h:126
Sym_tensor kk
Extrinsic curvature tensor
Definition compobj.h:153
virtual void del_deriv() const
Deletes all the derived quantities.
Definition compobj.C:155
void operator=(const Compobj &)
Assignment to another Compobj.
Definition compobj.C:175
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition compobj.C:239
Metric gamma
3-metric
Definition compobj.h:141
Scalar nn
Lapse function N .
Definition compobj.h:135
Vector beta
Shift vector .
Definition compobj.h:138
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition compobj.C:290
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition compobj.h:132
Base class for pure radial mappings.
Definition map.h:1536
Base class for coordinate mappings.
Definition map.h:670
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
virtual void sauve(FILE *) const
Save in a file.
Definition map.C:224
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
void sauve(FILE *fd, bool save_base=false) const
Saves into a file.
Definition mg3d.C:371
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
const Scalar & dsdt() const
Returns of *this .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Scalar & dsdr() const
Returns of *this .
void mult_sint()
Multiplication by .
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor_sym.C:372
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene prototypes.
Definition app_hor.h:64