bin_bhns_global.C

00001 /*
00002  *  Methods of class Bin_bhns to compute global quantities
00003  *
00004  *    (see file bin_bhns.h for documentation).
00005  *
00006  */
00007 
00008 /*
00009  *   Copyright (c) 2005-2007 Keisuke Taniguchi
00010  *
00011  *   This file is part of LORENE.
00012  *
00013  *   LORENE is free software; you can redistribute it and/or modify
00014  *   it under the terms of the GNU General Public License version 2
00015  *   as published by the Free Software Foundation.
00016  *
00017  *   LORENE is distributed in the hope that it will be useful,
00018  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00019  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00020  *   GNU General Public License for more details.
00021  *
00022  *   You should have received a copy of the GNU General Public License
00023  *   along with LORENE; if not, write to the Free Software
00024  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00025  *
00026  */
00027 
00028 char bin_bhns_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_global.C,v 1.2 2008/05/15 18:59:27 k_taniguchi Exp $" ;
00029 
00030 /*
00031  * $Id: bin_bhns_global.C,v 1.2 2008/05/15 18:59:27 k_taniguchi Exp $
00032  * $Log: bin_bhns_global.C,v $
00033  * Revision 1.2  2008/05/15 18:59:27  k_taniguchi
00034  * Introduction of new global quantities.
00035  *
00036  * Revision 1.1  2007/06/22 01:09:31  k_taniguchi
00037  * *** empty log message ***
00038  *
00039  *
00040  * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_global.C,v 1.2 2008/05/15 18:59:27 k_taniguchi Exp $
00041  *
00042  */
00043 
00044 // C++ headers
00045 //#include <>
00046 
00047 // C headers
00048 #include <math.h>
00049 
00050 // Lorene headers
00051 #include "bin_bhns.h"
00052 #include "blackhole.h"
00053 #include "unites.h"
00054 #include "utilitaires.h"
00055 #include "nbr_spx.h"
00056 
00057                //----------------------------//
00058                //          ADM mass          //
00059                //----------------------------//
00060 
00061 double Bin_bhns::mass_adm_bhns_surf() const {
00062 
00063     // Fundamental constants and units
00064     // -------------------------------
00065     using namespace Unites ;
00066 
00067     if (p_mass_adm_bhns_surf == 0x0) {   // a new computation is required
00068 
00069         double madm ;
00070 
00071     const Map& mp_bh = hole.get_mp() ;
00072         const Map& mp_ns = star.get_mp() ;
00073 
00074     Map_af mp_aff(mp_bh) ;
00075     Map_af mp_ns_aff(mp_ns) ;
00076 
00077     Scalar rr(mp_bh) ;
00078     rr = mp_bh.r ;
00079     rr.std_spectral_base() ;
00080 
00081     double mass = ggrav * hole.get_mass_bh() ;
00082 
00083     if (hole.is_kerrschild()) {
00084 
00085       cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
00086       abort() ;
00087 
00088     }
00089     else { // Isotropic coordinates with the maximal slicing
00090 
00091       //-------------------------------------
00092       //     Integration over the BH map
00093       //-------------------------------------
00094 
00095       // Sets C/M^2 for each case of the lapse boundary condition
00096       // --------------------------------------------------------
00097       double cc ;
00098 
00099       if (hole.has_bc_lapconf_nd()) {  // Neumann boundary condition
00100           if (hole.has_bc_lapconf_fs()) {  // First condition
00101               // d(\alpha \psi)/dr = 0
00102               // ---------------------
00103               cc = 2. * (sqrt(13.) - 1.) / 3. ;
00104           }
00105           else {  // Second condition
00106               // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
00107               // -----------------------------------------
00108               cc = 4. / 3. ;
00109           }
00110       }
00111       else {  // Dirichlet boundary condition
00112           if (hole.has_bc_lapconf_fs()) {  // First condition
00113               // (\alpha \psi) = 1/2
00114               // -------------------
00115           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00116           abort() ;
00117           }
00118           else {  // Second condition
00119               // (\alpha \psi) = 1/sqrt(2.) \psi_KS
00120               // ----------------------------------
00121           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00122           abort() ;
00123           //              cc = 2. * sqrt(2.) ;
00124           }
00125       }
00126 
00127       Scalar r_are(mp_bh) ;
00128       r_are = hole.r_coord(hole.has_bc_lapconf_nd(),
00129                    hole.has_bc_lapconf_fs()) ;
00130       r_are.std_spectral_base() ;
00131 
00132       // ADM mass by surface integral at infinity : dzpuis should be 2
00133       // ----------------------------------------
00134       const Scalar& confo_bh_auto_rs = hole.get_confo_auto_rs() ;
00135 
00136       Scalar lldconf_iso = confo_bh_auto_rs.dsdr() ;  // dzpuis = 2
00137       lldconf_iso.std_spectral_base() ;
00138 
00139       Scalar anoth(mp_bh) ;
00140       anoth = 0.5 * sqrt(r_are)
00141         * (sqrt(1. -2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
00142            - 1.) / rr ;
00143       anoth.std_spectral_base() ;
00144       anoth.annule_domain(0) ;
00145       anoth.raccord(1) ;
00146       anoth.inc_dzpuis(2) ;
00147 
00148       const Scalar& confo_ns_auto = star.get_confo_auto() ;
00149 
00150       Scalar lldconf_ns = confo_ns_auto.dsdr() ;  // dzpuis = 2
00151       lldconf_ns.std_spectral_base() ;
00152 
00153       madm =
00154         - 2.*(mp_aff.integrale_surface_infini(lldconf_iso+anoth))/qpig
00155         - 2.*(mp_ns_aff.integrale_surface_infini(lldconf_ns))/qpig ;
00156 
00157       cout << "ADM mass (surface) :   " << madm / msol << " [Mo]"
00158            << endl ;
00159 
00160     }
00161 
00162     p_mass_adm_bhns_surf = new double( madm ) ;
00163 
00164     }
00165 
00166     return *p_mass_adm_bhns_surf ;
00167 
00168 }
00169 
00170 
00171                //----------------------------//
00172                //          ADM mass          //
00173                //----------------------------//
00174 
00175 double Bin_bhns::mass_adm_bhns_vol() const {
00176 
00177     // Fundamental constants and units
00178     // -------------------------------
00179     using namespace Unites ;
00180 
00181     if (p_mass_adm_bhns_vol == 0x0) {   // a new computation is required
00182 
00183         double madm ;
00184     double integ_bh_s ;
00185     double integ_bh_v ;
00186     double integ_ns_v ;
00187 
00188     const Map& mp_bh = hole.get_mp() ;
00189         const Map& mp_ns = star.get_mp() ;
00190 
00191     double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
00192     Map_af mp_aff(mp_bh) ;
00193 
00194     Map_af mp_ns_aff(mp_ns) ;
00195 
00196     Scalar source_bh_surf(mp_bh) ;
00197     source_bh_surf.set_etat_qcq() ;
00198 
00199     Scalar source_bh_volm(mp_bh) ;
00200     source_bh_volm.set_etat_qcq() ;
00201 
00202     Scalar source_ns_volm(mp_ns) ;
00203     source_ns_volm.set_etat_qcq() ;
00204 
00205     Scalar rr(mp_bh) ;
00206     rr = mp_bh.r ;
00207     rr.std_spectral_base() ;
00208     Scalar st(mp_bh) ;
00209     st = mp_bh.sint ;
00210     st.std_spectral_base() ;
00211     Scalar ct(mp_bh) ;
00212     ct = mp_bh.cost ;
00213     ct.std_spectral_base() ;
00214     Scalar sp(mp_bh) ;
00215     sp = mp_bh.sinp ;
00216     sp.std_spectral_base() ;
00217     Scalar cp(mp_bh) ;
00218     cp = mp_bh.cosp ;
00219     cp.std_spectral_base() ;
00220 
00221     Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
00222     ll.set_etat_qcq() ;
00223     ll.set(1) = st % cp ;
00224     ll.set(2) = st % sp ;
00225     ll.set(3) = ct ;
00226     ll.std_spectral_base() ;
00227 
00228     const Vector& shift_auto_rs = hole.get_shift_auto_rs() ;
00229     const Vector& shift_comp = hole.get_shift_comp() ;
00230     const Tensor& dshift_comp = hole.get_d_shift_comp() ;
00231 
00232     Scalar divshift(mp_bh) ;  // dzpuis = 2
00233     divshift = shift_auto_rs(1).deriv(1) + shift_auto_rs(2).deriv(2)
00234       + shift_auto_rs(3).deriv(3) + dshift_comp(1,1)
00235       + dshift_comp(2,2) + dshift_comp(3,3) ;
00236     divshift.std_spectral_base() ;
00237 
00238     Scalar llshift(mp_bh) ;   // dzpuis = 0
00239     llshift = ll(1) % (shift_auto_rs(1) + shift_comp(1))
00240       + ll(2) % (shift_auto_rs(2) + shift_comp(2))
00241       + ll(3) % (shift_auto_rs(3) + shift_comp(3)) ;
00242     llshift.std_spectral_base() ;
00243 
00244     Scalar llshift_auto(mp_bh) ;   // dzpuis = 0
00245     llshift_auto = ll(1)%shift_auto_rs(1) + ll(2)%shift_auto_rs(2)
00246       + ll(3)%shift_auto_rs(3) ;
00247     llshift_auto.std_spectral_base() ;
00248 
00249     Scalar lldllsh = llshift_auto.dsdr()
00250       + ll(1) * (ll(1) % dshift_comp(1,1) + ll(2) % dshift_comp(1,2)
00251              + ll(3) % dshift_comp(1,3))
00252       + ll(2) * (ll(1) % dshift_comp(2,1) + ll(2) % dshift_comp(2,2)
00253              + ll(3) % dshift_comp(2,3))
00254       + ll(3) * (ll(1) % dshift_comp(3,1) + ll(2) % dshift_comp(3,2)
00255              + ll(3) % dshift_comp(3,3)) ;  // dzpuis = 2
00256     lldllsh.std_spectral_base() ;
00257 
00258     const Scalar& lapconf_bh = hole.get_lapconf_tot() ;
00259     const Scalar& lapconf_bh_auto_rs = hole.get_lapconf_auto_rs() ;
00260     const Scalar& lapconf_bh_comp = hole.get_lapconf_comp() ;
00261     const Scalar& confo_bh = hole.get_confo_tot() ;
00262     const Scalar& confo_bh_auto = hole.get_confo_auto() ;
00263     const Scalar& confo_bh_comp = hole.get_confo_comp() ;
00264     const Vector& dconfo_bh_comp = hole.get_d_confo_comp() ;
00265     const Scalar& taij_quad_tot_rs = hole.get_taij_quad_tot_rs() ;
00266     const Scalar& taij_quad_tot_rot = hole.get_taij_quad_tot_rot() ;
00267 
00268     const Scalar& taij_quad_auto_bh = hole.get_taij_quad_auto() ;
00269     const Scalar& taij_quad_comp_bh = hole.get_taij_quad_comp() ;
00270 
00271     const Scalar& confo_ns = star.get_confo_tot() ;
00272     const Scalar& confo_ns_auto = star.get_confo_auto() ;
00273     const Scalar& ener_euler = star.get_ener_euler() ;
00274 
00275     const Scalar& taij_quad_auto_ns = star.get_taij_quad_auto() ;
00276 
00277     Scalar lldconf = confo_bh_auto.dsdr() + ll(1)%dconfo_bh_comp(1)
00278       + ll(2)%dconfo_bh_comp(2) + ll(3)%dconfo_bh_comp(3) ;  // dzpuis = 2
00279     lldconf.std_spectral_base() ;
00280 
00281     double mass = ggrav * hole.get_mass_bh() ;
00282 
00283     if (hole.is_kerrschild()) {
00284 
00285       cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
00286       abort() ;
00287 
00288       /*
00289       Scalar lap_bh(mp_bh) ;
00290       lap_bh = 1./sqrt(1.+2.*mass/rr) ;
00291       lap_bh.std_spectral_base() ;
00292 
00293       Scalar lap_bh2(mp_bh) ;
00294       lap_bh2 = 1./(1.+2.*mass/rr) ;
00295       lap_bh2.std_spectral_base() ;
00296 
00297       Scalar omelld(mp_bh) ;
00298       omelld = omega * (ll(2) * (mp_bh.get_ori_x() - x_rot)
00299                 - ll(1) * (mp_bh.get_ori_y() - y_rot)) ;
00300       omelld.std_spectral_base() ;
00301 
00302       Scalar lldlldconf(mp_bh) ;  // dzpuis = 3
00303       lldlldconf = lldconf.dsdr() ;
00304       lldlldconf.std_spectral_base() ;
00305 
00306       //-------------------------------------
00307       //     Integration over the BH map
00308       //-------------------------------------
00309 
00310       // Surface integral <- dzpuis should be 0
00311       // --------------------------------------
00312       Scalar divshift_zero(divshift) ;
00313       divshift_zero.dec_dzpuis(2) ;
00314 
00315       Scalar lldllsh_zero(lldllsh) ;
00316       lldllsh_zero.dec_dzpuis(2) ;
00317 
00318       source_bh_surf = confo_bh
00319         * (1. - 2.*mass*lap_bh*confo_bh*confo_bh/lapse_bh/rr) / rr
00320         - pow(confo_bh, 3.)
00321         * ( divshift_zero - 3.*lldllsh_zero
00322         + 2. * lap_bh2 * mass * (llshift + omelld) / rr / rr
00323         + 4.*mass*lap_bh2*lap_bh*(1.+3.*mass/rr)
00324         *(lapse_bh_auto_rs + lapse_bh_comp)/rr/rr
00325         ) / 6. / lap_bh / lapse_bh ;
00326 
00327       source_bh_surf.std_spectral_base() ;
00328       source_bh_surf.annule_domain(0) ;
00329       source_bh_surf.raccord(1) ;
00330 
00331       integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
00332 
00333       // Volume integral <- dzpuis should be 4
00334       // -------------------------------------
00335       Scalar sou_bh1(mp_bh) ;
00336       sou_bh1 = 2.*lap_bh2*mass*lldlldconf/rr ;
00337       sou_bh1.std_spectral_base() ;
00338       sou_bh1.inc_dzpuis(1) ;
00339 
00340       Scalar sou_bh2(mp_bh) ;
00341       sou_bh2 = lap_bh2*lap_bh2*mass*(3.+8.*mass/rr)*lldconf/rr/rr ;
00342       sou_bh2.std_spectral_base() ;
00343       sou_bh2.inc_dzpuis(2) ;
00344 
00345       Scalar sou_bh3(mp_bh) ;
00346       sou_bh3 = pow(lap_bh2,3.)*mass*mass*confo_bh
00347         * ( (1.-lap_bh2/lapse_bh/lapse_bh)
00348         *(4.+12.*mass/rr+9.*mass*mass/rr/rr)*pow(confo_bh,4.)
00349         +3.*(1.+2.*mass/rr)*(1.-pow(confo_bh,4.)) )
00350         /3./pow(rr,4.) ;
00351       sou_bh3.std_spectral_base() ;
00352       sou_bh3.inc_dzpuis(4) ;
00353 
00354       source_bh_volm = 0.25 * (taij_quad_tot_rs + taij_quad_tot_rot)
00355         / pow(confo_bh,7.)
00356         - 2. * (sou_bh1 + sou_bh2 + sou_bh3) ;
00357 
00358       source_bh_volm.std_spectral_base() ;
00359       source_bh_volm.annule_domain(0) ;
00360 
00361       integ_bh_v = source_bh_volm.integrale() ;
00362 
00363       //-------------------------------------
00364       //     Integration over the NS map
00365       //-------------------------------------
00366 
00367       // Volume integral <- dzpuis should be 4
00368       // -------------------------------------
00369       source_ns_volm = pow(confo_ns, 5.) * ener_euler ;
00370       source_ns_volm.std_spectral_base() ;
00371       source_ns_volm.inc_dzpuis(4) ;
00372 
00373       integ_ns_v = source_ns_volm.integrale() ;
00374 
00375       cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
00376            << "  integ_bh_v : "
00377            << integ_bh_v/ qpig / msol
00378            << "  integ_ns_v : " << integ_ns_v/ msol << endl ;
00379 
00380       //------------------
00381       //     ADM mass
00382       //------------------
00383       madm = hole.get_mass_bh()
00384         + (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
00385 
00386       cout << "ADM mass :           " << madm / msol << " [Mo]"
00387            << endl ;
00388 
00389       // ADM mass by surface integral at infinity : dzpuis should be 2
00390       // ----------------------------------------
00391       double mmm = hole.get_mass_bh()
00392         - 2.*(mp_aff.integrale_surface_infini(lldconf))/qpig ;
00393 
00394       cout << "Another ADM mass :   " << mmm / msol << " [Mo]"
00395            << endl ;
00396       */
00397     }
00398     else { // Isotropic coordinates with the maximal slicing
00399 
00400       //-------------------------------------
00401       //     Integration over the BH map
00402       //-------------------------------------
00403 
00404       // Sets C/M^2 for each case of the lapse boundary condition
00405       // --------------------------------------------------------
00406       double cc ;
00407 
00408       if (hole.has_bc_lapconf_nd()) {  // Neumann boundary condition
00409           if (hole.has_bc_lapconf_fs()) {  // First condition
00410               // d(\alpha \psi)/dr = 0
00411               // ---------------------
00412               cc = 2. * (sqrt(13.) - 1.) / 3. ;
00413           }
00414           else {  // Second condition
00415               // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
00416               // -----------------------------------------
00417               cc = 4. / 3. ;
00418           }
00419       }
00420       else {  // Dirichlet boundary condition
00421           if (hole.has_bc_lapconf_fs()) {  // First condition
00422               // (\alpha \psi) = 1/2
00423               // -------------------
00424           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00425           abort() ;
00426           }
00427           else {  // Second condition
00428               // (\alpha \psi) = 1/sqrt(2.) \psi_KS
00429               // ----------------------------------
00430           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00431           abort() ;
00432           //              cc = 2. * sqrt(2.) ;
00433           }
00434       }
00435 
00436       Scalar r_are(mp_bh) ;
00437       r_are = hole.r_coord(hole.has_bc_lapconf_nd(),
00438                    hole.has_bc_lapconf_fs()) ;
00439       r_are.std_spectral_base() ;
00440 
00441       // Surface integral <- dzpuis should be 0
00442       // --------------------------------------
00443       Scalar divshift_zero(divshift) ;
00444       divshift_zero.dec_dzpuis(2) ;
00445 
00446       Scalar lldllsh_zero(lldllsh) ;
00447       lldllsh_zero.dec_dzpuis(2) ;
00448 
00449       source_bh_surf = confo_bh / rr
00450         - pow(confo_bh, 4.) * (divshift_zero - 3.*lldllsh_zero)
00451         /6./lapconf_bh
00452         - pow(confo_bh, 4.) * mass * mass * cc
00453         * sqrt(1. -2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
00454         / lapconf_bh / pow(r_are*rr,3.) ;
00455 
00456       source_bh_surf.std_spectral_base() ;
00457       source_bh_surf.annule_domain(0) ;
00458       source_bh_surf.raccord(1) ;
00459 
00460       integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
00461 
00462       // Volume integral <- dzpuis should be 4
00463       // -------------------------------------
00464       Scalar sou_bh1(mp_bh) ;
00465       sou_bh1 = 1.5 * pow(confo_bh,7.) * pow(mass*mass*cc,2.)
00466         * (1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
00467         / lapconf_bh / lapconf_bh / pow(r_are*rr,6.) ;
00468       sou_bh1.std_spectral_base() ;
00469       sou_bh1.inc_dzpuis(4) ;
00470 
00471       source_bh_volm = 0.25 * taij_quad_auto_bh / pow(confo_bh,7.)
00472         + 0.25 * taij_quad_comp_bh
00473         * (pow(confo_bh/(confo_bh_comp+0.5),7.)
00474            *pow((lapconf_bh_comp+0.5)/lapconf_bh,2.)
00475            - 1.) / pow(confo_bh_comp+0.5,7.)
00476         + sou_bh1 ;
00477 
00478       source_bh_volm.std_spectral_base() ;
00479       source_bh_volm.annule_domain(0) ;
00480 
00481       integ_bh_v = source_bh_volm.integrale() ;
00482 
00483       //-------------------------------------
00484       //     Integration over the NS map
00485       //-------------------------------------
00486 
00487       // Volume integral <- dzpuis should be 4
00488       // -------------------------------------
00489       Scalar sou_ns1(mp_ns) ;
00490       sou_ns1 = pow(confo_ns, 5.) * ener_euler ;
00491       sou_ns1.std_spectral_base() ;
00492       sou_ns1.inc_dzpuis(4) ;
00493 
00494       source_ns_volm = 0.25 * taij_quad_auto_ns
00495         / pow(confo_ns_auto+0.5, 7.) / qpig + sou_ns1 ;
00496 
00497       source_ns_volm.std_spectral_base() ;
00498 
00499       integ_ns_v = source_ns_volm.integrale() ;
00500 
00501       cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
00502            << "  integ_bh_v : "
00503            << integ_bh_v/ qpig / msol
00504            << "  integ_ns_v : " << integ_ns_v/ msol << endl ;
00505 
00506       //------------------
00507       //     ADM mass
00508       //------------------
00509       madm = (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
00510 
00511       cout << "ADM mass (volume) :    " << madm / msol << " [Mo]"
00512            << endl ;
00513 
00514     }
00515 
00516     p_mass_adm_bhns_vol = new double( madm ) ;
00517 
00518     }
00519 
00520     return *p_mass_adm_bhns_vol ;
00521 
00522 }
00523 
00524 
00525 
00526                //------------------------------//
00527                //          Komar mass          //
00528                //------------------------------//
00529 
00530 double Bin_bhns::mass_kom_bhns_surf() const {
00531 
00532     // Fundamental constants and units
00533     // -------------------------------
00534     using namespace Unites ;
00535 
00536     if (p_mass_kom_bhns_surf == 0x0) {   // a new computation is required
00537 
00538         double mkom ;
00539 
00540     const Map& mp_bh = hole.get_mp() ;
00541         const Map& mp_ns = star.get_mp() ;
00542 
00543     Map_af mp_aff(mp_bh) ;
00544     Map_af mp_ns_aff(mp_ns) ;
00545 
00546     Scalar rr(mp_bh) ;
00547     rr = mp_bh.r ;
00548     rr.std_spectral_base() ;
00549 
00550     double mass = ggrav * hole.get_mass_bh() ;
00551 
00552     if (hole.is_kerrschild()) {
00553 
00554       cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
00555       abort() ;
00556 
00557     }
00558     else { // Isotropic coordinates with the maximal slicing
00559 
00560       //-------------------------------------
00561       //     Integration over the BH map
00562       //-------------------------------------
00563 
00564       // Sets C/M^2 for each case of the lapse boundary condition
00565       // --------------------------------------------------------
00566       double cc ;
00567 
00568       if (hole.has_bc_lapconf_nd()) {  // Neumann boundary condition
00569           if (hole.has_bc_lapconf_fs()) {  // First condition
00570               // d(\alpha \psi)/dr = 0
00571               // ---------------------
00572               cc = 2. * (sqrt(13.) - 1.) / 3. ;
00573           }
00574           else {  // Second condition
00575               // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
00576               // -----------------------------------------
00577               cc = 4. / 3. ;
00578           }
00579       }
00580       else {  // Dirichlet boundary condition
00581           if (hole.has_bc_lapconf_fs()) {  // First condition
00582               // (\alpha \psi) = 1/2
00583               // -------------------
00584           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00585           abort() ;
00586           }
00587           else {  // Second condition
00588               // (\alpha \psi) = 1/sqrt(2.) \psi_KS
00589               // ----------------------------------
00590           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00591           abort() ;
00592           //              cc = 2. * sqrt(2.) ;
00593           }
00594       }
00595 
00596       Scalar r_are(mp_bh) ;
00597       r_are = hole.r_coord(hole.has_bc_lapconf_nd(),
00598                    hole.has_bc_lapconf_fs()) ;
00599       r_are.std_spectral_base() ;
00600 
00601       // Komar mass by surface integral at infinity : dzpuis should be 2
00602       // ------------------------------------------
00603       const Scalar& lapconf_bh_auto_rs = hole.get_lapconf_auto_rs() ;
00604       const Scalar& confo_bh_auto_rs = hole.get_confo_auto_rs() ;
00605 
00606       Scalar lldlap_bh(mp_bh) ;
00607       lldlap_bh = lapconf_bh_auto_rs.dsdr()
00608         - confo_bh_auto_rs.dsdr() ; // dzpuis = 2
00609       lldlap_bh.std_spectral_base() ;
00610 
00611       Scalar anoth(mp_bh) ;
00612       anoth = sqrt(r_are) * (1. - 1.5 *cc*cc*pow(mass/r_are/rr,4.)
00613                  - sqrt(1. - 2.*mass/r_are/rr
00614                     + cc*cc*pow(mass/r_are/rr,4.))) / rr ;
00615       anoth.std_spectral_base() ;
00616       anoth.annule_domain(0) ;
00617       anoth.raccord(1) ;
00618       anoth.inc_dzpuis(2) ;
00619 
00620       const Scalar& lapconf_ns_auto = star.get_lapconf_auto() ;
00621       const Scalar& confo_ns_auto = star.get_confo_auto() ;
00622 
00623       Scalar lldlap_ns(mp_ns) ;
00624       lldlap_ns = lapconf_ns_auto.dsdr() - confo_ns_auto.dsdr() ;
00625       lldlap_ns.std_spectral_base() ; // dzpuis = 2
00626 
00627       mkom =
00628         (mp_aff.integrale_surface_infini(lldlap_bh+anoth))/qpig
00629         + (mp_ns_aff.integrale_surface_infini(lldlap_ns))/qpig ;
00630 
00631       cout << "Komar mass (surface) : " << mkom / msol << " [Mo]"
00632            << endl ;
00633 
00634     }
00635 
00636     p_mass_kom_bhns_surf = new double( mkom ) ;
00637 
00638     }
00639 
00640     return *p_mass_kom_bhns_surf ;
00641 
00642 }
00643 
00644 
00645 
00646                //------------------------------//
00647                //          Komar mass          //
00648                //------------------------------//
00649 
00650 double Bin_bhns::mass_kom_bhns_vol() const {
00651 
00652     // Fundamental constants and units
00653     // -------------------------------
00654     using namespace Unites ;
00655 
00656     if (p_mass_kom_bhns_vol == 0x0) {   // a new computation is required
00657 
00658         double mkom ;
00659     double integ_bh_s ;
00660     double integ_bh_v ;
00661     double integ_ns_v ;
00662 
00663     const Map& mp_bh = hole.get_mp() ;
00664         const Map& mp_ns = star.get_mp() ;
00665 
00666     double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
00667     Map_af mp_aff(mp_bh) ;
00668 
00669     Map_af mp_ns_aff(mp_ns) ;
00670 
00671     Scalar source_bh_surf(mp_bh) ;
00672     source_bh_surf.set_etat_qcq() ;
00673 
00674     Scalar source_bh_volm(mp_bh) ;
00675     source_bh_volm.set_etat_qcq() ;
00676 
00677     Scalar source_ns_volm(mp_ns) ;
00678     source_ns_volm.set_etat_qcq() ;
00679 
00680     Scalar rr(mp_bh) ;
00681     rr = mp_bh.r ;
00682     rr.std_spectral_base() ;
00683     Scalar st(mp_bh) ;
00684     st = mp_bh.sint ;
00685     st.std_spectral_base() ;
00686     Scalar ct(mp_bh) ;
00687     ct = mp_bh.cost ;
00688     ct.std_spectral_base() ;
00689     Scalar sp(mp_bh) ;
00690     sp = mp_bh.sinp ;
00691     sp.std_spectral_base() ;
00692     Scalar cp(mp_bh) ;
00693     cp = mp_bh.cosp ;
00694     cp.std_spectral_base() ;
00695 
00696     Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
00697     ll.set_etat_qcq() ;
00698     ll.set(1) = st % cp ;
00699     ll.set(2) = st % sp ;
00700     ll.set(3) = ct ;
00701     ll.std_spectral_base() ;
00702 
00703     const Scalar& lapconf_bh = hole.get_lapconf_tot() ;
00704     const Scalar& lapconf_bh_auto_rs = hole.get_lapconf_auto_rs() ;
00705     const Scalar& lapconf_bh_comp = hole.get_lapconf_comp() ;
00706     const Vector& dlapconf_bh_comp = hole.get_d_lapconf_comp() ;
00707     const Scalar& confo_bh = hole.get_confo_tot() ;
00708     const Scalar& confo_bh_auto_rs = hole.get_confo_auto_rs() ;
00709     const Scalar& confo_bh_comp = hole.get_confo_comp() ;
00710     const Vector& dconfo_bh_comp = hole.get_d_confo_comp() ;
00711     const Scalar& taij_quad_tot_rs = hole.get_taij_quad_tot_rs() ;
00712     const Scalar& taij_quad_tot_rot = hole.get_taij_quad_tot_rot() ;
00713 
00714     const Scalar& taij_quad_auto_bh = hole.get_taij_quad_auto() ;
00715     const Scalar& taij_quad_comp_bh = hole.get_taij_quad_comp() ;
00716 
00717     const Vector& shift_auto_rs = hole.get_shift_auto_rs() ;
00718     const Vector& shift_comp = hole.get_shift_comp() ;
00719     const Tensor& dshift_comp = hole.get_d_shift_comp() ;
00720 
00721     Scalar divshift(mp_bh) ;  // dzpuis = 2
00722     divshift = shift_auto_rs(1).deriv(1) + shift_auto_rs(2).deriv(2)
00723       + shift_auto_rs(3).deriv(3) + dshift_comp(1,1)
00724       + dshift_comp(2,2) + dshift_comp(3,3) ;
00725     divshift.std_spectral_base() ;
00726 
00727     Scalar llshift_auto(mp_bh) ;   // dzpuis = 0
00728         llshift_auto = ll(1)%shift_auto_rs(1) + ll(2)%shift_auto_rs(2)
00729           + ll(3)%shift_auto_rs(3) ;
00730         llshift_auto.std_spectral_base() ;
00731 
00732     Scalar lldllsh = llshift_auto.dsdr()
00733           + ll(1) * (ll(1) % dshift_comp(1,1) + ll(2) % dshift_comp(1,2)
00734                      + ll(3) % dshift_comp(1,3))
00735           + ll(2) * (ll(1) % dshift_comp(2,1) + ll(2) % dshift_comp(2,2)
00736                      + ll(3) % dshift_comp(2,3))
00737           + ll(3) * (ll(1) % dshift_comp(3,1) + ll(2) % dshift_comp(3,2)
00738                      + ll(3) % dshift_comp(3,3)) ;  // dzpuis = 2
00739         lldllsh.std_spectral_base() ;
00740 
00741     const Scalar& lapconf_ns = star.get_lapconf_tot() ;
00742     const Scalar& ener_euler = star.get_ener_euler() ;
00743     const Scalar& s_euler = star.get_s_euler() ;
00744 
00745     const Scalar& confo_ns = star.get_confo_tot() ;
00746     const Scalar& lapconf_ns_auto = star.get_lapconf_auto() ;
00747     const Scalar& confo_ns_auto = star.get_confo_auto() ;
00748     const Vector& dconfo_ns_comp = star.get_d_confo_comp() ;
00749     const Scalar& taij_quad_auto_ns = star.get_taij_quad_auto() ;
00750 
00751     const Vector& dlapconf_bh_auto_rs = hole.get_d_lapconf_auto_rs() ;
00752     /*
00753     Vector dlc_bh_auto_rs(mp_bh, COV, mp_bh.get_bvect_cart()) ;
00754     dlc_bh_auto_rs.set_etat_qcq() ;
00755     for (int i=1; i<=3; i++) {
00756       dlc_bh_auto_rs.set(i) = lapconf_bh_auto_rs.deriv(i) ;
00757     }
00758     dlc_bh_auto_rs.std_spectral_base() ;
00759     */
00760 
00761     const Vector& dconfo_bh_auto_rs = hole.get_d_confo_auto_rs() ;
00762     /*
00763     Vector dc_bh_auto_rs(mp_bh, COV, mp_bh.get_bvect_cart()) ;
00764     dc_bh_auto_rs.set_etat_qcq() ;
00765     for (int i=1; i<=3; i++) {
00766       dc_bh_auto_rs.set(i) = confo_bh_auto_rs.deriv(i) ;
00767     }
00768     dc_bh_auto_rs.std_spectral_base() ;
00769     */
00770 
00771     const Vector& dlapconf_ns_auto = star.get_d_lapconf_auto() ;
00772     /*
00773     Vector dlc_ns_auto(mp_ns, COV, mp_ns.get_bvect_cart()) ;
00774     dlc_ns_auto.set_etat_qcq() ;
00775     for (int i=1; i<=3; i++) {
00776       dlc_ns_auto.set(i) = lapconf_ns_auto.deriv(i) ;
00777     }
00778     dlc_ns_auto.std_spectral_base() ;
00779     */
00780 
00781     const Vector& dconfo_ns_auto = star.get_d_confo_auto() ;
00782     /*
00783     Vector dc_ns_auto(mp_ns, COV, mp_ns.get_bvect_cart()) ;
00784     dc_ns_auto.set_etat_qcq() ;
00785     for (int i=1; i<=3; i++) {
00786       dc_ns_auto.set(i) = confo_ns_auto.deriv(i) ;
00787     }
00788     dc_ns_auto.std_spectral_base() ;
00789     */
00790     double mass = ggrav * hole.get_mass_bh() ;
00791 
00792     if (hole.is_kerrschild()) {
00793 
00794       cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
00795       abort() ;
00796       /*
00797       const Vector& shift_auto_rs = hole.get_shift_auto_rs() ;
00798       const Vector& shift_comp = hole.get_shift_comp() ;
00799 
00800       Scalar llshift(mp_bh) ;  // dzpuis = 0
00801       llshift = ll(1) % (shift_auto_rs(1) + shift_comp(1))
00802         + ll(2) % (shift_auto_rs(2) + shift_comp(2))
00803         + ll(3) % (shift_auto_rs(3) + shift_comp(3)) ;
00804       llshift.std_spectral_base() ;
00805 
00806       Scalar lldlldlap(mp_bh) ;  // dzpuis = 3
00807       lldlldlap = lldlap.dsdr() ;
00808       lldlldlap.std_spectral_base() ;
00809 
00810       Scalar lap_bh2(mp_bh) ;
00811       lap_bh2 = 1./(1.+2.*mass/rr) ;
00812       lap_bh2.std_spectral_base() ;
00813 
00814       Scalar omelld(mp_bh) ;
00815       omelld = omega * (ll(2) * (mp_bh.get_ori_x() - x_rot)
00816                 - ll(1) * (mp_bh.get_ori_y() - y_rot)) ;
00817       omelld.std_spectral_base() ;
00818 
00819       //-------------------------------------
00820       //     Integration over the BH map
00821       //-------------------------------------
00822 
00823       // Surface integral <- dzpuis should be 0
00824       // --------------------------------------
00825       source_bh_surf = lldlap ;
00826 
00827       source_bh_surf.std_spectral_base() ;
00828       source_bh_surf.annule_domain(0) ;
00829       source_bh_surf.raccord(1) ;
00830       source_bh_surf.dec_dzpuis(2) ;
00831 
00832       integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
00833 
00834       // Volume integral <- dzpuis should be 4
00835       // -------------------------------------
00836       Scalar sou_bh1(mp_bh) ;
00837       sou_bh1 = -2. * pow(lap_bh2,2.5) * mass * lldlconf / rr / rr ;
00838       sou_bh1.std_spectral_base() ;
00839       sou_bh1.inc_dzpuis(2) ;
00840 
00841       Scalar sou_bh2(mp_bh) ;
00842       sou_bh2 = 4.*mass*mass*pow(lap_bh2,3.)*(1.+3.*mass/rr)
00843         *(1.+3.*mass/rr)*(lapse_bh_auto_rs+lapse_bh_comp)
00844         *pow(confo_bh,4.)/3./pow(rr,4.) ;
00845       sou_bh2.std_spectral_base() ;
00846       sou_bh2.inc_dzpuis(4) ;
00847 
00848       Scalar sou_bh3(mp_bh) ;
00849       sou_bh3 = - 2.*mass*pow(lap_bh2,2.5)
00850         *(2.+10.*mass/rr+9.*mass*mass/rr/rr)
00851         *pow(confo_bh,4.)*(llshift+omelld)/pow(rr,3.) ;
00852       sou_bh3.std_spectral_base() ;
00853       sou_bh3.inc_dzpuis(4) ;
00854 
00855       Scalar sou_bh4(mp_bh) ;
00856       sou_bh4 = 2. * lap_bh2 * mass * lldlldlap / rr ;
00857       sou_bh4.std_spectral_base() ;
00858       sou_bh4.inc_dzpuis(1) ;
00859 
00860       Scalar sou_bh5(mp_bh) ;
00861       sou_bh5 = lap_bh2*lap_bh2*mass*(3.+8.*mass/rr)*lldlap/rr/rr ;
00862       sou_bh5.std_spectral_base() ;
00863       sou_bh5.inc_dzpuis(2) ;
00864 
00865       Scalar sou_bh6(mp_bh) ;
00866       sou_bh6 = 4.*pow(lap_bh2,3.5)*mass*mass
00867         *(2.*(sqrt(lap_bh2)/lapse_bh - 1.)*pow(confo_bh,4.)
00868           *(4.+12.*mass/rr+9.*mass*mass/rr/rr) + 3.*(pow(confo_bh,4.)-1.))
00869         /3./pow(rr,4.) ;
00870       sou_bh6.std_spectral_base() ;
00871       sou_bh6.inc_dzpuis(4) ;
00872 
00873       source_bh_volm = lapse_bh * (taij_quad_tot_rs + taij_quad_tot_rot)
00874         / pow(confo_bh,8.)
00875         - 2. * dlapdlcf + 4. * lap_bh2 * mass * lldlap * lldlconf / rr
00876         + sou_bh1 + sou_bh2 + sou_bh3 + sou_bh4 + sou_bh5 + sou_bh6 ;
00877 
00878       source_bh_volm.std_spectral_base() ;
00879       source_bh_volm.annule_domain(0) ;
00880 
00881       integ_bh_v = source_bh_volm.integrale() ;
00882 
00883       //-------------------------------------
00884       //     Integration over the NS map
00885       //-------------------------------------
00886 
00887       // Volume integral <- dzpuis should be 4
00888       // -------------------------------------
00889       source_ns_volm = lapse_ns * psi4_ns * (ener_euler + s_euler) ;
00890       source_ns_volm.std_spectral_base() ;
00891       source_ns_volm.inc_dzpuis(4) ;
00892 
00893       integ_ns_v = source_ns_volm.integrale() ;
00894 
00895       cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
00896            << "  integ_bh_v : "
00897            << integ_bh_v/ qpig / msol
00898            << "  integ_ns_v : " << integ_ns_v/ msol << endl ;
00899 
00900       //--------------------
00901       //     Komar mass
00902       //--------------------
00903       mkom = hole.get_mass_bh()
00904         + (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
00905 
00906       cout << "Komar mass :         " << mkom / msol << " [Mo]"
00907            << endl ;
00908 
00909       // Komar mass by surface integral at infinity : dzpuis should be 2
00910       // ------------------------------------------
00911       double mmm = hole.get_mass_bh()
00912         + (mp_aff.integrale_surface_infini(lldlap))/qpig ;
00913 
00914       cout << "Another Komar mass : " << mmm / msol << " [Mo]"  << endl ;
00915       */
00916     }
00917     else { // Isotropic coordinates with the maximal slicing
00918 
00919       //-------------------------------------
00920       //     Integration over the BH map
00921       //-------------------------------------
00922 
00923       // Sets C/M^2 for each case of the lapse boundary condition
00924       // --------------------------------------------------------
00925       double cc ;
00926 
00927       if (hole.has_bc_lapconf_nd()) {  // Neumann boundary condition
00928           if (hole.has_bc_lapconf_fs()) {  // First condition
00929               // d(\alpha \psi)/dr = 0
00930               // ---------------------
00931               cc = 2. * (sqrt(13.) - 1.) / 3. ;
00932           }
00933           else {  // Second condition
00934               // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
00935               // -----------------------------------------
00936               cc = 4. / 3. ;
00937           }
00938       }
00939       else {  // Dirichlet boundary condition
00940           if (hole.has_bc_lapconf_fs()) {  // First condition
00941               // (\alpha \psi) = 1/2
00942               // -------------------
00943           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00944           abort() ;
00945           }
00946           else {  // Second condition
00947               // (\alpha \psi) = 1/sqrt(2.) \psi_KS
00948               // ----------------------------------
00949           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
00950           abort() ;
00951           //              cc = 2. * sqrt(2.) ;
00952           }
00953       }
00954 
00955       Scalar r_are(mp_bh) ;
00956       r_are = hole.r_coord(hole.has_bc_lapconf_nd(),
00957                    hole.has_bc_lapconf_fs()) ;
00958       r_are.std_spectral_base() ;
00959 
00960       Scalar lldlapconf_is(mp_bh) ;
00961       lldlapconf_is = ll(1)%dlapconf_bh_auto_rs(1)
00962         + ll(2)%dlapconf_bh_auto_rs(2) + ll(3)%dlapconf_bh_auto_rs(3)
00963         + ll(1)%dlapconf_bh_comp(1) + ll(2)%dlapconf_bh_comp(2)
00964         + ll(3)%dlapconf_bh_comp(3) ;
00965       // dzpuis = 2
00966       lldlapconf_is.std_spectral_base() ;
00967 
00968       // Surface integral <- dzpuis should be 0
00969       // --------------------------------------
00970       Scalar divshift_zero(divshift) ;
00971       divshift_zero.dec_dzpuis(2) ;
00972 
00973       Scalar lldllsh_zero(lldllsh) ;
00974       lldllsh_zero.dec_dzpuis(2) ;
00975 
00976       Scalar sou_bh_s_psi(mp_bh) ;
00977       sou_bh_s_psi = 0.5 * confo_bh / rr
00978         - pow(confo_bh, 4.) * (divshift_zero - 3.*lldllsh_zero)
00979         / 12. / lapconf_bh
00980         - 0.5 * pow(confo_bh, 4.) * mass * mass * cc
00981         * sqrt(1. -2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
00982         / lapconf_bh / pow(r_are*rr,3.) ;
00983 
00984       sou_bh_s_psi.std_spectral_base() ;
00985       sou_bh_s_psi.annule_domain(0) ;
00986       sou_bh_s_psi.raccord(1) ;
00987 
00988       if (hole.has_bc_lapconf_nd()) {  // Neumann boundary condition
00989           if (hole.has_bc_lapconf_fs()) {  // First condition
00990 
00991         source_bh_surf = sou_bh_s_psi ;
00992 
00993         source_bh_surf.std_spectral_base() ;
00994         source_bh_surf.annule_domain(0) ;
00995         source_bh_surf.raccord(1) ;
00996 
00997           }
00998           else {
00999 
01000         Scalar sou_bh_s1(mp_bh) ;
01001         sou_bh_s1 = 0.5 * lapconf_bh / rr ;
01002         sou_bh_s1.std_spectral_base() ;
01003         sou_bh_s1.annule_domain(0) ;
01004         sou_bh_s1.raccord(1) ;
01005 
01006         source_bh_surf = sou_bh_s1 + sou_bh_s_psi ;
01007 
01008         source_bh_surf.std_spectral_base() ;
01009         source_bh_surf.annule_domain(0) ;
01010         source_bh_surf.raccord(1) ;
01011 
01012           }
01013       }
01014       else {
01015 
01016           Scalar sou_bh_s1(mp_bh) ;
01017           sou_bh_s1 = lldlapconf_is ;
01018           sou_bh_s1.std_spectral_base() ;
01019           sou_bh_s1.dec_dzpuis(2) ;
01020 
01021           Scalar sou_bh_s2(mp_bh) ;
01022           sou_bh_s2 = 0.5 * sqrt(r_are)
01023         * (1. - 3.*cc*cc*pow(mass/r_are/rr,4.)
01024            -sqrt(1. - 2.*mass/r_are/rr
01025              + cc*cc*pow(mass/r_are/rr,4.))) / rr ;
01026 
01027           sou_bh_s2.std_spectral_base() ;
01028 
01029           source_bh_surf = sou_bh_s1 + sou_bh_s2 + sou_bh_s_psi ;
01030 
01031           source_bh_surf.std_spectral_base() ;
01032           source_bh_surf.annule_domain(0) ;
01033           source_bh_surf.raccord(1) ;
01034 
01035       }
01036 
01037       integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
01038 
01039       // Volume integral <- dzpuis should be 4
01040       // -------------------------------------
01041       Scalar sou_bh1(mp_bh) ;
01042       sou_bh1 = 0.75 * pow(mass*mass*cc,2.)
01043         * (1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
01044         * (7.*pow(confo_bh,6.)/lapconf_bh
01045            + pow(confo_bh,7.)/lapconf_bh/lapconf_bh)
01046         / pow(r_are*rr,6.) ;
01047 
01048       sou_bh1.std_spectral_base() ;
01049       sou_bh1.inc_dzpuis(4) ;
01050 
01051       Scalar sou_bh2(mp_bh) ;
01052       sou_bh2 = 0.125 * (7.*lapconf_bh/pow(confo_bh,8.)
01053                  + 1./pow(confo_bh,7.)) * taij_quad_auto_bh ;
01054 
01055       sou_bh2.std_spectral_base() ;
01056 
01057       Scalar sou_bh3(mp_bh) ;
01058       sou_bh3 = 0.125 * (7.*((lapconf_bh_comp+0.5)/lapconf_bh
01059                  * pow(confo_bh/(confo_bh_comp+0.5),6.) - 1.)
01060                  * (lapconf_bh_comp+0.5)
01061                  / pow(confo_bh_comp+0.5,8.)
01062                  + (pow(confo_bh/(confo_bh_comp+0.5),7.)
01063                 *pow((lapconf_bh_comp+0.5)/lapconf_bh,2.)
01064                 - 1.) / pow(confo_bh_comp+0.5,7))
01065         * taij_quad_comp_bh ;
01066 
01067       sou_bh3.std_spectral_base() ;
01068 
01069       source_bh_volm = sou_bh1 + sou_bh2 + sou_bh3 ;
01070       source_bh_volm.std_spectral_base() ;
01071       source_bh_volm.annule_domain(0) ;
01072 
01073       integ_bh_v = source_bh_volm.integrale() ;
01074 
01075       //-------------------------------------
01076       //     Integration over the NS map
01077       //-------------------------------------
01078 
01079       // Volume integral <- dzpuis should be 4
01080       // -------------------------------------
01081       Scalar sou_ns1(mp_ns) ;
01082       sou_ns1 = 0.5 * pow(confo_ns,4.) * (lapconf_ns + confo_ns)
01083         * ener_euler + lapconf_ns * pow(confo_ns,4.) * s_euler ;
01084       sou_ns1.std_spectral_base() ;
01085       sou_ns1.inc_dzpuis(4) ;
01086 
01087       Scalar sou_ns2(mp_ns) ;
01088       sou_ns2 = 0.125 * (7.*(lapconf_ns_auto+0.5)/(confo_ns_auto+0.5)
01089                  + 1.) * taij_quad_auto_ns
01090         / pow(confo_ns_auto+0.5, 7.) / qpig ;
01091       sou_ns2.std_spectral_base() ;
01092 
01093       source_ns_volm = sou_ns1 + sou_ns2 ;
01094       source_ns_volm.std_spectral_base() ;
01095 
01096       integ_ns_v = source_ns_volm.integrale() ;
01097 
01098       cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
01099            << "  integ_bh_v : "
01100            << integ_bh_v/ qpig / msol
01101            << "  integ_ns_v : " << integ_ns_v/ msol << endl ;
01102 
01103       //--------------------
01104       //     Komar mass
01105       //--------------------
01106       mkom = (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
01107 
01108       cout << "Komar mass (volume) :  " << mkom / msol << " [Mo]"
01109            << endl ;
01110 
01111     }
01112 
01113     p_mass_kom_bhns_vol = new double( mkom ) ;
01114 
01115     }
01116 
01117     return *p_mass_kom_bhns_vol ;
01118 
01119 }
01120 
01121 
01122                //-----------------------------------------//
01123                //          Total linear momentum          //
01124                //-----------------------------------------//
01125 
01126 const Tbl& Bin_bhns::line_mom_bhns() const {
01127 
01128     // Fundamental constants and units
01129     // -------------------------------
01130     using namespace Unites ;
01131 
01132     if (p_line_mom_bhns == 0x0) {   // a new computation is required
01133 
01134         p_line_mom_bhns = new Tbl(3) ;
01135     p_line_mom_bhns->annule_hard() ;  // fills the double array with zeros
01136 
01137     if (hole.is_kerrschild()) {
01138 
01139       // Construction of a new grid and a new affine mapping
01140       // ---------------------------------------------------
01141       int nz = (hole.get_mp()).get_mg()->get_nzone() ;
01142       double* bornes = new double[nz+1] ;
01143       double radius = separ + 2. ;
01144 
01145       for (int i=nz-1; i>0; i--) {
01146         bornes[i] = radius ;
01147         radius /= 2. ;
01148       }
01149       bornes[0] = 0. ;
01150       bornes[nz] = __infinity ;
01151 
01152       Map_af mp_aff(*((hole.get_mp()).get_mg()), bornes) ;
01153       mp_aff.set_ori(0.,0.,0.) ;
01154 
01155       delete [] bornes ;
01156 
01157       // Construction of the extrinsic curvature
01158       // ---------------------------------------
01159       Vector shift_bh(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01160       shift_bh.set_etat_qcq() ;
01161 
01162       Vector shift_ns(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01163       shift_ns.set_etat_qcq() ;
01164 
01165       shift_bh.set(1).import(hole.get_shift_auto()(1)) ;
01166       shift_bh.set(2).import(hole.get_shift_auto()(2)) ;
01167       shift_bh.set(3).import(hole.get_shift_auto()(3)) ;
01168 
01169       shift_ns.set(1).import(star.get_shift_auto()(1)) ;
01170       shift_ns.set(2).import(star.get_shift_auto()(2)) ;
01171       shift_ns.set(3).import(star.get_shift_auto()(3)) ;
01172 
01173       Vector shift_tot = shift_bh + shift_ns ;
01174       shift_tot.std_spectral_base() ;
01175       shift_tot.annule(0, nz-2) ;
01176 
01177       Scalar stc(mp_aff) ;
01178       stc = mp_aff.sint ;
01179       stc.std_spectral_base() ;
01180       Scalar ctc(mp_aff) ;
01181       ctc = mp_aff.cost ;
01182       ctc.std_spectral_base() ;
01183       Scalar spc(mp_aff) ;
01184       spc = mp_aff.sinp ;
01185       spc.std_spectral_base() ;
01186       Scalar cpc(mp_aff) ;
01187       cpc = mp_aff.cosp ;
01188       cpc.std_spectral_base() ;
01189 
01190       Vector lc(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01191       lc.set_etat_qcq() ;
01192       lc.set(1) = stc * cpc ;
01193       lc.set(2) = stc * spc ;
01194       lc.set(3) = ctc ;
01195       lc.std_spectral_base() ;
01196 
01197       Vector kijlj(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01198       kijlj.set_etat_qcq() ;
01199 
01200       Scalar rr(mp_aff) ;
01201       rr = mp_aff.r ;
01202       rr.std_spectral_base() ;
01203 
01204       Scalar xbhsr(mp_aff) ;
01205       xbhsr = (hole.get_mp()).get_ori_x() / rr ;
01206       xbhsr.std_spectral_base() ;
01207 
01208       Scalar ybhsr(mp_aff) ;
01209       ybhsr = (hole.get_mp()).get_ori_y() / rr ;
01210       ybhsr.std_spectral_base() ;
01211 
01212       Scalar rl(mp_aff) ;
01213       rl = sqrt(1. - 2.*xbhsr*lc(1) - 2.*ybhsr*lc(2)
01214             + xbhsr*xbhsr + ybhsr*ybhsr) ;
01215       rl.std_spectral_base() ;
01216 
01217       Vector ll(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01218       ll.set_etat_qcq() ;
01219       ll.set(1) = (lc(1) - xbhsr) / rl ;
01220       ll.set(2) = (lc(2) - ybhsr)/ rl ;
01221       ll.set(3) = lc(3) / rl ;
01222       ll.std_spectral_base() ;
01223 
01224       Scalar lcll(mp_aff) ;
01225       lcll = lc(1)*ll(1) + lc(2)*ll(2) + lc(3)*ll(3) ;
01226       lcll.std_spectral_base() ;
01227 
01228       Scalar divshift(mp_aff) ;
01229       divshift = shift_tot(1).deriv(1) + shift_tot(2).deriv(2)
01230         + shift_tot(3).deriv(3) ;
01231       divshift.std_spectral_base() ;
01232 
01233       Scalar llshift(mp_aff) ;
01234       llshift = ll(1)*shift_tot(1) + ll(2)*shift_tot(2)
01235         + ll(3)*shift_tot(3) ;
01236       llshift.std_spectral_base() ;
01237 
01238       Scalar lcshift(mp_aff) ;
01239       lcshift = lc(1)*shift_tot(1) + lc(2)*shift_tot(2)
01240         + lc(3)*shift_tot(3) ;
01241       lcshift.std_spectral_base() ;
01242 
01243       Vector lcdshft(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01244       for (int i=1; i<=3; i++) {
01245         lcdshft.set(i) = lc(1)*(shift_tot(1).deriv(i))
01246           + lc(2)*(shift_tot(2).deriv(i))
01247           + lc(3)*(shift_tot(3).deriv(i)) ;
01248       }
01249       lcdshft.std_spectral_base() ;
01250 
01251       Vector dshift(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01252       for (int i=1; i<=3; i++) {
01253           dshift.set(i) = shift_tot(i).dsdr() ;
01254       }
01255       dshift.std_spectral_base() ;
01256 
01257       Vector lldshft(mp_aff, CON, mp_aff.get_bvect_cart()) ;
01258       for (int i=1; i<=3; i++) {
01259         lldshft.set(i) = ll(1)*(shift_tot(i).deriv(1))
01260           + ll(2)*(shift_tot(i).deriv(2))
01261           + ll(3)*(shift_tot(i).deriv(3)) ;
01262       }
01263       lldshft.std_spectral_base() ;
01264 
01265       double mass = ggrav * hole.get_mass_bh() ;
01266 
01267       Scalar lap_bh2(mp_aff) ;
01268       lap_bh2 = 1./(1.+2.*mass/rl/rr) ;
01269       lap_bh2.std_spectral_base() ;
01270 
01271       Scalar lap2hbh(mp_aff) ;
01272       lap2hbh = lap_bh2 * mass / rl / rr ;
01273       lap2hbh.std_spectral_base() ;
01274 
01275       Scalar omexsr(mp_aff) ;
01276       omexsr = omega * ((hole.get_mp()).get_ori_x() - x_rot)
01277         / rl / rr ;
01278       omexsr.std_spectral_base() ;
01279 
01280       Scalar omeysr(mp_aff) ;
01281       omeysr = omega * ((hole.get_mp()).get_ori_y() - y_rot)
01282         / rl / rr ;
01283       omeysr.std_spectral_base() ;
01284 
01285       Scalar kk(mp_aff) ;
01286       kk = 4.*sqrt(lap_bh2)*lap2hbh*(1.+3.*mass/rl/rr)/3./rl/rr ;
01287       kk.std_spectral_base() ;
01288 
01289       //-----------------------------------------------------------
01290       //     Surface integral at infinity : dzpuis should be 2
01291       //-----------------------------------------------------------
01292 
01293       // Source for X component
01294       // ----------------------
01295       Scalar kij_x1(mp_aff) ;
01296       kij_x1 = omexsr*ll(1)*lc(2) - omeysr*(ll(1)*lc(1)+lcll)
01297         + lcll*shift_tot(1)/rl/rr
01298         + ll(1)*lcshift/rl/rr
01299         + (lc(1)-lap_bh2*(9.+14.*mass/rl/rr)*ll(1)*lcll)
01300         *(llshift/rl/rr + omexsr*ll(2) - omeysr*ll(1))/3. ;
01301       kij_x1.std_spectral_base() ;
01302       kij_x1.inc_dzpuis(2) ;
01303 
01304       Scalar kij_x2(mp_aff) ;
01305       kij_x2 = kk * (lc(1) - 2.*lap2hbh*ll(1)*lcll) ;
01306       kij_x2.std_spectral_base() ;
01307       kij_x2.inc_dzpuis(2) ;
01308 
01309       kijlj.set(1) = lcdshft(1) + dshift(1) - 2.*lc(1)*divshift/3.
01310         + 2.*lap2hbh * (-ll(1)*(ll(1)*lcdshft(1) + ll(2)*lcdshft(2)
01311                     + ll(3)*lcdshft(3))
01312                 - lcll*lldshft(1)
01313                 + 2.*ll(1)*lcll*divshift/3.
01314                 + kij_x1)
01315         + kij_x2 ;
01316 
01317       // Source for Y component
01318       // ----------------------
01319       Scalar kij_y1(mp_aff) ;
01320       kij_y1 = omexsr*(ll(2)*lc(2)+lcll) - omeysr*ll(2)*lc(1)
01321         + lcll*shift_tot(2)/rl/rr
01322         + ll(2)*lcshift/rl/rr
01323         + (lc(2)-lap_bh2*(9.+14.*mass/rl/rr)*ll(2)*lcll)
01324         *(llshift/rl/rr + omexsr*ll(2) - omeysr*ll(1))/3. ;
01325       kij_y1.std_spectral_base() ;
01326       kij_y1.inc_dzpuis(2) ;
01327 
01328       Scalar kij_y2(mp_aff) ;
01329       kij_y2 = kk * (lc(2) - 2.*lap2hbh*ll(2)*lcll) ;
01330       kij_y2.std_spectral_base() ;
01331       kij_y2.inc_dzpuis(2) ;
01332 
01333       kijlj.set(2) = lcdshft(2) + dshift(2) - 2.*lc(2)*divshift/3.
01334         + 2.*lap2hbh * (-ll(2)*(ll(1)*lcdshft(1) + ll(2)*lcdshft(2)
01335                     + ll(3)*lcdshft(3))
01336                 - lcll*lldshft(2)
01337                 + 2.*ll(2)*lcll*divshift/3.
01338                 + kij_y1)
01339         + kij_y2 ;
01340 
01341       // Source for Z component
01342       // ----------------------
01343       Scalar kij_z1(mp_aff) ;
01344       kij_z1 = omexsr*ll(3)*lc(2) - omeysr*ll(3)*lc(1)
01345         + lcll*shift_tot(3)/rl/rr
01346         + ll(3)*lcshift/rl/rr
01347         + (lc(3)-lap_bh2*(9.+14.*mass/rl/rr)*ll(3)*lcll)
01348         *(llshift/rl/rr + omexsr*ll(2) - omeysr*ll(1))/3. ;
01349       kij_z1.std_spectral_base() ;
01350       kij_z1.inc_dzpuis(2) ;
01351 
01352       Scalar kij_z2(mp_aff) ;
01353       kij_z2 = kk * (lc(3) - 2.*lap2hbh*ll(3)*lcll) ;
01354       kij_z2.std_spectral_base() ;
01355       kij_z2.inc_dzpuis(2) ;
01356 
01357       kijlj.set(3) = lcdshft(3) + dshift(3) - 2.*lc(3)*divshift/3.
01358         + 2.*lap2hbh * (-ll(3)*(ll(1)*lcdshft(1) + ll(2)*lcdshft(2)
01359                     + ll(3)*lcdshft(3))
01360                 - lcll*lldshft(3)
01361                 + 2.*ll(3)*lcll*divshift/3.
01362                 + kij_z1)
01363         + kij_z2 ;
01364 
01365       kijlj.std_spectral_base() ;
01366 
01367       // X component  dzpuis should be 2
01368       // -----------
01369       double lm_x = mp_aff.integrale_surface_infini(kijlj(1)) ;
01370       p_line_mom_bhns->set(0) = 0.25 * lm_x / qpig ;
01371 
01372       // Y component
01373       // -----------
01374       double lm_y = mp_aff.integrale_surface_infini(kijlj(2)) ;
01375       p_line_mom_bhns->set(1) = 0.25 * lm_y / qpig ;
01376 
01377       // Z component
01378       // -----------
01379       double lm_z = mp_aff.integrale_surface_infini(kijlj(3)) ;
01380       p_line_mom_bhns->set(2) = 0.25 * lm_z / qpig ;
01381 
01382     }
01383     else {  // Isotropic coordinates with the maximal slicing
01384 
01385       /*
01386       // Method by the sourface integral at infinity
01387           // -------------------------------------------
01388 
01389       const Map& mp_bh = hole.get_mp() ;
01390       Map_af mp_aff(mp_bh) ;
01391 
01392       Scalar st(mp_bh) ;
01393       st = mp_bh.sint ;
01394       st.std_spectral_base() ;
01395       Scalar ct(mp_bh) ;
01396       ct = mp_bh.cost ;
01397       ct.std_spectral_base() ;
01398       Scalar sp(mp_bh) ;
01399       sp = mp_bh.sinp ;
01400       sp.std_spectral_base() ;
01401       Scalar cp(mp_bh) ;
01402       cp = mp_bh.cosp ;
01403       cp.std_spectral_base() ;
01404 
01405       Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01406       ll.set_etat_qcq() ;
01407       ll.set(1) = st * cp ;
01408       ll.set(2) = st * sp ;
01409       ll.set(3) = ct ;
01410       ll.std_spectral_base() ;
01411 
01412       const Scalar& confo_bh = hole.get_confo_tot() ;
01413       const Sym_tensor& taij_tot_rs = hole.get_taij_tot_rs() ;
01414 
01415       Vector kijlj(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01416       kijlj.set_etat_qcq() ;
01417 
01418       for (int i=1; i<=3; i++) {
01419         kijlj.set(i) = (taij_tot_rs(i,1)%ll(1)
01420                 + taij_tot_rs(i,2)%ll(2)
01421                 + taij_tot_rs(i,3)%ll(3)) / pow(confo_bh,10.) ;
01422       }
01423 
01424       kijlj.std_spectral_base() ;
01425 
01426       // X component
01427       // -----------
01428       double lm_x = mp_aff.integrale_surface_infini(kijlj(1)) ;
01429       p_line_mom_bhns->set(0) = 0.5 * lm_x / qpig ;
01430 
01431       // Y component
01432       // -----------
01433       double lm_y = mp_aff.integrale_surface_infini(kijlj(2)) ;
01434       p_line_mom_bhns->set(1) = 0.5 * lm_y / qpig ;
01435 
01436       // Z component
01437       // -----------
01438       double lm_z = mp_aff.integrale_surface_infini(kijlj(3)) ;
01439       p_line_mom_bhns->set(2) = 0.5 * lm_z / qpig ;
01440       */
01441 
01442       // Method by the volume integral and the surface integral
01443       // at the BH horizon
01444       // ------------------------------------------------------
01445 
01446       const Map& mp_bh = hole.get_mp() ;
01447       Map_af mp_aff(mp_bh) ;
01448       const Map& mp_ns = star.get_mp() ;
01449 
01450       const Sym_tensor& taij = hole.get_taij_tot() ;
01451       const Scalar& confo_ns = star.get_confo_tot() ;
01452       const Scalar& ee = star.get_ener_euler() ;
01453       const Scalar& pp = star.get_press() ;
01454       const Vector& uu = star.get_u_euler() ;
01455 
01456       double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
01457 
01458       Scalar st(mp_bh) ;
01459       st = mp_bh.sint ;
01460       st.std_spectral_base() ;
01461       Scalar ct(mp_bh) ;
01462       ct = mp_bh.cost ;
01463       ct.std_spectral_base() ;
01464       Scalar sp(mp_bh) ;
01465       sp = mp_bh.sinp ;
01466       sp.std_spectral_base() ;
01467       Scalar cp(mp_bh) ;
01468       cp = mp_bh.cosp ;
01469       cp.std_spectral_base() ;
01470 
01471       Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01472       ll.set_etat_qcq() ;
01473       ll.set(1) = st % cp ;
01474       ll.set(2) = st % sp ;
01475       ll.set(3) = ct ;
01476       ll.std_spectral_base() ;
01477 
01478       // X component
01479       // -----------
01480 
01481       //-------------------------------------
01482       //     Integration over the BH map
01483       //-------------------------------------
01484 
01485       // Surface integral <- dzpuis should be 0
01486       // --------------------------------------
01487       Scalar sou_bh_sx(mp_bh) ;
01488       sou_bh_sx = taij(1,1) * ll(1) + taij(1,2) * ll(2)
01489         + taij(1,3) * ll(3) ;
01490       sou_bh_sx.std_spectral_base() ;
01491       sou_bh_sx.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01492 
01493       sou_bh_sx.annule_domain(0) ;
01494       sou_bh_sx.raccord(1) ;
01495 
01496       double integ_bh_s_x = mp_aff.integrale_surface(sou_bh_sx,
01497                              radius_ah) ;
01498 
01499       //-------------------------------------
01500       //     Integration over the NS map
01501       //-------------------------------------
01502 
01503       // Volume integral <- dzpuis should be 4
01504       // -------------------------------------
01505       Scalar sou_ns_vx(mp_ns) ;
01506 
01507       sou_ns_vx = pow(confo_ns, 10.) * (ee + pp) * uu(1) ;
01508 
01509       sou_ns_vx.std_spectral_base() ;
01510       sou_ns_vx.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
01511 
01512       double integ_ns_v_x = sou_ns_vx.integrale() ;
01513 
01514       p_line_mom_bhns->set(0) = integ_ns_v_x + 0.5*integ_bh_s_x/qpig ;
01515 
01516       // Y component
01517       // -----------
01518 
01519       //-------------------------------------
01520       //     Integration over the BH map
01521       //-------------------------------------
01522 
01523       // Surface integral <- dzpuis should be 0
01524       // --------------------------------------
01525       Scalar sou_bh_sy(mp_bh) ;
01526       sou_bh_sy = taij(2,1) * ll(1) + taij(2,2) * ll(2)
01527         + taij(2,3) * ll(3) ;
01528       sou_bh_sy.std_spectral_base() ;
01529       sou_bh_sy.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01530 
01531       sou_bh_sy.annule_domain(0) ;
01532       sou_bh_sy.raccord(1) ;
01533 
01534       double integ_bh_s_y = mp_aff.integrale_surface(sou_bh_sy,
01535                              radius_ah) ;
01536 
01537       //-------------------------------------
01538       //     Integration over the NS map
01539       //-------------------------------------
01540 
01541       // Volume integral <- dzpuis should be 4
01542       // -------------------------------------
01543       Scalar sou_ns_vy(mp_ns) ;
01544 
01545       sou_ns_vy = pow(confo_ns, 10.) * (ee + pp) * uu(2) ;
01546 
01547       sou_ns_vy.std_spectral_base() ;
01548       sou_ns_vy.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
01549 
01550       double integ_ns_v_y = sou_ns_vy.integrale() ;
01551 
01552       p_line_mom_bhns->set(1) = integ_ns_v_y + 0.5*integ_bh_s_y/qpig ;
01553 
01554 
01555       // Z component
01556       // -----------
01557 
01558       //-------------------------------------
01559       //     Integration over the BH map
01560       //-------------------------------------
01561 
01562       // Surface integral <- dzpuis should be 0
01563       // --------------------------------------
01564       Scalar sou_bh_sz(mp_bh) ;
01565       sou_bh_sz = taij(3,1) * ll(1) + taij(3,2) * ll(2)
01566         + taij(3,3) * ll(3) ;
01567       sou_bh_sz.std_spectral_base() ;
01568       sou_bh_sz.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01569 
01570       sou_bh_sz.annule_domain(0) ;
01571       sou_bh_sz.raccord(1) ;
01572 
01573       double integ_bh_s_z = mp_aff.integrale_surface(sou_bh_sz,
01574                              radius_ah) ;
01575 
01576       //-------------------------------------
01577       //     Integration over the NS map
01578       //-------------------------------------
01579 
01580       // Volume integral <- dzpuis should be 4
01581       // -------------------------------------
01582       Scalar sou_ns_vz(mp_ns) ;
01583 
01584       sou_ns_vz = pow(confo_ns, 10.) * (ee + pp) * uu(3) ;
01585 
01586       sou_ns_vz.std_spectral_base() ;
01587       sou_ns_vz.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
01588 
01589       double integ_ns_v_z = sou_ns_vz.integrale() ;
01590 
01591       p_line_mom_bhns->set(2) = integ_ns_v_z + 0.5*integ_bh_s_z/qpig ;
01592 
01593     }
01594 
01595     }
01596 
01597     return *p_line_mom_bhns ;
01598 
01599 }
01600 
01601 
01602                //------------------------------------------//
01603                //          Total angular momentum          //
01604                //------------------------------------------//
01605 
01606 const Tbl& Bin_bhns::angu_mom_bhns() const {
01607 
01608     // Fundamental constants and units
01609     // -------------------------------
01610     using namespace Unites ;
01611 
01612     if (p_angu_mom_bhns == 0x0) {   // a new computation is required
01613 
01614         p_angu_mom_bhns = new Tbl(3) ;
01615     p_angu_mom_bhns->annule_hard() ;  // fills the double array with zeros
01616 
01617     double integ_bh_s_x ;
01618     double integ_bh_s_y ;
01619     double integ_bh_s_z ;
01620     double integ_bh_v_x ;
01621     double integ_bh_v_y ;
01622     double integ_bh_v_z ;
01623     double integ_ns_v_x ;
01624     double integ_ns_v_y ;
01625     double integ_ns_v_z ;
01626 
01627     const Map& mp_bh = hole.get_mp() ;
01628     const Map& mp_ns = star.get_mp() ;
01629 
01630     double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
01631     Map_af mp_aff(mp_bh) ;
01632 
01633     Scalar source_bh_surf(mp_bh) ;
01634     source_bh_surf.set_etat_qcq() ;
01635 
01636     Scalar source_bh_volm(mp_bh) ;
01637     source_bh_volm.set_etat_qcq() ;
01638 
01639     Scalar source_ns_volm(mp_ns) ;
01640     source_ns_volm.set_etat_qcq() ;
01641 
01642     Scalar rr(mp_bh) ;
01643     rr = mp_bh.r ;
01644     rr.std_spectral_base() ;
01645 
01646     Scalar st(mp_bh) ;
01647     st = mp_bh.sint ;
01648     st.std_spectral_base() ;
01649     Scalar ct(mp_bh) ;
01650     ct = mp_bh.cost ;
01651     ct.std_spectral_base() ;
01652     Scalar sp(mp_bh) ;
01653     sp = mp_bh.sinp ;
01654     sp.std_spectral_base() ;
01655     Scalar cp(mp_bh) ;
01656     cp = mp_bh.cosp ;
01657     cp.std_spectral_base() ;
01658 
01659     Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01660     ll.set_etat_qcq() ;
01661     ll.set(1) = st % cp ;
01662     ll.set(2) = st % sp ;
01663     ll.set(3) = ct ;
01664     ll.std_spectral_base() ;
01665 
01666     Scalar xx_bh(mp_bh) ;
01667     xx_bh = mp_bh.xa ;
01668     xx_bh.std_spectral_base() ;
01669     Scalar yy_bh(mp_bh) ;
01670     yy_bh = mp_bh.ya ;
01671     yy_bh.std_spectral_base() ;
01672     Scalar zz_bh(mp_bh) ;
01673     zz_bh = mp_bh.za ;
01674     zz_bh.std_spectral_base() ;
01675 
01676     Scalar xx_ns(mp_ns) ;
01677     xx_ns = mp_ns.xa ;
01678     xx_ns.std_spectral_base() ;
01679     Scalar yy_ns(mp_ns) ;
01680     yy_ns = mp_ns.ya ;
01681     yy_ns.std_spectral_base() ;
01682     Scalar zz_ns(mp_ns) ;
01683     zz_ns = mp_ns.za ;
01684     zz_ns.std_spectral_base() ;
01685 
01686     const Scalar& confo_bh = hole.get_confo_tot() ;
01687     const Sym_tensor& taij = hole.get_taij_tot() ;
01688     const Scalar& confo_ns = star.get_confo_tot() ;
01689     const Scalar& ee = star.get_ener_euler() ;
01690     const Scalar& pp = star.get_press() ;
01691     const Vector& uu = star.get_u_euler() ;
01692 
01693     Vector jbh_x(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01694     jbh_x.set_etat_qcq() ;
01695     for (int i=1; i<=3; i++)
01696       jbh_x.set(i) = yy_bh * taij(3,i) - zz_bh * taij(2,i) ;
01697 
01698     jbh_x.std_spectral_base() ;
01699 
01700     Vector jbh_y(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01701     jbh_y.set_etat_qcq() ;
01702     for (int i=1; i<=3; i++)
01703       jbh_y.set(i) = zz_bh * taij(1,i) - xx_bh * taij(3,i) ;
01704 
01705     jbh_y.std_spectral_base() ;
01706 
01707     Vector jbh_z(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01708     jbh_z.set_etat_qcq() ;
01709     for (int i=1; i<=3; i++)
01710       jbh_z.set(i) = xx_bh * taij(2,i) - yy_bh * taij(1,i) ;
01711 
01712     jbh_z.std_spectral_base() ;
01713 
01714     double mass = ggrav * hole.get_mass_bh() ;
01715 
01716     if (hole.is_kerrschild()) {
01717 
01718       double ori_bh = mp_bh.get_ori_x() ;
01719 
01720       Scalar lap_bh2(mp_bh) ;
01721       lap_bh2 = 1./(1.+2.*mass/rr) ;
01722       lap_bh2.std_spectral_base() ;
01723 
01724       Scalar lcnf(mp_bh) ;
01725       lcnf = log(confo_bh) ;
01726       lcnf.std_spectral_base() ;
01727 
01728       Vector jbhsr_x(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01729       jbhsr_x.set_etat_qcq() ;
01730       for (int i=1; i<=3; i++)
01731         jbhsr_x.set(i) = ll(2)*taij(3,i) - ll(3)*taij(2,i) ;
01732 
01733       jbhsr_x.std_spectral_base() ;
01734 
01735       Vector jbhsr_y(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01736       jbhsr_y.set_etat_qcq() ;
01737       for (int i=1; i<=3; i++)
01738         jbhsr_y.set(i) = ll(3)*taij(1,i) - (ll(1)+ori_bh/rr)*taij(3,i) ;
01739 
01740       jbhsr_y.std_spectral_base() ;
01741 
01742       Vector jbhsr_z(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01743       jbhsr_z.set_etat_qcq() ;
01744       for (int i=1; i<=3; i++)
01745         jbhsr_z.set(i) = (ll(1)+ori_bh/rr)*taij(2,i) - ll(2)*taij(1,i) ;
01746 
01747       jbhsr_z.std_spectral_base() ;
01748 
01749       Scalar tmp1(mp_bh) ;  // dzpuis = 0
01750       tmp1 = 2. * pow(lap_bh2,2.5) * mass * (1.+3.*mass/rr)
01751         * pow(confo_bh,6.) * ori_bh / 3. / rr / rr ;
01752       tmp1.std_spectral_base() ;
01753       tmp1.annule_domain(0) ;
01754 
01755       Scalar tmp2(mp_bh) ;  // dzpuis = 0
01756       tmp2 = 4. * sqrt(lap_bh2) * (1.+3.*mass/rr) * pow(confo_bh,6.) ;
01757       tmp2.std_spectral_base() ;
01758       tmp2.annule_domain(0) ;
01759 
01760       Scalar lltaij(mp_bh) ;  // dzpuis = 2
01761       lltaij = ll(1)*(ll(1)*taij(1,1)+ll(2)*taij(1,2)+ll(3)*taij(1,3))
01762         + ll(2)*(ll(1)*taij(2,1)+ll(2)*taij(2,2)+ll(3)*taij(2,3))
01763         + ll(3)*(ll(1)*taij(3,1)+ll(2)*taij(3,2)+ll(3)*taij(3,3)) ;
01764       lltaij.std_spectral_base() ;
01765       lltaij.dec_dzpuis(2) ; // dzpuis : 2 -> 0
01766 
01767       Scalar dlcnf(mp_bh) ;   // dzpuis = 2
01768       dlcnf = - 2. * lap_bh2 * tmp2 * mass * lcnf.dsdr() / rr ;
01769       dlcnf.std_spectral_base() ;
01770       dlcnf.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01771       dlcnf.annule_domain(0) ;
01772 
01773       Scalar tmp3(mp_bh) ;  // dzpuis = 0
01774       tmp3 = -2.*pow(lap_bh2,2.5)
01775         *(6.+32.*mass/rr+41.*mass*mass/rr/rr+24.*pow(mass,3.)/pow(rr,3.))
01776         *pow(confo_bh,6.)/3./rr
01777         + 3.* lltaij + dlcnf ;
01778       tmp3.std_spectral_base() ;
01779       tmp3.annule_domain(0) ;
01780 
01781       Scalar tmp4(mp_bh) ;  // dzpuis = 0
01782       tmp4 = lap_bh2 * mass / rr ;
01783       tmp4.std_spectral_base() ;
01784 
01785       //-------------//
01786       // X component //
01787       //-------------//
01788 
01789       //-------------------------------------
01790       //     Integration over the BH map
01791       //-------------------------------------
01792 
01793       // Surface integral <- dzpuis should be 0
01794       // --------------------------------------
01795       source_bh_surf = jbh_x(1)*ll(1) + jbh_x(2)*ll(2) + jbh_x(3)*ll(3) ;
01796 
01797       source_bh_surf.std_spectral_base() ;
01798       source_bh_surf.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01799       source_bh_surf.annule_domain(0) ;
01800       source_bh_surf.raccord(1) ;
01801 
01802       integ_bh_s_x = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
01803 
01804       // Volume integral <- dzpuis should be 4
01805       // -------------------------------------
01806       source_bh_volm = tmp4
01807         * ( jbhsr_x(1)*ll(1) + jbhsr_x(2)*ll(2) + jbhsr_x(3)*ll(3)
01808            + tmp2 * ( ll(2)*lcnf.dsdz() - ll(3)*lcnf.dsdy() ) ) ;
01809 
01810       source_bh_volm.std_spectral_base() ;
01811       source_bh_volm.inc_dzpuis(2) ;  // dzpuis : 2 -> 4
01812       source_bh_volm.annule_domain(0) ;
01813 
01814       integ_bh_v_x = source_bh_volm.integrale() ;
01815 
01816       //-------------------------------------
01817       //     Integration over the NS map
01818       //-------------------------------------
01819 
01820       // Volume integral <- dzpuis should be 4
01821       // -------------------------------------
01822       source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
01823         * (yy_ns*uu(3) - zz_ns*uu(2)) ;
01824 
01825       source_ns_volm.std_spectral_base() ;
01826       source_ns_volm.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
01827 
01828       integ_ns_v_x = source_ns_volm.integrale() ;
01829 
01830       p_angu_mom_bhns->set(0) = integ_ns_v_x
01831         + 0.5 * (integ_bh_s_x + integ_bh_v_x) / qpig ;
01832 
01833       //-------------//
01834       // Y component //
01835       //-------------//
01836 
01837       //-------------------------------------
01838       //     Integration over the BH map
01839       //-------------------------------------
01840 
01841       // Surface integral <- dzpuis should be 0
01842       // --------------------------------------
01843       Scalar jbh_y_ll(mp_bh) ;
01844       jbh_y_ll = jbh_y(1)*ll(1) + jbh_y(2)*ll(2) + jbh_y(3)*ll(3) ;
01845       jbh_y_ll.std_spectral_base() ;
01846       jbh_y_ll.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01847 
01848       source_bh_surf = jbh_y_ll - tmp1 * ll(3) ;
01849       source_bh_surf.std_spectral_base() ;
01850       source_bh_surf.annule_domain(0) ;
01851       source_bh_surf.raccord(1) ;
01852 
01853       integ_bh_s_y = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
01854 
01855       // Volume integral <- dzpuis should be 4
01856       // -------------------------------------
01857       Scalar tmp3_llz(mp_bh) ;
01858       tmp3_llz = tmp3 * ll(3) * ori_bh / rr ;
01859       tmp3_llz.std_spectral_base() ;
01860       tmp3_llz.inc_dzpuis(2) ;  // dzpuis : 0 -> 2
01861 
01862       source_bh_volm = tmp4
01863         * ( jbhsr_y(1)*ll(1) + jbhsr_y(2)*ll(2) + jbhsr_y(3)*ll(3)
01864         + tmp2 * ( ll(3)*lcnf.dsdx() - (ll(1)+ori_bh/rr)*lcnf.dsdz() )
01865         - tmp3_llz ) ;
01866 
01867       source_bh_volm.std_spectral_base() ;
01868       source_bh_volm.inc_dzpuis(2) ;  // dzpuis : 2 -> 4
01869       source_bh_volm.annule_domain(0) ;
01870 
01871       integ_bh_v_y = source_bh_volm.integrale() ;
01872 
01873       //-------------------------------------
01874       //     Integration over the NS map
01875       //-------------------------------------
01876 
01877       // Volume integral <- dzpuis should be 4
01878       // -------------------------------------
01879       source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
01880         * (zz_ns*uu(1) - xx_ns*uu(3)) ;
01881 
01882       source_ns_volm.std_spectral_base() ;
01883       source_ns_volm.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
01884 
01885       integ_ns_v_y = source_ns_volm.integrale() ;
01886 
01887       p_angu_mom_bhns->set(1) = integ_ns_v_y
01888         + 0.5 * (integ_bh_s_y + integ_bh_v_y) / qpig ;
01889 
01890       //-------------//
01891       // Z component //
01892       //-------------//
01893 
01894       //-------------------------------------
01895       //     Integration over the BH map
01896       //-------------------------------------
01897 
01898       // Surface integral <- dzpuis should be 0
01899       // --------------------------------------
01900       Scalar jbh_z_ll(mp_bh) ;
01901       jbh_z_ll = jbh_z(1)*ll(1) + jbh_z(2)*ll(2) + jbh_z(3)*ll(3) ;
01902       jbh_z_ll.std_spectral_base() ;
01903       jbh_z_ll.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
01904       source_bh_surf = jbh_z_ll + tmp1 * ll(2) ;
01905 
01906       source_bh_surf.std_spectral_base() ;
01907       source_bh_surf.annule_domain(0) ;
01908       source_bh_surf.raccord(1) ;
01909 
01910       integ_bh_s_z = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
01911 
01912       // Volume integral <- dzpuis should be 4
01913       // -------------------------------------
01914       Scalar tmp3_lly(mp_bh) ;
01915       tmp3_lly = tmp3 * ll(2) * ori_bh / rr ;
01916       tmp3_lly.std_spectral_base() ;
01917       tmp3_lly.inc_dzpuis(2) ;  // dzpuis : 0 -> 2
01918 
01919       source_bh_volm = tmp4
01920         * ( jbhsr_z(1)*ll(1) + jbhsr_z(2)*ll(2) + jbhsr_z(3)*ll(3)
01921         + tmp2 * ( (ll(1)+ori_bh/rr)*lcnf.dsdy() - ll(2)*lcnf.dsdx() )
01922         + tmp3_lly ) ;
01923 
01924       source_bh_volm.std_spectral_base() ;
01925       source_bh_volm.inc_dzpuis(2) ;  // dzpuis : 2 -> 4
01926       source_bh_volm.annule_domain(0) ;
01927 
01928       integ_bh_v_z = source_bh_volm.integrale() ;
01929 
01930       //-------------------------------------
01931       //     Integration over the NS map
01932       //-------------------------------------
01933 
01934       // Volume integral <- dzpuis should be 4
01935       // -------------------------------------
01936       source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
01937         * (xx_ns*uu(2) - yy_ns*uu(1)) ;
01938 
01939       source_ns_volm.std_spectral_base() ;
01940       source_ns_volm.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
01941 
01942       integ_ns_v_z = source_ns_volm.integrale() ;
01943 
01944       p_angu_mom_bhns->set(2) = integ_ns_v_z
01945         + 0.5 * (integ_bh_s_z + integ_bh_v_z) / qpig ;
01946 
01947     }
01948     else {  // Isotropic coordinates with the maximal slicing
01949 
01950       // Sets C/M^2 for each case of the lapse boundary condition
01951       // --------------------------------------------------------
01952       double cc ;
01953 
01954       if (hole.has_bc_lapconf_nd()) {  // Neumann boundary condition
01955           if (hole.has_bc_lapconf_fs()) {  // First condition
01956               // d(\alpha \psi)/dr = 0
01957               // ---------------------
01958               cc = 2. * (sqrt(13.) - 1.) / 3. ;
01959           }
01960           else {  // Second condition
01961               // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
01962               // -----------------------------------------
01963               cc = 4. / 3. ;
01964           }
01965       }
01966       else {  // Dirichlet boundary condition
01967           if (hole.has_bc_lapconf_fs()) {  // First condition
01968               // (\alpha \psi) = 1/2
01969               // -------------------
01970           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
01971           abort() ;
01972           }
01973           else {  // Second condition
01974               // (\alpha \psi) = 1/sqrt(2.) \psi_KS
01975               // ----------------------------------
01976           cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
01977           abort() ;
01978           //              cc = 2. * sqrt(2.) ;
01979           }
01980       }
01981 
01982       Scalar r_are(mp_bh) ;
01983       r_are = hole.r_coord(hole.has_bc_lapconf_nd(),
01984                    hole.has_bc_lapconf_fs()) ;
01985       r_are.std_spectral_base() ;
01986 
01987       const Scalar& lapconf_bh = hole.get_lapconf_tot() ;
01988       const Scalar& confo_bh = hole.get_confo_tot() ;
01989       const Sym_tensor& taij_rs = hole.get_taij_tot_rs() ;
01990 
01991       Vector jbh_rs_x(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01992       jbh_rs_x.set_etat_qcq() ;
01993       for (int i=1; i<=3; i++)
01994         jbh_rs_x.set(i) = yy_bh * taij_rs(3,i) - zz_bh * taij_rs(2,i) ;
01995 
01996       jbh_rs_x.std_spectral_base() ;
01997 
01998       Vector jbh_rs_y(mp_bh, CON, mp_bh.get_bvect_cart()) ;
01999       jbh_rs_y.set_etat_qcq() ;
02000       for (int i=1; i<=3; i++)
02001         jbh_rs_y.set(i) = zz_bh * taij_rs(1,i) - xx_bh * taij_rs(3,i) ;
02002 
02003       jbh_rs_y.std_spectral_base() ;
02004 
02005       Vector jbh_rs_z(mp_bh, CON, mp_bh.get_bvect_cart()) ;
02006       jbh_rs_z.set_etat_qcq() ;
02007       for (int i=1; i<=3; i++)
02008         jbh_rs_z.set(i) = xx_bh * taij_rs(2,i) - yy_bh * taij_rs(1,i) ;
02009 
02010       jbh_rs_z.std_spectral_base() ;
02011 
02012       Scalar tmp(mp_bh) ;
02013       tmp = - 2. * mass * mass * cc * pow(confo_bh,7.)
02014         * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
02015         / lapconf_bh / pow(r_are*rr,3.) ;
02016       tmp.std_spectral_base() ;
02017 
02018       //-------------//
02019       // X component //
02020       //-------------//
02021 
02022       //-------------------------------------
02023       //     Integration over the BH map
02024       //-------------------------------------
02025 
02026       // Surface integral <- dzpuis should be 0
02027       // --------------------------------------
02028       Scalar sou_bh_sx1(mp_bh) ;
02029       sou_bh_sx1 = jbh_rs_x(1)%ll(1) + jbh_rs_x(2)%ll(2)
02030         + jbh_rs_x(3)%ll(3) ;
02031       sou_bh_sx1.std_spectral_base() ;
02032       sou_bh_sx1.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
02033 
02034       Scalar sou_bh_sx2(mp_bh) ;
02035       sou_bh_sx2 = tmp * (yy_bh % ll(3) - zz_bh % ll(2)) ;
02036       sou_bh_sx2.std_spectral_base() ;
02037 
02038       source_bh_surf = sou_bh_sx1 + sou_bh_sx2 ;
02039 
02040       source_bh_surf.std_spectral_base() ;
02041       source_bh_surf.annule_domain(0) ;
02042       source_bh_surf.raccord(1) ;
02043 
02044       integ_bh_s_x = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
02045 
02046       //-------------------------------------
02047       //     Integration over the NS map
02048       //-------------------------------------
02049 
02050       // Volume integral <- dzpuis should be 4
02051       // -------------------------------------
02052       source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
02053         * (yy_ns*uu(3) - zz_ns*uu(2)) ;
02054 
02055       source_ns_volm.std_spectral_base() ;
02056       source_ns_volm.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
02057 
02058       integ_ns_v_x = source_ns_volm.integrale() ;
02059 
02060       p_angu_mom_bhns->set(0) = integ_ns_v_x + 0.5*integ_bh_s_x/qpig ;
02061 
02062 
02063       //-------------//
02064       // Y component //
02065       //-------------//
02066 
02067       //-------------------------------------
02068       //     Integration over the BH map
02069       //-------------------------------------
02070 
02071       // Surface integral <- dzpuis should be 0
02072       // --------------------------------------
02073       Scalar sou_bh_sy1(mp_bh) ;
02074       sou_bh_sy1 = jbh_rs_y(1)%ll(1) + jbh_rs_y(2)%ll(2)
02075         + jbh_rs_y(3)%ll(3) ;
02076       sou_bh_sy1.std_spectral_base() ;
02077       sou_bh_sy1.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
02078 
02079       Scalar sou_bh_sy2(mp_bh) ;
02080       sou_bh_sy2 = tmp * (zz_bh % ll(1) - xx_bh % ll(3)) ;
02081       sou_bh_sy2.std_spectral_base() ;
02082 
02083       source_bh_surf = sou_bh_sy1 + sou_bh_sy2 ;
02084 
02085       source_bh_surf.std_spectral_base() ;
02086       source_bh_surf.annule_domain(0) ;
02087       source_bh_surf.raccord(1) ;
02088 
02089       integ_bh_s_y = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
02090 
02091       //-------------------------------------
02092       //     Integration over the NS map
02093       //-------------------------------------
02094 
02095       // Volume integral <- dzpuis should be 4
02096       // -------------------------------------
02097       source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
02098         * (zz_ns*uu(1) - xx_ns*uu(3)) ;
02099 
02100       source_ns_volm.std_spectral_base() ;
02101       source_ns_volm.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
02102 
02103       integ_ns_v_y = source_ns_volm.integrale() ;
02104 
02105       p_angu_mom_bhns->set(1) = integ_ns_v_y + 0.5*integ_bh_s_y/qpig ;
02106 
02107 
02108       //-------------//
02109       // Z component //
02110       //-------------//
02111 
02112       //-------------------------------------
02113       //     Integration over the BH map
02114       //-------------------------------------
02115 
02116       // Surface integral <- dzpuis should be 0
02117       // --------------------------------------
02118       Scalar sou_bh_sz1(mp_bh) ;
02119       sou_bh_sz1 = jbh_rs_z(1)%ll(1) + jbh_rs_z(2)%ll(2)
02120         + jbh_rs_z(3)%ll(3) ;
02121       sou_bh_sz1.std_spectral_base() ;
02122       sou_bh_sz1.dec_dzpuis(2) ;  // dzpuis : 2 -> 0
02123 
02124       Scalar sou_bh_sz2(mp_bh) ;
02125       sou_bh_sz2 = tmp * (xx_bh % ll(2) - yy_bh % ll(1)) ;
02126       sou_bh_sz2.std_spectral_base() ;
02127 
02128       source_bh_surf = sou_bh_sz1 + sou_bh_sz2 ;
02129 
02130       source_bh_surf.std_spectral_base() ;
02131       source_bh_surf.annule_domain(0) ;
02132       source_bh_surf.raccord(1) ;
02133 
02134       integ_bh_s_z = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
02135 
02136       //-------------------------------------
02137       //     Integration over the NS map
02138       //-------------------------------------
02139 
02140       // Volume integral <- dzpuis should be 4
02141       // -------------------------------------
02142       source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
02143         * (xx_ns*uu(2) - yy_ns*uu(1)) ;
02144 
02145       source_ns_volm.std_spectral_base() ;
02146       source_ns_volm.inc_dzpuis(4) ;  // dzpuis : 0 -> 4
02147 
02148       integ_ns_v_z = source_ns_volm.integrale() ;
02149 
02150       p_angu_mom_bhns->set(2) = integ_ns_v_z + 0.5*integ_bh_s_z/qpig ;
02151 
02152       }
02153 
02154     }
02155 
02156     return *p_angu_mom_bhns ;
02157 
02158 }
02159 
02160 
02161                //-----------------------------------------------//
02162                //          Error on the virial theorem          //
02163                //-----------------------------------------------//
02164 
02165 double Bin_bhns::virial_bhns_surf() const {
02166 
02167     if (p_virial_bhns_surf == 0x0) {   // a new computation is required
02168 
02169       double virial = 1. - mass_kom_bhns_surf() / mass_adm_bhns_surf() ;
02170 
02171       p_virial_bhns_surf = new double( virial ) ;
02172 
02173     }
02174 
02175     return *p_virial_bhns_surf ;
02176 
02177 }
02178 
02179 
02180 double Bin_bhns::virial_bhns_vol() const {
02181 
02182     if (p_virial_bhns_vol == 0x0) {   // a new computation is required
02183 
02184       double virial = 1. - mass_kom_bhns_vol() / mass_adm_bhns_vol() ;
02185 
02186       p_virial_bhns_vol = new double( virial ) ;
02187 
02188     }
02189 
02190     return *p_virial_bhns_vol ;
02191 
02192 }
02193 
02194                //--------------------------------------------------//
02195                //          X coordinate of the barycenter          //
02196                //--------------------------------------------------//
02197 
02198 double Bin_bhns::xa_barycenter() const {
02199 
02200     using namespace Unites ;
02201 
02202     if (p_xa_barycenter == 0x0) {    // a new computation is required
02203 
02204         double mass_b = star.mass_b_bhns(hole.is_kerrschild(),
02205                      hole.get_mass_bh(), separ) ;
02206 
02207     const Map& mp = star.get_mp() ;
02208     Scalar xxa(mp) ;
02209     xxa = mp.xa ;   // Absolute X coordinate
02210     xxa.std_spectral_base() ;
02211 
02212     Scalar tmp(mp) ;
02213 
02214     if (hole.is_kerrschild()) {
02215 
02216         Scalar xx(mp) ;
02217         xx = mp.x ;
02218         xx.std_spectral_base() ;
02219         Scalar yy(mp) ;
02220         yy = mp.y ;
02221         yy.std_spectral_base() ;
02222         Scalar zz(mp) ;
02223         zz = mp.z ;
02224         zz.std_spectral_base() ;
02225 
02226         double yns = (star.get_mp()).get_ori_y() ;
02227 
02228         Scalar rbh(mp) ;
02229         rbh = sqrt( (xx+separ)*(xx+separ) + (yy+yns)*(yy+yns) + zz*zz ) ;
02230         rbh.std_spectral_base() ;
02231 
02232         double mass = ggrav * hole.get_mass_bh() ;
02233 
02234         tmp = sqrt( 1.+2.*mass/rbh ) ;
02235         tmp.std_spectral_base() ;
02236 
02237     }
02238     else { // Isotropic coordinates with the maximal slicing
02239 
02240         tmp = 1. ;
02241         tmp.std_spectral_base() ;
02242 
02243     }
02244 
02245     Scalar confo = star.get_confo_tot() ;
02246     confo.std_spectral_base() ;
02247 
02248     Scalar g_euler = star.get_gam_euler() ;
02249     g_euler.std_spectral_base() ;
02250 
02251     Scalar nbary = star.get_nbar() ;
02252     nbary.std_spectral_base() ;
02253 
02254     Scalar dens = pow(confo, 6.) * g_euler * nbary * tmp * xxa ;
02255     dens.std_spectral_base() ;
02256     dens.inc_dzpuis(4) ;
02257     assert(dens.get_dzpuis() == 4) ;
02258 
02259     double xa_bary = dens.integrale() / mass_b ;
02260 
02261     p_xa_barycenter = new double( xa_bary ) ;
02262 
02263     }
02264 
02265     return *p_xa_barycenter ;
02266 
02267 }
02268 
02269 
02270                //--------------------------------------------------//
02271                //          Y coordinate of the barycenter          //
02272                //--------------------------------------------------//
02273 
02274 double Bin_bhns::ya_barycenter() const {
02275 
02276     using namespace Unites ;
02277 
02278     if (p_ya_barycenter == 0x0) {    // a new computation is required
02279 
02280         double mass_b = star.mass_b_bhns(hole.is_kerrschild(),
02281                      hole.get_mass_bh(), separ) ;
02282 
02283     const Map& mp = star.get_mp() ;
02284     Scalar yya(mp) ;
02285     yya = mp.ya ;   // Absolute Y coordinate
02286     yya.std_spectral_base() ;
02287 
02288     Scalar tmp(mp) ;
02289 
02290     if (hole.is_kerrschild()) {
02291 
02292         Scalar xx(mp) ;
02293         xx = mp.x ;
02294         xx.std_spectral_base() ;
02295         Scalar yy(mp) ;
02296         yy = mp.y ;
02297         yy.std_spectral_base() ;
02298         Scalar zz(mp) ;
02299         zz = mp.z ;
02300         zz.std_spectral_base() ;
02301 
02302         double yns = (star.get_mp()).get_ori_y() ;
02303 
02304         Scalar rbh(mp) ;
02305         rbh = sqrt( (xx+separ)*(xx+separ) + (yy+yns)*(yy+yns) + zz*zz ) ;
02306         rbh.std_spectral_base() ;
02307 
02308         double mass = ggrav * hole.get_mass_bh() ;
02309 
02310         tmp = sqrt( 1.+2.*mass/rbh ) ;
02311         tmp.std_spectral_base() ;
02312 
02313     }
02314     else { // Isotropic coordinates with the maximal slicing
02315 
02316         tmp = 1. ;
02317         tmp.std_spectral_base() ;
02318 
02319     }
02320 
02321     Scalar confo = star.get_confo_tot() ;
02322     confo.std_spectral_base() ;
02323 
02324     Scalar g_euler = star.get_gam_euler() ;
02325     g_euler.std_spectral_base() ;
02326 
02327     Scalar nbary = star.get_nbar() ;
02328     nbary.std_spectral_base() ;
02329 
02330     Scalar dens = pow(confo, 6.) * g_euler * nbary * tmp * yya ;
02331     dens.std_spectral_base() ;
02332     dens.inc_dzpuis(4) ;
02333     assert(dens.get_dzpuis() == 4) ;
02334 
02335     double ya_bary = dens.integrale() / mass_b ;
02336 
02337     p_ya_barycenter = new double( ya_bary ) ;
02338 
02339     }
02340 
02341     return *p_ya_barycenter ;
02342 
02343 }
Generated on Tue Feb 7 01:35:14 2012 for LORENE by  doxygen 1.4.6