et_rot_bifluid.C

/*
 * Methods for two fluids rotating relativistic stars.
 *
 * See the file et_rot_bifluid.h for documentation
 *
 */

/*
 *   Copyright (c) 2001 Jerome Novak
 *
 *   This file is part of LORENE.
 *
 *   LORENE is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   LORENE is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with LORENE; if not, write to the Free Software
 *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


char et_rot_bifluid_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_bifluid.C,v 1.12 2004/03/25 10:29:04 j_novak Exp $" ;

/*
 * $Id: et_rot_bifluid.C,v 1.12 2004/03/25 10:29:04 j_novak Exp $
 * $Log: et_rot_bifluid.C,v $
 * Revision 1.12  2004/03/25 10:29:04  j_novak
 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
 *
 * Revision 1.11  2003/12/04 14:28:26  r_prix
 * allow for the case of "slow-rot-style" EOS inversion, in which we need to adapt
 * the inner domain to n_outer=0 instead of mu_outer=0 ...
 * (this should only be used for comparison to analytic slow-rot solution!)
 *
 * Revision 1.10  2003/11/20 14:01:26  r_prix
 * changed member names to better conform to Lorene coding standards:
 * J_euler -> j_euler, EpS_euler -> enerps_euler, Delta_car -> delta_car
 *
 * Revision 1.9  2003/11/18 18:38:11  r_prix
 * use of new member EpS_euler: matter sources in equilibrium() and global quantities
 * no longer distinguish Newtonian/relativistic, as all terms should have the right limit...
 *
 * Revision 1.8  2003/11/17 13:49:43  r_prix
 * - moved superluminal check into hydro_euler()
 * - removed some warnings
 *
 * Revision 1.7  2003/11/13 12:07:57  r_prix
 * *) changed xxx2 -> Delta_car
 * *) added (non 2-fluid specific!) members sphph_euler J_euler
 * *) more or less rewritten hydro_euler() to see if I understand it ;)
 *   - somewhat simplified and more adapted to the notation used in our notes/paper.
 *   - Main difference: u_euler is no longer used!!, the "output" instead
 *     consists of ener_euler, s_euler, sphph_euler and J_euler, which are
 *     the general 3+1 components for Tmunu.
 *
 * Revision 1.6  2003/09/17 08:27:50  j_novak
 * New methods: mass_b1() and mass_b2().
 *
 * Revision 1.5  2002/10/18 08:42:58  j_novak
 * Take into account the sign for uuu and uuu2
 *
 * Revision 1.4  2002/01/16 15:03:28  j_novak
 * *** empty log message ***
 *
 * Revision 1.3  2002/01/11 14:09:34  j_novak
 * Added newtonian version for 2-fluid stars
 *
 * Revision 1.2  2001/12/04 21:27:53  e_gourgoulhon
 *
 * All writing/reading to a binary file are now performed according to
 * the big endian convention, whatever the system is big endian or
 * small endian, thanks to the functions fwrite_be and fread_be
 *
 * Revision 1.1.1.1  2001/11/20 15:19:28  e_gourgoulhon
 * LORENE
 *
 * Revision 1.3  2001/08/28  16:04:22  novak
 * Use of new definition of relative velocity and new declarations for EOS
 *
 * Revision 1.2  2001/08/27 09:58:43  novak
 * *** empty log message ***
 *
 * Revision 1.1  2001/06/22 15:39:17  novak
 * Initial revision
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_bifluid.C,v 1.12 2004/03/25 10:29:04 j_novak Exp $
 *
 */
// Headers C
#include "math.h"

// Headers Lorene
#include "et_rot_bifluid.h"
#include "utilitaires.h"
#include "unites.h"     

                //--------------//
                // Constructors //
                //--------------//
