binstar.h

Go to the documentation of this file.

//
// Created by mario on 04/12/18.
//
 
#ifndef SEVN_BINSTAR_H
#define SEVN_BINSTAR_H
 
#include <BinaryProperty.h>
#include <Processes.h>
#include <star.h>
#include <cmath>
#include <utilities.h>
#include <string>
#include <sstream>
#include <params.h>
 
#include <sevnlog.h>
using sevnstd::SevnLogging;
 
#include <lookup_and_phases.h>
#include "Orbit.h"
 
using namespace Lookup;
 
 
class Binstar {
 
public:
    bool broken; 
    bool onesurvived;  
    bool empty;
    bool break_at_remnant;  
    bool break_at_broken; 
    bool repeatstep;  
    bool mix; 
    bool comenv; 
    bool is_swallowed[2]; 
    double last_Timestep; 
    //TODO the print options should be removed from the classes stars and binaries and handled by another class only devoted to that
    bool print_all_steps;  
    bool print_per_phase;  
    bool print_only_end;  
    bool print_events; 
    bool print_rlo; 
    //Check disable, some process can force the e or a to jump to some values
    //in that case we have to disable the check for adaptive time step.
    bool disable_e_check; 
    bool disable_a_check; 
    bool disable_DM_check; 
    bool disable_OmegaRem_NS_check; 
    bool disable_stellar_rotation_check; 
    bool force_tiny_dt; 
    //Other variables
    double CE_E_tomatch = utilities::NULL_DOUBLE; 
    //Guard
    bool evolution_step_completed;
 
 
    Binstar(){}
 
    //Proper constructors
    Binstar(IO *_io, std::vector<std::string> &params, size_t &_ID, unsigned long _rseed=0){
 
 
        //GI 03/02/2023: Based on a issue regarding data leak after failed initialisation, I implement
        //this solution to robustly clean the heap allocated variables after a failed object initilisation.
        //Indeed, in such cases, the object is never build, so the desctructor is not available.
        //So, I have implemented a private method default_destructor that is called from the destructor to clean the heap
        //but can be called also here. So we use a try-catch-rethrow schema to check if the initilisation
        //was successful, otherwise we clean the heap and then rethrow the error.
        //We use this schema also for the other two constructors.
 
        //Try the initialisation
        try {
            /***********************************
            * Some variable assignment
            ************************************/
            io = _io;
            ID = _ID;
            break_at_remnant = false;
            broken = onesurvived = empty = false;
            break_at_remnant = break_at_broken = false;
            init_params = params;
 
            print_all_steps = false;
            print_per_phase = false;
            print_events = false;
            print_rlo = false;
 
            last_Timestep = 0.0;
            repeatstep = false;
            mix = false;
            comenv = false;
            evolution_step_completed = false;
            is_swallowed[0] = is_swallowed[1] = false;
            disable_e_check = disable_a_check = disable_DM_check = disable_OmegaRem_NS_check = disable_stellar_rotation_check = force_tiny_dt = false;
            rseed = _rseed;
            /************************************/
 
 
            state.resize(io->printcolumns_binary.size());
 
            for (size_t i = 0; i < BinaryProperty::all.size(); i++) {
                property.push_back(BinaryProperty::all[i]->Instance()); //create instances for all properties
            }
 
            for (size_t i = 0; i < Process::all.size(); i++) {
                //std:cout<<"Procss "<< Process::all[i]->name() << " " << i <<std::endl;
                process.push_back(Process::all[i]->Instance(_io)); //create instances for all properties
            }
 
 
 
            /***********************************
            * Init Stars and binary params
            ************************************/
            init(init_params);
            /************************************/
        }
        //Catch any kind of error, dangerous in general, but we use it just to clean the heap and then
        //the error is rethrown and handled somewhere else.
        catch (...){ //Catch all errors
            default_destructor();  //Clean the heap associated to class members
            throw; //Rethrown the exact same  error
        }
 
 
    }
 
