LORENE
bhole_kij.C
1/*
2 * Copyright (c) 2001 Philippe Grandclement
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char bhole_kij_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_kij.C,v 1.6 2014/10/13 08:52:40 j_novak Exp $" ;
24
25/*
26 * $Id: bhole_kij.C,v 1.6 2014/10/13 08:52:40 j_novak Exp $
27 * $Log: bhole_kij.C,v $
28 * Revision 1.6 2014/10/13 08:52:40 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.5 2014/10/06 15:12:58 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.4 2005/08/29 15:10:14 p_grandclement
35 * Addition of things needed :
36 * 1) For BBH with different masses
37 * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
38 * WORKING YET !!!
39 *
40 * Revision 1.3 2004/05/27 07:10:22 p_grandclement
41 * Correction of some shadowed variables
42 *
43 * Revision 1.2 2002/10/16 14:36:33 j_novak
44 * Reorganization of #include instructions of standard C++, in order to
45 * use experimental version 3 of gcc.
46 *
47 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
48 * LORENE
49 *
50 * Revision 1.6 2001/05/07 09:12:02 phil
51 * *** empty log message ***
52 *
53 * Revision 1.5 2001/04/30 09:22:52 phil
54 * cas omega = 0
55 *
56 * Revision 1.4 2001/04/27 11:44:01 phil
57 * correction devant assurer la symetrie entre les deux trous
58 *
59 * Revision 1.3 2001/04/26 12:04:44 phil
60 * *** empty log message ***
61 *
62 * Revision 1.2 2001/04/25 15:54:30 phil
63 * corrections diverses rien de bien mechant
64 *
65 * Revision 1.1 2001/04/25 15:10:23 phil
66 * Initial revision
67 *
68 *
69 * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_kij.C,v 1.6 2014/10/13 08:52:40 j_novak Exp $
70 *
71 */
72
73
74//standard
75#include <cstdlib>
76#include <cmath>
77
78// Lorene
79#include "nbr_spx.h"
80#include "tenseur.h"
81#include "bhole.h"
82#include "proto.h"
83#include "utilitaires.h"
84#include "graphique.h"
85
86namespace Lorene {
87//calcul de kij total. (la regularisation ayant ete faite)
88void Bhole_binaire::fait_tkij () {
89
90 hole1.tkij_tot.set_etat_qcq() ;
91 hole2.tkij_tot.set_etat_qcq() ;
92
93 // On construit a_ij a partir du shift ...
94 // taij tot doit etre nul sur les deux horizons.
95 hole1.fait_taij_auto () ;
96 hole2.fait_taij_auto () ;
97
98 // On trouve les trucs du compagnon
99 hole1.taij_comp.set_etat_qcq() ;
100 hole2.taij_comp.set_etat_qcq() ;
101
102 Tenseur sans_dz_un (hole1.taij_auto) ;
103 sans_dz_un.dec2_dzpuis() ;
104 Tenseur sans_dz_deux (hole2.taij_auto) ;
105 sans_dz_deux.dec2_dzpuis() ;
106
107 // ON DOIT VERIFIER SI LES DEUX Aij sont alignes ou non !
108 // Les bases des deux vecteurs doivent etre alignees ou non alignees :
109 double orientation_un = hole1.taij_auto.get_mp()->get_rot_phi() ;
110 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
111 double orientation_deux = hole2.taij_auto.get_mp()->get_rot_phi() ;
112 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
113 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
114
115 // Les importations :
116 if (hole2.taij_auto.get_etat() == ETATZERO)
117 hole1.taij_comp.set_etat_zero() ;
118 else {
119 hole1.taij_comp.set(0, 0).import_asymy(sans_dz_deux(0, 0)) ;
120 hole1.taij_comp.set(0, 1).import_symy(sans_dz_deux(0, 1)) ;
121 hole1.taij_comp.set(0, 2).import_asymy(same_orient*sans_dz_deux(0, 2)) ;
122 hole1.taij_comp.set(1, 1).import_asymy(sans_dz_deux(1, 1)) ;
123 hole1.taij_comp.set(1, 2).import_symy(same_orient*sans_dz_deux(1, 2)) ;
124 hole1.taij_comp.set(2, 2).import_asymy(sans_dz_deux(2, 2)) ;
125 }
126
127 if (hole1.taij_auto.get_etat() == ETATZERO)
128 hole2.taij_comp.set_etat_zero() ;
129 else {
130 hole2.taij_comp.set(0, 0).import_asymy(sans_dz_un(0, 0)) ;
131 hole2.taij_comp.set(0, 1).import_symy(sans_dz_un(0, 1)) ;
132 hole2.taij_comp.set(0, 2).import_asymy(same_orient*sans_dz_un(0, 2)) ;
133 hole2.taij_comp.set(1, 1).import_asymy(sans_dz_un(1, 1)) ;
134 hole2.taij_comp.set(1, 2).import_symy(same_orient*sans_dz_un(1, 2)) ;
135 hole2.taij_comp.set(2, 2).import_asymy(sans_dz_un(2, 2)) ;
136 }
137
138 hole1.taij_comp.set_std_base() ;
139 hole2.taij_comp.set_std_base() ;
140 hole1.taij_comp.inc2_dzpuis() ;
141 hole2.taij_comp.inc2_dzpuis() ;
142
143 // Et enfin les trucs totaux...
144 hole1.taij_tot = hole1.taij_auto + hole1.taij_comp ;
145 hole2.taij_tot = hole2.taij_auto + hole2.taij_comp ;
146
147 if ((hole1.taij_tot.get_etat() == ETATZERO) &&
148 (hole2.taij_tot.get_etat() == ETATZERO)) {
149
150 hole1.tkij_tot.set_etat_zero() ;
151 hole1.tkij_auto.set_etat_zero() ;
152 hole2.tkij_tot.set_etat_zero() ;
153 hole2.tkij_auto.set_etat_zero() ;
154 }
155 else {
156 int nz_un = hole1.mp.get_mg()->get_nzone() ;
157 int nz_deux = hole2.mp.get_mg()->get_nzone() ;
158
159 Cmp ntot_un (hole1.n_tot()) ;
160 ntot_un = division_xpun (ntot_un, 0) ;
161 ntot_un.raccord(1) ;
162
163 Cmp ntot_deux (hole2.n_tot()) ;
164 ntot_deux = division_xpun (ntot_deux, 0) ;
165 ntot_deux.raccord(1) ;
166
167 // Boucle sur les composantes :
168 for (int lig = 0 ; lig<3 ; lig++)
169 for (int col = lig ; col<3 ; col++) {
170
171 // Le sens d orientation
172 int ind = 1 ;
173 if (lig !=2)
174 ind *= -1 ;
175 if (col != 2)
176 ind *= -1 ;
177 if (same_orient == 1)
178 ind = 1 ;
179
180 // Pres de H1 :
181 Cmp auxi_un (hole1.taij_tot(lig, col)/2.) ;
182 auxi_un.dec2_dzpuis() ;
183 auxi_un = division_xpun (auxi_un, 0) ;
184 auxi_un = auxi_un / ntot_un ;
185 auxi_un.raccord(1) ;
186
187 // Pres de H2 :
188 Cmp auxi_deux (hole2.taij_tot(lig, col)/2.) ;
189 auxi_deux.dec2_dzpuis() ;
190 auxi_deux = division_xpun (auxi_deux, 0) ;
191 auxi_deux = auxi_deux / ntot_deux ;
192 auxi_deux.raccord(1) ;
193
194 // copie :
195 Cmp copie_un (hole1.taij_tot(lig, col)) ;
196 copie_un.dec2_dzpuis() ;
197
198 Cmp copie_deux (hole2.taij_tot(lig, col)) ;
199 copie_deux.dec2_dzpuis() ;
200
201 // Double les rayons limites :
202 double lim_un = hole1.mp.get_alpha()[1] + hole1.mp.get_beta()[1] ;
203 double lim_deux = hole2.mp.get_alpha()[1] + hole2.mp.get_beta()[1] ;
204
205 Mtbl xabs_un (hole1.mp.xa) ;
206 Mtbl yabs_un (hole1.mp.ya) ;
207 Mtbl zabs_un (hole1.mp.za) ;
208
209 Mtbl xabs_deux (hole2.mp.xa) ;
210 Mtbl yabs_deux (hole2.mp.ya) ;
211 Mtbl zabs_deux (hole2.mp.za) ;
212
213 double xabs, yabs, zabs, air, theta, phi ;
214
215 // On boucle sur les autres zones :
216 for (int l=2 ; l<nz_un ; l++) {
217
218 int nr = hole1.mp.get_mg()->get_nr (l) ;
219
220 if (l==nz_un-1)
221 nr -- ;
222
223 int np = hole1.mp.get_mg()->get_np (l) ;
224 int nt = hole1.mp.get_mg()->get_nt (l) ;
225
226 for (int k=0 ; k<np ; k++)
227 for (int j=0 ; j<nt ; j++)
228 for (int i=0 ; i<nr ; i++) {
229
230 xabs = xabs_un (l, k, j, i) ;
231 yabs = yabs_un (l, k, j, i) ;
232 zabs = zabs_un (l, k, j, i) ;
233
234 // les coordonnees du point en 2 :
235 hole2.mp.convert_absolute
236 (xabs, yabs, zabs, air, theta, phi) ;
237
238 if (air >= lim_deux)
239 // On est loin du trou 2 : pas de pb :
240 auxi_un.set(l, k, j, i) =
241 copie_un(l, k, j, i) / ntot_un (l, k, j, i)/2. ;
242 else
243 // On est pres du trou deux :
244 auxi_un.set(l, k, j, i) =
245 ind * auxi_deux.val_point (air, theta, phi) ;
246 }
247
248 // Cas infini :
249 if (l==nz_un-1)
250 for (int k=0 ; k<np ; k++)
251 for (int j=0 ; j<nt ; j++)
252 auxi_un.set(nz_un-1, k, j, nr-1) = 0 ;
253 }
254
255 // Le second trou :
256 for (int l=2 ; l<nz_deux ; l++) {
257
258 int nr = hole2.mp.get_mg()->get_nr (l) ;
259
260 if (l==nz_deux-1)
261 nr -- ;
262
263 int np = hole2.mp.get_mg()->get_np (l) ;
264 int nt = hole2.mp.get_mg()->get_nt (l) ;
265
266 for (int k=0 ; k<np ; k++)
267 for (int j=0 ; j<nt ; j++)
268 for (int i=0 ; i<nr ; i++) {
269
270 xabs = xabs_deux (l, k, j, i) ;
271 yabs = yabs_deux (l, k, j, i) ;
272 zabs = zabs_deux (l, k, j, i) ;
273
274 // les coordonnees du point en 1 :
275 hole1.mp.convert_absolute
276 (xabs, yabs, zabs, air, theta, phi) ;
277
278 if (air >= lim_un)
279 // On est loin du trou 1 : pas de pb :
280 auxi_deux.set(l, k, j, i) =
281 copie_deux(l, k, j, i) / ntot_deux (l, k, j, i) /2.;
282 else
283 // On est pres du trou deux :
284 auxi_deux.set(l, k, j, i) =
285 ind * auxi_un.val_point (air, theta, phi) ;
286 }
287 // Cas infini :
288 if (l==nz_deux-1)
289 for (int k=0 ; k<np ; k++)
290 for (int j=0 ; j<nt ; j++)
291 auxi_deux.set(nz_deux-1, k, j, nr-1) = 0 ;
292 }
293
294 auxi_un.inc2_dzpuis() ;
295 auxi_deux.inc2_dzpuis() ;
296
297 hole1.tkij_tot.set(lig, col) = auxi_un ;
298 hole2.tkij_tot.set(lig, col) = auxi_deux ;
299 }
300
301 hole1.tkij_auto.set_etat_qcq() ;
302 hole2.tkij_auto.set_etat_qcq() ;
303
304 for (int lig=0 ; lig<3 ; lig++)
305 for (int col=lig ; col<3 ; col++) {
306 hole1.tkij_auto.set(lig, col) = hole1.tkij_tot(lig, col)*
307 hole1.decouple ;
308 hole2.tkij_auto.set(lig, col) = hole2.tkij_tot(lig, col)*
309 hole2.decouple ;
310 }
311 }
312}
313
314void Bhole_binaire::fait_decouple () {
315
316 int nz_un = hole1.mp.get_mg()->get_nzone() ;
317 int nz_deux = hole2.mp.get_mg()->get_nzone() ;
318
319 // On determine R_limite (pour le moment en tout cas...) :
320 double distance = hole1.mp.get_ori_x() - hole2.mp.get_ori_x() ;
321 double lim_un = distance/2. ;
322 double lim_deux = distance/2. ;
323 double int_un = distance/6. ;
324 double int_deux = distance/6. ;
325
326 // Les fonctions de base
327 Cmp fonction_f_un (hole1.mp) ;
328 fonction_f_un = 0.5*pow(
329 cos((hole1.mp.r-int_un)*M_PI/2./(lim_un-int_un)), 2.)+0.5 ;
330 fonction_f_un.std_base_scal();
331
332 Cmp fonction_g_un (hole1.mp) ;
333 fonction_g_un = 0.5*pow
334 (sin((hole1.mp.r-int_un)*M_PI/2./(lim_un-int_un)), 2.) ;
335 fonction_g_un.std_base_scal();
336
337 Cmp fonction_f_deux (hole2.mp) ;
338 fonction_f_deux = 0.5*pow(
339 cos((hole2.mp.r-int_deux)*M_PI/2./(lim_deux-int_deux)), 2.)+0.5 ;
340 fonction_f_deux.std_base_scal();
341
342 Cmp fonction_g_deux (hole2.mp) ;
343 fonction_g_deux = 0.5*pow(
344 sin((hole2.mp.r-int_deux)*M_PI/2./(lim_un-int_deux)), 2.) ;
345 fonction_g_deux.std_base_scal();
346
347 // Les fonctions totales :
348 Cmp decouple_un (hole1.mp) ;
349 decouple_un.allocate_all() ;
350 Cmp decouple_deux (hole2.mp) ;
351 decouple_deux.allocate_all() ;
352
353 Mtbl xabs_un (hole1.mp.xa) ;
354 Mtbl yabs_un (hole1.mp.ya) ;
355 Mtbl zabs_un (hole1.mp.za) ;
356
357 Mtbl xabs_deux (hole2.mp.xa) ;
358 Mtbl yabs_deux (hole2.mp.ya) ;
359 Mtbl zabs_deux (hole2.mp.za) ;
360
361 double xabs, yabs, zabs, air_un, air_deux, theta, phi ;
362
363 // On boucle sur les autres zones :
364 for (int l=0 ; l<nz_un ; l++) {
365 int nr = hole1.mp.get_mg()->get_nr (l) ;
366
367 if (l==nz_un-1)
368 nr -- ;
369
370 int np = hole1.mp.get_mg()->get_np (l) ;
371 int nt = hole1.mp.get_mg()->get_nt (l) ;
372
373 for (int k=0 ; k<np ; k++)
374 for (int j=0 ; j<nt ; j++)
375 for (int i=0 ; i<nr ; i++) {
376
377 xabs = xabs_un (l, k, j, i) ;
378 yabs = yabs_un (l, k, j, i) ;
379 zabs = zabs_un (l, k, j, i) ;
380
381 // les coordonnees du point :
382 hole1.mp.convert_absolute
383 (xabs, yabs, zabs, air_un, theta, phi) ;
384 hole2.mp.convert_absolute
385 (xabs, yabs, zabs, air_deux, theta, phi) ;
386
387 if (air_un <= lim_un)
388 if (air_un < int_un)
389 decouple_un.set(l, k, j, i) = 1 ;
390 else
391 // pres du trou un :
392 decouple_un.set(l, k, j, i) =
393 fonction_f_un (l, k, j, i) ;
394 else
395 if (air_deux <= lim_deux)
396 if (air_deux < int_deux)
397 decouple_un.set(l, k, j, i) = 0 ;
398 else
399 // On est pres du trou deux :
400 decouple_un.set(l, k, j, i) =
401 fonction_g_deux.val_point (air_deux, theta, phi) ;
402
403 else
404 // On est loin des deux trous :
405 decouple_un.set(l, k, j, i) = 0.5 ;
406 }
407
408 // Cas infini :
409 if (l==nz_un-1)
410 for (int k=0 ; k<np ; k++)
411 for (int j=0 ; j<nt ; j++)
412 decouple_un.set(nz_un-1, k, j, nr) = 0.5 ;
413 }
414
415 for (int l=0 ; l<nz_deux ; l++) {
416
417 int nr = hole2.mp.get_mg()->get_nr (l) ;
418
419 if (l==nz_deux-1)
420 nr -- ;
421
422 int np = hole2.mp.get_mg()->get_np (l) ;
423 int nt = hole2.mp.get_mg()->get_nt (l) ;
424
425 for (int k=0 ; k<np ; k++)
426 for (int j=0 ; j<nt ; j++)
427 for (int i=0 ; i<nr ; i++) {
428
429 xabs = xabs_deux (l, k, j, i) ;
430 yabs = yabs_deux (l, k, j, i) ;
431 zabs = zabs_deux (l, k, j, i) ;
432
433 // les coordonnees du point :
434 hole1.mp.convert_absolute
435 (xabs, yabs, zabs, air_un, theta, phi) ;
436 hole2.mp.convert_absolute
437 (xabs, yabs, zabs, air_deux, theta, phi) ;
438
439 if (air_deux <= lim_deux)
440 if (air_deux < int_deux)
441 decouple_deux.set(l, k, j, i) = 1 ;
442 else
443 // pres du trou deux :
444 decouple_deux.set(l, k, j, i) =
445 fonction_f_deux (l, k, j, i) ;
446 else
447 if (air_un <= lim_un)
448 if (air_un < int_un)
449 decouple_deux.set(l, k, j, i) = 0 ;
450 else
451 // On est pres du trou un :
452 decouple_deux.set(l, k, j, i) =
453 fonction_g_un.val_point (air_un, theta, phi) ;
454
455 else
456 // On est loin des deux trous :
457 decouple_deux.set(l, k, j, i) = 0.5 ;
458 }
459
460 // Cas infini :
461 if (l==nz_deux-1)
462 for (int k=0 ; k<np ; k++)
463 for (int j=0 ; j<nt ; j++)
464 decouple_deux.set(nz_deux-1, k, j, nr) = 0.5 ;
465 }
466
467 hole1.decouple = decouple_un ;
468 hole2.decouple = decouple_deux ;
469}
470}
Lorene::Bhole_binaire::hole1
Bhole hole1
Black hole one.
Definition bhole.h:762
Lorene::Bhole_binaire::hole2
Bhole hole2
Black hole two.
Definition bhole.h:763
Lorene::Bhole_binaire::fait_decouple
void fait_decouple()
Calculates {tt decouple} which is used to obtain tkij_auto by the formula : tkij_auto = decouple * tk...
Definition bhole_kij.C:314
Lorene::Bhole_binaire::fait_tkij
void fait_tkij()
Calculation af the extrinsic curvature tensor.
Definition bhole_kij.C:88
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::allocate_all
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition cmp.C:323
Lorene::Cmp::std_base_scal
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
Lorene::Cmp::inc2_dzpuis
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition cmp_r_manip.C:192
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Lorene::Cmp::dec2_dzpuis
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition cmp_r_manip.C:180
Lorene::Cmp::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point , by means of the spectral...
Definition cmp.C:732
Lorene::Cmp::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition cmp_raccord.C:170
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::dec2_dzpuis
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720