// Standard constructor
// --------------------
Et_rot_bifluid::Et_rot_bifluid(Map& mpi, int nzet_i, bool relat, const Eos_bifluid& eos_i):
  Etoile_rot(mpi, nzet_i, relat, *eos_i.trans2Eos()), 
  eos(eos_i),
  ent2(mpi),
  nbar2(mpi),
  sphph_euler(mpi),
  j_euler(mpi, 1, CON, mp.get_bvect_cart()), 
  enerps_euler(mpi),
  uuu2(mpi),
  gam_euler2(mpi),
  delta_car(mpi)
{
  // All the matter quantities are initialized to zero :
  nbar2 = 0 ;
  ent2 = 0 ; 
  sphph_euler = 0;
  j_euler = 0;
  enerps_euler = 0;

  gam_euler.set_std_base() ; 
  
  // Initialization to a static state : 
  omega2 = 0 ; 
  uuu2 = 0 ; 
  gam_euler2 = 1 ; 
  delta_car = 0 ;
  
  // Pointers of derived quantities initialized to zero : 
  set_der_0x0() ;
  
}

// Copy constructor
// ----------------

Et_rot_bifluid::Et_rot_bifluid(const Et_rot_bifluid& et):
  Etoile_rot(et), 
  eos(et.eos),
  ent2(et.ent2),
  nbar2(et.nbar2),
  sphph_euler(et.sphph_euler),
  j_euler(et.j_euler),
  enerps_euler(et.enerps_euler),
  uuu2(et.uuu2),
  gam_euler2(et.gam_euler2),
  delta_car(et.delta_car)
{
  omega2 = et.omega2 ; 
  
  // Pointers of derived quantities initialized to zero : 
  set_der_0x0() ;
}    


// Constructor from a file 
// ------------------------
Et_rot_bifluid::Et_rot_bifluid(Map& mpi, const Eos_bifluid& eos_i, FILE* fich):
  Etoile_rot(mpi, *eos_i.trans2Eos(), fich),
  eos(eos_i),
  ent2(mpi),
  nbar2(mpi),
  sphph_euler(mpi),
  j_euler(mpi, 1, CON, mp.get_bvect_cart()), 
  enerps_euler(mpi),
  uuu2(mpi),
  gam_euler2(mpi),
  delta_car(mpi)
{

  // Etoile parameters
  // -----------------
  // omega2 is read in the file:     
  fread_be(&omega2, sizeof(double), 1, fich) ;      
  
  
  // Read of the saved fields:
  // ------------------------
  
  Tenseur ent2_file(mp, fich) ; 
  ent2 = ent2_file ; 
        
  // All other fields are initialized to zero : 
  // ----------------------------------------
  uuu2 = 0 ;
  delta_car = 0 ;
  
  // Pointers of derived quantities initialized to zero 
  // --------------------------------------------------
  set_der_0x0() ;
  
}

                //------------//
                // Destructor //
                //------------//

Et_rot_bifluid::~Et_rot_bifluid(){

  del_deriv() ; 

}

        //----------------------------------//
        // Management of derived quantities //
        //----------------------------------//

void Et_rot_bifluid::del_deriv() const {

  Etoile_rot::del_deriv() ; 
  
  if (p_ray_eq2 != 0x0) delete p_ray_eq2 ; 
  if (p_ray_eq2_pis2 != 0x0) delete p_ray_eq2_pis2 ; 
  if (p_ray_eq2_pi != 0x0) delete p_ray_eq2_pi ; 
  if (p_ray_pole2 != 0x0) delete p_ray_pole2 ; 
  if (p_l_surf2 != 0x0) delete p_l_surf2 ; 
  if (p_xi_surf2 != 0x0) delete p_xi_surf2 ;
  if (p_r_circ2 != 0x0) delete p_r_circ2 ;
  if (p_aplat2 != 0x0) delete p_aplat2 ; 
  if (p_mass_b1 != 0x0) delete p_mass_b1 ;
  if (p_mass_b2 != 0x0) delete p_mass_b2 ;
  
  set_der_0x0() ; 
}               




void Et_rot_bifluid::set_der_0x0() const {

  Etoile_rot::set_der_0x0() ;
  
  p_ray_eq2 = 0x0 ;
  p_ray_eq2_pis2 = 0x0 ; 
  p_ray_eq2_pi = 0x0 ; 
  p_ray_pole2 = 0x0 ; 
  p_l_surf2 = 0x0 ; 
  p_xi_surf2 = 0x0 ; 
  p_r_circ2 = 0x0 ;
  p_aplat2 = 0x0 ;
  p_mass_b1 = 0x0;
  p_mass_b2 = 0x0;
}               

