Updates to SeBa to add functionalities to AMUSE (#977)

* adding functionality to SeBa; the binary type and mass transfer type are now retreavable from within Amuse

* adding option to restart SeBa simulations from within AMUSE, through new_advanced_particle. Need to specify the relative_mass. Optional parameters are stellar type, age, core_mass, CO_coremass and radius

---------

Co-authored-by: Silvia Toonen <s.toonen@bham.ac.uk>

src/amuse/community/seba/download.py

@@ -93,7 +93,7 @@ def new_option_parser():
     result = OptionParser()
     result.add_option(
         "--seba-version",
-        default='2f6e7f37a53167b4b0dcd6c723dff7b5ee1aecba',
+        default='94e9b1d6ba1466d288a12e3afaa1eba5bca6ddca',
         dest="seba_version",
         help="SeBa commit hash to download",
         type="string"

src/amuse/community/seba/interface.cc

@@ -1,4 +1,3 @@
-#include "node.h"
 #include "single_star.h"
 #include "main_sequence.h"
 #include "worker_code.h"
@@ -189,6 +188,51 @@ local int translate_stellar_type_to_int(stellar_type stp, const real mass) {
 }
 
 
+local stellar_type translate_int_to_stellar_type(int type_number) {
+
+  switch (type_number) {
+    case 19:
+      return Brown_Dwarf;
+    case 0:
+    case 1:
+      return Main_Sequence;
+    case 2:
+      return Hertzsprung_Gap;
+    case 3:
+      return Sub_Giant;
+    case 4:
+      return Horizontal_Branch;
+    case 5:
+    case 6:
+      return Super_Giant;
+//      return Hyper_Giant;
+    case 7:
+      return Helium_Star;
+    case 8:
+    case 9:
+      return Helium_Giant;
+    case 10:
+      return Helium_Dwarf;
+    case 11:
+      return Carbon_Dwarf;
+    case 12:
+      return Oxygen_Dwarf;
+    case 13:
+      return Neutron_Star;
+    case 14:
+      return Black_Hole;
+    case 15:
+      return Disintegrated;
+    case 17:
+      return Proto_Star;
+    case 18:
+      return Planet;
+    case -1:
+      return no_of_stellar_type;
+  }
+}
+
+
 
 local int translate_binary_type_to_int(binary_type btp) {
   switch (btp) {
@@ -419,6 +463,47 @@ int new_particle(int * index_of_the_star, double mass){
     return 0;
 }
 
+
+int new_advanced_particle(int * index_of_the_star, double mass,  double relative_mass, int type_number,  double age, double core_mass, double COcore_mass,  double radius){
+
+    if (relative_mass == 0) return new_particle(index_of_the_star, mass);
+    if (age < 0) return -1;
+
+    node * new_node = new node();
+    new_node->set_label(next_seba_id);
+    new_node->set_parent(seba_root);
+    new_node->set_mass(mass);
+    mapping_from_id_to_node[next_seba_id] = new_node;
+
+    if(seba_insertion_point == 0) {
+        seba_insertion_point = new_node;
+        seba_root->set_oldest_daughter(new_node);
+    } else {
+        seba_insertion_point->set_younger_sister(new_node);
+        new_node->set_elder_sister(seba_insertion_point);
+        seba_insertion_point = new_node;
+    }
+
+    stellar_type seba_stellar_type = translate_int_to_stellar_type(type_number);    
+
+    addstar(new_node, seba_time, seba_stellar_type, seba_metallicity, 0, false);
+    new_node->get_starbase()->set_time_offset(seba_time);
+    *index_of_the_star = next_seba_id;
+
+    next_seba_id++;
+
+    new_node->get_starbase()->set_relative_age(age);
+    new_node->get_starbase()->set_core_mass(core_mass);
+    new_node->get_starbase()->set_COcore_mass(COcore_mass);
+    new_node->get_starbase()->set_effective_radius(radius);
+
+
+    return 0;
+}
+
+
+
+
 int delete_star(int index_of_the_star){
 
     map<int, nodeptr>::iterator i = mapping_from_id_to_node.find(index_of_the_star);
@@ -728,6 +813,7 @@ int get_wind_mass_loss_rate(int index_of_the_star, double * wind_mass_loss_rate)
 
 
 
+
 int evolve_one_step(int index_of_the_star){
     int error_code = 0;
     int n_steps_per_phase = 10;
@@ -966,3 +1052,19 @@ int set_rotation_period(int index_of_the_star, double value){
     return error_code;
 }
 
+int get_binary_type(int index_of_the_star, double * value){
+    int error_code = 0;
+    node * seba_node = get_seba_node_from_index(index_of_the_star, &error_code);
+    if(error_code < 0) {return error_code;}
+    *value= seba_node->get_starbase()->get_bin_type();
+    return error_code;
+}
+
+int get_mass_transfer_type(int index_of_the_star, double * value){
+    int error_code = 0;
+    node * seba_node = get_seba_node_from_index(index_of_the_star, &error_code);
+    if(error_code < 0) {return error_code;}
+    *value= seba_node->get_starbase()->get_current_mass_transfer_type();
+    return error_code;
+}
+

src/amuse/community/seba/interface.py

@@ -45,6 +45,54 @@ def evolve_star():
         function.can_handle_array = True
         return function
 
+    @legacy_function
+    def new_advanced_particle():
+        """
+        Define a new star in the code. The star will start with the given mass.
+        """
+        function = LegacyFunctionSpecification()
+        function.can_handle_array = True
+        function.addParameter(
+            'index_of_the_star', dtype='int32', direction=function.OUT,
+            description=(
+                "The new index for the star. This index can be used to refer "
+                "to this star in other functions"
+            )
+        )
+        function.addParameter(
+            'mass', dtype='float64', direction=function.IN,
+            description="The initial mass of the star")
+        function.addParameter(
+            'relative_mass', dtype='float64', direction=function.IN, default=0,
+            description="The initial relative mass of the star")
+        function.addParameter(
+            'stellar_type', dtype='int32', direction=function.IN, default=0,
+            description="The initial stellar type of the star")
+        function.addParameter(
+            'age', dtype='float64', direction=function.IN, default=0,
+            description="The initial age of the star")
+        function.addParameter(
+            'core_mass', dtype='float64', direction=function.IN, default=0, 
+            description="The initial core mass of the star")
+        function.addParameter(
+            'COcore_mass', dtype='float64', direction=function.IN, default=0,
+            description="The initial CO core mass of the star")
+        function.addParameter(
+            'radius', dtype='float64', direction=function.IN, default=0,
+            description="The initial radius of the star")
+        # function.addParameter(
+        #     'age_tag', dtype='float64', direction=function.IN,
+        #     description="Starting age of the star *to be specified exactly*")
+        function.result_type = 'int32'
+        function.result_doc = """
+        0 - OK
+            New star was loaded and the index_of_the_star parameter set.
+        -1 - ERROR
+            New star could not be created.
+        """
+        return function
+
+
     @legacy_function
     def new_binary():
         """
@@ -145,6 +193,47 @@ def get_semi_major_axis():
             A binary with the given index was not found.
         """
         return function
+
+    @legacy_function
+    def get_binary_type():
+        """
+        Retrieve the current binary type of the binary star.
+        """
+        function = LegacyFunctionSpecification()
+        function.can_handle_array = True
+        function.addParameter('index_of_the_star', dtype='int32', direction=function.IN
+            , description="The index of the star to get the value of")
+        function.addParameter('value', dtype='float64', direction=function.OUT
+            , description="The current binary type.")
+        function.result_type = 'int32'
+        function.result_doc = """
+        0 - OK
+            The value has been set.
+        -1 - ERROR
+            A binary with the given index was not found.
+        """
+        return function
+
+    @legacy_function
+    def get_mass_transfer_type():
+        """
+        Retrieve the current mass transfer type of the binary star.
+        """
+        function = LegacyFunctionSpecification()
+        function.can_handle_array = True
+        function.addParameter('index_of_the_star', dtype='int32', direction=function.IN
+            , description="The index of the star to get the value of")
+        function.addParameter('value', dtype='float64', direction=function.OUT
+            , description="The current mass transfer type.")
+        function.result_type = 'int32'
+        function.result_doc = """
+        0 - OK
+            The value has been set.
+        -1 - ERROR
+            A binary with the given index was not found.
+        """
+        return function
+
 
     @legacy_function
     def get_core_mass():
@@ -251,8 +340,6 @@ def merge_with_other_star():
         """
         return function
 
-
-
     @legacy_function
     def refresh_memory():
         """
@@ -875,6 +962,11 @@ def define_methods(self, handler):
             (units.Myr,),
             (handler.ERROR_CODE,)
         )
+        handler.add_method(
+            "new_advanced_particle",
+            (units.MSun, units.MSun, units.stellar_type, units.Myr, units.MSun, units.MSun, units.RSun),
+            (handler.INDEX, handler.ERROR_CODE)
+        )        
         handler.add_method(
             "new_binary",
             (units.RSun, handler.NO_UNIT, handler.LINK('particles'), handler.LINK('particles')),
@@ -895,6 +987,16 @@ def define_methods(self, handler):
             (handler.INDEX,),
             (units.RSun, handler.ERROR_CODE,)
         )
+        handler.add_method(
+            "get_binary_type",
+            (handler.INDEX,),
+            (handler.NO_UNIT, handler.ERROR_CODE,)
+        )
+        handler.add_method(
+            "get_mass_transfer_type",
+            (handler.INDEX,),
+            (handler.NO_UNIT, handler.ERROR_CODE,)
+        )        
         handler.add_method(
             "get_core_mass",
             (handler.INDEX,),
@@ -1075,14 +1177,16 @@ def define_parameters(self, handler):
     def define_particle_sets(self, handler):
 
         handler.define_set('particles', 'index_of_the_star')
-        handler.set_new('particles', 'new_particle')
+#        handler.set_new('particles', 'new_particle')
+        handler.set_new('particles', 'new_advanced_particle')        
         handler.set_delete('particles', 'delete_star')
-
+
+
         handler.add_getter('particles', 'get_radius', names = ('radius',))
         handler.add_getter('particles', 'get_stellar_type', names = ('stellar_type',))
         handler.add_getter('particles', 'get_mass', names = ('mass',))
         handler.add_getter('particles', 'get_core_mass', names = ('core_mass',))
-        handler.add_getter('particles', 'get_COcore_mass', names = ('CO_core_mass',))
+        handler.add_getter('particles', 'get_COcore_mass', names = ('COcore_mass',))
         handler.add_getter('particles', 'get_envelope_mass', names = ('envelope_mass',))
         handler.add_getter('particles', 'get_core_radius', names = ('core_radius',))
         handler.add_getter('particles', 'get_age', names = ('age',))
@@ -1091,6 +1195,7 @@ def define_particle_sets(self, handler):
         handler.add_getter('particles', 'get_luminosity', names = ('luminosity',))
         handler.add_getter('particles', 'get_temperature', names = ('temperature',))
         handler.add_getter('particles', 'get_natal_kick_velocity', names = ('natal_kick_x','natal_kick_y','natal_kick_z'))
+
         handler.add_getter('particles', 'get_convective_envelope_mass', names = ('convective_envelope_mass',))
         handler.add_getter('particles', 'get_convective_envelope_radius', names = ('convective_envelope_radius',))
         handler.add_getter('particles', 'get_gyration_radius', names = ('gyration_radius',))
@@ -1120,13 +1225,17 @@ def define_particle_sets(self, handler):
         handler.set_delete('binaries', 'delete_binary')
 
         handler.add_getter('binaries', 'get_semi_major_axis', names = ('semi_major_axis',))
+        handler.add_setter('binaries', 'set_semi_major_axis', names = ('semi_major_axis',))
         handler.add_getter('binaries', 'get_eccentricity', names = ('eccentricity',))
+        handler.add_setter('binaries', 'set_eccentricity', names = ('eccentricity',))
+
         handler.add_getter('binaries', 'get_mass', names = ('mass',))
         handler.add_getter('binaries', 'get_time_step', names = ('time_step',))
         handler.add_getter('binaries', 'get_age', names = ('age',))
         handler.add_getter("binaries", 'get_children_of_binary')
-        handler.add_setter('binaries', 'set_semi_major_axis', names = ('semi_major_axis',))
-        handler.add_setter('binaries', 'set_eccentricity', names = ('eccentricity',))
+        handler.add_getter('binaries', 'get_binary_type', names = ('binary_type',))
+        handler.add_getter('binaries', 'get_mass_transfer_type', names = ('mass_transfer_type',))
+
         handler.add_method('binaries', 'merge_the_binary')
 
 
@@ -1135,6 +1244,11 @@ def define_state(self, handler):
         se.StellarEvolution.define_state(self, handler)
 
         self.stopping_conditions.define_state(handler)
+
+        handler.add_method('EDIT', 'new_advanced_particle')
+        handler.add_method('UPDATE', 'new_advanced_particle')
+        handler.add_transition('RUN', 'UPDATE', 'new_advanced_particle', False)
+