binhor_hh.C
/*
* Copyright (c) 2004-2005 Francois Limousin
* Jose-Luis Jaramillo
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: binhor_hh.C,v 1.5 2014/10/13 08:52:42 j_novak Exp $
* $Log: binhor_hh.C,v $
* Revision 1.5 2014/10/13 08:52:42 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2014/10/06 15:13:01 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.3 2008/01/09 14:28:58 jl_jaramillo
* Improved the construction of hh1 and hh2
*
* Revision 1.2 2007/08/22 16:10:35 f_limousin
* Correction of many errors in binhor_hh.C
*
* Revision 1.1 2007/04/18 14:27:19 f_limousin
* First version
*
*
*
* $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_hh.C,v 1.5 2014/10/13 08:52:42 j_novak Exp $
*
*/
//standard
#include <cstdlib>
// Lorene
#include "tensor.h"
#include "cmp.h"
#include "isol_hor.h"
#include "graphique.h"
#include "utilitaires.h"
namespace Lorene {
Sym_tensor Bin_hor::hh_Samaya_hole1 () {
//========
// Grid 1
//========
int nz1 = hole1.mp.get_mg()->get_nzone() ;
int nz2 = hole2.mp.get_mg()->get_nzone() ;
// General coordinate values
const Coord& xx_1 = hole1.mp.x ;
const Coord& yy_1 = hole1.mp.y ;
const Coord& zz_1 = hole1.mp.z ;
//========
// Grid 2
//========
// General coordinate values
const Coord& xx_2 = hole2.mp.x ;
const Coord& yy_2 = hole2.mp.y ;
const Coord& zz_2 = hole2.mp.z ;
//===================================
// Definition of the relevant vectors
//===================================
// Coordinate vector from hole 1 in the grid 1: nn1
//--------------------------------------------------
Vector rr1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
rr1.set(1) = xx_1 ;
rr1.set(2) = yy_1 ;
rr1.set(3) = zz_1 ;
rr1.std_spectral_base() ;
// Norm r1
Scalar r1 (hole1.mp) ;
r1 = hole1.mp.r ;
r1.std_spectral_base() ;
Scalar temp1 (r1) ;
temp1.raccord(1) ;
r1.set_domain(0) = temp1.domain(0) ;
// Unitary vector
Vector nn1 (rr1);
nn1 = nn1/r1 ;
for (int i=0; i<hole1.mp.get_mg()->get_nr(nz1-1); i++)
for (int j=0; j<hole1.mp.get_mg()->get_nt(nz1-1); j++)
for (int k=0; k<hole1.mp.get_mg()->get_np(nz1-1); k++)
for (int ind=1; ind<=3; ind++){
nn1.set(ind).set_grid_point(nz1-1,k,j,i) = nn1(ind).val_grid_point(1,k,j,0) ;
}
//cout << "nn1(1)" << endl << nn1(1) << endl ;
//des_profile(nn1(1), 0., 20., M_PI/2, 0.);
//des_profile(nn1(2), 0., 20., M_PI/2, M_PI/2);
//des_profile(nn1(3), 0., 20., 0., 0.);
//cout << "nn1(1)" << endl << nn1(1) << endl ;
//cout << "nn1(2)" << endl << nn1(2) << endl ;
//cout << "nn1(3)" << endl << nn1(3) << endl ;
//cout << "r1" << endl << r1 << endl ;
// Coordinate vector from hole 2 in the grid 2: nn2_2
//-----------------------------------------------------
Vector rr2_2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
rr2_2.set(1) = xx_2 ;
rr2_2.set(2) = yy_2 ;
rr2_2.set(3) = zz_2 ;
rr2_2.std_spectral_base() ;
// Norm r2_g2
Scalar r2_2 (hole2.mp) ;
r2_2 = hole2.mp.r ;
r2_2.std_spectral_base() ;
Scalar temp2 (r2_2) ;
temp2.raccord(1) ;
r2_2.set_domain(0) = temp2.domain(0) ;
// Unitary vector
Vector nn2_2 (rr2_2);
nn2_2 = nn2_2/r2_2 ;
for (int i=0; i<hole2.mp.get_mg()->get_nr(nz2-1); i++)
for (int j=0; j<hole2.mp.get_mg()->get_nt(nz2-1); j++)
for (int k=0; k<hole2.mp.get_mg()->get_np(nz2-1); k++)
for (int ind=1; ind<=3; ind++){
nn2_2.set(ind).set_grid_point(nz2-1,k,j,i) = nn2_2(ind).val_grid_point(1,k,j,0) ;
}
Scalar unsr1 (hole1.mp) ;
unsr1 = 1./hole1.mp.r ;
unsr1.std_spectral_base() ;
unsr1.raccord(1) ;
/*
Scalar unsr1_2 (hole2.mp) ;
unsr1_2.set_etat_qcq() ;
unsr1_2.import(unsr1) ;
unsr1_2.set_spectral_va().set_base(unsr1.get_spectral_va().get_base()) ;
Scalar r2sr1_2 (hole2.mp) ;
r2sr1_2 = r2_2*unsr1_2 ;
r2sr1_2.set_outer_boundary(nz2-1, 1.) ;
des_meridian(r2sr1_2, 0., 20., "r2sr1_2", 10) ;
arrete() ;
des_profile(r2sr1_2, 0., 20., M_PI/2, M_PI) ;
des_profile(r2sr1_2, 0., 20., M_PI/2, 0) ;
Scalar r2sr1 (hole1.mp) ;
r2sr1.set_etat_qcq() ;
r2sr1.import(r2sr1_2) ;
r2sr1.set_spectral_va().set_base(r2sr1_2.get_spectral_va().get_base()) ;
des_meridian(r2sr1, 0., 20., "r2sr1", 11) ;
arrete() ;
des_profile(r2sr1, 0., 20., M_PI/2, M_PI) ;
des_profile(r2sr1, 0., 20., M_PI/2, 0) ;
*/
// Coordinate vector from hole 2 in the grid 1: nn2
//-----------------------------------------------------
Vector nn2 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
nn2_2.change_triad(hole1.mp.get_bvect_cart()) ;
nn2.set_etat_qcq() ;
for (int i=1 ; i<=3 ; i++){
nn2.set(i).import(nn2_2(i)) ;
nn2.set(i).set_spectral_va().set_base(nn2_2(i).get_spectral_va().get_base()) ;
}
// r2/r1
// -----
Scalar unsr2_2 (hole2.mp) ;
unsr2_2 = 1./hole2.mp.r ;
unsr2_2.std_spectral_base() ;
unsr2_2.raccord(1) ;
Scalar unsr2 (hole1.mp) ;
unsr2.set_etat_qcq() ;
unsr2.import(unsr2_2) ;
unsr2.set_spectral_va().set_base(unsr2_2.get_spectral_va().get_base()) ;
Scalar r1sr2 (unsr2*r1) ;
r1sr2.set_outer_boundary(nz1-1, 1.) ;
Scalar r2sr1 (1./unsr2*unsr1) ;
r2sr1.set_outer_boundary(nz1-1, 1.) ;
/*
des_meridian(r2sr1, 0., 20., "r2sr1", 14) ;
arrete() ;
des_profile(r2sr1, 0., 20., M_PI/2, M_PI) ;
des_profile(r2sr1, 0., 20., M_PI/2, 0) ;
des_meridian(1./r2sr1, 0., 20., "1./r2sr1", 12) ;
arrete() ;
des_profile(1./r2sr1, 0., 20., M_PI/2, M_PI) ;
des_profile(1./r2sr1, 0., 20., M_PI/2, 0) ;
des_meridian(1./r1, 0., 20., "1./r1", 13) ;
arrete() ;
des_profile(1./r1, 0., 20., M_PI/2, M_PI) ;
des_profile(1./r1, 0., 20., M_PI/2, 0) ;
des_meridian(1./(r1*r2sr1), 0., 20., "1./r1*r2sr1", 14) ;
arrete() ;
des_profile(1./(r1*r2sr1), 0., 20., M_PI/2, M_PI) ;
des_profile(1./(r1*r2sr1), 0., 20., M_PI/2, 0) ;
*/
// Coordinate vector from hole 1 to hole 2 in the grid 1: nn12
//----------------------------------------------------------------
// Warning! Valid only in the symmetric case (for the general case it would
// necessary to construct this whole function as a Bin_hor function
Vector rr12 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
rr12.set(1) = hole1.mp.get_ori_x() - hole2.mp.get_ori_x() ;
rr12.set(2) = hole1.mp.get_ori_y() - hole2.mp.get_ori_y() ;
rr12.set(3) = hole1.mp.get_ori_z() - hole2.mp.get_ori_z() ;
rr12.std_spectral_base() ;
//Norm r12
Scalar r12 (hole1.mp) ;
r12 = sqrt( rr12(1)*rr12(1) + rr12(2)*rr12(2) + rr12(3)*rr12(3)) ;
r12.std_spectral_base() ;
// Unitary vector
Vector nn12 ( rr12 );
nn12 = nn12/ r12 ;
Scalar f_delta (hole1.mp) ;
Scalar f_delta_zec (hole1.mp) ;
Scalar f_1_1 (hole1.mp) ;
Scalar f_1_1_zec (hole1.mp) ;
Scalar f_1_12 (hole1.mp) ;
Scalar f_1_12_zec (hole1.mp) ;
Scalar f_12_12 (hole1.mp) ;
Scalar f_1_2 (hole1.mp) ;
f_delta.set_etat_qcq() ;
f_delta_zec.set_etat_qcq() ;
f_1_1.set_etat_qcq() ;
f_1_1_zec.set_etat_qcq() ;
f_1_12.set_etat_qcq() ;
f_1_12_zec.set_etat_qcq() ;
f_12_12.set_etat_qcq() ;
f_1_2.set_etat_qcq() ;
// Function exp(-(r-r_0)^2/sigma^2)
// --------------------------------
double r0 = hole1.mp.val_r(nz1-2, 1, 0, 0) ;
double sigma = 1.*r0 ;
Scalar rr (hole1.mp) ;
rr = hole1.mp.r ;
Scalar fexp (hole1.mp) ;
fexp = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
for (int ii=0; ii<nz1-1; ii++)
fexp.set_domain(ii) = 1. ;
fexp.set_outer_boundary(nz1-1, 0) ;
fexp.std_spectral_base() ;
// Conformal metric
//=================
// tilde{gamma}- \delta = m_1*m_2* ( f_delta \delta_{ij}
// + f_1_1 nn1*nn1 + f_1_12 nn1*nn12
// + f_12_12 nn12*nn12
// + f_2_2 nn2*nn2 + f_2_12 nn2*nn12
// + f_1_2 nn1*nn2
// f_delta
//--------
f_delta = -5.*r1/(8.*r12*r12*r12) - 15./(8.*r1*r12) +
5.*r1*r1*unsr2/(8.*r12*r12*r12) + 1./(r1+1./unsr2+r12)/(r1+1./unsr2+r12)*
(1 + r1/r12 + r12/r1 - r1sr2 - r1*r1sr2/r12 + r12*r12*unsr2/(2*r1)) +
1./(r1+1./unsr2+r12)*(-7./r1 + 2./r12) ;
f_delta.annule_domain(nz1-1) ;
f_delta_zec = - 15./(8.*r1*r12) + 1./(r1+1./unsr2+r12)/(r1+1./unsr2+r12)*
(1 + r1/r12 + r12/r1 - r1sr2 - r1*r1sr2/r12 + r12*r12*unsr2/(2*r1)) +
1./(r1+1./unsr2+r12)*(-7./r1 + 2./r12) ;
f_delta_zec += fexp*(-5.*r1/(8.*r12*r12*r12)+5.*r1*r1*unsr2/(8.*r12*r12*r12)) ;
f_delta_zec.set_outer_boundary(nz1-1, 0.) ;
for (int i=0 ;i<nz1-1 ; i++){
f_delta_zec.annule_domain(i) ;
}
f_delta = f_delta + f_delta_zec ;
/*
des_meridian(f_delta, 0., 20., "f_delta", 10) ;
arrete() ;
des_profile(f_delta, 0., 20., M_PI/2, M_PI) ;
des_profile(f_delta, 0., 20., M_PI/2, 0) ;
des_profile(f_delta, 0., 20., 0, M_PI) ;
des_coupe_z (f_delta, 0., 2) ;
des_coupe_z (f_delta, 0., 3) ;
des_coupe_z (f_delta, 0., 4) ;
des_coupe_z (f_delta, 0., 5) ;
*/
// f_1_1
//------
f_1_1 = r1/(8.*r12*r12*r12) + 11./(8.*r1*r12) -
1./(8.*r1*unsr2*unsr2*r12*r12*r12) + 7./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) +
7./r1/(r1+1./unsr2+r12) ;
f_1_1.annule_domain(nz1-1) ;
f_1_1_zec = 11./(8.*r1*r12) + 7./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) +
7./r1/(r1+1./unsr2+r12) ;
f_1_1_zec += fexp*(r1/(8.*r12*r12*r12)-1./(8.*r1*unsr2*unsr2*r12*r12*r12)) ;
f_1_1_zec.set_outer_boundary(nz1-1, 0.) ;
for (int i=0 ; i<nz1-1 ; i++){
f_1_1_zec.annule_domain(i) ;
}
f_1_1 = f_1_1 + f_1_1_zec ;
/*
des_meridian(f_1_1, 0., 20., "f_1_1", 14) ;
arrete() ;
des_profile(f_1_1, 0., 20., M_PI/2, M_PI) ;
des_profile(f_1_1, 0., 20., M_PI/2, 0) ;
des_profile(f_1_1, 0., 20., 0, M_PI) ;
des_coupe_z (f_1_1, 0., 2) ;
des_coupe_z (f_1_1, 0., 3) ;
des_coupe_z (f_1_1, 0., 4) ;
des_coupe_z (f_1_1, 0., 5) ;
*/
// f_1_12
//------
f_1_12 = - 7./(2*r12*r12) + 8./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) ;
f_1_12.annule_domain(nz1-1) ;
f_1_12_zec = 8./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) ;
f_1_12_zec += fexp*(- 7./(2*r12*r12)) ;
f_1_12_zec.set_outer_boundary(nz1-1, 0.) ;
for (int i=0 ; i<nz1-1 ; i++){
f_1_12_zec.annule_domain(i) ;
}
f_1_12 = f_1_12 + f_1_12_zec ;
/*
des_meridian(f_1_12, 0., 40., "f_1_12", 15) ;
arrete() ;
des_profile(f_1_12, 0., 20., M_PI/2, M_PI) ;
des_profile(f_1_12, 0., 20., M_PI/2, 0) ;
des_profile(f_1_12, 0., 20., 0, M_PI) ;
des_coupe_z (f_1_12, 0., 2) ;
des_coupe_z (f_1_12, 0., 3) ;
des_coupe_z (f_1_12, 0., 4) ;
des_coupe_z (f_1_12, 0., 5) ;
*/
// f_12_12
//-------
f_12_12 = (-4./(r1+1./unsr2+r12)/(r1+1./unsr2+r12) - // facteur 2 ???????
4./r12/(r1+1./unsr2+r12)) ;
f_12_12.set_outer_boundary(nz1-1, 0.) ;
/*
des_meridian(f_12_12, 0., 40., "f_12_12", 15) ;
arrete() ;
des_profile(f_12_12, 0., 20., M_PI/2, M_PI) ;
des_profile(f_12_12, 0., 20., M_PI/2, 0) ;
des_profile(f_12_12, 0., 20., 0, M_PI) ;
des_coupe_z (f_12_12, 0., 2) ;
des_coupe_z (f_12_12, 0., 3) ;
des_coupe_z (f_12_12, 0., 4) ;
des_coupe_z (f_12_12, 0., 5) ;
*/
// f_1_2
//-------
f_1_2 = 11./(r1+1./unsr2+r12)/(r1+1./unsr2+r12); // facteur 2 ???????
f_1_2.set_outer_boundary(nz1-1, 0.) ;
/*
des_meridian(f_1_2, 0., 40., "f_1_1", 15) ;
arrete() ;
des_profile(f_1_2, 0., 20., M_PI/2, M_PI) ;
des_profile(f_1_2, 0., 20., M_PI/2, 0) ;
des_profile(f_1_2, 0., 20., 0, M_PI) ;
des_coupe_z (f_1_2, 0., 2) ;
des_coupe_z (f_1_2, 0., 3) ;
des_coupe_z (f_1_2, 0., 4) ;
des_coupe_z (f_1_2, 0., 5) ;
*/
// First part of the correction metric (needed to be complemented by the (1 <-> 2) term
Sym_tensor hh_temp(hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
for (int i=1 ; i<= 3 ; i++){
for (int j=i ; j<= 3 ; j++){
hh_temp.set(i,j) = f_delta * hole1.ff.con()(i,j) + f_1_1 * nn1(i)*nn1(j)
+ f_1_12 * 0.5 *(nn1(i) * nn12(j) + nn1(j) * nn12(i))
+ f_12_12 * nn12(i)*nn12(j)
+ f_1_2 * 0.5*(nn1(i)*nn2(j) + nn1(j)*nn2(i) ) ;
}
}
/*
des_meridian(hh_temp, 0., 20., "hh_temp", 25) ;
arrete() ;
for (int i=1 ; i<= 3 ; i++)
for (int j=i ; j<= 3 ; j++){
des_profile(hh_temp(i,j), 0., 20., M_PI/2, M_PI) ;
des_profile(hh_temp(i,j), 0., 20., M_PI/2, 0) ;
des_profile(hh_temp(i,j), 0., 20., 0, M_PI) ;
des_coupe_z (hh_temp(i,j), 0., 5) ;
}
*/
return hh_temp ;
}
Sym_tensor Bin_hor::hh_Samaya_hole2() {
//========
// Grid 1
//========
int nz1 = hole1.mp.get_mg()->get_nzone() ;
int nz2 = hole2.mp.get_mg()->get_nzone() ;
// General coordinate values
const Coord& xx_1 = hole1.mp.x ;
const Coord& yy_1 = hole1.mp.y ;
const Coord& zz_1 = hole1.mp.z ;
//========
// Grid 2
//========
// General coordinate values
const Coord& xx_2 = hole2.mp.x ;
const Coord& yy_2 = hole2.mp.y ;
const Coord& zz_2 = hole2.mp.z ;
//===================================
// Definition of the relevant vectors
//===================================
// Coordinate vector from hole 2 in the grid 2: nn2
//--------------------------------------------------
Vector rr2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
rr2.set(1) = xx_2 ;
rr2.set(2) = yy_2 ;
rr2.set(3) = zz_2 ;
rr2.std_spectral_base() ;
// Norm r2
Scalar r2 (hole2.mp) ;
r2 = hole1.mp.r ;
r2.std_spectral_base() ;
Scalar temp2 (r2) ;
temp2.raccord(1) ;
r2.set_domain(0) = temp2.domain(0) ;
// Unitary vector
Vector nn2 (rr2);
nn2 = nn2/r2 ;
for (int i=0; i<hole2.mp.get_mg()->get_nr(nz2-1); i++)
for (int j=0; j<hole2.mp.get_mg()->get_nt(nz2-1); j++)
for (int k=0; k<hole2.mp.get_mg()->get_np(nz2-1); k++)
for (int ind=1; ind<=3; ind++){
nn2.set(ind).set_grid_point(nz2-1,k,j,i) = nn2(ind).val_grid_point(1,k,j,0) ;
}
// Coordinate vector from hole 1 in the grid 1: nn1_1
//-----------------------------------------------------
Vector rr1_1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
rr1_1.set(1) = xx_1 ;
rr1_1.set(2) = yy_1 ;
rr1_1.set(3) = zz_1 ;
rr1_1.std_spectral_base() ;
// Norm r1_g1
Scalar r1_1 (hole1.mp) ;
r1_1 = hole1.mp.r ;
r1_1.std_spectral_base() ;
Scalar temp1 (r1_1) ;
temp1.raccord(1) ;
r1_1.set_domain(0) = temp1.domain(0) ;
// Unitary vector
Vector nn1_1 (rr1_1);
nn1_1 = nn1_1/r1_1 ;
for (int i=0; i<hole1.mp.get_mg()->get_nr(nz1-1); i++)
for (int j=0; j<hole1.mp.get_mg()->get_nt(nz1-1); j++)
for (int k=0; k<hole1.mp.get_mg()->get_np(nz1-1); k++)
for (int ind=1; ind<=3; ind++){
nn1_1.set(ind).set_grid_point(nz1-1,k,j,i) = nn1_1(ind).val_grid_point(1,k,j,0) ;
}
Scalar unsr2 (hole2.mp) ;
unsr2 = 1./hole2.mp.r ;
unsr2.std_spectral_base() ;
unsr2.raccord(1) ;
// Coordinate vector from hole 1 in the grid 2: nn1
//-----------------------------------------------------
Vector nn1 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
nn1_1.change_triad(hole2.mp.get_bvect_cart()) ;
nn1.set_etat_qcq() ;
for (int i=1 ; i<=3 ; i++){
nn1.set(i).import(nn1_1(i)) ;
nn1.set(i).set_spectral_va().set_base(nn1_1(i).get_spectral_va().get_base()) ;
}
// r1/r2
// -----
Scalar unsr1_1 (hole1.mp) ;
unsr1_1 = 1./hole1.mp.r ;
unsr1_1.std_spectral_base() ;
unsr1_1.raccord(1) ;
Scalar unsr1 (hole2.mp) ;
unsr1.set_etat_qcq() ;
unsr1.import(unsr1_1) ;
unsr1.set_spectral_va().set_base(unsr1_1.get_spectral_va().get_base()) ;
Scalar r2sr1 (unsr1*r2) ;
r2sr1.set_outer_boundary(nz2-1, 1.) ;
Scalar r1sr2 (1./unsr1*unsr2) ;
r1sr2.set_outer_boundary(nz2-1, 1.) ;
// Coordinate vector from hole 2 to hole 1 in the grid 2: nn21
//----------------------------------------------------------------
// Warning! Valid only in the symmetric case (for the general case it would
// necessary to construct this whole function as a Bin_hor function
Vector rr21 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
rr21.set(1) = hole2.mp.get_ori_x() - hole1.mp.get_ori_x() ;
rr21.set(2) = hole2.mp.get_ori_y() - hole1.mp.get_ori_y() ;
rr21.set(3) = hole2.mp.get_ori_z() - hole1.mp.get_ori_z() ;
rr21.std_spectral_base() ;
//Norm r21
Scalar r21 (hole2.mp) ;
r21 = sqrt( rr21(1)*rr21(1) + rr21(2)*rr21(2) + rr21(3)*rr21(3)) ;
r21.std_spectral_base() ;
// Unitary vector
Vector nn21 ( rr21 );
nn21 = nn21/ r21 ;
Scalar f_delta (hole2.mp) ;
Scalar f_delta_zec (hole2.mp) ;
Scalar f_2_2 (hole2.mp) ;
Scalar f_2_2_zec (hole2.mp) ;
Scalar f_2_21 (hole2.mp) ;
Scalar f_2_21_zec (hole2.mp) ;
Scalar f_21_21 (hole2.mp) ;
Scalar f_2_1 (hole2.mp) ;
f_delta.set_etat_qcq() ;
f_delta_zec.set_etat_qcq() ;
f_2_2.set_etat_qcq() ;
f_2_2_zec.set_etat_qcq() ;
f_2_21.set_etat_qcq() ;
f_2_21_zec.set_etat_qcq() ;
f_21_21.set_etat_qcq() ;
f_2_1.set_etat_qcq() ;
// Function exp(-(r-r_0)^2/sigma^2)
// --------------------------------
double r0 = hole2.mp.val_r(nz2-2, 1, 0, 0) ;
double sigma = 1.*r0 ;
Scalar rr (hole2.mp) ;
rr = hole2.mp.r ;
Scalar fexp (hole2.mp) ;
fexp = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
for (int ii=0; ii<nz2-1; ii++)
fexp.set_domain(ii) = 1. ;
fexp.set_outer_boundary(nz2-1, 0) ;
fexp.std_spectral_base() ;
// Conformal metric
//=================
// tilde{gamma}- \delta = m_1*m_2* ( f_delta \delta_{ij}
// + f_2_2 nn2*nn2 + f_2_21 nn2*nn21
// + f_21_21 nn21*nn21
// + f_1_1 nn1*nn1 + f_1_21 nn1*nn21
// + f_2_1 nn2*nn1
// f_delta
//--------
f_delta = -5.*r2/(8.*r21*r21*r21) - 15./(8.*r2*r21) +
5.*r2*r2*unsr1/(8.*r21*r21*r21) + 1./(r2+1./unsr1+r21)/(r2+1./unsr1+r21)*
(1 + r2/r21 + r21/r2 - r2sr1 - r2*r2sr1/r21 + r21*r21*unsr1/(2*r2)) +
1./(r2+1./unsr1+r21)*(-7./r2 + 2./r21) ;
f_delta.annule_domain(nz2-1) ;
f_delta_zec = - 15./(8.*r2*r21) + 1./(r2+1./unsr1+r21)/(r2+1./unsr1+r21)*
(1 + r2/r21 + r21/r2 - r2sr1 - r2*r2sr1/r21 + r21*r21*unsr1/(2*r2)) +
1./(r2+1./unsr1+r21)*(-7./r2 + 2./r21) ;
f_delta_zec += fexp*(-5.*r2/(8.*r21*r21*r21)+5.*r2*r2*unsr1/(8.*r21*r21*r21)) ;
f_delta_zec.set_outer_boundary(nz2-1, 0.) ;
for (int i=0 ;i<nz2-1 ; i++){
f_delta_zec.annule_domain(i) ;
}
f_delta = f_delta + f_delta_zec ;
/*
des_meridian(f_delta, 0., 20., "f_delta", 10) ;
arrete() ;
des_profile(f_delta, 0., 20., M_PI/2, M_PI) ;
des_profile(f_delta, 0., 20., M_PI/2, 0) ;
des_profile(f_delta, 0., 20., 0, M_PI) ;
des_coupe_z (f_delta, 0., 2) ;
des_coupe_z (f_delta, 0., 3) ;
des_coupe_z (f_delta, 0., 4) ;
des_coupe_z (f_delta, 0., 5) ;
*/
// f_2_2
//------
f_2_2 = r2/(8.*r21*r21*r21) + 11./(8.*r2*r21) -
1./(8.*r2*unsr1*unsr1*r21*r21*r21) + 7./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) +
7./r2/(r2+1./unsr1+r21) ;
f_2_2.annule_domain(nz2-1) ;
f_2_2_zec = 11./(8.*r2*r21) + 7./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) +
7./r2/(r2+1./unsr1+r21) ;
f_2_2_zec += fexp*(r2/(8.*r21*r21*r21)-1./(8.*r2*unsr1*unsr1*r21*r21*r21)) ;
f_2_2_zec.set_outer_boundary(nz2-1, 0.) ;
for (int i=0 ; i<nz2-1 ; i++){
f_2_2_zec.annule_domain(i) ;
}
f_2_2 = f_2_2 + f_2_2_zec ;
/*
des_meridian(f_2_2, 0., 20., "f_2_2", 14) ;
arrete() ;
des_profile(f_2_2, 0., 20., M_PI/2, M_PI) ;
des_profile(f_2_2, 0., 20., M_PI/2, 0) ;
des_profile(f_2_2, 0., 20., 0, M_PI) ;
des_coupe_z (f_2_2, 0., 2) ;
des_coupe_z (f_2_2, 0., 3) ;
des_coupe_z (f_2_2, 0., 4) ;
des_coupe_z (f_2_2, 0., 5) ;
*/
// f_2_21
//------
f_2_21 = - 7./(2*r21*r21) + 8./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) ;
f_2_21.annule_domain(nz2-1) ;
f_2_21_zec = 8./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) ;
f_2_21_zec += fexp*(- 7./(2*r21*r21)) ;
f_2_21_zec.set_outer_boundary(nz2-1, 0.) ;
for (int i=0 ; i<nz2-1 ; i++){
f_2_21_zec.annule_domain(i) ;
}
f_2_21 = f_2_21 + f_2_21_zec ;
/*
des_meridian(f_2_21, 0., 40., "f_2_21", 15) ;
arrete() ;
des_profile(f_2_21, 0., 20., M_PI/2, M_PI) ;
des_profile(f_2_21, 0., 20., M_PI/2, 0) ;
des_profile(f_2_21, 0., 20., 0, M_PI) ;
des_coupe_z (f_2_21, 0., 2) ;
des_coupe_z (f_2_21, 0., 3) ;
des_coupe_z (f_2_21, 0., 4) ;
des_coupe_z (f_2_21, 0., 5) ;
*/
// f_21_21
//-------
f_21_21 = (-4./(r2+1./unsr1+r21)/(r2+1./unsr1+r21) - // facteur 2 ???????
4./r21/(r2+1./unsr1+r21)) ;
f_21_21.set_outer_boundary(nz2-1, 0.) ;
/*
des_meridian(f_21_21, 0., 40., "f_21_21", 15) ;
arrete() ;
des_profile(f_21_21, 0., 20., M_PI/2, M_PI) ;
des_profile(f_21_21, 0., 20., M_PI/2, 0) ;
des_profile(f_21_21, 0., 20., 0, M_PI) ;
des_coupe_z (f_21_21, 0., 2) ;
des_coupe_z (f_21_21, 0., 3) ;
des_coupe_z (f_21_21, 0., 4) ;
des_coupe_z (f_21_21, 0., 5) ;
*/
// f_2_1
//-------
f_2_1 = 11./(r2+1./unsr1+r21)/(r2+1./unsr1+r21); // facteur 2 ???????
f_2_1.set_outer_boundary(nz2-1, 0.) ;
/*
des_meridian(f_2_1, 0., 40., "f_2_1", 15) ;
arrete() ;
des_profile(f_2_1, 0., 20., M_PI/2, M_PI) ;
des_profile(f_2_1, 0., 20., M_PI/2, 0) ;
des_profile(f_2_1, 0., 20., 0, M_PI) ;
des_coupe_z (f_2_1, 0., 2) ;
des_coupe_z (f_2_1, 0., 3) ;
des_coupe_z (f_2_1, 0., 4) ;
des_coupe_z (f_2_1, 0., 5) ;
*/
// First part of the correction metric (needed to be complemented by the (1 <-> 2) term
Sym_tensor hh_temp(hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
for (int i=1 ; i<= 3 ; i++){
for (int j=i ; j<= 3 ; j++){
hh_temp.set(i,j) = f_delta * hole2.ff.con()(i,j) + f_2_2 * nn2(i)*nn2(j)
- f_2_21 * 0.5 *(nn2(i) * nn21(j) + nn2(j) * nn21(i))
+ f_21_21 * nn21(i)*nn21(j)
+ f_2_1 * 0.5*(nn2(i)*nn1(j) + nn2(j)*nn1(i) ) ;
}
}
/*
des_meridian(hh_temp, 0., 20., "hh_temp", 25) ;
arrete() ;
for (int i=1 ; i<= 3 ; i++)
for (int j=i ; j<= 3 ; j++){
des_profile(hh_temp(i,j), 0., 20., M_PI/2, M_PI) ;
des_profile(hh_temp(i,j), 0., 20., M_PI/2, 0) ;
des_profile(hh_temp(i,j), 0., 20., 0, M_PI) ;
des_coupe_z (hh_temp(i,j), 0., 5) ;
}
*/
return hh_temp ;
}
void Bin_hor::set_hh_Samaya() {
Sym_tensor hh1 ( hh_Samaya_hole1() ) ;
Sym_tensor hh2 ( hh_Samaya_hole2() ) ;
// Definition of the surface
// -------------------------
Cmp surface_un (hole1.mp) ;
surface_un = pow(hole1.mp.r, 2.)-pow(hole1.get_radius(), 2.) ;
surface_un.annule(hole1.mp.get_mg()->get_nzone()-1) ;
surface_un.std_base_scal() ;
Cmp surface_deux (hole2.mp) ;
surface_deux = pow(hole2.mp.r, 2.)-pow(hole2.get_radius(), 2.) ;
surface_deux.annule(hole1.mp.get_mg()->get_nzone()-1) ;
surface_deux.std_base_scal() ;
/*
double ta = 12 ;
for (int i=1 ; i<= 3 ; i++)
for (int j=i ; j<= 3 ; j++){
Cmp dessin_un (hh1(i,j)) ;
dessin_un.annule(0) ;
Cmp dessin_deux (hh2(i,j)) ;
dessin_deux.annule(0) ;
des_coupe_bin_z (dessin_un, dessin_deux, 0,
-ta, ta, -ta, ta, "hh(1,1)", &surface_un, &surface_deux,
false, 15, 300, 300) ;
}
*/
// Importation
// ----------------
Sym_tensor hh2_1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
hh2_1.set_etat_qcq() ;
Sym_tensor hh1_2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
hh1_2.set_etat_qcq() ;
/*
Scalar temp (hh1(1,1)) ;
temp.annule_domain(0) ;
des_profile(temp, 0., 4., M_PI/2, M_PI) ;
des_profile(temp, 0., 4., M_PI/2, 0) ;
des_profile(temp, 0., 4., 0, M_PI) ;
des_coupe_z (temp, 0., 5) ;
temp.raccord(1) ;
des_profile(temp, 0., 4., M_PI/2, M_PI) ;
des_profile(temp, 0., 4., M_PI/2, 0) ;
des_profile(temp, 0., 4., 0, M_PI) ;
des_coupe_z (temp, 0., 5) ;
*/
/*
for (int i=1 ; i<= 3 ; i++)
for (int j=i ; j<= 3 ; j++){
des_profile(hh1(i,j), 0., 20., M_PI/2, M_PI) ;
des_profile(hh1(i,j), 0., 20., M_PI/2, 0) ;
des_profile(hh1(i,j), 0., 20., 0, M_PI) ;
des_coupe_z (hh1(i,j), 0., 5) ;
}
*/
for (int i=1 ; i<=3 ; i++)
for (int j=i ; j<=3 ; j++){
hh1.set(i,j).raccord(1) ;
hh2.set(i,j).raccord(1) ;
}
/*
for (int i=1 ; i<= 3 ; i++)
for (int j=i ; j<= 3 ; j++){
des_profile(hh1(i,j), 0., 20., M_PI/2, M_PI) ;
des_profile(hh1(i,j), 0., 20., M_PI/2, 0) ;
des_profile(hh1(i,j), 0., 20., 0, M_PI) ;
des_coupe_z (hh1(i,j), 0., 5) ;
}
*/
hh2.change_triad(hole1.mp.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++){
for (int j=i ; j<=3 ; j++){
hh2_1.set(i,j).import(hh2(i,j)) ;
hh2_1.set(i,j).set_spectral_va().set_base(hh2(i,j).get_spectral_va().get_base()) ;
}
}
hh2.change_triad(hole2.mp.get_bvect_cart()) ;
hh1.change_triad(hole2.mp.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++){
for (int j=i ; j<=3 ; j++){
hh1_2.set(i,j).import(hh1(i,j)) ;
hh1_2.set(i,j).set_spectral_va().set_base(hh1(i,j).get_spectral_va().get_base()) ;
}
}
hh1.change_triad(hole1.mp.get_bvect_cart()) ;
double m1, m2 ;
m1 = pow(hole1.area_hor()/(16.*M_PI) + hole1.ang_mom_hor()/hole1.radius, 0.5) ;
m2 = pow(hole2.area_hor()/(16.*M_PI) + hole2.ang_mom_hor()/hole2.radius, 0.5) ;
hh1 = hh1 + hh2_1 ;
hh2 = hh2 + hh1_2 ;
cout << hole1.mp.r << endl ;
cout << hole1.mp.phi << endl ;
cout << hole1.mp.tet << endl ;
//des_meridian(hh1, 0., 20., "hh1 cart", 20) ;
for (int i=1 ; i<= 3 ; i++)
for (int j=i ; j<= 3 ; j++){
// des_profile(hh1(i,j), 0., 20., M_PI/2, M_PI) ;
//des_profile(hh1(i,j), 0., 20., M_PI/2, 0) ;
//des_profile(hh1(i,j), 0., 20., 0, M_PI) ;
des_coupe_z (hh1(i,j), 0., 5) ;
}
hh1.change_triad(hole1.mp.get_bvect_spher()) ;
hh2.change_triad(hole2.mp.get_bvect_spher()) ;
hole1.hh = m1*m2* hh1 ;
hole2.hh = m1*m2* hh2 ;
Metric tgam_1 ( hole1.ff.con() + hh1 ) ;
Metric tgam_2 ( hole2.ff.con() + hh2 ) ;
hole1.tgam = tgam_1 ;
hole2.tgam = tgam_2 ;
des_meridian(hh1, 0., 20., "hh1", 0) ;
}
}
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