void Et_rot_bifluid::del_hydro_euler() {

  Etoile_rot::del_hydro_euler() ; 
  sphph_euler.set_etat_nondef();
  j_euler.set_etat_nondef();
  enerps_euler.set_etat_nondef();
  uuu2.set_etat_nondef();
  gam_euler2.set_etat_nondef() ; 
  delta_car.set_etat_nondef();
  
  del_deriv() ; 
}               

// Assignment of the enthalpy field
// --------------------------------

void Et_rot_bifluid::set_enthalpies(const Cmp& ent_i, const Cmp& ent2_i) {
    
  ent = ent_i ; 
  ent2 = ent2_i ;
    
  // Update of (nbar, ener, press) :
  equation_of_state() ; 
    
  // The derived quantities are obsolete:
  del_deriv() ; 
    
}


                //--------------//
                //  Assignment  //
                //--------------//

// Assignment to another Et_rot_bifluid
// --------------------------------
void Et_rot_bifluid::operator=(const Et_rot_bifluid& et) {

    // Assignment of quantities common to all the derived classes of Etoile
    Etoile_rot::operator=(et) ;     

    assert( &(et.eos) == &eos ) ;       // Same EOS
    // Assignement of proper quantities of class Et_rot_bifluid
    omega2 = et.omega2 ; 

    ent2 = et.ent2 ;
    nbar2 = et.nbar2 ;
    sphph_euler = et.sphph_euler;
    j_euler = et.j_euler;
    enerps_euler = et.enerps_euler;
    uuu2 = et.uuu2 ;
    gam_euler2 = et.gam_euler2 ;
    delta_car = et.delta_car ;

    
    del_deriv() ;  // Deletes all derived quantities

}   

                //--------------//
                //    Outputs   //
                //--------------//

// Save in a file
// --------------
void Et_rot_bifluid::sauve(FILE* fich) const {
    
    Etoile_rot::sauve(fich) ; 
    
    fwrite_be(&omega2, sizeof(double), 1, fich) ;       
    
    ent2.sauve(fich) ; 
    
    
}

// Printing
// --------

ostream& Et_rot_bifluid::operator>>(ostream& ost) const {
    
  using namespace Unites ;

    Etoile::operator>>(ost) ; 
    
    ost << endl ; 
    ost << "Bifluid rotating star" << endl ; 
    ost << "-------------" << endl ; 
    
    double freq = omega / (2.*M_PI) ;  
    ost << "Omega1 : " << omega * f_unit 
        << " rad/s     f : " << freq * f_unit << " Hz" << endl ; 
    ost << "Rotation period 1: " << 1000. / (freq * f_unit) << " ms"
        << endl ;
       
    double freq2 = omega2 / (2.*M_PI) ;  
    ost << "Omega2 : " << omega2 * f_unit 
        << " rad/s     f : " << freq2 * f_unit << " Hz" << endl ; 
    ost << "Rotation period 2: " << 1000. / (freq2 * f_unit) << " ms"
        << endl ;
       
    double nphi_c = nphi()(0, 0, 0, 0) ;
    if ( (omega==0) && (nphi_c==0) ) {
        ost << "Central N^phi :               " << nphi_c << endl ;
    }
    else{
        ost << "Central N^phi/Omega :    " << nphi_c / omega << endl ;
    }
    if (omega2!=0) 
      ost << "Central N^phi/Omega2 :     " << nphi_c / omega2 << endl ;
    
    ost << "Error on the virial identity GRV2 : " << endl ; 
    ost << "GRV2 = " << grv2() << endl ; 
    ost << "Error on the virial identity GRV3 : " << endl ; 
    double xgrv3 = grv3(&ost) ; 
    ost << "GRV3 = " << xgrv3 << endl ; 

    ost << "Circumferential equatorial radius R_circ :     " 
    << r_circ()/km << " km" << endl ;  
    ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km" 
     << endl ;  
    ost << "Flattening r_pole/r_eq :        " << aplat() << endl ; 
    ost << "Circumferential equatorial radius R_circ2 :     " 
    << r_circ2()/km << " km" << endl ;  
    ost << "Coordinate equatorial radius r_eq2 : " << ray_eq2()/km << " km" 
     << endl ;  
    ost << "Flattening r_pole2/r_eq2 :        " << aplat2() << endl ; 

    int lsurf = nzet - 1; 
    int nt = mp.get_mg()->get_nt(lsurf) ; 
    int nr = mp.get_mg()->get_nr(lsurf) ; 
    ost << "Equatorial value of the velocity U:         " 
     << uuu()(nzet-1, 0, nt-1, nr-1) << " c" << endl ; 
    ost << "Equatorial value of the velocity U2:         " 
     << uuu2()(nzet-1, 0, nt-1, nr-1) << " c" << endl ; 
    ost << "Redshift at the equator (forward) : " << z_eqf() << endl ; 
    ost << "Redshift at the equator (backward): " << z_eqb() << endl ; 
    ost << "Redshift at the pole              : " << z_pole() << endl ; 


    ost << "Central value of log(N)        : " 
    << logn()(0, 0, 0, 0) << endl ; 

    ost << "Central value of dzeta=log(AN) : " 
    << dzeta()(0, 0, 0, 0) << endl ; 

    ost << "Central value of B^2 : " << b_car()(0,0,0,0) <<  endl ; 

    Tbl diff_a_b = diffrel( a_car(), b_car() ) ;
    ost << 
    "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
       << endl ;
    for (int l=0; l<diff_a_b.get_taille(); l++) {
        ost << diff_a_b(l) << "  " ;
    }
    ost << endl ;       

    return ost ;
    
}