    inline double getp(const size_t &id) const {return property[id]->get(this);}
    inline double getp_0(const size_t &id) const {return property[id]->get_0(this);}
    inline  std::string get_name() const {return name;}
    inline size_t get_ID() const {return ID;}
    inline Star *getstar(const size_t &id)  {return star[id];}
    inline Process *getprocess(const size_t &id) {return process[id];}
    inline const std::vector<Process*>& getprocesses() const {return process;}
    inline const std::vector<double> & getstate() {return state;}
    inline double get_last_Timestep() const {return last_Timestep;}
 
    inline unsigned  int get_id_more_evolved_star() const{
        //If one of the star is empty return the other one
        if (star[0]->isempty and !star[1]->isempty)
            return 1;
        else if (star[1]->isempty)
            return 0;
        //If both are remnant the older is the more compact following the order HeWD,COWD,ONEWD,NS,BH
        else if (star[0]->isremnant and star[1]->isremnant and star[0]->getp(RemnantType::ID)!=star[1]->getp(RemnantType::ID))
            return star[0]->getp(RemnantType::ID)>star[1]->getp(RemnantType::ID) ? 0 : 1;
        //If both are remnant of the same type, the more evolved is the one with the largest mass
        else if (star[0]->isremnant and star[1]->isremnant)
            return star[0]->getp(Mass::ID)>star[1]->getp(Mass::ID) ? 0 : 1;
        //In the stars have the same phase (not remnant handled before), return the more evolved is the one with the larger plife
        else if (star[0]->getp(Phase::ID)==star[1]->getp(Phase::ID))
            return star[0]->plife()>=star[1]->plife() ? 0 : 1;
        //Finally, if the stars are not empty, not remnant and don't have the same phase
        //the more evolved is the one with the largest phase.
        else
            return star[0]->getp(Phase::ID)>=star[1]->getp(Phase::ID) ? 0 : 1;
    }
 
    inline Star *get_more_evolved_star() {
 
        unsigned int _id=get_id_more_evolved_star();
 
        return getstar(_id);
    }
 
    double Radx(size_t starID);
 
 
 
    inline double get_svpar_num(std::string name) { return io->svpar.get_num(name);};
    inline std::string get_svpar_str(std::string name) { return io->svpar.get_str(name);};
    inline bool get_svpar_bool(std::string name) { return io->svpar.get_bool(name);};
    inline unsigned long get_rseed(){
        if (rseed==0)
            svlog.critical("rseed has not been initialised",__FILE__,__LINE__,sevnstd::bse_error());
        return rseed;
    }
 
