Processes.h

Go to the documentation of this file.

//
// Created by spera on 13/02/19.
//
 
#ifndef SEVN_PROCESSES_H
#define SEVN_PROCESSES_H
 
 
#include <iostream>
#include <vector>
#include <string>
#include <map>
#include <random>
#include <property.h>
#include <BinaryProperty.h>
#include <Orbit.h>
#include <IO.h>
#include <lookup_and_phases.h>
 
#include <sevnlog.h>
using sevnstd::SevnLogging;
 
class Star;
class Binstar;
class Orbital_change;
 
class Process {
 
public:
    Process(){
        size++;
        VB.resize(BinaryProperty::all.size());
        VS.resize(2);
 
        VS[0].resize(Property::all.size());
        VS[1].resize(Property::all.size());
 
        set_V_to_0();
    }
 
    virtual ~ Process();
 
    static std::vector<Process*> all;
    typedef std::map<std::string,size_t> _PrintMap;
    static _PrintMap PrintMap;
 
    virtual Process * Instance(_UNUSED IO *_io){ return nullptr; }
    virtual inline std::string name(){return "Property (generic)";}
    virtual int evolve(_UNUSED Binstar *binstar){return EXIT_SUCCESS;}
    virtual int special_evolve(_UNUSED Binstar *binstar){return 0;} //0 Means system is not broken
 
 
    double  get(const size_t &id){ return VB[id]; }
    double  get_var(const size_t &starID, const size_t &propID){ return VS[starID][propID]; }
 
    std::string & get_msg() {return message;}
 
    inline void restore(){set_V_to_0();}
 
    //Special functions to  modify the VB matrix
    int modify_SemimajorDV_by_a_factor(double factor){
        VB[Semimajor::ID]*=factor;
        return EXIT_SUCCESS;
    }
    int modify_EccentricityDV_by_a_factor(double factor){
        VB[Eccentricity::ID]*=factor;
        return EXIT_SUCCESS;
    }
 
    bool is_mass_transfer_happening(){
        if (VS[0][Mass::ID]==0 && VS[1][Mass::ID]==0)
            return false;
        else
            return true;
    }
 
    virtual bool is_process_ongoing() const {return orb_change->is_process_ongoing();};
 
    inline void set_event(double code){event_code=code;}
    inline double get_event(){return event_code;}
    inline void reset_event(){event_code=-1;}
 
    inline void special_evolution_alarm_switch_on(){special_evolution_alarm=true;}
    inline void special_evolution_alarm_switch_off(){special_evolution_alarm=false;}
    inline bool is_special_evolution_alarm_on() const {return special_evolution_alarm;}
 
 
    virtual double estimate_accreted_mass(_UNUSED  double DM,  _UNUSED Star *donor, _UNUSED Star *accretor, _UNUSED Binstar *binstar) const {return 0.0;}
 
private:
    static size_t size;
    std::vector<double> VB; 
    std::vector<std::vector<double>> VS; 
    std::string message;
    double event_code=-1; //
    bool special_evolution_alarm=false; 
    //std::string log_message_core
 
protected:
 
    //Notice the pointer are set to nullptr because now the Collision process don't initialise them
    //so when the Process destructor is called it could be possible that it will try to delete a fake memory
    //address. With orb_change=nullptr this will not happen and we can avoid SEGFAULT ERRORS
    Star *donor= nullptr;
    Star *accretor= nullptr;
    SevnLogging svlog;
    Orbital_change *orb_change= nullptr;
 
    void Register(Process *_p, size_t *id, const std::string &_name){
        Process::all.push_back(_p);
        *id = Process::size - 1;
        Process::PrintMap.insert(std::pair<std::string,size_t>(_name, *id));
        svlog.debug("Binary process "+ name() + " registered" + " (Nproperties: " + utilities::n2s(all.size(),__FILE__,__LINE__) + ")");
 
    }
    std::uniform_real_distribution<double> _uniform_real;
 
 
    void set(const size_t &id, const double &value){ VB[id] = value; }
    void set_var(const size_t &starID, const size_t &propID, const double &value){ VS[starID][propID] = value; }
    void set_msg(const std::string &str) {message = str;}
 
    inline void set_V_to_0(){
 
        for(size_t i = 0; i < VB.size(); i++)
            VB[i] = 0.0;
 
        for(size_t i = 0; i < 2; i++){
            for(size_t k = 0; k < VS[i].size(); k++){
                VS[i][k] = 0.0;
            }
        }
    }
 
};
 
 
class MaccretionProcess : public  Process{
public:
 
    virtual double NS_DBmag_accretion  (_UNUSED Binstar *b, _UNUSED Star *s, _UNUSED double DM) const;
    virtual double NS_DOmegaRem_accretion  (_UNUSED Binstar *b, _UNUSED Star *s, _UNUSED double DM) const;
 
protected:
 
    virtual int accrete_mass(_UNUSED Binstar* binstar){return EXIT_SUCCESS;}
 