void Et_rot_bifluid::partial_display(ostream& ost) const {
    
  using namespace Unites ;

    double freq = omega / (2.*M_PI) ;  
    ost << "Omega : " << omega * f_unit 
        << " rad/s     f : " << freq * f_unit << " Hz" << endl ; 
    ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
     << endl ;
    ost << endl << "Central enthalpy : " << ent()(0,0,0,0) << " c^2" << endl ; 
    ost << "Central proper baryon density : " << nbar()(0,0,0,0) 
    << " x 0.1 fm^-3" << endl ; 
    double freq2 = omega2 / (2.*M_PI) ;  
    ost << "Omega2 : " << omega2 * f_unit 
        << " rad/s     f : " << freq2 * f_unit << " Hz" << endl ; 
    ost << "Rotation period 2: " << 1000. / (freq2 * f_unit) << " ms"
     << endl ;
    ost << endl << "Central enthalpy 2: " << ent2()(0,0,0,0) << " c^2" << endl ; 
    ost << "Central proper baryon density 2: " << nbar2()(0,0,0,0) 
    << " x 0.1 fm^-3" << endl ; 
   ost << "Central proper energy density : " << ener()(0,0,0,0) 
    << " rho_nuc c^2" << endl ; 
    ost << "Central pressure : " << press()(0,0,0,0) 
    << " rho_nuc c^2" << endl ; 

    ost << "Central value of log(N)        : " 
    << logn()(0, 0, 0, 0) << endl ; 
    ost << "Central lapse N :      " << nnn()(0,0,0,0) <<  endl ; 
    ost << "Central value of dzeta=log(AN) : " 
    << dzeta()(0, 0, 0, 0) << endl ; 
    ost << "Central value of A^2 : " << a_car()(0,0,0,0) <<  endl ; 
    ost << "Central value of B^2 : " << b_car()(0,0,0,0) <<  endl ; 

    double nphi_c = nphi()(0, 0, 0, 0) ;
    if ( (omega==0) && (nphi_c==0) ) {
        ost << "Central N^phi :               " << nphi_c << endl ;
    }
    else{
        ost << "Central N^phi/Omega :         " << nphi_c / omega << endl ;
    }
        

    int lsurf = nzet - 1; 
    int nt = mp.get_mg()->get_nt(lsurf) ; 
    int nr = mp.get_mg()->get_nr(lsurf) ; 
    ost << "Equatorial value of the velocity U:         " 
     << uuu()(nzet-1, 0, nt-1, nr-1) << " c" << endl ; 

    ost << "Equatorial value of the velocity U2:         " 
     << uuu2()(nzet-1, 0, nt-1, nr-1) << " c" << endl ; 

    ost << endl 
    << "Coordinate equatorial radius r_eq =    " 
    << ray_eq()/km << " km" << endl ;  
    ost << "Flattening r_pole/r_eq :        " << aplat() << endl ; 
    ost << endl 
    << "Coordinate equatorial radius r_eq2 =    " 
    << ray_eq2()/km << " km" << endl ;  
    ost << "Flattening r_pole2/r_eq2 :        " << aplat2() << endl ; 

}

