27 char bin_ns_bh_kij_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_kij.C,v 1.11 2014/10/13 08:52:43 j_novak Exp $" ;
82 #include "bin_ns_bh.h" 84 #include "utilitaires.h" 86 #include "graphique.h" 109 for (
int i=0 ; i<3 ; i++)
110 for (
int j=i ; j<3 ; j++)
112 for (
int i=0 ; i<3 ; i++)
121 int nz_bh = hole.mp.get_mg()->get_nzone() ;
122 int nz_ns = star.mp.get_mg()->get_nzone() ;
125 double distance = star.mp.get_ori_x()-hole.mp.get_ori_x() ;
126 double lim_ns = distance/2. ;
127 double lim_bh = distance/2. ;
128 double int_ns = lim_ns/3. ;
129 double int_bh = lim_bh/3. ;
132 Cmp decouple_ns (star.mp) ;
134 Cmp decouple_bh (hole.mp) ;
137 Mtbl xabs_ns (star.mp.xa) ;
138 Mtbl yabs_ns (star.mp.ya) ;
139 Mtbl zabs_ns (star.mp.za) ;
141 Mtbl xabs_bh (hole.mp.xa) ;
142 Mtbl yabs_bh (hole.mp.ya) ;
143 Mtbl zabs_bh (hole.mp.za) ;
145 double xabs, yabs, zabs, air_ns, air_bh, theta, phi ;
148 for (
int l=0 ; l<nz_ns ; l++) {
149 int nr = star.mp.get_mg()->get_nr (l) ;
154 int np = star.mp.get_mg()->get_np (l) ;
155 int nt = star.mp.get_mg()->get_nt (l) ;
157 for (
int k=0 ; k<np ; k++)
158 for (
int j=0 ; j<nt ; j++)
159 for (
int i=0 ; i<nr ; i++) {
161 xabs = xabs_ns (l, k, j, i) ;
162 yabs = yabs_ns (l, k, j, i) ;
163 zabs = zabs_ns (l, k, j, i) ;
166 star.mp.convert_absolute
167 (xabs, yabs, zabs, air_ns, theta, phi) ;
168 hole.mp.convert_absolute
169 (xabs, yabs, zabs, air_bh, theta, phi) ;
171 if (air_ns <= lim_ns)
173 decouple_ns.
set(l, k, j, i) = 1 ;
175 decouple_ns.
set(l, k, j, i) =
176 0.5*
pow(
cos((air_ns-int_ns)*M_PI/2./(lim_ns-int_ns)), 2.)+0.5
179 if (air_bh <= lim_bh)
181 decouple_ns.
set(l, k, j, i) = 0 ;
183 decouple_ns.
set(l, k, j, i) = 0.5*
184 pow(
sin((air_bh-int_bh)*M_PI/2./(lim_bh-int_bh)), 2.)
188 decouple_ns.
set(l, k, j, i) = 0.5 ;
193 for (
int k=0 ; k<np ; k++)
194 for (
int j=0 ; j<nt ; j++)
195 decouple_ns.
set(nz_ns-1, k, j, nr) = 0.5 ;
198 for (
int l=0 ; l<nz_bh ; l++) {
199 int nr = hole.mp.get_mg()->get_nr (l) ;
204 int np = hole.mp.get_mg()->get_np (l) ;
205 int nt = hole.mp.get_mg()->get_nt (l) ;
207 for (
int k=0 ; k<np ; k++)
208 for (
int j=0 ; j<nt ; j++)
209 for (
int i=0 ; i<nr ; i++) {
211 xabs = xabs_bh (l, k, j, i) ;
212 yabs = yabs_bh (l, k, j, i) ;
213 zabs = zabs_bh (l, k, j, i) ;
216 star.mp.convert_absolute
217 (xabs, yabs, zabs, air_ns, theta, phi) ;
218 hole.mp.convert_absolute
219 (xabs, yabs, zabs, air_bh, theta, phi) ;
221 if (air_bh <= lim_bh)
223 decouple_bh.
set(l, k, j, i) = 1 ;
225 decouple_bh.
set(l, k, j, i) = 0.5*
226 pow(
cos((air_bh-int_bh)*M_PI/2./(lim_bh-int_bh)), 2.)+0.5 ;
228 if (air_ns <= lim_ns)
230 decouple_bh.
set(l, k, j, i) = 0 ;
232 decouple_bh.
set(l, k, j, i) = 0.5*
233 pow(
sin((air_ns-int_ns)*M_PI/2./(lim_ns-int_ns)), 2.) ;
237 decouple_bh.
set(l, k, j, i) = 0.5 ;
242 for (
int k=0 ; k<np ; k++)
243 for (
int j=0 ; j<nt ; j++)
244 decouple_bh.
set(nz_bh-1, k, j, nr) = 0.5 ;
250 star.decouple = decouple_ns ;
251 hole.decouple = decouple_bh ;
261 double norme_hole = 0 ;
262 double norme_star = 0 ;
264 for (
int i=0 ; i<3 ; i++) {
265 norme_hole +=
max(
norme(hole.get_shift_auto()(i))) ;
266 norme_star +=
max(
norme(star.get_shift_auto()(i))) ;
270 bool zero_shift_hole = (norme_hole <1e-14) ?
true :
false ;
272 bool zero_shift_star = (norme_star <1e-14) ?
true :
false ;
274 assert (zero_shift_hole == zero_shift_star) ;
276 if (zero_shift_star ==
true) {
277 hole.tkij_tot.set_etat_zero() ;
278 hole.tkij_auto.set_etat_zero() ;
280 hole.taij_tot.set_etat_zero() ;
281 hole.taij_auto.set_etat_zero() ;
282 hole.taij_comp.set_etat_zero() ;
284 star.tkij_auto.set_etat_zero() ;
285 star.tkij_comp.set_etat_zero() ;
286 star.akcar_comp.set_etat_zero() ;
287 star.akcar_auto.set_etat_zero() ;
292 int nnt = hole.mp.get_mg()->get_nt(1) ;
293 int nnp = hole.mp.get_mg()->get_np(1) ;
295 for (
int k=0; k<nnp; k++)
296 for (
int j=0; j<nnt; j++){
297 if (fabs((hole.n_auto+hole.n_comp)()(1, k, j , 0)) < 1e-2){
305 if (bound_nn != 0 || lim_nn != 0){
307 hole.fait_taij_auto () ;
308 star.fait_taij_auto() ;
311 hole.taij_comp.set_etat_qcq() ;
322 if (hole.taij_auto.get_etat() == ETATZERO)
325 ns_taij_comp.
set(0, 0).
import(copie_bh(0, 0)) ;
326 ns_taij_comp.
set(0, 1).
import(copie_bh(0, 1)) ;
327 ns_taij_comp.
set(0, 2).
import(copie_bh(0, 2)) ;
328 ns_taij_comp.
set(1, 1).
import(copie_bh(1, 1)) ;
329 ns_taij_comp.
set(1, 2).
import(copie_bh(1, 2)) ;
330 ns_taij_comp.
set(2, 2).
import(copie_bh(2, 2)) ;
335 if (star.taij_auto.get_etat() == ETATZERO)
336 hole.taij_comp.set_etat_zero() ;
338 hole.taij_comp.set(0, 0).import(copie_ns(0, 0)) ;
339 hole.taij_comp.set(0, 1).import(copie_ns(0, 1)) ;
340 hole.taij_comp.set(0, 2).import(copie_ns(0, 2)) ;
341 hole.taij_comp.set(1, 1).import(copie_ns(1, 1)) ;
342 hole.taij_comp.set(1, 2).import(copie_ns(1, 2)) ;
343 hole.taij_comp.set(2, 2).import(copie_ns(2, 2)) ;
344 hole.taij_comp.set_triad(*copie_ns.
get_triad()) ;
345 hole.taij_comp.change_triad (hole.mp.get_bvect_cart()) ;
348 hole.taij_comp.set_std_base() ;
350 hole.taij_comp.inc2_dzpuis() ;
354 hole.taij_tot = hole.taij_auto + hole.taij_comp ;
355 Tenseur_sym ns_taij_tot (star.taij_auto + ns_taij_comp) ;
356 star.taij_tot = ns_taij_tot ;
360 star.tkij_auto.set_etat_qcq() ;
361 star.tkij_comp.set_etat_qcq() ;
362 hole.tkij_tot.set_etat_qcq() ;
363 hole.tkij_auto.set_etat_qcq() ;
365 for (
int i = 0 ; i<3 ; i++)
366 for (
int j = i ; j<3 ; j++) {
367 star.tkij_tot.set(i,j) = 0.5*star.taij_tot(i,j)/star.nnn() ;
370 hole.tkij_tot.set(i,j) = 0.5*hole.taij_tot(i,j)/hole.n_tot() ;
374 for (
int lig=0 ; lig<3 ; lig++)
375 for (
int col=lig ; col<3 ; col++) {
376 star.tkij_auto.set(lig, col) = star.tkij_tot(lig, col)*
378 star.tkij_comp.set(lig, col) = star.tkij_tot(lig, col)*
380 hole.tkij_auto.set(lig, col) = hole.tkij_tot(lig, col)*
383 star.tkij_auto.set_std_base() ;
384 star.tkij_comp.set_std_base() ;
385 hole.tkij_auto.set_std_base() ;
389 hole.tkij_tot.set_etat_qcq() ;
390 star.tkij_tot.set_etat_qcq() ;
394 hole.fait_taij_auto () ;
395 star.fait_taij_auto() ;
398 hole.taij_comp.set_etat_qcq() ;
409 if (hole.taij_auto.get_etat() == ETATZERO)
422 if (star.taij_auto.get_etat() == ETATZERO)
423 hole.taij_comp.set_etat_zero() ;
425 hole.taij_comp.set(0, 0).import_asymy(copie_ns(0, 0)) ;
426 hole.taij_comp.set(0, 1).import_symy(copie_ns(0, 1)) ;
427 hole.taij_comp.set(0, 2).import_asymy(copie_ns(0, 2)) ;
428 hole.taij_comp.set(1, 1).import_asymy(copie_ns(1, 1)) ;
429 hole.taij_comp.set(1, 2).import_symy(copie_ns(1, 2)) ;
430 hole.taij_comp.set(2, 2).import_asymy(copie_ns(2, 2)) ;
431 hole.taij_comp.set_triad(*copie_ns.
get_triad()) ;
432 hole.taij_comp.change_triad (hole.mp.get_bvect_cart()) ;
435 hole.taij_comp.set_std_base() ;
437 hole.taij_comp.inc2_dzpuis() ;
441 hole.taij_tot = hole.taij_auto + hole.taij_comp ;
442 Tenseur_sym ns_taij_tot (star.taij_auto + ns_taij_comp) ;
443 star.taij_tot = ns_taij_tot ;
446 Cmp ntot_ns (star.get_nnn()()) ;
448 Cmp ntot_bh (hole.n_tot()) ;
450 ntot_bh = division_xpun (ntot_bh, 0) ;
456 for (
int lig = 0 ; lig<3 ; lig++)
457 for (
int col = lig ; col<3 ; col++) {
460 Cmp auxi_bh (hole.taij_tot(lig, col)/2.) ;
462 auxi_bh = division_xpun (auxi_bh, 0) ;
465 auxi_bh = auxi_bh / ntot_bh ;
469 Cmp auxi_ns (star.mp) ;
473 Cmp copie_ns_bis (ns_taij_tot(lig, col)) ;
477 double lim_bh = hole.mp.get_alpha()[1] + hole.mp.get_beta()[1] ;
479 Mtbl xabs_ns (star.mp.xa) ;
480 Mtbl yabs_ns (star.mp.ya) ;
481 Mtbl zabs_ns (star.mp.za) ;
483 double xabs, yabs, zabs, air, theta, phi ;
486 for (
int l=0 ; l<nz_ns ; l++) {
488 int nr = star.mp.get_mg()->get_nr (l) ;
493 int np = star.mp.get_mg()->get_np (l) ;
494 int nt = star.mp.get_mg()->get_nt (l) ;
496 for (
int k=0 ; k<np ; k++)
497 for (
int j=0 ; j<nt ; j++)
498 for (
int i=0 ; i<nr ; i++) {
500 xabs = xabs_ns (l, k, j, i) ;
501 yabs = yabs_ns (l, k, j, i) ;
502 zabs = zabs_ns (l, k, j, i) ;
505 hole.mp.convert_absolute
506 (xabs, yabs, zabs, air, theta, phi) ;
510 auxi_ns.
set(l, k, j, i) =
511 copie_ns_bis(l, k, j, i) / ntot_ns (l, k, j, i)/2. ;
514 auxi_ns.
set(l, k, j, i) = auxi_bh.
val_point (air, theta, phi) ;
519 for (
int k=0 ; k<np ; k++)
520 for (
int j=0 ; j<nt ; j++)
521 auxi_ns.
set(nz_ns-1, k, j, nr) = 0 ;
525 star.tkij_tot.
set(lig, col) = auxi_ns ;
526 hole.tkij_tot.
set(lig, col) = auxi_bh ;
529 star.tkij_tot.set_std_base() ;
530 hole.tkij_tot.set_std_base() ;
533 hole.tkij_tot.inc2_dzpuis() ;
540 star.tkij_auto.set_etat_qcq() ;
541 star.tkij_comp.set_etat_qcq() ;
542 hole.tkij_auto.set_etat_qcq() ;
544 for (
int lig=0 ; lig<3 ; lig++)
545 for (
int col=lig ; col<3 ; col++) {
546 star.tkij_auto.set(lig, col) = star.tkij_tot(lig, col)*
548 star.tkij_comp.set(lig, col) = star.tkij_tot(lig, col)*
550 hole.tkij_auto.set(lig, col) = hole.tkij_tot(lig, col)*
553 star.tkij_auto.set_std_base() ;
554 star.tkij_comp.set_std_base() ;
555 hole.tkij_auto.set_std_base() ;
560 star.akcar_auto.set_etat_qcq() ;
561 star.akcar_auto.set() = 0 ;
563 for (
int i=0; i<3; i++)
564 for (
int j=0; j<3; j++)
565 star.akcar_auto.set() += star.tkij_auto(i, j) % star.tkij_auto(i, j) ;
567 star.akcar_auto.set_std_base() ;
568 star.akcar_auto = star.a_car % star.akcar_auto ;
570 star.akcar_comp.set_etat_qcq() ;
571 star.akcar_comp.set() = 0 ;
573 for (
int i=0; i<3; i++)
574 for (
int j=0; j<3; j++)
575 star.akcar_comp.set() += star.tkij_auto(i, j) % star.tkij_comp(i, j) ;
577 star.akcar_comp.set_std_base() ;
578 star.akcar_comp = star.a_car % star.akcar_comp ;
const Map *const mp
Reference mapping.
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void dec2_dzpuis()
dzpuis -= 2 ;
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Cmp cos(const Cmp &)
Cosine.
void inc2_dzpuis()
dzpuis += 2 ;
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
void fait_taij_auto()
Computes (LB)^{ij} auto.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
void fait_decouple()
Function used to compute the { decouple} functions for both the NS and the BH.
Map & mp
Mapping associated with the star.
int get_nzone() const
Returns the number of domains.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
int get_etat() const
Returns the logical state.
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Cmp pow(const Cmp &, int)
Power .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
void import_symy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Tbl & set(int l)
Read/write of the value in a given domain.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Tenseur_sym taij_auto
Part of the extrinsic curvature tensor $ A^{ij} = 2 N K^{ij}$ generated by { shift_auto}.
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
void fait_tkij(int bound_nn=-1, double lim_nn=0)
Computation of the extrinsic curvature tensor for both { star} and { bhole}.
Cmp sin(const Cmp &)
Sine.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void import(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
void import_asymy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
double val_point(double r, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point , by means of the spectral...
Tensor handling *** DEPRECATED : use class Tensor instead ***.
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)