    //TODO Move to the src file
    virtual void evolve() {
 
        evolution_step_completed=false;
 
        svlog.debug(std::to_string(property[BTimestep::ID]->get()));
        svlog.debug(std::to_string(property[BTimestep::ID]->get_0()));
        //utilities::wait(__PRETTY_FUNCTION__);
 
 
        //double twait=4.5668373678e+01;
        double twait=1E30;
        if (getp(BWorldtime::ID)>twait){
            svlog.debug("S0 T= " + utilities::n2s(star[0]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[0]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[0]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("S1 T= " + utilities::n2s(star[1]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[1]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[1]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("B T= " + utilities::n2s(getp(BTimestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(getp_0(BTimestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(getp(BWorldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            utilities::wait("Before synch_dt");
        }
 
 
 
        /*************************************
        * Synchronise proposed timesteps
        *************************************/
        synchronise_dt_star();
 
 
 
        if (getp(BWorldtime::ID)>twait){
            svlog.debug("S0 T= " + utilities::n2s(star[0]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[0]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[0]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("S1 T= " + utilities::n2s(star[1]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[1]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[1]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("B T= " + utilities::n2s(getp(BTimestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(getp_0(BTimestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(getp(BWorldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            utilities::wait("After synch_dt");
        }
 
        /*************************************
        * Single stellar evolution
        *************************************/
        //At this point timesteps are in sync
        star[0]->evolve();
        //timestep of star 0 may have changed at this point
        star[1]->evolve();
        //timestep of star 1 may have changed at this point
 
 
 
        //std::cout<<" 0: "<<star[0]->getp_0(Timestep::ID)<<" "<<star[0]->getp(Timestep::ID)<<std::endl;
        //std::cout<<" 1: "<<star[1]->getp_0(Timestep::ID)<<" "<<star[1]->getp(Timestep::ID)<<std::endl;
 
        if (getp(BWorldtime::ID)>twait){
            svlog.debug("S0 T= " + utilities::n2s(star[0]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[0]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[0]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("S1 T= " + utilities::n2s(star[1]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[1]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[1]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("B T= " + utilities::n2s(getp(BTimestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(getp_0(BTimestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(getp(BWorldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.pdebug("S0 type",star[0]->getp(Phase::ID),__FILE__,__LINE__);
            utilities::wait("After evolve star");
        }
 
        //Check if the the two stars have been evolved with the same timestep.
        //If not repeat the evolution of the star with the largest used timestep (V0).
        //Plus update the new timestep of the binary.
        check_and_sync_sse();
 
 
 
        if (getp(BWorldtime::ID)>twait){
            svlog.debug("S0 T= " + utilities::n2s(star[0]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[0]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[0]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("S1 T= " + utilities::n2s(star[1]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[1]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[1]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("B T= " + utilities::n2s(getp(BTimestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(getp_0(BTimestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(getp(BWorldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            utilities::wait("After check and synch");
        }
 
 
        /*************************************
        * Binary stellar evolution
        *************************************/
        //Here the timestep of the binary is equal to the one just done by the single stars (s->getp_0(Timestep::ID)).
        //So it is equal to one defined in synchronise_dt_star if the stars have not repeated the evoluion in evolve.
        //In that case it is equal to the actual timestep of the single stellar evolution
        evolve_binary();
        //NOT NEEDED ANYMORE
        //binary_evolution_status = evolve_binary();
        //Now if the the Binary properties have not been uptated, update in any case the BTimestep to be syncrhonised with stars Timestepgit
        //if (binary_evolution_status==utilities::BIN_EV_NOT_DONE){
         //   std::cout<<" WE "<<std::endl;
            //property[BTimestep::ID]->evolve(this);
        //    std::cout<<" WE2 "<<std::endl;
        //}
 
        if (getp(BWorldtime::ID)>twait){
            svlog.debug("S0 T= " + utilities::n2s(star[0]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[0]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[0]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("S1 T= " + utilities::n2s(star[1]->getp(Timestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(star[1]->getp_0(Timestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(star[1]->getp(Worldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            svlog.debug("B T= " + utilities::n2s(getp(BTimestep::ID),__FILE__,__LINE__) + " T0=" + utilities::n2s(getp_0(BTimestep::ID),__FILE__,__LINE__) + " WT=" + utilities::n2s(getp(BWorldtime::ID),__FILE__,__LINE__),__FILE__,__LINE__);
            utilities::wait("After evolv binary");
        }
 
 
 
        //Check if the binary evolution has re-evolved the system with a smaller timestep
        //In that case just re-evolve the stars with this timestep and again the binary.
        //NB: This check assumes that the stars are already synchronised (check_and_sync_sse) and that
        //the binary has been evolved with the actual time step of the stellar evolution. Therefore,
        // the new proposed binary dt can be only smaller than the stellar dt. If it is larger
        check_and_sync_bse();
        //property[BTimestep::ID]->resynch(1e30);
        //At this point stars and binary processes have been evolved with the same timestep
 
        evolution_step_done();
 
 
        /****************************************************************
        * FINAL STUFF. All the updates that need to be done after a complete evolution step (all the repetitions already done) go here
        * Last checks, change tracks and check mass limits for remnants
        * Plus,  evolve the derived property so that they are updated considering processes and change of tracks
        *****************************************************************/
        //TODO The points below can be completely skipped if the binary is broken
        Star * this_star = star[0];
        Star * other_star = star[1];
 
        for (int i=0; i<2; i++){
            //First of all check if it has to explode
            //TODO this is a temporary solution, we have to think about how to handle SNIa
            if (this_star->is_exploding_as_SNI()){
                    //TODO Can we remove the difference between explode_as_SNI(binstar *b) and explode_as_SNI(), i.e. between Star::set_empty and Star::set_empty_in_bse
                    this_star->explode_as_SNI(this); //Use explode_as_SNI(Binstar* b) so that timestep_0 is reset to the current value
                    set_broken();
                    set_onesurvived();
                    //TODO, Here we are using a arbitrary process just to set the event, we should have a proper way to set the event without calling the preoceese
                    process[RocheLobe::ID]->set_event((double)Lookup::EventsList::SNIa);
            }
            check_and_set_QHE(this_star); //Check if we have to flag QHE evolution
            this_star->find_new_track(); //Check if we have to change tracks
            this_star->check_remnant_mass_limits(); //Check if we have to performe a Crem transition (e.g. from NS to BH)
            //If the Mass is extremely small set it to empty.
            if (this_star->getp(Mass::ID)<1e-3){this_star->set_empty();}
 
            //Update derived properties (basic properties could have been changes during binary processes)
            //TODO Now this check is made also inside binary_evolution. I have to check  if we can remove this call from here
            update_derived_properties_star(this_star);
 
            //update Maximum CO and minimum HE (only after the CO core has been formed)
            if (get_svpar_bool("ev_set_maxCO") and std::isnan(this_star->get_MCO_max()) and this_star->getp(Phase::ID) >= TerminalCoreHeBurning and !this_star->amiremnant()){
                this_star->set_MCO_max();
            }
            if (get_svpar_bool("ev_set_minHE") and std::isnan(this_star->get_MHE_min()) and this_star->getp(Phase::ID) >= TerminalCoreHeBurning and !this_star->amiremnant()){
                this_star->set_MHE_min();
            }
            utilities::swap_stars(this_star,other_star);
        }
 
        //Very last thing Binary event special evolve
        property[BEvent::ID]->special_evolve(this);
 
    }
 
    //TODO tO save space maybe we can flush the allstates data to the output file after a given number of stored elements.
    //TODO use a binary property NextOutput
    inline void recordstate_w_timeupdate() {
 
        evolution_guard(__FILE__,__LINE__);
 
        //TODO check that they are in sync
        star[0]->recordstate_w_timeupdate();
        star[1]->recordstate_w_timeupdate();
 
        recordstate_binary();
 
    }
    inline void recordstate() {
 
        //TODO check that they are in sync
        star[0]->recordstate();
        star[1]->recordstate();
 
        recordstate_binary();
    }
    inline void recordstate_binary(){
 
        for(size_t i = 0; i < io->printcolumns_binary.size(); i++){
            //Concerning the Timestep, put the used one not the proposed one
            if ((size_t)io->printIDs_binary[i]==BTimestep::ID)
                state[i] = getp_0((size_t)io->printIDs_binary[i]);
                //state[i] = star[0]->getp_0(Timestep::ID);
            else
                state[i] = getp((size_t)io->printIDs_binary[i]);
 
        }
 
 
        combinedstate = star[0]->getstate();
        combinedstate.insert(combinedstate.end(), star[1]->getstate().begin(), star[1]->getstate().end());
        combinedstate.insert(combinedstate.end(), state.begin(), state.end());
 
        allstates.push_back(combinedstate);
 
        //TODO The single recordstate of the two stars are now stored in combinedstate. Can we reset and free the memory of star.allstates?
        for (auto & cc: combinedstate)
            svlog.debug("Combined state: "+utilities::n2s(cc,__FILE__,__LINE__));
    }
 
    inline void set_broken(){
        broken = true;
        check_and_set_broken();
        //std::string w = utilities::log_print("BROKEN", this);
        //print_to_log(w);
    }
    inline void set_empty(){onesurvived=false; empty=true;}
    inline void set_onesurvived(){
        //set_onesurvived is triggered when a star in the system becomes empty:
        //if oneservived is already true, the binasry system is now empty
        //if empty the binary is empty and onesurvived has to be false.
        if (empty)
            onesurvived = false;
        //if true and not empty, not the binary system becomes empty and onesurvived is reset to false
        if (onesurvived){
            onesurvived = false;
            empty = true;
        }
        //If binary is not empty and onesurvived was false set onesurvived to true
        else
            onesurvived= true;
    }
    inline bool breaktrigger() const {
 
 
        //if (star[0]->isempty and star[1]->isempty) return true; //the system is empty
        if ( (broken or onesurvived) and break_at_broken) return true; //system is broken or just one object in the system
        else if (star[0]->isremnant and star[1]->isremnant and break_at_remnant) return true;
        else if ( (!star[0]->isremnant or !star[1]->isremnant)  and break_at_remnant ){
            return false;
        }
        else if(!break_at_remnant){
            return getp(BWorldtime::ID)>=get_tf() or  utilities::areEqual(getp(BWorldtime::ID),get_tf());
        }
        else{
            SevnLogging svlog_tmp;
            svlog_tmp.critical("Star breaktrigger check failed (Maybe Worldtimes are not synchronised?)",__FILE__,__LINE__);
        }
 
 
        return false;
 
    }
 
    inline bool is_process_ongoing(const size_t &id) const {return process[id]->is_process_ongoing();}
 
 
    inline bool processes_alarm(const size_t &id) const {return process[id]->is_special_evolution_alarm_on();}
 
    bool process_alarm_any() const;
 
    inline void reset_all_processes_alarms(){for (auto& proc : process) proc->special_evolution_alarm_switch_off();};
 
 
    inline void reset_events(){
        for (auto& proc : process)
            proc->reset_event();
    }
    inline double get_event_from_processes(){
 
        bool RLOB_triggered=false;
        bool Merger_triggered=false;
        bool Collision_triggered=false;
        bool CE_triggered=false;
 
        double event=-1.0;
        for (const auto& proc : process){
            double this_event=proc->get_event();
            if (this_event==EventsList::RLOBegin) RLOB_triggered=true;
            else  if (this_event==EventsList::Merger) Merger_triggered=true;
            else if (this_event==EventsList::Collision) Collision_triggered=true;
            else if (this_event==EventsList::CE) CE_triggered=true;
            else if (this_event==EventsList::CE_Merger){CE_triggered=true;Merger_triggered=true;}
            if (this_event>event){
                event=this_event;
            }
        }
 
        //Check composite events
        if (Collision_triggered and CE_triggered and Merger_triggered) event=EventsList::Collision_CE_Merger;
        else if (Collision_triggered and CE_triggered) event=EventsList::Collision_CE;
        else if (Collision_triggered and Merger_triggered) event=EventsList::Collision_Merger;
        else if (RLOB_triggered and CE_triggered and Merger_triggered) event=EventsList::RLOB_CE_Merger;
        else if (RLOB_triggered and CE_triggered) event=EventsList::RLOB_CE;
        else if (RLOB_triggered and Merger_triggered) event=EventsList::RLOB_Merger;
        else if (CE_triggered and Merger_triggered) event=EventsList::CE_Merger;
 
        return event;
    }
 
 
    //TODO, outoput stuff should be removed from the star and binary classes and handeld by a separate class
    bool isoutputtime();
 
    inline bool printall() {return print_all_steps;}
    inline bool notprint() {return print_only_end;}
    void print(){
        io->print_log();
        io->print_evolved_summary(name, rseed, ID);
        io->print_formatted_output(allstates, name, rseed, ID,true);
    }
    void print_failed(bool include_in_output=false){
        io->print_log();
        io->print_failed_summary(name, rseed, ID);
        if (include_in_output)
            io->print_formatted_output(allstates, name, rseed, ID,true);
    }
 
    void  inline print_to_log(std::string& message){
        io->log_put(message);
    }
 
    void init(const std::vector<std::string> &params);
 
    inline std::vector<double>  get_zams() {
 
        std::vector<double> zams{star[0]->get_zams(), star[1]->get_zams()};
        return zams;
    };
 
    inline std::vector<double>  get_Z() {
 
        std::vector<double> Z{star[0]->get_Z(), star[1]->get_Z()};
        return Z;
    };
 
    inline void set_tf(const double a, const char* file, const int line) {
 
        if(std::isnan(a) || std::isinf(a))
            svlog.critical("Tf set to INF or NAN", file, line);
        else if(a < 0.0)
            svlog.critical("Tf cannot be negative", file, line);
        else
            tf = a;
    }
 
    inline void set_dtout(const double a, const char* file, const int line) {
 
        if(std::isnan(a) || std::isinf(a))
            svlog.critical("dtout set to INF or NAN", file, line);
        else if(a <0.0)
            svlog.critical("dtout cannot be negative", file, line);
        else
            dtout = a;
    }
 
 
    inline void evolution_step_done(){
        last_Timestep = getp_0(BTimestep::ID); //V0 in BTimestep contains the Timestep of the step currently done, V is the time predicted for the new Timestep
        evolution_step_completed=true;
    }
    inline void evolution_guard(const char *file_input = nullptr, int line_input = -1){
        if (!evolution_step_completed)
            svlog.critical("The function recordstate_w_timeupdate can be called only after an evolution step",file_input,line_input,sevnstd::bse_error());
    }
    inline double  get_tf() const { return tf;};
 
    //TODO We have to take into account that the two stars can have different phases times. For now we just disable the print_per_phase option
    inline double get_dtout()  {
        //if(print_per_phase)
        //    return dtout_phase[(size_t)properties[Phase::ID]->get()];   //TODO it can be or fixed or print all time steps or print all timesteps diveded by 2 (print every 2), or divided by 3 (print every 3) and so on...
        //else
        //    return dtout;
        return dtout;
    }
 
    inline std::string get_id_name() const {
        std::string mss = "ID: "+utilities::n2s(ID,__FILE__,__LINE__) + " Name: "+utilities::n2s(name,__FILE__,__LINE__);
        return mss;
    }
 
    utilities::sn_explosion check_accretion_on_compact(size_t donorID, size_t accretorID, double DMaccreted);
 
    void check_accretion_on_compact(){
 
        Star *donor    = star[0];
        Star *accretor = star[1];
        double DM_accretor;
 
        for (size_t i=0; i<2; i++){
 
            //Consider all the process, NB DM can be negative
            //Cycle over all the processes to get the total DM variation.
            DM_accretor=0;
            for (auto &proc : process){
                DM_accretor += proc->get_var(accretor->get_ID(), Mass::ID);
            }
 
            //Check only if it is really accreting
            if (DM_accretor>0)
                check_accretion_on_compact(donor->get_ID(), accretor->get_ID(), DM_accretor);
 
 
            utilities::swap_stars(donor, accretor);
        }
    }
 
    bool check_and_set_broken(){
        //If broken or onesurvived or empty set properties to broken  and return true
        if (broken or onesurvived or empty) {
            for (auto &prop: property)
                prop->set_broken(this);
            return true;
        }
 
        return false;
    }
 
 
    void check_and_set_bavera_xspin();
    void check_if_applied_bavera_xspin();
    double SetXspinBavera(const double, const double) const; // Computes the Xspin following Bavera et al. (2021). Used only for the second BH spin in the case of a binary undergoing WR+BH --> BH+BH.
 
 
 
protected:
 
    double tf; 
    double dtout; 
    std::vector<double> state;  
    std::vector<double> combinedstate; 
    std::vector<std::vector<double>> allstates; 
 
 
    void synchronise_dt_star();
 
    void check_and_sync_sse();
 
 
    void check_and_sync_bse();
 
 
    void check_nakedHe_or_nakedCO_after_binary_evolution();
 
    void check_AngMomSpin_after_binary_evolution();
 
    void limit_and_correct_mass_transfer_for_donor_from_binary(Star *donor, Star *accretor,
                                                               std::size_t mass_limiting_property_id);
 
    bool check_and_set_QHE(Star *accretor);
 
    void update_derived_properties_star(Star *s);
 
 
    inline void reset_evolution_flags(){mix=comenv=is_swallowed[0]=is_swallowed[1]=false;}
 
 
protected:
    SevnLogging svlog;  
private:
 
    vector<Process*> process;  
    vector<BinaryProperty*> property; 
    Star *star[2]={nullptr, nullptr}; 
    IO *io; 
    size_t ID; 
    std::vector<std::string> init_params; 
    std::string name;  
    unsigned long rseed; 
 
    utilities::bse_evolution evolve_binary();
 
    inline void resynch(const double &dt){
 
        for (auto &prop : property) prop->restore();
        property[BTimestep::ID]->resynch(dt);
 
    }
 
    void default_destructor();
 
    void call_stars_constructor(std::vector<std::string> &params_star1, std::vector<std::string> &params_star2);
 
    void init_stars(const std::vector<std::string> &params);
 
    void init_stars_legacy(const std::vector<std::string> &params);
 
    inline void set_initial_semimajor(double a_ini){
        if (a_ini<=0){
            svlog.critical("Error on binary initialisation (" + get_id_name() +"): Semimajor is " + utilities::n2s(a_ini, __FILE__, __LINE__) + ": out of range (Semimajor>0)", __FILE__, __LINE__,sevnstd::sevnio_error());
        }
        property[Semimajor::ID]->init(a_ini);
    }
 
    inline void set_initial_semimajor(std::string a_ini){
 
        double a_inid;
        try{
            a_inid = utilities::s2n<double>(a_ini, __FILE__, __LINE__);
        }
        catch(sevnstd::sevnio_error& e){
            svlog.critical("Impossible to transform the given initial Semimajor axis to a number ("+
            get_id_name()+"): its value is "+a_ini,__FILE__,__LINE__,sevnstd::sevnio_error());
        }
        set_initial_semimajor(a_inid);
    }
 
 
    inline void set_initial_eccentricity(double ecc_ini){
        if (ecc_ini<0 or ecc_ini>=1){
            svlog.critical("Error on binary initialisation (" + get_id_name() +"): Eccentricity is " + utilities::n2s(ecc_ini, __FILE__, __LINE__) + ": out of range (0<=Eccentricity<1)", __FILE__, __LINE__,sevnstd::sevnio_error());
        }
        property[Eccentricity::ID]->init(ecc_ini);
    }
 
    inline void set_initial_eccentricity(std::string ecc_ini){
 
        double ecc_inid;
        try{
            ecc_inid = utilities::s2n<double>(ecc_ini, __FILE__, __LINE__);
        }
        catch(sevnstd::sevnio_error& e){
            svlog.critical("Impossible to transform the given initial Eccentricity  to a number ("+
                           get_id_name()+"): its value is "+ecc_ini,__FILE__,__LINE__,sevnstd::sevnio_error());
        }
        set_initial_eccentricity(ecc_inid);
    }
 
 
    inline void init_binary_properties(const std::vector<std::string> &params){
 
        int inchoice=inputmapbin.at(io->binput_mode).first;
 
        std::string a,ecc;
 
        if (inchoice==InputBinaryOption::_new){
            a=params[10];
            ecc=params[11];
        }
        else if (inchoice==InputBinaryOption::_legacy){
            a=params[6];
            ecc=params[7];
        }
 
        //Init semimajor
        set_initial_semimajor(a);
 
        //Init eccentricity
        set_initial_eccentricity(ecc);
 
    }
 
    void  init_other_params();
 
    void set_break_at_broken(std::string &tf){
        if (tf=="broken"){
            break_at_broken=true;
            tf="14000"; //Set time to Hubble time
        }
        else{
            break_at_broken=false;
        }
    }
 
    inline void set_rseed(const unsigned long a, const char* file, const int line){
 
        if(std::isnan(a) || std::isinf(a))
            svlog.critical("Rseed set to INF or NAN", file, line);
        else if (a<=0)
            svlog.critical("Rseed cannot be negative or 0", file, line);
        else
            rseed = a;
    }
 
    void inline check_init_param(const std::vector<std::string> &params){
 
 
        int par_expected_size = inputmapbin.at(io->binput_mode).second; //Number of expected parameters given the binput_mode stored in io.
        //If rseed is provided we have an extra parameter that is the seed.
        if (io->rseed_provided())
            par_expected_size+=1;
        if ((int)params.size()!=par_expected_size){
            std::string par_exp=utilities::n2s(par_expected_size, __FILE__, __LINE__);
            std::string inparam_size=utilities::n2s(params.size(), __FILE__, __LINE__);
            std::string this_ID=utilities::n2s(ID, __FILE__, __LINE__);
            svlog.critical("Number of Binary params needs to be " +par_exp+", it is instead " + inparam_size + " (ID_bin: " + this_ID +"). Maybe you are using inputs containing the ran"
                                                                                                                                       "dom seed, but rseed parameter is false." , __FILE__, __LINE__);
        }
    }
 
public:
 
    ~Binstar();
 
};
 
 
 
#endif //SEVN_BINSTAR_H