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