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