//
//   Equation of state
//

void Et_rot_bifluid::equation_of_state() {

  Cmp ent_eos = ent() ;
  Cmp ent2_eos = ent2() ;
  Tenseur rel_vel(delta_car) ;

  if (nzet > 1) {
    // Slight rescale of the enthalpy field in case of 2 domains inside the
    //  star
  
    if (nzet > 2) {
      cout << "Et_rot_bifluid::equation_of_state: not ready yet for nzet > 2 !" << endl ;       
    }
    
    double epsilon = 1.e-12 ;
  
    const Mg3d* mg = mp.get_mg() ;
    int nz = mg->get_nzone() ;
    
    Mtbl xi(mg) ;
    xi.set_etat_qcq() ;
    for (int l=0; l<nz; l++) {
      xi.t[l]->set_etat_qcq() ;
      for (int k=0; k<mg->get_np(l); k++) {
    for (int j=0; j<mg->get_nt(l); j++) {
      for (int i=0; i<mg->get_nr(l); i++) {
        xi.set(l,k,j,i) =
          mg->get_grille3d(l)->x[i] ;
      }
    }
      }
      
    }
  
    Cmp fact_ent(mp) ;
    fact_ent.allocate_all() ;
    
    fact_ent.set(0) = 1 + epsilon * xi(0) * xi(0) ;
    fact_ent.set(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
  
    for (int l=nzet; l<nz; l++) {
      fact_ent.set(l) = 1 ;
    }

    ent_eos = fact_ent * ent_eos ;
    ent2_eos = fact_ent * ent2_eos ;
    ent_eos.std_base_scal() ;
    ent2_eos.std_base_scal() ;

  } // if nzet > 1
  
  
  // Call to the EOS
  nbar.set_etat_qcq() ;
  nbar2.set_etat_qcq() ;
  ener.set_etat_qcq() ;
  press.set_etat_qcq() ;

  eos.calcule_tout(ent_eos, ent2_eos, rel_vel(), nbar.set(), nbar2.set(), 
           ener.set(), press.set(), nzet)  ; 
  
  // Set the bases for spectral expansion 
  nbar.set_std_base() ; 
  nbar2.set_std_base() ; 
  ener.set_std_base() ; 
  press.set_std_base() ; 
  
  // The derived quantities are obsolete
  del_deriv() ; 
  
}

//
// Computation of hydro quantities
//

void Et_rot_bifluid::hydro_euler(){

  const Mg3d* mg = mp.get_mg(); 
  int nz = mg->get_nzone() ; 
  int nzm1 = nz - 1 ; 

    // RP: I prefer to use the 3-vector j_euler instead of u_euler
    // for better physical "encapsulation" 
    // (i.e. --> use same form of Poisson-equations for all etoile sub-classes!)
    u_euler.set_etat_nondef(); // make sure it's not used

    // (Norm of) Euler-velocity of the first fluid
    //------------------------------
    uuu.set_etat_qcq();

    uuu.set() = bbb() * (omega - nphi() ) / nnn();
    uuu.annule(nzm1) ; 

    // gosh, we have to exclude the thing being zero here... :(
    if( uuu.get_etat() != ETATZERO )
      {
    (uuu.set()).std_base_scal() ;
    (uuu.set()).mult_rsint();
      }
    uuu.set_std_base();
    

    // (Norm of) Euler-velocity of the second fluid
    //----------------------------------------
    uuu2.set_etat_qcq();
    
    uuu2.set() = bbb() * (omega2 - nphi() ) / nnn();
    uuu2.annule(nzm1) ; 

    if( uuu2.get_etat() != ETATZERO )
      {
    (uuu2.set()).std_base_scal();
    (uuu2.set()).mult_rsint();
      }
    uuu2.set_std_base();

    // Sanity check:
    // Is one of the new velocities larger than c in the equatorial plane ?
    //----------------------------------------

    // Index of the point at phi=0, theta=pi/2 at the surface of the star:
    int l_b = nzet - 1 ; 
    int j_b = mg->get_nt(l_b) - 1 ; 

    bool superlum = false ; 
    if (relativistic) {
      for (int l=0; l<nzet; l++) {
    for (int i=0; i<mg->get_nr(l); i++) {
      
      double u1 = uuu()(l, 0, j_b, i) ; 
      double u2 = uuu2()(l, 0, j_b, i) ;
      if ((u1 >= 1.) || (u2>=1.)) {     // superluminal velocity
        superlum = true ; 
        cout << "U > c  for l, i : " << l << "  " << i 
         << "   U1 = " << u1 << endl ;
        cout << "   U2 = " << u2 << endl ;
      }
    }
      }
      if ( superlum ) {
    cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ; 
    abort() ;
      }
    }


    Tenseur uuu_car = uuu * uuu;
    Tenseur uuu2_car = uuu2 * uuu2;

    // Lorentz factors
    // --------------
    Tenseur gam_car = 1.0 / (1.0 - unsurc2 * uuu_car) ; 
    gam_euler = sqrt(gam_car) ; 
    gam_euler.set_std_base() ;  // sets the standard spectral bases for a scalar field

    Tenseur gam2_car = 1.0 /(1.0 - unsurc2 * uuu2_car) ;
    gam_euler2 = sqrt(gam2_car) ;
    gam_euler2.set_std_base() ;

    // Update of "relative velocity" squared: $\Delta^2$
    // ---------------------------
    delta_car = (uuu - uuu2)*(uuu - uuu2) / ( (1 - unsurc2* uuu*uuu2) *(1 - unsurc2* uuu*uuu2 ) ) ;
    delta_car.set_std_base() ;

    // some auxiliary EOS variables
    Tenseur Ann(ent) ;
    Ann = gam_car()*nbar()*nbar() * eos.get_Knn(nbar(), nbar2(), delta_car(), nzet);
    Tenseur Anp(ent) ;
    Anp = gam_euler()*gam_euler2()* nbar()*nbar2()* eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
    Tenseur App(ent) ;
    App = gam2_car()*nbar2()*nbar2() *eos.get_Kpp(nbar(), nbar2(), delta_car(), nzet);
  

    //  Energy density E with respect to the Eulerian observer
    //------------------------------------
    // use of ener_euler is deprecated, because it's useless in Newtonian limit!
    ener_euler.set_etat_nondef(); // make sure, it's not used

    Tenseur E_euler(mp);
    E_euler =  Ann + 2*Anp + App - press ;
    E_euler.set_std_base() ; 
    

    // S^phi_phi component of stress-tensor S^i_j
    //------------------------------------
    sphph_euler = press() + Ann()*uuu_car() + 2*Anp()*uuu()*uuu2() + App()* uuu2_car();
    sphph_euler.set_std_base();


    // Trace of the stress tensor with respect to the Eulerian observer
    //------------------------------------
    s_euler = 2*press() + sphph_euler();
    s_euler.set_std_base() ; 

    // The combination enerps_euler := (E + S_i^i) which has Newtonian limit -> rho
    if (relativistic)
      enerps_euler = E_euler + s_euler;
    else
      enerps_euler = eos.get_m1() * nbar() + eos.get_m2() * nbar2();


    // the (flat-space) angular-momentum 3-vector j_euler^i
    //-----------------------------------
    Tenseur Jph(mp);   // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
    Jph = Ann*uuu + Anp*(uuu + uuu2) + App*uuu2 ;

    j_euler.set_etat_qcq();

    j_euler.set(0) = 0;     // r tetrad component
    j_euler.set(1) = 0;     // theta tetrad component
    j_euler.set(2) = Jph()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
    j_euler.set_triad (mp.get_bvect_spher());
    j_euler.set_std_base();

    // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
    j_euler.change_triad( mp.get_bvect_cart() ) ;   // Triad = Cartesian triad

    if( (j_euler(0).get_etat() == ETATZERO)&&(j_euler(1).get_etat() == ETATZERO)&&(j_euler(2).get_etat()==ETATZERO))
      j_euler = 0;

    // The derived quantities are obsolete
    // -----------------------------------
    del_deriv() ;                

} // hydro_euler()

---