diff --git a/README.md b/README.md index 530c20b..b41f8cc 100644 --- a/README.md +++ b/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. diff --git a/dstar/SeBa.C b/dstar/SeBa.C index f4890ee..4461fba 100644 --- a/dstar/SeBa.C +++ b/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,7 +356,7 @@ int main(int argc, char ** argv) { int input_seed=0; char seedlog[64]; - char paramlog[120]; + char paramlog[400]; //check_help(); @@ -266,7 +364,7 @@ int main(int argc, char ** argv) { 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!"); diff --git a/dstar/init/initdouble.C b/dstar/init/initdouble.C index 395057c..0aac480 100644 --- a/dstar/init/initdouble.C +++ b/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: "<set_velocity(velocity); @@ -1987,7 +1987,7 @@ void double_star::merge_elements(star* consumer, // dump("binev.data", false); cerr << "Merger is: "<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()"<set_first_contact(false); if (is_binary_component()){ - get_binary()->dump("SeBa.data", true); + get_binary()->dump(cnsts.output_file, true); } else { } @@ -238,7 +238,7 @@ void single_star::post_constructor() { if (is_binary_component() && get_binary()->roche_radius(this) > radius) { - get_binary()->dump("SeBa.data", true); + get_binary()->dump(cnsts.output_file, true); } if (remnant()) { diff --git a/std/gethist.C b/std/gethist.C index 06ea82f..36c06f6 100644 --- a/std/gethist.C +++ b/std/gethist.C @@ -5,7 +5,7 @@ * Institute for Advanced Study, Princeton, NJ, USA * version 2: Dec 1992 Piet Hut -- adopted to the new C++-based starlab *............................................................................. - * non-local function: + * non-local function: * gethist *............................................................................. * A string is returned which contains the date at which a command is @@ -45,12 +45,12 @@ local bool contains_white(char * s) /*----------------------------------------------------------------------------- * eff_length -- effective length, needed to stored a command line argument: * one for each character, and two extra for the surrounding - * "" quotes, in case the string contains one or more + * "" quotes, in case the string contains one or more * internal white spaces. * example: * "hi" --> hi * "hi there" --> "hi there" - * + * * eff_length("hi") = 2 * eff_length("hi there") = 10 *----------------------------------------------------------------------------- @@ -67,7 +67,7 @@ local int eff_length(char * s) * DATE_TOT_LENGTH -- total length, in characters, of the date plus the * date separator; this is the length reserved in total * for the date information by gethist(). - * DATE_STR_LENGTH -- The fixed standard length of the date containing + * DATE_STR_LENGTH -- The fixed standard length of the date containing * string, as returned by the UNIX operating system, * apart from the last NULL character; * DATE_SEPARATOR -- This string can be changed whenever desired: @@ -96,7 +96,7 @@ local int eff_length(char * s) * the command of the calling program was given, and the full * command line. * note: - * simply printing successive command line arguments does not + * simply printing successive command line arguments does not * provide the full information, since the following two commands * gobble this thing and gobble "this thing" * would both be printed as gobble this thing . @@ -109,7 +109,7 @@ char *gethist(int argc, char ** argv) { int i, j, k; int n; - long int clock; /* will contain time in seconds since 1970 */ + time_t clock; /* will contain time in seconds since 1970 */ char *hist_string; /* will point to newly allocated string */ /* * determine length of new character array, and allocate memory for it: @@ -130,7 +130,7 @@ char *gethist(int argc, char ** argv) n += user_length; hist_string = new char[n]; - + if (hist_string == NULL) { cerr << "gethist: not enough memory left for hist_string\n"; @@ -140,8 +140,9 @@ char *gethist(int argc, char ** argv) * write the date: */ strcpy(hist_string, INTRO_PART); - clock = (long int)time(0); /* see time(2), UNIX manual */ - strcpy(hist_string + INTRO_LENGTH, ctime((time_t*)&clock)); + clock = time(0); /* see time(2), UNIX manual */ + const char *date = ctime(&clock); + strcpy(hist_string + INTRO_LENGTH, date ? date :"unknown date \n"); /* see ctime(3C), UNIX manual */ strcpy(hist_string + INTRO_LENGTH + DATE_STR_LENGTH, DATE_SEPARATOR); strcpy(hist_string + INTRO_LENGTH + DATE_STR_LENGTH + DATE_SEP_LENGTH,