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 char 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 
62 namespace Lorene {
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()) ;
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 
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()) ;
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 
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 
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  }
873 
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  }
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 
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 }
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
void set_hh_Samaya()
Calculation of the Post-Newtonian correction to .
Definition: binhor_hh.C:786
Sym_tensor hh_Samaya_hole2()
Calculation of the hole2 part of the Post-Newtonian correction to .
Definition: binhor_hh.C:453
Sym_tensor hh_Samaya_hole1()
Calculation of the hole1 part of the Post-Newtonian correction to .
Definition: binhor_hh.C:63
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Definition: map_af_radius.C:96
double get_ori_z() const
Returns the z coordinate of the origin.
Definition: map.h:772
Coord y
y coordinate centered on the grid
Definition: map.h:727
Coord r
r coordinate centered on the grid
Definition: map.h:718
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:770
Coord tet
coordinate centered on the grid
Definition: map.h:719
Coord z
z coordinate centered on the grid
Definition: map.h:728
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
Coord x
x coordinate centered on the grid
Definition: map.h:726
Coord phi
coordinate centered on the grid
Definition: map.h:720
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
Metric for tensor calculation.
Definition: metric.h:90
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
double & 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
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition: scalar.h:625
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
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition: scalar_import.C:68
double radius
Radius of the horizon in LORENE's units.
Definition: isol_hor.h:906
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
double area_hor() const
Area of the horizon.
Definition: single_param.C:88
double ang_mom_hor() const
Angular momentum (modulo)
Definition: single_param.C:107
double get_radius() const
Returns the radius of the horizon.
Definition: isol_hor.h:1046
Sym_tensor hh
Deviation metric.
Definition: isol_hor.h:983
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
Tensor field of valence 1.
Definition: vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
void des_coupe_z(const Scalar &uu, double z0, int nzdes, const char *title=0x0, const Scalar *defsurf=0x0, double zoom=1.2, bool draw_bound=true, int ncour=15, int nx=100, int ny=100)
Draws isocontour lines of a Scalar in a plane Z=constant.
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Lorene prototypes.
Definition: app_hor.h:64