eos_compose_fit.h

/*
 *  Definition of Lorene class Eos_compose_fit
 */

/*
 *   Copyright (c) 2022 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
 *
 */


#ifndef __EOS_COMPOSE_FIT_H_
#define __EOS_COMPOSE_FIT_H_

/*
 * $Id: eos_compose_fit.h,v 1.6 2023/07/12 15:39:30 j_novak Exp $
 * $Log: eos_compose_fit.h,v $
 * Revision 1.6  2023/07/12 15:39:30  j_novak
 * Reversed order of matching the EoS fit: it now starts from the crust (polytrope) and adjusts the kappa of this crust so as to have the right pressure at the highest density in the original table.
 *
 * Revision 1.4  2023/01/27 16:10:35  j_novak
 * A polytrope (class Eos_poly) is used for low and high enthalpies.
 *
 * Revision 1.3  2022/07/21 12:33:51  j_novak
 * Improved comments
 *
 * Revision 1.2  2022/07/20 12:59:43  j_novak
 * Added methods for saving into files and constructor from a file
 *
 * Revision 1.1  2022/04/15 13:39:24  j_novak
 * New class Eos_compose_fit to generate fitted EoSs from CompOSE tables.
 *
 *
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/eos_compose_fit.h,v 1.6 2023/07/12 15:39:30 j_novak Exp $
 *
 */

// Standard C++
#include "eos.h"
#include "scalar.h"

// Lorene classes
namespace Lorene {


            //------------------------------------//
            //       class Eos_compose_fit        //
            //------------------------------------//


class Eos_compose_fit : public Eos {

  // Data
  //--------
  
 protected:
  string tablename ;
        
  int n_coefs ;

  double nb_min ;

  double nb_mid ;

  double nb_max ;

  double hmin, hmax ;

  const Eos_poly* p_eos_low ;

  const Eos_poly* p_eos_high ;

  const Mg3d* mg ;
  
  const Map_af* mp ;
        
  Scalar* log_p ;
        
  Scalar* log_e ;

  Scalar* log_nb ;
        
  Scalar* log_cs2 ;

  
  // Constructors - Destructor
  // -------------------------
 public:
  
  Eos_compose_fit(const string& param_file) ;   

    
 protected:
  Eos_compose_fit(FILE* ) ;
  
  Eos_compose_fit(ifstream&) ;
  
 private:   
  Eos_compose_fit(const Eos_compose_fit& ) ;    

  friend Eos* Eos::eos_from_file(FILE* ) ;
  friend Eos* Eos::eos_from_file(ifstream& ) ;
  
 public:
  virtual ~Eos_compose_fit() ;          

  // Miscellaneous
  // -------------

 public :
  virtual bool operator==(const Eos& ) const ;
  
  virtual bool operator!=(const Eos& ) const ;
  
  virtual int identify() const ;

protected:

  void read_and_compute(ifstream&) ;

  void read_compose_data(int& nbp, Tbl*& logh, Tbl*& logp, Tbl*& loge,
             Tbl*& lognb, Tbl*& gam1) ;

  void adiabatic_index_fit(int i_min, int i_max, const Tbl& logh_read,
               const Tbl& logp_read, const Tbl& loge_read,
               const Tbl& lognb_read, const Tbl& gam1_read) ;

  double integrate_equations(double kappa, double mean_gam, const Scalar& gamma1,
                 const Scalar& log_ent, const Scalar& ent, Scalar& YYY) ;
  
  // Outputs
  // -------
 public:

  virtual void sauve(FILE* ) const ;    

  void write_lorene_table(const string&, int nlines = 200) const ; 
  
 protected:
  
  virtual ostream& operator>>(ostream &) const ;    

  // Accessors
  //-------------

  const string& get_tablename() const {return tablename ;} ;

  double get_nbmin() const { return nb_min ;} ;  
  double get_nbmax() const { return nb_max ;} ;  
  
    // Computational functions
    // -----------------------

    public:
        virtual double nbar_ent_p(double ent, const Param* par=0x0) const ;      

        virtual double ener_ent_p(double ent, const Param* par=0x0) const ;

        virtual double press_ent_p(double ent, const Param* par=0x0) const ;

        virtual double der_nbar_ent_p(double ent, const Param* par=0x0) const ; 
       
        virtual double der_ener_ent_p(double ent, const Param* par=0x0) const ; 
       
        virtual double der_press_ent_p(double ent, const Param* par=0x0) const ; 

        virtual double der_press_nbar_p(double ent, const Param* par=0x0) const ; 
        
    virtual double csound_square_ent_p(double, const Param* par=0x0) const ;
  
};



}
#endif