webassembly & command line options

README.md

@@ -109,7 +109,7 @@ cd run1; ../SeBa -R -n 100000 -N 0      -s 1001 > log.txt 2>&1; cd ..
 cd run2; ../SeBa -R -n 100000 -N 100000 -s 1002 > log.txt 2>&1; cd ..
 cat run1/SeBa.data run2/SeBa.data > all.data
 ```
-The binaries in `all.data` then have the IDs 0 to 199999, each exactly once.
+The binaries in `all.data` then have the IDs 0 to 199999, each exactly once. Instead of separate directories, you can also give each part its own output file with `-O`.
 
 ## Options
 
@@ -124,8 +124,10 @@ The binaries in `all.data` then have the IDs 0 to 199999, each exactly once.
 | `-T` or `-t` | end time of the simulation in Myr, the same for all binaries | 13500 | > 0 |
 | `-z` | metallicity, the same for all binaries; to vary the metallicity, run several simulations | 0.02 | 0.0001 – 0.03 |
 | `-I` | read the initial conditions from a file (columns a e M m z); binaries outside the allowed range are skipped | | |
+| `-O` | name of the output file | `SeBa.data` | |
 | `-D` | stop the evolution of a binary when it merges or is disrupted | off | |
 | `-S` | stop the evolution of a binary when one of the stars becomes a remnant | off | |
+| `-K`, `-k` | stellar type of the primary and the secondary at the start, by name (for example `main_sequence`) | main_sequence | |
 
 ### Initial parameters (with `-R`)
 
@@ -155,12 +157,58 @@ Instead of the numbers, the distributions can also be given by name with the cap
 
 Note on the minimum separation: SeBa makes sure that the two stars do not touch at the start. The actual minimum is the larger of the value given with `-a` and the separation at which the two zero-age main-sequence stars touch at pericentre. So with `-a 0` (the default), the actual minimum is set by the sizes of the stars.
 
+### Model parameters
+
+These parameters of the binary evolution can only be set with a word, not with a letter, for example `./SeBa -R -n 1000 --alpha 1 --ce-method 3 --SN_kick hobbs`. Without them SeBa uses its standard model. The values used are printed at the start of the log.
+
+| Option | Meaning | Default | Valid range |
+|---|---|---|---|
+| `--alpha` | common-envelope efficiency α | 4.0 | > 0 |
+| `--lambda` | envelope binding-energy parameter λ | 0.5 | > 0 |
+| `--gamma` | γ of the γ common-envelope formalism | 1.75 | > 0 |
+| `--beta` | specific angular momentum loss β | 2.5 | |
+| `--ce-method` | common-envelope prescription: 1 = α-γ, 2 = γ-γ, 3 = α-α | 1 | 1 – 3 |
+| `--jloss-method` | angular momentum loss during mass transfer: 1 = SeBa standard (with β), 2 = specific angular momentum of the orbit, 3 = of the donor (Jeans mode), 4 = of the accretor (isotropic re-emission) | 1 | 1 – 4 |
+| `--bh-mass-method` | black hole mass: 1 = SeBa (based on Fryer & Kalogera 2001), 2 = Fryer et al. (2012) delayed, 3 = Fryer et al. (2012) rapid, 4 = Fryer et al. (2012) StarTrack | 2 | 1 – 4 |
+| `--SN_kick` | supernova kick distribution, by name or number: 0 = standard (Verbunt et al. 2017), 1 = `none`, 2 = `maxwellian`, 3 = `paczynski`, 4 = `hobbs` (Hobbs et al. 2005), 5 = `arzoumanian` (Arzoumanian et al. 2002), 6 = `verbunt` (Verbunt et al. 2017), 7 = `delta` (every kick the same) | 0 | 0 – 7 |
+| `--SN_kick_dispersion` | velocity dispersion of the kick in km/s; used by `maxwellian` and `paczynski`, and the kick itself for `delta` | 600 | ≥ 0 |
+
+### Long option names
+
+Most letter options also have a long name, which makes scripts easier to read. These two commands do the same:
+```
+./SeBa -R -n 1000 -s 5 -O pop.data
+./SeBa --random --number 1000 --seed 5 --output pop.data
+```
+
+| Long name | Letter |
+|---|---|
+| `--random` | `-R` |
+| `--number` | `-n` |
+| `--first-id` | `-N` |
+| `--seed` | `-s` |
+| `--time` | `-T` |
+| `--metallicity` | `-z` |
+| `--input`, `--output` | `-I`, `-O` |
+| `--stop-at-merger`, `--stop-at-remnant` | `-D`, `-S` |
+| `--type1`, `--type2` | `-K`, `-k` |
+| `--mass1` or `--max-mass`, `--mass2` or `--min-mass` | `-M`, `-m` |
+| `--mass-function`, `--mass-function-name`, `--mass-exponent` | `-f`, `-F`, `-x` |
+| `--min-q`, `--max-q` | `-q`, `-Q` |
+| `--q-distribution`, `--q-distribution-name`, `--q-exponent` | `-p`, `-P`, `-w` |
+| `--separation` or `--min-separation`, `--max-separation` | `-a`, `-A` |
+| `--separation-distribution`, `--separation-distribution-name`, `--separation-exponent` | `-g`, `-G`, `-y` |
+| `--eccentricity` or `--min-eccentricity`, `--max-eccentricity` | `-e`, `-E` |
+| `--eccentricity-distribution`, `--eccentricity-distribution-name`, `--eccentricity-exponent` | `-u`, `-U`, `-v` |
+
+A long name is followed by its value after a space (`--seed 5`); the form `--seed=5` is not understood.
+
 
 ## Understanding the SeBa output
 
 ### The output file
 
-When you run SeBa, the output is written to a file called `SeBa.data` in the current directory. If this file already exists, SeBa does not overwrite it but adds the new output at the end. So remove or rename an old `SeBa.data` before a new run, unless you want to add to it.
+When you run SeBa, the output is written to a file called `SeBa.data` in the current directory. If this file already exists, SeBa does not overwrite it but adds the new output at the end. So remove or rename an old `SeBa.data` before a new run, unless you want to add to it. With `-O name` SeBa writes to the file `name` instead (with the same rule: an existing file is added to).
 
 SeBa adds the evolutionary history of each binary in the SeBa.data file. Every line represents a moment in the evolution of the binary when something interesting happened, for example one of the stars transitions from the main-sequence to the hertzsprung gap, or mass transfer starts or stops. So the lines are written at these events, not at fixed time steps. The meaning of the columns is defined below. The first column represents a unique identifier for each binary.
 

dstar/SeBa.C

@@ -194,11 +194,109 @@ local bool  evolve_binary(dyn * bi,
 
 int main(int argc, char ** argv) {
 
+    // Model parameters: these are set only with a word, e.g. "--alpha 2",
+    // and have no one-letter option. They are read here and taken out of
+    // argv, so pgetopt never sees them.
+    struct { const char *name; real *value; int *method; } word_options[] = {
+        {"--alpha", &cnsts.alpha_ce, NULL},    // common-envelope efficiency
+        {"--lambda", &cnsts.lambda_ce, NULL},  // envelope binding energy
+        {"--beta", &cnsts.beta, NULL},         // specific angular momentum loss
+        {"--gamma", &cnsts.gamma, NULL},       // gamma common-envelope formalism
+        {"--ce-method", NULL, &cnsts.ce_method},       // 1: alpha-gamma, 2: gamma-gamma, 3: alpha-alpha
+        {"--jloss-method", NULL, &cnsts.jloss_method}, // 1: SeBa default, 2: orbit, 3: donor, 4: accretor
+        {"--bh-mass-method", NULL, &cnsts.bh_mass_method}, // 1: SeBa (Fryer 2001), 2: Fryer 2012 delayed, 3: rapid, 4: Startrack
+        {"--SN_kick_dispersion", &cnsts.v_disp, NULL},    // km/s, for the maxwellian, paczynski and delta kicks
+    };
+    // --SN_kick: the supernova kick distribution, by name or by number
+    // (the number is the value in enum super_nova_kick_distribution, constants.h)
+    const char *kick_names[] = {"standard", "none", "maxwellian", "paczynski",
+                                "hobbs", "arzoumanian", "verbunt", "delta"};
+    const int n_kick = sizeof(kick_names) / sizeof(kick_names[0]);
+    const int n_word = sizeof(word_options) / sizeof(word_options[0]);
+    int kept = 1;
+    for (int i = 1; i < argc; i++) {
+        if (strcmp(argv[i], "--SN_kick") == 0) {
+            if (i + 1 >= argc) {
+                cerr << "error: --SN_kick needs a value" << endl;
+                exit(1);
+            }
+            const char *v = argv[++i];
+            int kick = -1;
+            for (int n = 0; n < n_kick; n++)
+                if (strcmp(v, kick_names[n]) == 0) kick = n;
+            if (kick < 0 && v[0] >= '0' && v[0] <= '9') kick = atoi(v);
+            if (kick < 0 || kick >= n_kick) {
+                cerr << "error: --SN_kick must be one of 0-7 or none, maxwellian, paczynski,"
+                     << " hobbs, arzoumanian, verbunt, delta" << endl;
+                exit(1);
+            }
+            cnsts.pk = (super_nova_kick_distribution) kick;
+            continue;
+        }
+        int k = 0;
+        while (k < n_word && strcmp(argv[i], word_options[k].name) != 0) k++;
+        if (k == n_word) {             // not a model parameter: keep it
+            argv[kept++] = argv[i];
+            continue;
+        }
+        if (i + 1 >= argc) {
+            cerr << "error: " << argv[i] << " needs a value" << endl;
+            exit(1);
+        }
+        if (word_options[k].value) *word_options[k].value = atof(argv[++i]);
+        else *word_options[k].method = atoi(argv[++i]);
+    }
+    argv[kept] = NULL;
+    argc = kept;
+    if (cnsts.ce_method < 1 || cnsts.ce_method > 3)
+        err_exit("--ce-method must be 1, 2 or 3");
+    if (cnsts.jloss_method < 1 || cnsts.jloss_method > 4)
+        err_exit("--jloss-method must be 1, 2, 3 or 4");
+    if (cnsts.bh_mass_method < 1 || cnsts.bh_mass_method > 4)
+        err_exit("--bh-mass-method must be 1, 2, 3 or 4");
+    if (cnsts.v_disp < 0)
+        err_exit("--SN_kick_dispersion must be 0 or larger");
+
+    // Long option names: each --name is translated into its one-letter
+    // option before pgetopt reads the command line. The letters keep working.
+    // Only the form "--alpha 2" is understood, not "--alpha=2".
+    static const char *long_options[][2] = {
+        {"--random", "-R"}, {"--number", "-n"}, {"--first-id", "-N"},
+        {"--seed", "-s"}, {"--time", "-T"}, {"--metallicity", "-z"},
+        {"--mass1", "-M"}, {"--max-mass", "-M"},
+        {"--mass2", "-m"}, {"--min-mass", "-m"},
+        {"--separation", "-a"}, {"--min-separation", "-a"},
+        {"--max-separation", "-A"},
+        {"--eccentricity", "-e"}, {"--min-eccentricity", "-e"},
+        {"--max-eccentricity", "-E"},
+        {"--min-q", "-q"}, {"--max-q", "-Q"},
+        {"--mass-function", "-f"}, {"--mass-function-name", "-F"},
+        {"--mass-exponent", "-x"},
+        {"--q-distribution", "-p"}, {"--q-distribution-name", "-P"},
+        {"--q-exponent", "-w"},
+        {"--separation-distribution", "-g"},
+        {"--separation-distribution-name", "-G"},
+        {"--separation-exponent", "-y"},
+        {"--eccentricity-distribution", "-u"},
+        {"--eccentricity-distribution-name", "-U"},
+        {"--eccentricity-exponent", "-v"},
+        {"--input", "-I"}, {"--output", "-O"},
+        {"--stop-at-merger", "-D"}, {"--stop-at-remnant", "-S"},
+        {"--type1", "-K"}, {"--type2", "-k"},
+    };
+    const int n_long = sizeof(long_options) / sizeof(long_options[0]);
+    for (int i = 1; i < argc; i++)
+        for (int k = 0; k < n_long; k++)
+            if (strcmp(argv[i], long_options[k][0]) == 0) {
+                argv[i] = (char *) long_options[k][1];
+                break;
+            }
+
     bool e_flag = false;
     bool R_flag = false;
     bool F_flag = false;
     bool I_flag = false;
-    bool O_flag = false;//doesn't work fully as "SeBa.data" is written in the code multiple times
+    bool O_flag = false;
     bool P_flag = false;
     bool U_flag = false;
     bool G_flag = false;
@@ -258,15 +356,15 @@ int main(int argc, char ** argv) {
 
     int input_seed=0;
     char seedlog[64];
-    char paramlog[120];
+    char paramlog[400];
 
     //check_help();
 
 //    seba_counters* new_seba_counters = new seba_counters;
 
     extern char *poptarg;
     int c;
-    const char *param_string = "n:N:RDSM:m:x:F:f:A:a:y:G:g:E:e:v:U:u:Q:q:T:t:I:O:w:P:p:n:s:z:C:c:";
+    const char *param_string = "n:N:RDSM:m:x:F:f:A:a:y:G:g:E:e:v:U:u:Q:q:T:t:I:O:w:P:p:n:s:z:K:k:";
 
     while ((c = pgetopt(argc, argv, param_string)) != -1)
 	switch(c) {
@@ -321,6 +419,7 @@ int main(int argc, char ** argv) {
 		      break;
             case 'O': O_flag = true;
 		      output_filename = poptarg;
+		      cnsts.output_file = poptarg;
 		      break;
             case 'w': q_exp = atof(poptarg);
 		      break;
@@ -337,26 +436,30 @@ int main(int argc, char ** argv) {
 		      break;
            case 'z': metal = atof(poptarg);
                 break;
-           case 'C': strcpy(star_type_string, poptarg);
+           case 'K': strcpy(star_type_string, poptarg);
 	           primary_type = extract_stellar_type_string(star_type_string);
                 break;
-           case 'c': strcpy(star_type_string, poptarg);
+           case 'k': strcpy(star_type_string, poptarg);
 	           secondary_type = extract_stellar_type_string(star_type_string);
                 break;
             case '?': params_to_usage(cerr, argv[0], param_string);
+		      cerr << "Model parameters (only as a word): --alpha --lambda --beta --gamma"
+		           << " --ce-method --jloss-method --bh-mass-method --SN_kick --SN_kick_dispersion" << endl;
 		      exit(1);
 	}
 
     int actual_seed = srandinter(input_seed);
     cerr << "random number generator seed = " << actual_seed << endl;
     sprintf(paramlog,
-	   "   alpha  = %3.1f\n   lambda = %3.1f\n   beta   = %3.1f\n   gamma  = %4.2f\n   CE_method = %d\n   Jloss_method = %d \n",
+	   "   alpha  = %3.1f\n   lambda = %3.1f\n   beta   = %3.1f\n   gamma  = %4.2f\n   CE_method = %d\n   Jloss_method = %d \n   BH_mass_method = %d\n   SN_kick = %s (%d), dispersion = %.0f km/s\n",
 	    cnsts.parameters(common_envelope_efficiency),
 	    cnsts.parameters(envelope_binding_energy),
 	    cnsts.parameters(specific_angular_momentum_loss),
 	    cnsts.parameters(dynamic_mass_transfer_gamma),
             cnsts.use_common_envelope_method(),
-            cnsts.use_jloss_method()
+            cnsts.use_jloss_method(),
+            cnsts.use_black_hole_mass_method(),
+            kick_names[cnsts.pk], (int) cnsts.pk, cnsts.v_disp
           );
 
     if (n <= 0) err_exit("mknodes: N > 0 required!");

dstar/init/initdouble.C

@@ -689,7 +689,7 @@ cerr << "addstar(node... called from adddouble(node ..." << endl;
                          type);
 
 	    // added spz:8Aug2002
-	    new_double->dump("SeBa.data", true);
+	    new_double->dump(cnsts.output_file, true);
 
 	    if (REPORT_ADD_DOUBLE) {
 	      cerr<<"double: "<<new_double<<endl;

dstar/starclass/double_star.C

@@ -869,7 +869,7 @@ void double_star::perform_mass_transfer(const real dt,
 // (GN+SilT Mar  2 2011) new dumping regime
     if (!first_contact) {
 
-        dump("SeBa.data", true);
+        dump(cnsts.output_file, true);
         first_contact = true;
     }
 
@@ -1693,7 +1693,7 @@ void double_star::double_spiral_in() {
 
        // (GN+SilT Mar  2 2011) new dumping regime
        bin_type = Double_Spiral_In;
-       dump("SeBa.data", true);
+       dump(cnsts.output_file, true);
        current_mass_transfer_type = Unknown;
 
 
@@ -1850,7 +1850,7 @@ void double_star::spiral_in(star* larger,
 
       // (GN+SilT Mar  2 2011) new dumping regime
       bin_type = Spiral_In;
-      dump("SeBa.data", true);
+      dump(cnsts.output_file, true);
       current_mass_transfer_type = Unknown;
 
 
@@ -1967,7 +1967,7 @@ void double_star::merge_elements(star* consumer,
 //cerr<<"bin: "<<identity<<endl;
   bin_type = Merged;
 //  dump(cerr, false);
-  dump("SeBa.data", true);
+  dump(cnsts.output_file, true);
 
   if (!get_use_hdyn()) {
     consumer->set_velocity(velocity);
@@ -1987,7 +1987,7 @@ void double_star::merge_elements(star* consumer,
 //  dump("binev.data", false);
 cerr << "Merger is: "<<endl;
   dump(cerr, false);
-  dump("SeBa.data", true);
+  dump(cnsts.output_file, true);
 
   //get_seba_counters()->mergers++;
 }
@@ -2233,7 +2233,7 @@ void double_star::contact_binary(real dt) {
     if (!first_contact) {
 
         bin_type = Contact;
-        dump("SeBa.data", true);
+        dump(cnsts.output_file, true);
         first_contact = true;
     }
 
@@ -2277,7 +2277,7 @@ void double_star::angular_momentum_envelope_ejection(star* larger,
 //  cerr << "double_star::angular_momentum_envelope_ejection()"<<endl;
   // (GN+SilT Mar  2 2011) new dumping regime
   bin_type = Common_Envelope;
-  dump("SeBa.data", true);
+  dump(cnsts.output_file, true);
   current_mass_transfer_type = Unknown;
 
 

include/star/constants.h

@@ -204,14 +204,20 @@ enum dynamics_update_parameter {stellar_mass_update_limit,
 extern
 class stellar_evolution_constants {  // Easy to have a name for compiling.
   public:
-  stellar_evolution_constants() {
-    pk = internally_decided_velocity_kick;
-    v_disp = 600;
-  }
+  stellar_evolution_constants();    // sets the default values, see constants.C
 
   super_nova_kick_distribution pk;
   real v_disp;
 
+  char * output_file;
+  real alpha_ce;
+  real lambda_ce;
+  real beta;
+  real gamma;
+  int ce_method;
+  int jloss_method;
+  int bh_mass_method;
+
   real mathematics(mathematical_constant);
   real physics(physics_constants);