et_rot_mag.h

/*
 *  Definition of Lorene class Et_rot_mag
 *
 */

/*
 *   Copyright (c) 2002 Emmanuel Marcq
 *   Copyright (c) 2002,2013 Jerome Novak
 *   Copyright (c) 2013 Deberati Chatterjee
 *
 *   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
 *
 */


#ifndef __ET_ROT_MAG_H_ 
#define __ET_ROT_MAG_H_ 

/*
 * $Id: et_rot_mag.h,v 1.26 2015/06/12 12:38:24 j_novak Exp $
 * $Log: et_rot_mag.h,v $
 * Revision 1.26  2015/06/12 12:38:24  j_novak
 * Implementation of the corrected formula for the quadrupole momentum.
 *
 * Revision 1.25  2014/10/21 09:23:53  j_novak
 * Addition of global functions mass_g(), angu_mom(), grv2/3() and mom_quad().
 *
 * Revision 1.24  2014/10/13 08:52:34  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.23  2014/05/27 12:32:28  j_novak
 * Added possibility to converge to a given magnetic moment.
 *
 * Revision 1.22  2014/05/13 10:06:12  j_novak
 * Change of magnetic units, to make the Lorene unit system coherent. Magnetic field is now expressed in Lorene units. Improvement on the comments on units.
 *
 * Revision 1.21  2014/04/29 13:46:06  j_novak
 * Addition of switches 'use_B_in_eos' and 'include_magnetisation' to control the model.
 *
 * Revision 1.20  2014/04/28 12:48:12  j_novak
 * Minor modifications.
 *
 * Revision 1.19  2013/12/13 16:36:51  j_novak
 * Addition and computation of magnetisation terms in the Einstein equations.
 *
 * Revision 1.18  2013/11/25 13:52:11  j_novak
 * New class Et_magnetisation to include magnetization terms in the stress energy tensor.
 *
 * Revision 1.17  2012/08/12 17:48:36  p_cerda
 * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
 *
 * Revision 1.16  2011/10/06 14:55:36  j_novak
 * equation_of_state() is now virtual to be able to call to the magnetized
 * Eos_mag.
 *
 * Revision 1.15  2005/06/02 11:35:27  j_novak
 * Added members for sving to a file and reading from it.
 *
 * Revision 1.14  2004/03/22 13:12:41  j_novak
 * Modification of comments to use doxygen instead of doc++
 *
 * Revision 1.13  2002/10/11 11:47:35  j_novak
 * Et_rot_mag::MHD_comput is now virtual.
 * Use of standard constructor for Tenseur mtmp in Et_rot_mag::equilibrium_mag
 *
 * Revision 1.12  2002/10/09 07:54:29  j_novak
 * Et_rot_bifluid and Et_rot_mag inheritate virtually from Etoile_rot
 *
 * Revision 1.11  2002/08/02 15:07:41  j_novak
 * Member function determinant has been added to the class Metrique.
 * A better handling of spectral bases is now implemented for the class Tenseur.
 *
 * Revision 1.10  2002/06/05 15:15:59  j_novak
 * The case of non-adapted mapping is treated.
 * parmag.d and parrot.d have been merged.
 *
 * Revision 1.9  2002/06/03 13:23:16  j_novak
 * The case when the mapping is not adapted is now treated
 *
 * Revision 1.8  2002/06/03 13:00:45  e_marcq
 *
 * conduc parameter read in parmag.d
 *
 * Revision 1.7  2002/05/30 16:06:30  j_novak
 * added the right et_rot_mag.h
 *
 * Revision 1.6  2002/05/20 08:27:59  j_novak
 * *** empty log message ***
 *
 * Revision 1.5  2002/05/17 15:08:01  e_marcq
 *
 * Rotation progressive plug-in, units corrected, Q and a_j new member data
 *
 * Revision 1.4  2002/05/15 09:53:59  j_novak
 * First operational version
 *
 * Revision 1.3  2002/05/14 13:38:36  e_marcq
 *
 *
 * Unit update, new outputs
 *
 * Revision 1.1  2002/05/10 09:26:51  j_novak
 * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/et_rot_mag.h,v 1.26 2015/06/12 12:38:24 j_novak Exp $
 *
 */

// Headers Lorene

#include "etoile.h"
#include "tensor.h"

namespace Lorene {

// Local prototype (for determining the surface)
Cmp prolonge_c1(const Cmp& uu, const int nzet) ;

class Et_rot_mag : public Etoile_rot {
  
  // Data : 
  // -----
 protected:

  Cmp A_t ; 
  Cmp A_phi; 

  Cmp B_phi; 
  Cmp j_t; 
  Cmp j_phi; 

  Tenseur E_em; 

  Tenseur Jp_em; 

  Tenseur Srr_em;

  Tenseur Spp_em; 

  double Q ;
  double a_j ; 
  int conduc ; 

  // Constructors - Destructor
  // -------------------------
 public:

  Et_rot_mag(Map& mp_i, int nzet_i, bool relat, const Eos& eos_i, 
         const int cond); 


  Et_rot_mag(const Et_rot_mag& ) ;       


  Et_rot_mag(Map& mp_i, const Eos& eos_i, FILE* fich, int withbphi=0) ;     
        

  virtual ~Et_rot_mag() ;           
 

  // Memory management
  // -----------------
 protected:

  virtual void del_deriv() const ; 
    
  virtual void set_der_0x0() const ; 

  virtual void del_hydro_euler() ; 
    

  // Mutators / assignment
  // ---------------------
 public:

  void operator=(const Et_rot_mag&) ;   

  /* /\** Computes the proper baryon and energy density, as well as */
  /*  *  pressure from the enthalpy. */
  /*  *\/ */
  /* virtual void equation_of_state() ;  */
    
