LORENE
binhor_hh.C
1/*
2 * Copyright (c) 2004-2005 Francois Limousin
3 * Jose-Luis Jaramillo
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24
25
26/*
27 * $Id: binhor_hh.C,v 1.6 2016/12/05 16:17:46 j_novak Exp $
28 * $Log: binhor_hh.C,v $
29 * Revision 1.6 2016/12/05 16:17:46 j_novak
30 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31 *
32 * Revision 1.5 2014/10/13 08:52:42 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.4 2014/10/06 15:13:01 j_novak
36 * Modified #include directives to use c++ syntax.
37 *
38 * Revision 1.3 2008/01/09 14:28:58 jl_jaramillo
39 * Improved the construction of hh1 and hh2
40 *
41 * Revision 1.2 2007/08/22 16:10:35 f_limousin
42 * Correction of many errors in binhor_hh.C
43 *
44 * Revision 1.1 2007/04/18 14:27:19 f_limousin
45 * First version
46 *
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_hh.C,v 1.6 2016/12/05 16:17:46 j_novak Exp $
50 *
51 */
52
53//standard
54#include <cstdlib>
55
56// Lorene
57#include "tensor.h"
58#include "cmp.h"
59#include "isol_hor.h"
60#include "graphique.h"
61#include "utilitaires.h"
62
63
64
65namespace Lorene {
66Sym_tensor Bin_hor::hh_Samaya_hole1 () {
67
68 //========
69 // Grid 1
70 //========
71 int nz1 = hole1.mp.get_mg()->get_nzone() ;
72 int nz2 = hole2.mp.get_mg()->get_nzone() ;
73
74 // General coordinate values
75 const Coord& xx_1 = hole1.mp.x ;
76 const Coord& yy_1 = hole1.mp.y ;
77 const Coord& zz_1 = hole1.mp.z ;
78
79 //========
80 // Grid 2
81 //========
82
83 // General coordinate values
84 const Coord& xx_2 = hole2.mp.x ;
85 const Coord& yy_2 = hole2.mp.y ;
86 const Coord& zz_2 = hole2.mp.z ;
87
88 //===================================
89 // Definition of the relevant vectors
90 //===================================
91
92 // Coordinate vector from hole 1 in the grid 1: nn1
93 //--------------------------------------------------
94 Vector rr1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
95 rr1.set(1) = xx_1 ;
96 rr1.set(2) = yy_1 ;
97 rr1.set(3) = zz_1 ;
98 rr1.std_spectral_base() ;
99
100 // Norm r1
101 Scalar r1 (hole1.mp) ;
102 r1 = hole1.mp.r ;
103 r1.std_spectral_base() ;
104 Scalar temp1 (r1) ;
105 temp1.raccord(1) ;
106 r1.set_domain(0) = temp1.domain(0) ;
107
108 // Unitary vector
109 Vector nn1 (rr1);
110 nn1 = nn1/r1 ;
111
112 for (int i=0; i<hole1.mp.get_mg()->get_nr(nz1-1); i++)
113 for (int j=0; j<hole1.mp.get_mg()->get_nt(nz1-1); j++)
114 for (int k=0; k<hole1.mp.get_mg()->get_np(nz1-1); k++)
115 for (int ind=1; ind<=3; ind++){
116 nn1.set(ind).set_grid_point(nz1-1,k,j,i) = nn1(ind).val_grid_point(1,k,j,0) ;
117 }
118
119
120 //cout << "nn1(1)" << endl << nn1(1) << endl ;
121 //des_profile(nn1(1), 0., 20., M_PI/2, 0.);
122 //des_profile(nn1(2), 0., 20., M_PI/2, M_PI/2);
123 //des_profile(nn1(3), 0., 20., 0., 0.);
124 //cout << "nn1(1)" << endl << nn1(1) << endl ;
125 //cout << "nn1(2)" << endl << nn1(2) << endl ;
126 //cout << "nn1(3)" << endl << nn1(3) << endl ;
127 //cout << "r1" << endl << r1 << endl ;
128
129
130 // Coordinate vector from hole 2 in the grid 2: nn2_2
131 //-----------------------------------------------------
132 Vector rr2_2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
133 rr2_2.set(1) = xx_2 ;
134 rr2_2.set(2) = yy_2 ;
135 rr2_2.set(3) = zz_2 ;
136 rr2_2.std_spectral_base() ;
137
138 // Norm r2_g2
139 Scalar r2_2 (hole2.mp) ;
140 r2_2 = hole2.mp.r ;
141 r2_2.std_spectral_base() ;
142 Scalar temp2 (r2_2) ;
143 temp2.raccord(1) ;
144 r2_2.set_domain(0) = temp2.domain(0) ;
145
146 // Unitary vector
147 Vector nn2_2 (rr2_2);
148 nn2_2 = nn2_2/r2_2 ;
149
150 for (int i=0; i<hole2.mp.get_mg()->get_nr(nz2-1); i++)
151 for (int j=0; j<hole2.mp.get_mg()->get_nt(nz2-1); j++)
152 for (int k=0; k<hole2.mp.get_mg()->get_np(nz2-1); k++)
153 for (int ind=1; ind<=3; ind++){
154 nn2_2.set(ind).set_grid_point(nz2-1,k,j,i) = nn2_2(ind).val_grid_point(1,k,j,0) ;
155 }
156
157 Scalar unsr1 (hole1.mp) ;
158 unsr1 = 1./hole1.mp.r ;
159 unsr1.std_spectral_base() ;
160 unsr1.raccord(1) ;
161
162 /*
163 Scalar unsr1_2 (hole2.mp) ;
164 unsr1_2.set_etat_qcq() ;
165 unsr1_2.import(unsr1) ;
166 unsr1_2.set_spectral_va().set_base(unsr1.get_spectral_va().get_base()) ;
167
168 Scalar r2sr1_2 (hole2.mp) ;
169 r2sr1_2 = r2_2*unsr1_2 ;
170 r2sr1_2.set_outer_boundary(nz2-1, 1.) ;
171
172 des_meridian(r2sr1_2, 0., 20., "r2sr1_2", 10) ;
173 arrete() ;
174 des_profile(r2sr1_2, 0., 20., M_PI/2, M_PI) ;
175 des_profile(r2sr1_2, 0., 20., M_PI/2, 0) ;
176
177 Scalar r2sr1 (hole1.mp) ;
178 r2sr1.set_etat_qcq() ;
179 r2sr1.import(r2sr1_2) ;
180 r2sr1.set_spectral_va().set_base(r2sr1_2.get_spectral_va().get_base()) ;
181
182 des_meridian(r2sr1, 0., 20., "r2sr1", 11) ;
183 arrete() ;
184 des_profile(r2sr1, 0., 20., M_PI/2, M_PI) ;
185 des_profile(r2sr1, 0., 20., M_PI/2, 0) ;
186
187 */
188
189
190 // Coordinate vector from hole 2 in the grid 1: nn2
191 //-----------------------------------------------------
192 Vector nn2 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
193 nn2_2.change_triad(hole1.mp.get_bvect_cart()) ;
194 nn2.set_etat_qcq() ;
195 for (int i=1 ; i<=3 ; i++){
196 nn2.set(i).import(nn2_2(i)) ;
197 nn2.set(i).set_spectral_va().set_base(nn2_2(i).get_spectral_va().get_base()) ;
198 }
199
200 // r2/r1
201 // -----
202 Scalar unsr2_2 (hole2.mp) ;
203 unsr2_2 = 1./hole2.mp.r ;
204 unsr2_2.std_spectral_base() ;
205 unsr2_2.raccord(1) ;
206
207 Scalar unsr2 (hole1.mp) ;
208 unsr2.set_etat_qcq() ;
209 unsr2.import(unsr2_2) ;
210 unsr2.set_spectral_va().set_base(unsr2_2.get_spectral_va().get_base()) ;
211
212 Scalar r1sr2 (unsr2*r1) ;
213 r1sr2.set_outer_boundary(nz1-1, 1.) ;
214
215 Scalar r2sr1 (1./unsr2*unsr1) ;
216 r2sr1.set_outer_boundary(nz1-1, 1.) ;
217 /*
218 des_meridian(r2sr1, 0., 20., "r2sr1", 14) ;
219 arrete() ;
220 des_profile(r2sr1, 0., 20., M_PI/2, M_PI) ;
221 des_profile(r2sr1, 0., 20., M_PI/2, 0) ;
222
223 des_meridian(1./r2sr1, 0., 20., "1./r2sr1", 12) ;
224 arrete() ;
225 des_profile(1./r2sr1, 0., 20., M_PI/2, M_PI) ;
226 des_profile(1./r2sr1, 0., 20., M_PI/2, 0) ;
227
228 des_meridian(1./r1, 0., 20., "1./r1", 13) ;
229 arrete() ;
230 des_profile(1./r1, 0., 20., M_PI/2, M_PI) ;
231 des_profile(1./r1, 0., 20., M_PI/2, 0) ;
232
233 des_meridian(1./(r1*r2sr1), 0., 20., "1./r1*r2sr1", 14) ;
234 arrete() ;
235 des_profile(1./(r1*r2sr1), 0., 20., M_PI/2, M_PI) ;
236 des_profile(1./(r1*r2sr1), 0., 20., M_PI/2, 0) ;
237 */
238
239 // Coordinate vector from hole 1 to hole 2 in the grid 1: nn12
240 //----------------------------------------------------------------
241 // Warning! Valid only in the symmetric case (for the general case it would
242 // necessary to construct this whole function as a Bin_hor function
243 Vector rr12 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
244 rr12.set(1) = hole1.mp.get_ori_x() - hole2.mp.get_ori_x() ;
245 rr12.set(2) = hole1.mp.get_ori_y() - hole2.mp.get_ori_y() ;
246 rr12.set(3) = hole1.mp.get_ori_z() - hole2.mp.get_ori_z() ;
247 rr12.std_spectral_base() ;
248
249 //Norm r12
250 Scalar r12 (hole1.mp) ;
251 r12 = sqrt( rr12(1)*rr12(1) + rr12(2)*rr12(2) + rr12(3)*rr12(3)) ;
252 r12.std_spectral_base() ;
253
254 // Unitary vector
255 Vector nn12 ( rr12 );
256 nn12 = nn12/ r12 ;
257
258
259 Scalar f_delta (hole1.mp) ;
260 Scalar f_delta_zec (hole1.mp) ;
261 Scalar f_1_1 (hole1.mp) ;
262 Scalar f_1_1_zec (hole1.mp) ;
263 Scalar f_1_12 (hole1.mp) ;
264 Scalar f_1_12_zec (hole1.mp) ;
265 Scalar f_12_12 (hole1.mp) ;
266 Scalar f_1_2 (hole1.mp) ;
267
268 f_delta.set_etat_qcq() ;
269 f_delta_zec.set_etat_qcq() ;
270 f_1_1.set_etat_qcq() ;
271 f_1_1_zec.set_etat_qcq() ;
272 f_1_12.set_etat_qcq() ;
273 f_1_12_zec.set_etat_qcq() ;
274 f_12_12.set_etat_qcq() ;
275 f_1_2.set_etat_qcq() ;
276
277 // Function exp(-(r-r_0)^2/sigma^2)
278 // --------------------------------
279
280 double r0 = hole1.mp.val_r(nz1-2, 1, 0, 0) ;
281 double sigma = 1.*r0 ;
282
283 Scalar rr (hole1.mp) ;
284 rr = hole1.mp.r ;
285
286 Scalar fexp (hole1.mp) ;
287 fexp = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
288 for (int ii=0; ii<nz1-1; ii++)
289 fexp.set_domain(ii) = 1. ;
290 fexp.set_outer_boundary(nz1-1, 0) ;
291 fexp.std_spectral_base() ;
292
293 // Conformal metric
294 //=================
295
296 // tilde{gamma}- \delta = m_1*m_2* ( f_delta \delta_{ij}
297 // + f_1_1 nn1*nn1 + f_1_12 nn1*nn12
298 // + f_12_12 nn12*nn12
299 // + f_2_2 nn2*nn2 + f_2_12 nn2*nn12
300 // + f_1_2 nn1*nn2
301
302 // f_delta
303 //--------
304 f_delta = -5.*r1/(8.*r12*r12*r12) - 15./(8.*r1*r12) +
305 5.*r1*r1*unsr2/(8.*r12*r12*r12) + 1./(r1+1./unsr2+r12)/(r1+1./unsr2+r12)*
306 (1 + r1/r12 + r12/r1 - r1sr2 - r1*r1sr2/r12 + r12*r12*unsr2/(2*r1)) +
307 1./(r1+1./unsr2+r12)*(-7./r1 + 2./r12) ;
308
309 f_delta.annule_domain(nz1-1) ;
310
311 f_delta_zec = - 15./(8.*r1*r12) + 1./(r1+1./unsr2+r12)/(r1+1./unsr2+r12)*
312 (1 + r1/r12 + r12/r1 - r1sr2 - r1*r1sr2/r12 + r12*r12*unsr2/(2*r1)) +
313 1./(r1+1./unsr2+r12)*(-7./r1 + 2./r12) ;
314 f_delta_zec += fexp*(-5.*r1/(8.*r12*r12*r12)+5.*r1*r1*unsr2/(8.*r12*r12*r12)) ;
315
316 f_delta_zec.set_outer_boundary(nz1-1, 0.) ;
317 for (int i=0 ;i<nz1-1 ; i++){
318 f_delta_zec.annule_domain(i) ;
319 }
320
321 f_delta = f_delta + f_delta_zec ;
322
323 /*
324 des_meridian(f_delta, 0., 20., "f_delta", 10) ;
325 arrete() ;
326 des_profile(f_delta, 0., 20., M_PI/2, M_PI) ;
327 des_profile(f_delta, 0., 20., M_PI/2, 0) ;
328 des_profile(f_delta, 0., 20., 0, M_PI) ;
329 des_coupe_z (f_delta, 0., 2) ;
330 des_coupe_z (f_delta, 0., 3) ;
331 des_coupe_z (f_delta, 0., 4) ;
332 des_coupe_z (f_delta, 0., 5) ;
333 */
334
335 // f_1_1
336 //------
337 f_1_1 = r1/(8.*r12*r12*r12) + 11./(8.*r1*r12) -
338 1./(8.*r1*unsr2*unsr2*r12*r12*r12) + 7./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) +
339 7./r1/(r1+1./unsr2+r12) ;
340 f_1_1.annule_domain(nz1-1) ;
341
342 f_1_1_zec = 11./(8.*r1*r12) + 7./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) +
343 7./r1/(r1+1./unsr2+r12) ;
344 f_1_1_zec += fexp*(r1/(8.*r12*r12*r12)-1./(8.*r1*unsr2*unsr2*r12*r12*r12)) ;
345 f_1_1_zec.set_outer_boundary(nz1-1, 0.) ;
346
347 for (int i=0 ; i<nz1-1 ; i++){
348 f_1_1_zec.annule_domain(i) ;
349 }
350
351 f_1_1 = f_1_1 + f_1_1_zec ;
352
353 /*
354 des_meridian(f_1_1, 0., 20., "f_1_1", 14) ;
355 arrete() ;
356 des_profile(f_1_1, 0., 20., M_PI/2, M_PI) ;
357 des_profile(f_1_1, 0., 20., M_PI/2, 0) ;
358 des_profile(f_1_1, 0., 20., 0, M_PI) ;
359 des_coupe_z (f_1_1, 0., 2) ;
360 des_coupe_z (f_1_1, 0., 3) ;
361 des_coupe_z (f_1_1, 0., 4) ;
362 des_coupe_z (f_1_1, 0., 5) ;
363 */
364
365 // f_1_12
366 //------
367 f_1_12 = - 7./(2*r12*r12) + 8./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) ;
368 f_1_12.annule_domain(nz1-1) ;
369
370 f_1_12_zec = 8./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) ;
371 f_1_12_zec += fexp*(- 7./(2*r12*r12)) ;
372 f_1_12_zec.set_outer_boundary(nz1-1, 0.) ;
373
374 for (int i=0 ; i<nz1-1 ; i++){
375 f_1_12_zec.annule_domain(i) ;
376 }
377
378 f_1_12 = f_1_12 + f_1_12_zec ;
379
380 /*
381 des_meridian(f_1_12, 0., 40., "f_1_12", 15) ;
382 arrete() ;
383 des_profile(f_1_12, 0., 20., M_PI/2, M_PI) ;
384 des_profile(f_1_12, 0., 20., M_PI/2, 0) ;
385 des_profile(f_1_12, 0., 20., 0, M_PI) ;
386 des_coupe_z (f_1_12, 0., 2) ;
387 des_coupe_z (f_1_12, 0., 3) ;
388 des_coupe_z (f_1_12, 0., 4) ;
389 des_coupe_z (f_1_12, 0., 5) ;
390 */
391
392 // f_12_12
393 //-------
394 f_12_12 = (-4./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) - // facteur 2 ???????
395 4./r12/(r1+1./unsr2+r12)) ;
396 f_12_12.set_outer_boundary(nz1-1, 0.) ;
397
398 /*
399 des_meridian(f_12_12, 0., 40., "f_12_12", 15) ;
400 arrete() ;
401 des_profile(f_12_12, 0., 20., M_PI/2, M_PI) ;
402 des_profile(f_12_12, 0., 20., M_PI/2, 0) ;
403 des_profile(f_12_12, 0., 20., 0, M_PI) ;
404 des_coupe_z (f_12_12, 0., 2) ;
405 des_coupe_z (f_12_12, 0., 3) ;
406 des_coupe_z (f_12_12, 0., 4) ;
407 des_coupe_z (f_12_12, 0., 5) ;
408 */
409
410 // f_1_2
411 //-------
412 f_1_2 = 11./(r1+1./unsr2+r12)/(r1+1./unsr2+r12); // facteur 2 ???????
413 f_1_2.set_outer_boundary(nz1-1, 0.) ;
414
415 /*
416 des_meridian(f_1_2, 0., 40., "f_1_1", 15) ;
417 arrete() ;
418 des_profile(f_1_2, 0., 20., M_PI/2, M_PI) ;
419 des_profile(f_1_2, 0., 20., M_PI/2, 0) ;
420 des_profile(f_1_2, 0., 20., 0, M_PI) ;
421 des_coupe_z (f_1_2, 0., 2) ;
422 des_coupe_z (f_1_2, 0., 3) ;
423 des_coupe_z (f_1_2, 0., 4) ;
424 des_coupe_z (f_1_2, 0., 5) ;
425 */
426
427 // First part of the correction metric (needed to be complemented by the (1 <-> 2) term
428
429 Sym_tensor hh_temp(hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
430
431 for (int i=1 ; i<= 3 ; i++){
432 for (int j=i ; j<= 3 ; j++){
433 hh_temp.set(i,j) = f_delta * hole1.ff.con()(i,j) + f_1_1 * nn1(i)*nn1(j)
434 + f_1_12 * 0.5 *(nn1(i) * nn12(j) + nn1(j) * nn12(i))
435 + f_12_12 * nn12(i)*nn12(j)
436 + f_1_2 * 0.5*(nn1(i)*nn2(j) + nn1(j)*nn2(i) ) ;
437 }
438 }
439 /*
440 des_meridian(hh_temp, 0., 20., "hh_temp", 25) ;
441 arrete() ;
442 for (int i=1 ; i<= 3 ; i++)
443 for (int j=i ; j<= 3 ; j++){
444 des_profile(hh_temp(i,j), 0., 20., M_PI/2, M_PI) ;
445 des_profile(hh_temp(i,j), 0., 20., M_PI/2, 0) ;
446 des_profile(hh_temp(i,j), 0., 20., 0, M_PI) ;
447 des_coupe_z (hh_temp(i,j), 0., 5) ;
448 }
449 */
450
451 return hh_temp ;
452
453}
454
455
456Sym_tensor Bin_hor::hh_Samaya_hole2() {
457
458
459 //========
460 // Grid 1
461 //========
462 int nz1 = hole1.mp.get_mg()->get_nzone() ;
463 int nz2 = hole2.mp.get_mg()->get_nzone() ;
464
465 // General coordinate values
466 const Coord& xx_1 = hole1.mp.x ;
467 const Coord& yy_1 = hole1.mp.y ;
468 const Coord& zz_1 = hole1.mp.z ;
469
470 //========
471 // Grid 2
472 //========
473
474 // General coordinate values
475 const Coord& xx_2 = hole2.mp.x ;
476 const Coord& yy_2 = hole2.mp.y ;
477 const Coord& zz_2 = hole2.mp.z ;
478
479
480 //===================================
481 // Definition of the relevant vectors
482 //===================================
483
484 // Coordinate vector from hole 2 in the grid 2: nn2
485 //--------------------------------------------------
486 Vector rr2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
487 rr2.set(1) = xx_2 ;
488 rr2.set(2) = yy_2 ;
489 rr2.set(3) = zz_2 ;
490 rr2.std_spectral_base() ;
491
492 // Norm r2
493 Scalar r2 (hole2.mp) ;
494 r2 = hole1.mp.r ;
495 r2.std_spectral_base() ;
496 Scalar temp2 (r2) ;
497 temp2.raccord(1) ;
498 r2.set_domain(0) = temp2.domain(0) ;
499
500 // Unitary vector
501 Vector nn2 (rr2);
502 nn2 = nn2/r2 ;
503
504 for (int i=0; i<hole2.mp.get_mg()->get_nr(nz2-1); i++)
505 for (int j=0; j<hole2.mp.get_mg()->get_nt(nz2-1); j++)
506 for (int k=0; k<hole2.mp.get_mg()->get_np(nz2-1); k++)
507 for (int ind=1; ind<=3; ind++){
508 nn2.set(ind).set_grid_point(nz2-1,k,j,i) = nn2(ind).val_grid_point(1,k,j,0) ;
509 }
510
511 // Coordinate vector from hole 1 in the grid 1: nn1_1
512 //-----------------------------------------------------
513 Vector rr1_1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
514 rr1_1.set(1) = xx_1 ;
515 rr1_1.set(2) = yy_1 ;
516 rr1_1.set(3) = zz_1 ;
517 rr1_1.std_spectral_base() ;
518
519 // Norm r1_g1
520 Scalar r1_1 (hole1.mp) ;
521 r1_1 = hole1.mp.r ;
522 r1_1.std_spectral_base() ;
523 Scalar temp1 (r1_1) ;
524 temp1.raccord(1) ;
525 r1_1.set_domain(0) = temp1.domain(0) ;
526
527 // Unitary vector
528 Vector nn1_1 (rr1_1);
529 nn1_1 = nn1_1/r1_1 ;
530
531 for (int i=0; i<hole1.mp.get_mg()->get_nr(nz1-1); i++)
532 for (int j=0; j<hole1.mp.get_mg()->get_nt(nz1-1); j++)
533 for (int k=0; k<hole1.mp.get_mg()->get_np(nz1-1); k++)
534 for (int ind=1; ind<=3; ind++){
535 nn1_1.set(ind).set_grid_point(nz1-1,k,j,i) = nn1_1(ind).val_grid_point(1,k,j,0) ;
536 }
537
538 Scalar unsr2 (hole2.mp) ;
539 unsr2 = 1./hole2.mp.r ;
540 unsr2.std_spectral_base() ;
541 unsr2.raccord(1) ;
542
543 // Coordinate vector from hole 1 in the grid 2: nn1
544 //-----------------------------------------------------
545 Vector nn1 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
546 nn1_1.change_triad(hole2.mp.get_bvect_cart()) ;
547 nn1.set_etat_qcq() ;
548 for (int i=1 ; i<=3 ; i++){
549 nn1.set(i).import(nn1_1(i)) ;
550 nn1.set(i).set_spectral_va().set_base(nn1_1(i).get_spectral_va().get_base()) ;
551 }
552
553 // r1/r2
554 // -----
555 Scalar unsr1_1 (hole1.mp) ;
556 unsr1_1 = 1./hole1.mp.r ;
557 unsr1_1.std_spectral_base() ;
558 unsr1_1.raccord(1) ;
559
560 Scalar unsr1 (hole2.mp) ;
561 unsr1.set_etat_qcq() ;
562 unsr1.import(unsr1_1) ;
563 unsr1.set_spectral_va().set_base(unsr1_1.get_spectral_va().get_base()) ;
564
565 Scalar r2sr1 (unsr1*r2) ;
566 r2sr1.set_outer_boundary(nz2-1, 1.) ;
567
568 Scalar r1sr2 (1./unsr1*unsr2) ;
569 r1sr2.set_outer_boundary(nz2-1, 1.) ;
570
571 // Coordinate vector from hole 2 to hole 1 in the grid 2: nn21
572 //----------------------------------------------------------------
573 // Warning! Valid only in the symmetric case (for the general case it would
574 // necessary to construct this whole function as a Bin_hor function
575 Vector rr21 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
576 rr21.set(1) = hole2.mp.get_ori_x() - hole1.mp.get_ori_x() ;
577 rr21.set(2) = hole2.mp.get_ori_y() - hole1.mp.get_ori_y() ;
578 rr21.set(3) = hole2.mp.get_ori_z() - hole1.mp.get_ori_z() ;
579 rr21.std_spectral_base() ;
580
581 //Norm r21
582 Scalar r21 (hole2.mp) ;
583 r21 = sqrt( rr21(1)*rr21(1) + rr21(2)*rr21(2) + rr21(3)*rr21(3)) ;
584 r21.std_spectral_base() ;
585
586 // Unitary vector
587 Vector nn21 ( rr21 );
588 nn21 = nn21/ r21 ;
589
590
591 Scalar f_delta (hole2.mp) ;
592 Scalar f_delta_zec (hole2.mp) ;
593 Scalar f_2_2 (hole2.mp) ;
594 Scalar f_2_2_zec (hole2.mp) ;
595 Scalar f_2_21 (hole2.mp) ;
596 Scalar f_2_21_zec (hole2.mp) ;
597 Scalar f_21_21 (hole2.mp) ;
598 Scalar f_2_1 (hole2.mp) ;
599
600 f_delta.set_etat_qcq() ;
601 f_delta_zec.set_etat_qcq() ;
602 f_2_2.set_etat_qcq() ;
603 f_2_2_zec.set_etat_qcq() ;
604 f_2_21.set_etat_qcq() ;
605 f_2_21_zec.set_etat_qcq() ;
606 f_21_21.set_etat_qcq() ;
607 f_2_1.set_etat_qcq() ;
608
609 // Function exp(-(r-r_0)^2/sigma^2)
610 // --------------------------------
611
612 double r0 = hole2.mp.val_r(nz2-2, 1, 0, 0) ;
613 double sigma = 1.*r0 ;
614
615 Scalar rr (hole2.mp) ;
616 rr = hole2.mp.r ;
617
618 Scalar fexp (hole2.mp) ;
619 fexp = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
620 for (int ii=0; ii<nz2-1; ii++)
621 fexp.set_domain(ii) = 1. ;
622 fexp.set_outer_boundary(nz2-1, 0) ;
623 fexp.std_spectral_base() ;
624
625 // Conformal metric
626 //=================
627
628 // tilde{gamma}- \delta = m_1*m_2* ( f_delta \delta_{ij}
629 // + f_2_2 nn2*nn2 + f_2_21 nn2*nn21
630 // + f_21_21 nn21*nn21
631 // + f_1_1 nn1*nn1 + f_1_21 nn1*nn21
632 // + f_2_1 nn2*nn1
633
634 // f_delta
635 //--------
636 f_delta = -5.*r2/(8.*r21*r21*r21) - 15./(8.*r2*r21) +
637 5.*r2*r2*unsr1/(8.*r21*r21*r21) + 1./(r2+1./unsr1+r21)/(r2+1./unsr1+r21)*
638 (1 + r2/r21 + r21/r2 - r2sr1 - r2*r2sr1/r21 + r21*r21*unsr1/(2*r2)) +
639 1./(r2+1./unsr1+r21)*(-7./r2 + 2./r21) ;
640
641 f_delta.annule_domain(nz2-1) ;
642
643 f_delta_zec = - 15./(8.*r2*r21) + 1./(r2+1./unsr1+r21)/(r2+1./unsr1+r21)*
644 (1 + r2/r21 + r21/r2 - r2sr1 - r2*r2sr1/r21 + r21*r21*unsr1/(2*r2)) +
645 1./(r2+1./unsr1+r21)*(-7./r2 + 2./r21) ;
646 f_delta_zec += fexp*(-5.*r2/(8.*r21*r21*r21)+5.*r2*r2*unsr1/(8.*r21*r21*r21)) ;
647
648 f_delta_zec.set_outer_boundary(nz2-1, 0.) ;
649 for (int i=0 ;i<nz2-1 ; i++){
650 f_delta_zec.annule_domain(i) ;
651 }
652
653 f_delta = f_delta + f_delta_zec ;
654
655 /*
656 des_meridian(f_delta, 0., 20., "f_delta", 10) ;
657 arrete() ;
658 des_profile(f_delta, 0., 20., M_PI/2, M_PI) ;
659 des_profile(f_delta, 0., 20., M_PI/2, 0) ;
660 des_profile(f_delta, 0., 20., 0, M_PI) ;
661 des_coupe_z (f_delta, 0., 2) ;
662 des_coupe_z (f_delta, 0., 3) ;
663 des_coupe_z (f_delta, 0., 4) ;
664 des_coupe_z (f_delta, 0., 5) ;
665 */
666
667 // f_2_2
668 //------
669 f_2_2 = r2/(8.*r21*r21*r21) + 11./(8.*r2*r21) -
670 1./(8.*r2*unsr1*unsr1*r21*r21*r21) + 7./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) +
671 7./r2/(r2+1./unsr1+r21) ;
672 f_2_2.annule_domain(nz2-1) ;
673
674 f_2_2_zec = 11./(8.*r2*r21) + 7./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) +
675 7./r2/(r2+1./unsr1+r21) ;
676 f_2_2_zec += fexp*(r2/(8.*r21*r21*r21)-1./(8.*r2*unsr1*unsr1*r21*r21*r21)) ;
677 f_2_2_zec.set_outer_boundary(nz2-1, 0.) ;
678
679 for (int i=0 ; i<nz2-1 ; i++){
680 f_2_2_zec.annule_domain(i) ;
681 }
682
683 f_2_2 = f_2_2 + f_2_2_zec ;
684
685 /*
686 des_meridian(f_2_2, 0., 20., "f_2_2", 14) ;
687 arrete() ;
688 des_profile(f_2_2, 0., 20., M_PI/2, M_PI) ;
689 des_profile(f_2_2, 0., 20., M_PI/2, 0) ;
690 des_profile(f_2_2, 0., 20., 0, M_PI) ;
691 des_coupe_z (f_2_2, 0., 2) ;
692 des_coupe_z (f_2_2, 0., 3) ;
693 des_coupe_z (f_2_2, 0., 4) ;
694 des_coupe_z (f_2_2, 0., 5) ;
695 */
696
697 // f_2_21
698 //------
699 f_2_21 = - 7./(2*r21*r21) + 8./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) ;
700 f_2_21.annule_domain(nz2-1) ;
701
702 f_2_21_zec = 8./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) ;
703 f_2_21_zec += fexp*(- 7./(2*r21*r21)) ;
704 f_2_21_zec.set_outer_boundary(nz2-1, 0.) ;
705
706 for (int i=0 ; i<nz2-1 ; i++){
707 f_2_21_zec.annule_domain(i) ;
708 }
709
710 f_2_21 = f_2_21 + f_2_21_zec ;
711
712 /*
713 des_meridian(f_2_21, 0., 40., "f_2_21", 15) ;
714 arrete() ;
715 des_profile(f_2_21, 0., 20., M_PI/2, M_PI) ;
716 des_profile(f_2_21, 0., 20., M_PI/2, 0) ;
717 des_profile(f_2_21, 0., 20., 0, M_PI) ;
718 des_coupe_z (f_2_21, 0., 2) ;
719 des_coupe_z (f_2_21, 0., 3) ;
720 des_coupe_z (f_2_21, 0., 4) ;
721 des_coupe_z (f_2_21, 0., 5) ;
722 */
723
724 // f_21_21
725 //-------
726 f_21_21 = (-4./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) - // facteur 2 ???????
727 4./r21/(r2+1./unsr1+r21)) ;
728 f_21_21.set_outer_boundary(nz2-1, 0.) ;
729
730 /*
731 des_meridian(f_21_21, 0., 40., "f_21_21", 15) ;
732 arrete() ;
733 des_profile(f_21_21, 0., 20., M_PI/2, M_PI) ;
734 des_profile(f_21_21, 0., 20., M_PI/2, 0) ;
735 des_profile(f_21_21, 0., 20., 0, M_PI) ;
736 des_coupe_z (f_21_21, 0., 2) ;
737 des_coupe_z (f_21_21, 0., 3) ;
738 des_coupe_z (f_21_21, 0., 4) ;
739 des_coupe_z (f_21_21, 0., 5) ;
740 */
741
742 // f_2_1
743 //-------
744 f_2_1 = 11./(r2+1./unsr1+r21)/(r2+1./unsr1+r21); // facteur 2 ???????
745 f_2_1.set_outer_boundary(nz2-1, 0.) ;
746
747 /*
748 des_meridian(f_2_1, 0., 40., "f_2_1", 15) ;
749 arrete() ;
750 des_profile(f_2_1, 0., 20., M_PI/2, M_PI) ;
751 des_profile(f_2_1, 0., 20., M_PI/2, 0) ;
752 des_profile(f_2_1, 0., 20., 0, M_PI) ;
753 des_coupe_z (f_2_1, 0., 2) ;
754 des_coupe_z (f_2_1, 0., 3) ;
755 des_coupe_z (f_2_1, 0., 4) ;
756 des_coupe_z (f_2_1, 0., 5) ;
757 */
758
759 // First part of the correction metric (needed to be complemented by the (1 <-> 2) term
760
761 Sym_tensor hh_temp(hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
762
763 for (int i=1 ; i<= 3 ; i++){
764 for (int j=i ; j<= 3 ; j++){
765 hh_temp.set(i,j) = f_delta * hole2.ff.con()(i,j) + f_2_2 * nn2(i)*nn2(j)
766 - f_2_21 * 0.5 *(nn2(i) * nn21(j) + nn2(j) * nn21(i))
767 + f_21_21 * nn21(i)*nn21(j)
768 + f_2_1 * 0.5*(nn2(i)*nn1(j) + nn2(j)*nn1(i) ) ;
769 }
770 }
771 /*
772 des_meridian(hh_temp, 0., 20., "hh_temp", 25) ;
773 arrete() ;
774 for (int i=1 ; i<= 3 ; i++)
775 for (int j=i ; j<= 3 ; j++){
776 des_profile(hh_temp(i,j), 0., 20., M_PI/2, M_PI) ;
777 des_profile(hh_temp(i,j), 0., 20., M_PI/2, 0) ;
778 des_profile(hh_temp(i,j), 0., 20., 0, M_PI) ;
779 des_coupe_z (hh_temp(i,j), 0., 5) ;
780 }
781 */
782
783 return hh_temp ;
784
785
786
787}
788
789void Bin_hor::set_hh_Samaya() {
790
791 Sym_tensor hh1 ( hh_Samaya_hole1() ) ;
792 Sym_tensor hh2 ( hh_Samaya_hole2() ) ;
793
794 // Definition of the surface
795 // -------------------------
796
797 Cmp surface_un (hole1.mp) ;
798 surface_un = pow(hole1.mp.r, 2.)-pow(hole1.get_radius(), 2.) ;
799 surface_un.annule(hole1.mp.get_mg()->get_nzone()-1) ;
800 surface_un.std_base_scal() ;
801
802 Cmp surface_deux (hole2.mp) ;
803 surface_deux = pow(hole2.mp.r, 2.)-pow(hole2.get_radius(), 2.) ;
804 surface_deux.annule(hole1.mp.get_mg()->get_nzone()-1) ;
805 surface_deux.std_base_scal() ;
806 /*
807 double ta = 12 ;
808 for (int i=1 ; i<= 3 ; i++)
809 for (int j=i ; j<= 3 ; j++){
810 Cmp dessin_un (hh1(i,j)) ;
811 dessin_un.annule(0) ;
812
813 Cmp dessin_deux (hh2(i,j)) ;
814 dessin_deux.annule(0) ;
815
816 des_coupe_bin_z (dessin_un, dessin_deux, 0,
817 -ta, ta, -ta, ta, "hh(1,1)", &surface_un, &surface_deux,
818 false, 15, 300, 300) ;
819 }
820 */
821
822 // Importation
823 // ----------------
824
825 Sym_tensor hh2_1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
826 hh2_1.set_etat_qcq() ;
827 Sym_tensor hh1_2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
828 hh1_2.set_etat_qcq() ;
829
830 /*
831 Scalar temp (hh1(1,1)) ;
832 temp.annule_domain(0) ;
833 des_profile(temp, 0., 4., M_PI/2, M_PI) ;
834 des_profile(temp, 0., 4., M_PI/2, 0) ;
835 des_profile(temp, 0., 4., 0, M_PI) ;
836 des_coupe_z (temp, 0., 5) ;
837 temp.raccord(1) ;
838 des_profile(temp, 0., 4., M_PI/2, M_PI) ;
839 des_profile(temp, 0., 4., M_PI/2, 0) ;
840 des_profile(temp, 0., 4., 0, M_PI) ;
841 des_coupe_z (temp, 0., 5) ;
842 */
843
844 /*
845 for (int i=1 ; i<= 3 ; i++)
846 for (int j=i ; j<= 3 ; j++){
847 des_profile(hh1(i,j), 0., 20., M_PI/2, M_PI) ;
848 des_profile(hh1(i,j), 0., 20., M_PI/2, 0) ;
849 des_profile(hh1(i,j), 0., 20., 0, M_PI) ;
850 des_coupe_z (hh1(i,j), 0., 5) ;
851 }
852 */
853 for (int i=1 ; i<=3 ; i++)
854 for (int j=i ; j<=3 ; j++){
855 hh1.set(i,j).raccord(1) ;
856 hh2.set(i,j).raccord(1) ;
857 }
858 /*
859 for (int i=1 ; i<= 3 ; i++)
860 for (int j=i ; j<= 3 ; j++){
861 des_profile(hh1(i,j), 0., 20., M_PI/2, M_PI) ;
862 des_profile(hh1(i,j), 0., 20., M_PI/2, 0) ;
863 des_profile(hh1(i,j), 0., 20., 0, M_PI) ;
864 des_coupe_z (hh1(i,j), 0., 5) ;
865 }
866 */
867
868 hh2.change_triad(hole1.mp.get_bvect_cart()) ;
869 for (int i=1 ; i<=3 ; i++){
870 for (int j=i ; j<=3 ; j++){
871 hh2_1.set(i,j).import(hh2(i,j)) ;
872 hh2_1.set(i,j).set_spectral_va().set_base(hh2(i,j).get_spectral_va().get_base()) ;
873 }
874 }
875 hh2.change_triad(hole2.mp.get_bvect_cart()) ;
876
877 hh1.change_triad(hole2.mp.get_bvect_cart()) ;
878 for (int i=1 ; i<=3 ; i++){
879 for (int j=i ; j<=3 ; j++){
880 hh1_2.set(i,j).import(hh1(i,j)) ;
881 hh1_2.set(i,j).set_spectral_va().set_base(hh1(i,j).get_spectral_va().get_base()) ;
882 }
883 }
884 hh1.change_triad(hole1.mp.get_bvect_cart()) ;
885
886 double m1, m2 ;
887 m1 = pow(hole1.area_hor()/(16.*M_PI) + hole1.ang_mom_hor()/hole1.radius, 0.5) ;
888 m2 = pow(hole2.area_hor()/(16.*M_PI) + hole2.ang_mom_hor()/hole2.radius, 0.5) ;
889
890
891 hh1 = hh1 + hh2_1 ;
892 hh2 = hh2 + hh1_2 ;
893
894 cout << hole1.mp.r << endl ;
895 cout << hole1.mp.phi << endl ;
896 cout << hole1.mp.tet << endl ;
897
898
899 //des_meridian(hh1, 0., 20., "hh1 cart", 20) ;
900 for (int i=1 ; i<= 3 ; i++)
901 for (int j=i ; j<= 3 ; j++){
902 // des_profile(hh1(i,j), 0., 20., M_PI/2, M_PI) ;
903 //des_profile(hh1(i,j), 0., 20., M_PI/2, 0) ;
904 //des_profile(hh1(i,j), 0., 20., 0, M_PI) ;
905 des_coupe_z (hh1(i,j), 0., 5) ;
906 }
907
908 hh1.change_triad(hole1.mp.get_bvect_spher()) ;
909 hh2.change_triad(hole2.mp.get_bvect_spher()) ;
910
911 hole1.hh = m1*m2* hh1 ;
912 hole2.hh = m1*m2* hh2 ;
913
914 Metric tgam_1 ( hole1.ff.con() + hh1 ) ;
915 Metric tgam_2 ( hole2.ff.con() + hh2 ) ;
916
917 hole1.tgam = tgam_1 ;
918 hole2.tgam = tgam_2 ;
919
920
921 des_meridian(hh1, 0., 20., "hh1", 0) ;
922
923
924}
925}
Lorene::Bin_hor::hole1
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Lorene::Bin_hor::hole2
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
Lorene::Bin_hor::set_hh_Samaya
void set_hh_Samaya()
Calculation of the Post-Newtonian correction to .
Definition binhor_hh.C:789
Lorene::Bin_hor::hh_Samaya_hole2
Sym_tensor hh_Samaya_hole2()
Calculation of the hole2 part of the Post-Newtonian correction to .
Definition binhor_hh.C:456
Lorene::Bin_hor::hh_Samaya_hole1
Sym_tensor hh_Samaya_hole1()
Calculation of the hole1 part of the Post-Newtonian correction to .
Definition binhor_hh.C:66
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
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::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:351
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Metric
Metric for tensor calculation.
Definition metric.h:90
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition scalar_raccord.C:65
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:359
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:790
Lorene::Scalar::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:621
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:610
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:607
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:690
Lorene::Scalar::domain
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition scalar.h:631
Lorene::Scalar::set_outer_boundary
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition scalar_manip.C:321
Lorene::Scalar::import
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition scalar_import.C:71
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:226
Lorene::Valeur::get_base
const Base_val & get_base() const
Return the bases for spectral expansions (member base ).
Definition valeur.h:490
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base ).
Definition valeur.C:813
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:319
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition vector_change_triad.C:78
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:299
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::Tensor::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tensor_change_triad.C:80
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:675
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:490
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:663
Lorene
Lorene prototypes.
Definition app_hor.h:67