    virtual void handle_NS_massaccretion(_UNUSED Binstar *b, _UNUSED Star *s, _UNUSED double DM);
 
 
};
 
class RocheLobe : public MaccretionProcess{
 
public:
 
    //The RocheLobe in a eccentric orbit
    //https://iopscience.iop.org/article/10.1086/513736/pdf
    RocheLobe(_UNUSED IO *_io= nullptr, bool reg = true) {
        if (reg) {
            Register(this, &ID, name());
        }
 
        if (_io!= nullptr){
            std::string orb_change_name=Lookup::rlmap_name.at(Lookup::rlmap.at(_io->RL_mode));
            orb_change = Orbital_change_RL::Instance(orb_change_name);
            if(orb_change== nullptr) svlog.critical("RL mode " + _io->RL_mode + " not available", __FILE__, __LINE__);
            //TODO io should be passed directly to Instance
            orb_change->set_options(_io);
        }
 
    }
 
    static size_t ID;
    static RocheLobe _rochelobe;
 
    RocheLobe *Instance(IO *_io) override {
        return (new RocheLobe(_io, false));
    }
 
    inline std::string name() override { return "RocheLobe"; }
 
    int evolve(Binstar *binstar) override;
 
    int special_evolve(Binstar *binstar) override;
 
    std::string log_message_start(Binstar *binstar, double q, double qcrit, bool swap_donor_accretor=false);
 
    std::string log_message_end(Binstar *binstar, bool swap_donor_accretor=false);
 
    static std::string log_message_swallowed(Binstar *binstar, Star *swallowed, Star *other,double DM_accreted);
    std::string log_message_swallowed(Binstar *binstar, double DM_accreted);
 
protected:
 
    bool RLO_last_step=false;
 
    std::string _log_message(Binstar *binstar, bool swap_donor_accretor=false);
 
    inline int reset_DM_global(){DM_global_0=DM_global_1=0.0; return EXIT_SUCCESS;}
 
    bool check_and_handle_swapping_condition(Binstar *binstar);
 
    int reset_dMcumul_RLO_in_stars(Binstar *binstar);
 
private:
 
    double DM_global_0=0.0, DM_global_1=0.0;
    double id_donor=-1.;
    bool RLOswapped=false; 
};
 
class CommonEnvelope : public MaccretionProcess{
 
public:
    CommonEnvelope(_UNUSED IO *_io= nullptr, bool reg = true) {
        if (reg) {
            Register(this, &ID, name());
        }
 
        if (_io!= nullptr){
            std::string orb_change_name=Lookup::cemap_name.at(Lookup::cemap.at(_io->CE_mode));
            orb_change = Orbital_change_CE::Instance(orb_change_name);
            if(orb_change== nullptr) svlog.critical("CE mode " + _io->CE_mode + " not available", __FILE__, __LINE__);
        }
 
 
    }
 
    static size_t ID;
    static CommonEnvelope _commonenvelope;
 
    CommonEnvelope *Instance(IO *_io) override {
        return (new CommonEnvelope(_io, false));
    }
 
    inline std::string name() override { return "CommonEnvelope"; }
 
    int special_evolve(Binstar *binstar) override;
 
    static std::string log_mess(Binstar *binstar, Star *primary, Star *secondary);
 
    std::string log_message_swallowed(Binstar *binstar);
 
 
protected:
 
    int init(Binstar *binstar);
 
    int whoisprimary(Binstar *binstar);
 
    double final_mass_after_coalescence();
 
    std::string main_coalesce(Binstar *binstar);
    std::string coalesce(Binstar *binstar);
    std::string coalesce_with_bindingEnergy(Binstar *binstar);
    int lose_the_envelope(Binstar *binstar);
 
    inline double hurley_rm_x(double Z){
 
        double zeta  = std::log(Z/0.02);
 
        double ret = 0.30406 + 0.0805*zeta + 0.0897*zeta*zeta + 0.0878*zeta*zeta*zeta + 0.0222*zeta*zeta*zeta*zeta;
 
 
        return ret;
 
    }
 
    double hurley_final_mass(const double Ebindf,  const double toll=1e-3, const double maxit=100){
 
        //Check
        if (!isinisialised)
            svlog.critical("CE process has not been initialised",__FILE__,__LINE__,sevnstd::ce_error());
 
        double XX = 1 + hurley_rm_x(primary->get_Z());
 
 
        double MC22=secondary->getp(MHE::ID);
        double MC1=M_core_primary;
        double M1=primary->getp(Mass::ID);
        double M2=secondary->getp(Mass::ID);
        double DELY, DERI;
        int Niter = 0;
 
        double Eratio = Ebindf/Ebind_ini;
 
        double CONST = (pow(M1+M2,XX))*(M1-MC1+M2-MC22)*Eratio;
 
        double MF=std::max(MC1+MC22, (M1+M2)*pow((Eratio),(1/XX)));
        DELY = 10*toll*MF;
 
 
        while(std::abs(DELY/MF)>toll && Niter<maxit){
 
            DELY = pow(MF,XX)*(MF-MC1-MC22) - CONST;
 
            DERI = pow(MF,(XX-1))*((1+XX)*MF - XX*(MC1+MC22));
 
            MF -= DELY/DERI;
 
            Niter++;
 
 
 
        }
 
 
 
        return MF;
 
    }
 
 
 
private:
 