  // Accessors
  // ---------
 public:
  bool is_conduct() const {return (conduc==1) ;} ;
  const Cmp& get_At() const {return A_t ; } ; 
  const Cmp& get_Aphi() const {return A_phi ;} ;

  const Cmp& get_Bphi() const {return B_phi ;} ;
  const Cmp& get_jt() const {return j_t ; } ;
  const Cmp& get_jphi() const {return j_phi ;} ;
  const Tenseur& get_Eem() const {return E_em ; } ;

  const Tenseur& get_Jpem() const {return Jp_em ;} ;

  const Tenseur& get_Srrem() const {return Srr_em ; } ;

  const Tenseur& get_Sppem() const {return Spp_em ;} ;

  double get_Q() const {return Q ;} ;
  double get_a_j() const {return a_j ;} ;

  // Outputs
  // -------
 public:
  virtual void sauve(FILE* ) const ;        
   
  virtual ostream& operator>>(ostream& ) const ;    

  // Global quantities
  // -----------------
 public:
  Tenseur Elec() const ; 
  Tenseur Magn() const ; 
  virtual void MHD_comput() ; 
  virtual double mass_g() const ;       
  virtual double angu_mom() const ;  
  virtual double grv2() const ; 
  virtual double tsw() const ; 
  double MagMom() const ; 
  double Q_comput() const; 

  double Q_int() const; 

  double GyroMag() const ; 

  virtual double grv3(ostream* ost = 0x0) const ;   

  virtual double mom_quad_old() const ;

  // Computational routines
  // ----------------------
 public: 
  void magnet_comput(const int adapt_flag,
                 Cmp (*f_j)(const Cmp& x, const double),
              Param& par_poisson_At, Param& par_poisson_Avect) ;


  virtual void magnet_comput_plus(const int adapt_flag, const int initial_j,
      const Tbl an_j,
      Cmp (*f_j)(const Cmp& x, const Tbl),
      const Tbl bn_j,
      Cmp (*g_j)(const Cmp& x, const Tbl),
      Cmp (*N_j)(const Cmp& x, const Tbl),
      Param& par_poisson_At, Param& par_poisson_Avect) ;
    
  void equilibrium_mag(double ent_c, double omega0, double fact_omega, 
               int nzadapt, const Tbl& ent_limit, const Itbl& icontrol, 
               const Tbl& control, double mbar_wanted, double aexp_mass, 
               Tbl& diff, const double Q0, const double a_j0, 
               Cmp (*f_j)(const Cmp& x, const double), 
               Cmp (*M_j)(const Cmp& x,const double));


      void equilibrium_mag_plus( const Itbl& icontrol, const Tbl& control, 
               Tbl& diff,
                       const int initial_j,
                       const Tbl an_j, 
               Cmp (*f_j)(const Cmp& x, const Tbl), 
               Cmp (*M_j)(const Cmp& x,const Tbl),
                       const Tbl bn_j, 
               Cmp (*g_j)(const Cmp& x, const Tbl), 
               Cmp (*N_j)(const Cmp& x,const Tbl),
                       const double relax_mag);
};

class Et_magnetisation : public Et_rot_mag {

  // Data : 
  // -----
 protected:

  bool use_B_in_eos ; 

  bool include_magnetisation ;

  Scalar xmag ; 

  Scalar E_I; 

  Vector J_I; 

  Sym_tensor Sij_I;

  // Constructors - Destructor
  // -------------------------
 public:

  Et_magnetisation(Map& mp_i, int nzet_i, bool relat, const Eos& eos_i,
           bool include_mag=true, bool use_B = true); 


  Et_magnetisation(const Et_magnetisation& ) ;       


  Et_magnetisation(Map& mp_i, const Eos& eos_i, FILE* fich) ;       

  virtual ~Et_magnetisation() ;         
 
  // Mutators / assignment
  // ---------------------
 public:

  void operator=(const Et_magnetisation&) ; 

  virtual void equation_of_state() ; 

  // Accessors
  // ---------
 public:
  bool B_in_eos() const {return use_B_in_eos;};

  bool use_magnetisation() const {return include_magnetisation;} ;

  const Scalar& get_magnetisation() const {return xmag;} ;

  const Scalar& get_E_I() const {return E_I;} ;

  const Vector& get_J_I() const {return J_I;} ;

  const Sym_tensor& get_Sij_I() const {return Sij_I;} ; 
  
  // Outputs
  // -------
 public:
  virtual void sauve(FILE* ) const ;        
   
  virtual ostream& operator>>(ostream& ) const ;    

  // Global quantities
  // -----------------
 public:
  virtual double mass_g() const ;       
  virtual double angu_mom() const ;  
  virtual double grv2() const ; 

  virtual double grv3(ostream* ost = 0x0) const ;   

  virtual double mom_quad_old() const ;

  virtual double mom_quad_Bo() const ;

  // Computational routines
  // ----------------------
 public: 
  virtual void magnet_comput(const int adapt_flag,
                 Cmp (*f_j)(const Cmp& x, const double),
              Param& par_poisson_At, Param& par_poisson_Avect) ;

  virtual void MHD_comput() ; 

  void equilibrium_mag(double ent_c, double omega0, double fact_omega, 
               int nzadapt, const Tbl& ent_limit, const Itbl& icontrol, 
               const Tbl& control, double mbar_wanted, 
               double magmom_wanted, double aexp_mass, 
               Tbl& diff, double Q0, double a_j0, 
               Cmp (*f_j)(const Cmp& x, const double), 
               Cmp (*M_j)(const Cmp& x,const double));

};

}
#endif