    //Common variables
    Star * primary;
    Star * secondary;
 
    //CE lambda, alpha formalism
    double lambda; //Concentraction factor to define the envelope binding energy
    double alpha; //Fraction of orbital energy used to expel the envelope
 
    //Orbital
    double a_fin, ecc_fin;
    //Mass and radii
    double M_env_primary, M_core_primary;
    double R_core_primary;
    double M_env_secondary, M_core_secondary;
    double R_core_secondary;
    double RL_primary_final, RL_secondary_final; //Roche Lobe of primary and secondary after CE
    //Energies
    double Ebind_ini; // Initial envelope binding energy
    double Eorb_ini, Eorb_fin; //Initial and final orbital energy of the cores
 
    //Initilisation check
    bool isinisialised=false;
};
 
class SNKicks : public Process{
 
public:
    SNKicks(_UNUSED IO *_io= nullptr, bool reg = true) {
        if (reg) {
            Register(this, &ID, name());
        }
 
        if (_io!= nullptr){
            std::string orb_change_name=Lookup::snkmap_name.at(Lookup::snkmap.at(_io->SNK_mode));
            orb_change = Orbital_change_SNKicks::Instance(orb_change_name);
            if(orb_change== nullptr) svlog.critical("RL mode " + _io->RL_mode + " not available", __FILE__, __LINE__);
        }
 
    }
 
    static size_t ID;
    static SNKicks _snkicks;
 
    SNKicks *Instance(IO *_io) override {
        return (new SNKicks(_io, false));
    }
 
    inline std::string name() override { return "SNKick"; }
 
    int special_evolve(Binstar *binstar) override;
 
    static std::string log_message(Binstar *binstar,double a_fin,double e_fin, double cos_nu, double vcom);
 
 
};
 
class GWrad : public Process{
 
public:
    GWrad(_UNUSED IO *_io= nullptr, bool reg = true) {
        if (reg) {
            Register(this, &ID, name());
        }
 
        if (_io!= nullptr){
            std::string orb_change_name=Lookup::gwmap_name.at(Lookup::gwmap.at(_io->GW_mode));
            orb_change = Orbital_change_GW::Instance(orb_change_name);
            if(orb_change== nullptr) svlog.critical("GW mode " + _io->GW_mode + " not available", __FILE__, __LINE__);
        }
 
 
    }
 
    static size_t ID;
    static GWrad _gwrad;
 
    GWrad *Instance(IO *_io) override {
        return (new GWrad(_io, false));
    }
 
    inline std::string name() override { return "GWrad"; }
 
    int evolve(Binstar *binstar) override;
 
private:
 
    bool first_time_GW = true; //This true until the GW orbital energy loss starts to be estimated (see inside evolve)
 
};
 
class Tides : public Process{
 
public:
    Tides(IO * _io= nullptr, bool reg = true) {
        if (reg) {
            Register(this, &ID, name());
        }
 
        if (_io!= nullptr){
            std::string orb_change_name=Lookup::tidesmap_name.at(Lookup::tidesmap.at(_io->tides_mode));
            orb_change = Orbital_change_Tides::Instance(orb_change_name);
            if(orb_change== nullptr) svlog.critical("Tides mode " + _io->tides_mode + " not available", __FILE__, __LINE__);
        }
    }
 
    static size_t ID;
    static Tides _tides;
 
    Tides *Instance(IO * _io) override {
        return (new Tides(_io,false));
    }
 
    inline std::string name() override { return "Tides"; }
 
    int evolve(Binstar *binstar) override;
};
 
class Mix : public  Process{
 
public:
 
    Mix(_UNUSED IO *_io= nullptr, bool reg = true) {
        if (reg) {
            Register(this, &ID, name());
        }
 
        if (_io!= nullptr){
            std::string orb_change_name=Lookup::mixmap_name.at(Lookup::mixmap.at(_io->mix_mode));
            orb_change = Orbital_change_Mix::Instance(orb_change_name);
            if(orb_change== nullptr) svlog.critical("Mix mode " + _io->mix_mode + " not available", __FILE__, __LINE__);
        }
 
 
    }
 
    static size_t ID;
    static Mix _mix;
 
    Mix *Instance(IO *_io) override {
        return (new Mix(_io, false));
    }
 
    inline std::string name() override { return "Mix"; }
 
    int special_evolve(Binstar *binstar) override;
 
    static std::string log_message(Binstar *binstar, Star *accretor, Star *donor);
 
 
};
 
 
 
 
 
 
 
 
#endif //SEVN_PROCESSES_H