Loading [MathJax]/extensions/tex2jax.js
PeTar
N-body code for collisional gravitational systems
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
bse_interface.h
Go to the documentation of this file.
1 #pragma once
2 #include <iostream>
3 #include <iomanip>
4 #include <cassert>
5 #include <cstdio>
6 #include <string>
7 #include <getopt.h>
8 #include "../src/io.hpp"
9 
51 #if (defined BSEBBF) || (defined BSEEMP)
52 extern "C" {
53  // The COMMON variables in BSE-base codes should be declared here
54  // Then, PeTar can read and modify these variables
55  // A struct in C corresponds to a COMMON block in Fortran.
56  // The suffix '_' is added to the name of COMMON block in C.
57 
58  extern struct{
59  double neta;
60  double bwind;
61  double hewind;
62  } value1_;
63 
64  extern struct{
65  double alpha;
66  double lambda;
67  } value2_;
68 
69  extern struct{
70  int idum;
71  } value3_;
72 
73  extern struct{
74  double sigma;
75  double mxns;
76  int bhflag;
77  } value4_;
78 
79  extern struct{
80  double beta;
81  double xi;
82  double bhwacc;
83  double epsnov;
84  double eddfac;
85  double gamma;
86  } value5_;
87 
88  extern struct{
89  int ceflag;
90  int tflag;
91  int ifflag;
92  int nsflag;
93  int wdflag;
94  } flags_;
95 
96  extern struct{
97  int psflag;
98  int kmech;
99  int ecflag;
100  } flags2_;
101 
102  extern struct{
103  double pts1;
104  double pts2;
105  double pts3;
106  } points_;
107 
108  extern struct{
109  int idum2;
110  int iy;
111  int ir[32];
112  } rand3_;
113 
114 
115  // The Fortran function name + '_' is the C version. All arguments should be in pointer type.
116 
118 #ifdef BSEEMP
119  void zcnsts_(double* z, double* zpars, int* trackmode);
120 #else
121  void zcnsts_(double* z, double* zpars);
122 #endif
123 
125  void instar_();
126 
128  void evolv1_(int* kw, double* mass, double* mt, double* r, double* lum, double* mc, double* rc, double* menv, double* renv, double* ospin,
129  double* epoch, double* tm, double* tphys, double* tphysf, double* dtp, double* z, double* zpars, double* vs);
130 
132  void evolv2_(int* kw, double* mass, double* mt, double* r, double* lum, double* mc, double* rc, double* menv, double* renv, double* ospin,
133  double* epoch, double* tm, double* tphys, double* tphysf, double* dtp, double* z, double* zpars,
134  double* period, double* ecc, double* bse_event, double* vkick);
135 
136  void star_(int* kw, double* mass, double* mt, double* tm, double* tn, double* tscls, double* lums, double* GB, double* zpars);
137 
138  void deltat_(int* kw, double* age, double* tm, double* tn, double* tscls, double* dt, double* dtr);
139 
140  void trdot_(int* kw, double* m0, double* mt, double* r, double* mc, double* rc, double* age, double* dt, double* dtr, double* zpars);
141 
142  void trflow_(int* kw, double* m0, double* mt, double* r, double* mc, double* rc, double* age, double* dt, double* semi, double* ecc, double* zpars);
143 
144  void mix_(double* m0, double* mt, double* age, int* kw, double* zpars);
145 
146  //void comenv_(double* m01, double* m1, double* mc1, double* aj1, double* jspin1, int* kw1,
147  // double* m02, double* m2, double* mc2, double* aj2, double* jspin2, int* kw2,
148  // double* zpars, double* ecc, double* sep, double* jorb, double* vkick1, double* vkick2, int* coel);
149 
150  void merge_(int* kw, double* m0, double* mt, double* r, double* mc, double* rc, double* menv, double* renv, double* ospin, double* age, double* semi, double* ecc, double* vkick, double* zpars);
151 
152  void printconst_();
153 }
154 
155 #elif MOBSE
156 extern "C" {
157  extern struct{
158  double neta;
159  double bwind;
160  double hewind;
161  } value1_;
162 
163  extern struct{
164  double alpha;
165  double lambda;
166  } value2_;
167 
168  extern struct{
169  int idum;
170  } value3_;
171 
172  extern struct{
173  double sigma1;
174  double sigma2;
175  double mxns;
176  int bhflag;
177  } value4_;
178 
179  extern struct{
180  double beta;
181  double xi;
182  double bhwacc;
183  double epsnov;
184  double eddfac;
185  double gamma;
186  } value5_;
187 
188  extern struct{
189  int ceflag;
190  int tflag;
191  int ifflag;
192  int nsflag;
193  int wdflag;
194  int piflag;
195  } flags_;
196 
197 // extern struct{
198 // int psflag; ///> PPSN condition
199 // int kmech; ///> kick mechanism
200 // int ecflag; ///> ECS switcher
201 // } flags2_;
202 
203  extern struct{
204  double pts1;
205  double pts2;
206  double pts3;
207  } points_;
208 
209  extern struct{
210  int idum2;
211  int iy;
212  int ir[32];
213  } rand3_;
214 
216  void zcnsts_(double* z, double* zpars);
217 
219  void instar_();
220 
222  void evolv1_(int* kw, double* mass, double* mt, double* r, double* lum, double* mc, double* rc, double* menv, double* renv, double* ospin,
223  double* epoch, double* tm, double* tphys, double* tphysf, double* dtp, double* z, double* zpars, double* vkick);
224 
226  void evolv2_(int* kw, double* mass, double* mt, double* r, double* lum, double* mc, double* rc, double* menv, double* renv, double* ospin,
227  double* epoch, double* tm, double* tphys, double* tphysf, double* dtp, double* z, double* zpars,
228  double* period, double* ecc, double* bse_event, double* vkick);
229 
230  void star_(int* kw, double* mass, double* mt, double* tm, double* tn, double* tscls, double* lums, double* GB, double* zpars);
231 
232  void deltat_(int* kw, double* age, double* tm, double* tn, double* tscls, double* dt, double* dtr);
233 
234  void trdot_(int* kw, double* m0, double* mt, double* r, double* mc, double* rc, double* age, double* dt, double* dtr, double* zpars);
235 
236  void trflow_(int* kw, double* m0, double* mt, double* r, double* mc, double* rc, double* age, double* dt, double* semi, double* ecc, double* zpars);
237 
238  //void mix_(double* m0, double* mt, double* age, int* kw, double* zpars);
239 
240  void mix_(double* m0, double* mt, double* age, int* kw, double* zpars, int* krol);
241 
242  //void comenv_(double* m01, double* m1, double* mc1, double* aj1, double* jspin1, int* kw1,
243  // double* m02, double* m2, double* mc2, double* aj2, double* jspin2, int* kw2,
244  // double* zpars, double* ecc, double* sep, double* jorb, double* vkick1, double* vkick2, int* coel);
245 
246  void merge_(int* kw, double* m0, double* mt, double* r, double* mc, double* rc, double* menv, double* renv, double* ospin, double* age, double* semi, double* ecc, double* vkick, double* zpars);
247 
248  void printconst_();
249 }
250 #endif
251 
253 
257 struct StarParameter{
258  long long int kw;
259  double m0;
260  double mt;
261  double r;
262  double mc;
263  double rc;
264  double ospin;
265  double epoch;
266  double tphys;
267  double lum;
268 
270 
276  void initial(double _mass, int _kw=1, double _ospin=0.0, double _epoch=0.0) {
277  kw = _kw;
278  m0 = _mass;
279  mt = _mass;
280  r = 0.0;
281  mc = 0.0;
282  rc = 0.0;
283  ospin = _ospin;
284  epoch = _epoch;
285  tphys = _epoch;
286  }
287 
289 
291  void writeAscii(FILE* fp) const{
292  fprintf(fp, "%lld %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e %26.17e ",
293  this->kw, this->m0, this->mt, this->r, this->mc, this->rc, this->ospin, this->epoch, this->tphys, this->lum);
294  }
295 
297 
299  void readAscii(FILE* fp) {
300  int rcount=fscanf(fp, "%lld %lf %lf %lf %lf %lf %lf %lf %lf %lf ",
301  &this->kw, &this->m0, &this->mt, &this->r, &this->mc, &this->rc, &this->ospin, &this->epoch, &this->tphys, & this->lum);
302  if(rcount<10) {
303  std::cerr<<"Error: Data reading fails! requiring data number is 10, only obtain "<<rcount<<".\n";
304  abort();
305  }
306  }
307 
309  void print(std::ostream & fout) const{
310  fout<<" type= "<<kw
311  <<" m0[M*]= "<<m0
312  <<" m[M*]= "<<mt
313  <<" rad[R*]= "<<r
314  <<" mc[M*]= "<<mc
315  <<" rc[M*]= "<<rc
316  <<" spin= "<<ospin
317  <<" epoch= "<<epoch
318  <<" t[myr]= "<<tphys
319  <<" lum[L*]= "<<lum;
320  }
321 
323 
327  static void printColumnTitle(std::ostream & _fout, const int _width=20) {
328  _fout<<std::setw(_width)<<"s_type"
329  <<std::setw(_width)<<"s_mass0[M*]"
330  <<std::setw(_width)<<"s_mass[M*]"
331  <<std::setw(_width)<<"s_rad[R*]"
332  <<std::setw(_width)<<"s_mcore[M*]"
333  <<std::setw(_width)<<"s_rcore[R*]"
334  <<std::setw(_width)<<"s_spin"
335  <<std::setw(_width)<<"s_epoch[Myr]"
336  <<std::setw(_width)<<"s_time[Myr]"
337  <<std::setw(_width)<<"s_lum[L*]";
338  }
339 
341 
345  void printColumn(std::ostream & _fout, const int _width=20) const{
346  _fout<<std::setw(_width)<<kw
347  <<std::setw(_width)<<m0
348  <<std::setw(_width)<<mt
349  <<std::setw(_width)<<r
350  <<std::setw(_width)<<mc
351  <<std::setw(_width)<<rc
352  <<std::setw(_width)<<ospin
353  <<std::setw(_width)<<epoch
354  <<std::setw(_width)<<tphys
355  <<std::setw(_width)<<lum;
356  }
357 
359 
364  static int printTitleWithMeaning(std::ostream & _fout, const int _counter=0, const int _offset=0) {
365  int counter = _counter;
366  counter++;
367  _fout<<std::setw(_offset)<<" "<<counter<<". s_type: SSE/BSE based code stellar type\n";
368  counter++;
369  _fout<<std::setw(_offset)<<" "<<counter<<". s_mass0: initial mass at each evolution stage [Msun]\n";
370  counter++;
371  _fout<<std::setw(_offset)<<" "<<counter<<". s_mass: current mass [Msun]\n";
372  counter++;
373  _fout<<std::setw(_offset)<<" "<<counter<<". s_rad: stellar radius [Rsun]\n";
374  counter++;
375  _fout<<std::setw(_offset)<<" "<<counter<<". s_mcore: stellar core mass [Msun]\n";
376  counter++;
377  _fout<<std::setw(_offset)<<" "<<counter<<". s_rcore: stellar core radius [Rsun]\n";
378  counter++;
379  _fout<<std::setw(_offset)<<" "<<counter<<". s_spin: stellar rotation\n";
380  counter++;
381  _fout<<std::setw(_offset)<<" "<<counter<<". s_epoch: time offset at each evolution stage [Myr]\n";
382  counter++;
383  _fout<<std::setw(_offset)<<" "<<counter<<". s_time: physical time [Myr]\n";
384  counter++;
385  _fout<<std::setw(_offset)<<" "<<counter<<". s_lum: luminosity [Lsun]\n";
386  return counter;
387  }
388 };
389 
391 
393 struct StarParameterOut{
394  long long int kw0;
395  double dtmiss;
396  double menv;
397  double renv;
398  double tm;
399  double vkick[4];
400  double dm;
401 
402  StarParameterOut(): kw0(0), dtmiss(0.0), menv(0.0), renv(0.0), tm(0.0), vkick{0.0}, dm(0.0) {}
403 
405 
409  static void printColumnTitle(std::ostream & _fout, const int _width=20) {
410  _fout<<std::setw(_width)<<"type0"
411  <<std::setw(_width)<<"m_CE[M*]"
412  <<std::setw(_width)<<"r_CE[R*]"
413  <<std::setw(_width)<<"t_MS[Myr]"
414  <<std::setw(_width)<<"vk.x[km/s]"
415  <<std::setw(_width)<<"vk.y[km/s]"
416  <<std::setw(_width)<<"vk.z[km/s]"
417  <<std::setw(_width)<<"vk[km/s]"
418  <<std::setw(_width)<<"m_loss[M*]";
419  }
420 
422 
426  void printColumn(std::ostream & _fout, const int _width=20) const{
427  _fout<<std::setw(_width)<<kw0
428  <<std::setw(_width)<<menv
429  <<std::setw(_width)<<renv
430  <<std::setw(_width)<<tm
431  <<std::setw(_width)<<vkick[0]
432  <<std::setw(_width)<<vkick[1]
433  <<std::setw(_width)<<vkick[2]
434  <<std::setw(_width)<<vkick[3]
435  <<std::setw(_width)<<dm;
436  }
437 
439  void print(std::ostream & fout) const{
440  fout<<" type0= "<<kw0
441  <<" menv[M*]= "<<menv
442  <<" renv[R*]= "<<renv
443  <<" t_MS[Myr]= "<<tm
444  <<" vk.x[km/s]= "<<vkick[0]
445  <<" vk.y[km/s]= "<<vkick[1]
446  <<" vk.z[km/s]= "<<vkick[2]
447  <<" vk[km/s]= "<<vkick[3]
448  <<" m_loss[M*]= "<<dm;
449  }
450 
451 };
452 
454 
458 static double EstimateRocheRadiusOverSemi(double& _q) {
459  double p = std::pow(_q,1.0/3.0);
460  double p2= p*p;
461  double rad_semi = 0.49*p2/(0.6*p2 + std::log(1.0+p));
462  return rad_semi;
463 }
464 
466 
477 static double EstimateGRTimescale(StarParameter& _star1, StarParameter& _star2, double& _semi, double& _ecc) {
478  double ecc2 = _ecc*_ecc;
479  double omecc2 = 1.0 - ecc2;
480  double sqome2 = std::sqrt(omecc2);
481  double sqome5 = std::pow(sqome2,5.0);
482  double semi2 = _semi*_semi;
483  // (32/5)(G^{7/2}/c^5 is ~8.3x10^{-10} in the unit of Msun, Rsun, and year.
484  double djgr = 8.315e-10*_star1.mt*_star2.mt*(_star1.mt+_star2.mt)/(semi2*semi2)*(1.0+0.875*ecc2)/sqome5;
485  double dtr = 1.0e-9/djgr*(1-_ecc); // in Myr, with 0.001 coefficent and ecc correction factor
486  return dtr;
487 }
488 
489 
491 
493 class BinaryEvent{
494 public:
495  // Tanikawa's BH model
496  double record[20][9];
497  //double record[20][81];
498  //
499 
501  void recordEvent(const StarParameter& _p1, const StarParameter& _p2, const double _semi, const double _ecc, const int _type, const int _index) {
502  record[0][_index] = std::min(_p1.tphys, _p2.tphys);
503  record[1][_index] = _p1.mt;
504  record[2][_index] = _p2.mt;
505  record[3][_index] = _p1.kw;
506  record[4][_index] = _p2.kw;
507  record[5][_index] = _semi;
508  record[6][_index] = _ecc;
509  double q = _p1.mt/_p2.mt;
510  double rl1 = EstimateRocheRadiusOverSemi(q);
511  q = 1.0/q;
512  double rl2 = EstimateRocheRadiusOverSemi(q);
513  record[7][_index] = _p1.r/(rl1*_semi);
514  record[8][_index] = _p2.r/(rl2*_semi);
515  record[9][_index] = _type;
516  record[10][_index] = _p1.lum;
517  record[11][_index] = _p2.lum;
518  record[12][_index] = _p1.r;
519  record[13][_index] = _p2.r;
520  record[14][_index] = _p1.mc;
521  record[15][_index] = _p2.mc;
522  record[16][_index] = _p1.rc;
523  record[17][_index] = _p2.rc;
524  record[18][_index] = _p1.ospin;
525  record[19][_index] = _p2.ospin;
526  }
527 
529  void setEventIndexEnd(const int index) {
530  record[9][index] = -1;
531  }
532 
534  int getEventNMax() const {
535  return 8;
536  }
537 
539  int getEventIndexInit() const {
540  return 8;
541  }
542 
544  int getType(const int index) const {
545  return int(record[9][index]);
546  }
547 
549  void print(std::ostream & fout, const int index) const{
550  fout<<" t[Myr]= "<<record[0][index]
551  <<" m1[M*]= "<<record[1][index]
552  <<" m2[M*]= "<<record[2][index]
553  <<" type1= "<<int(record[3][index])
554  <<" type2= "<<int(record[4][index])
555  <<" semi[R*]= "<<record[5][index]
556  <<" ecc= "<<record[6][index]
557  <<" rad1[Ro]= "<<record[7][index]
558  <<" rad2[Ro]= "<<record[8][index]
559  <<" lum1[L*]= "<<record[10][index]
560  <<" lum2[L*]= "<<record[11][index]
561  <<" rad1[R*]= "<<record[12][index]
562  <<" rad2[R*]= "<<record[13][index]
563  <<" mc1[M*]= "<<record[14][index]
564  <<" mc2[M*]= "<<record[15][index]
565  <<" rc1[R*]= "<<record[16][index]
566  <<" rc2[R*]= "<<record[17][index]
567  <<" ospin1= "<<record[18][index]
568  <<" ospin2= "<<record[19][index];
569  }
570 
572 
576  static void printColumnTitle(std::ostream & _fout, const int _width=20) {
577  _fout<<std::setw(_width)<<"t[Myr]"
578  <<std::setw(_width)<<"m1[M*]"
579  <<std::setw(_width)<<"m2[M*]"
580  <<std::setw(_width)<<"type1"
581  <<std::setw(_width)<<"type2"
582  <<std::setw(_width)<<"semi[R*]"
583  <<std::setw(_width)<<"ecc"
584  <<std::setw(_width)<<"rad1[Ro]"
585  <<std::setw(_width)<<"rad2[Ro]"
586  <<std::setw(_width)<<"btype"
587  <<std::setw(_width)<<"lum1[L*]"
588  <<std::setw(_width)<<"lum2[L*]"
589  <<std::setw(_width)<<"rad1[R*]"
590  <<std::setw(_width)<<"rad2[R*]"
591  <<std::setw(_width)<<"mc1[M*]"
592  <<std::setw(_width)<<"mc2[M*]"
593  <<std::setw(_width)<<"rc1[R*]"
594  <<std::setw(_width)<<"rc2[R*]"
595  <<std::setw(_width)<<"ospin1"
596  <<std::setw(_width)<<"ospin2";
597 
598  }
599 
601 
606  void printColumn(std::ostream & _fout, const int _index, const int _width=20) const{
607  for (int i=0; i<3; i++) _fout<<std::setw(_width)<<record[i][_index];
608  for (int i=3; i<5; i++) _fout<<std::setw(_width)<<int(record[i][_index]);
609  for (int i=5; i<9; i++) _fout<<std::setw(_width)<<record[i][_index];
610  _fout<<std::setw(_width)<<int(record[9][_index]);
611  for (int i=10; i<20; i++) _fout<<std::setw(_width)<<record[i][_index];
612  }
613 };
614 
616 
619 class IOParamsBSE{
620 public:
621  IOParamsContainer input_par_store;
622  IOParams<double> neta;
623  IOParams<double> bwind;
624  IOParams<double> hewind;
625  IOParams<double> alpha;
626  IOParams<double> lambda;
627  IOParams<double> beta;
628  IOParams<double> xi;
629  IOParams<double> bhwacc;
630  IOParams<double> epsnov;
631  IOParams<double> eddfac;
632  IOParams<double> gamma;
633  //IOParams<double> mxns;
634 #if (defined BSEBBF) || (defined BSEEMP)
635  IOParams<double> sigma;
636 #elif MOBSE
637  IOParams<double> sigma1;
638  IOParams<double> sigma2;
639 #endif
640  IOParams<long long int> ceflag;
641  IOParams<long long int> tflag;
642  //IOParams<long long int> ifflag;
643  IOParams<long long int> wdflag;
644  IOParams<long long int> bhflag;
645  IOParams<long long int> nsflag;
646 #if (defined BSEBBF) || (defined BSEEMP)
647  IOParams<long long int> psflag;
648  IOParams<long long int> kmech;
649  IOParams<long long int> ecflag;
650 #elif MOBSE
651  IOParams<long long int> piflag;
652 #endif
653  IOParams<double> pts1;
654  IOParams<double> pts2;
655  IOParams<double> pts3;
656 #ifdef BSEEMP
657  IOParams<long long int> trackmode;
658 #endif
659  IOParams<long long int> idum;
660  IOParams<double> tscale;
661  IOParams<double> rscale;
662  IOParams<double> mscale;
663  IOParams<double> vscale;
664  IOParams<double> z;
665 
666  bool print_flag;
667 
668 #if (defined BSEBBF) || (defined BSEEMP)
669  IOParamsBSE(): input_par_store(),
670  neta (input_par_store, 0.5, "bse-neta", "Reimers mass-loss coefficent [neta*4x10^-13]"),
671  bwind (input_par_store, 0.0, "bse-bwind", "Binary enhanced mass loss parameter; inactive for single"),
672  hewind(input_par_store, 1.0, "bse-hewind","Helium star mass loss factor"),
673  //mxns (input_par_store, 1.0, "bse-mxns", "Helium star mass loss factor"),
674  alpha (input_par_store, 3.0, "bse-alpha", "Common-envelope efficiency parameter"),
675  lambda(input_par_store, 0.5, "bse-lambda", "Binding energy factor for common envelope evolution"),
676  beta (input_par_store, 0.125, "bse-beta", "wind velocity factor: proportional to vwind**2"),
677  xi (input_par_store, 1.0, "bse-xi", "wind accretion efficiency factor"),
678  bhwacc(input_par_store, 1.5, "bse-bhwacc", "Bondi-Hoyle wind accretion factor"),
679  epsnov(input_par_store, 0.001, "bse-epsnov", "The fraction of accreted matter retained in nova eruption"),
680  eddfac(input_par_store, 1.0, "bse-eddfac", "Eddington limit factor for mass transfer"),
681  gamma (input_par_store, -1.0, "bse-gamma", "Angular momentum factor for mass lost during Roche"),
682  sigma (input_par_store, 265.0, "bse-sigma", "Dispersion in the Maxwellian for the SN kick speed [km/s]"),
683  ceflag(input_par_store, 0, "bse-ceflag", "if =3, activates de Kool common-envelope model"),
684  tflag (input_par_store, 1, "bse-tflag", "if >0, activates tidal circularisation"),
685  //ifflag(input_par_store, 2, "bse-ifflag", "if > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800"),
686  wdflag(input_par_store, 1, "bse-wdflag", "if >0, uses WD IFMR of HPE, 1995, MNRAS, 272, 800"),
687  bhflag(input_par_store, 2, "bse-bhflag", "BH kick option: 0: no kick; 1: same as NS; 2: scaled by fallback"),
688  nsflag(input_par_store, 3, "bse-nsflag", "NS/BH foramtion options: 0: original SSE; 1: Belczynski (2002); 2: Belczynski (2008); 3: Fryer (2012) rapid SN; 4: Fryer (2012) delayed SN; 5: Eldridge & Tout (2004)"),
689  psflag(input_par_store, 1, "bse-psflag", "PPSN condition (Belczynski 2016): 0: no PPSN; 1: strong; (Leung 2019): 2: moderate; 3: weak"),
690  kmech (input_par_store, 1, "bse-kmech", "Kick mechanism: 1: standard momentum-conserving; 2: convection-asymmetry-driven; 3: collapse-asymmerty-driven; 4: neutrino driven"),
691  ecflag(input_par_store, 1, "bse-ecflag", "if >0, ECS is switched on"),
692  pts1 (input_par_store, 0.05, "bse-pts1", "time step of MS"),
693  pts2 (input_par_store, 0.01, "bse-pts2", "time step of GB, CHeB, AGB, HeGB"),
694  pts3 (input_par_store, 0.02, "bse-pts3", "time step of HG, HeMS"),
695 #ifdef BSEEMP
696  trackmode(input_par_store, 2, "bse-trackmode", "star evolution option, need to make a soft link to the data directory in bse-interface/bseEmp/emptrack/: 1: L model (larger overshoot; directory name: ffbonn); 2: M model (smaller overshoot; directory name: ffgeneva); See details in Appendix A of Tanikawa et al. (2022, ApJ, 926, 83)"),
697 #endif
698  idum (input_par_store, 1234, "bse-idum", "random number seed used by the kick routine"),
699  tscale(input_par_store, 1.0, "bse-tscale", "Time scale factor from input data unit (IN) to Myr (time[Myr]=time[IN]*tscale)"),
700  rscale(input_par_store, 1.0, "bse-rscale", "Radius scale factor from input data unit (IN) to Rsun (r[Rsun]=r[IN]*rscale)"),
701  mscale(input_par_store, 1.0, "bse-mscale", "Mass scale factor from input data unit (IN) to Msun (m[Msun]=m[IN]*mscale)"),
702  vscale(input_par_store, 1.0, "bse-vscale", "Velocity scale factor from input data unit(IN) to km/s (v[km/s]=v[IN]*vscale)"),
703 #ifdef BSEEMP
704  z (input_par_store, 0.001, "bse-metallicity", "Metallicity Z, ranging from 0 to 0.03; when Z<0.0001, using EMP track, please make a symbolic link in the simulation directory to the track directory according to the --bse-trackmode option"),
705 #else
706  z (input_par_store, 0.001, "bse-metallicity", "Metallicity Z, ranging from 0.0001 to 0.03"),
707 #endif
708  print_flag(false) {}
709 #elif MOBSE
710  IOParamsBSE(): input_par_store(),
711  neta (input_par_store, 0.5, "mobse-neta", "Reimers mass-loss coefficent [neta*4x10^-13]"),
712  bwind (input_par_store, 0.0, "mobse-wind", "Binary enhanced mass loss parameter; inactive for single"),
713  hewind(input_par_store, 1.0, "mobse-hewind", "Helium star mass loss factor"),
714  //mxns (input_par_store, 1.0, "Helium star mass loss factor"),
715  alpha (input_par_store, 3.0, "mobse-alpha", "Common-envelope efficiency parameter"),
716  lambda(input_par_store, 0.1, "mobse-lambda", "Binding energy factor for common envelope evolution"),
717  beta (input_par_store, 0.125, "mobse-beta", "wind velocity factor: proportional to vwind**2"),
718  xi (input_par_store, 1.0, "mobse-xi", "wind accretion efficiency factor"),
719  bhwacc(input_par_store, 1.5, "mobse-bwacc", "Bondi-Hoyle wind accretion factor"),
720  epsnov(input_par_store, 0.001, "mobse-epsnov", "The fraction of accreted matter retained in nova eruption"),
721  eddfac(input_par_store, 1.0, "mobse-eddfac", "Eddington limit factor for mass transfer"),
722  gamma (input_par_store, -1.0, "mobse-gamma", "Angular momentum factor for mass lost during Roche"),
723  sigma1 (input_par_store, 265.0, "mobse-sigma1", "Dispersion in the Maxwellian for the CCSN kick speed [km/s]"),
724  sigma2 (input_par_store, 265.0, "mobse-sigma2", "Dispersion in the Maxwellian for the ECSN kick speed [km/s]"),
725  ceflag(input_par_store, 0, "mobse-cflag", "if =3, activates de Kool common-envelope model"),
726  tflag (input_par_store, 1, "mobse-tflag", "if >0, activates tidal circularisation"),
727  //ifflag(input_par_store, 2, "if > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800"),
728  wdflag(input_par_store, 1, "mobse-wdflag", "if >0, uses WD IFMR of HPE, 1995, MNRAS, 272, 800"),
729  bhflag(input_par_store, 3, "mobse-bhflag", "BH kick option: 0: no kick; 1: same as NS; 2: scaled by fallback; 3: Giacobbo&Mapelli (2020)"),
730  nsflag(input_par_store, 3, "mobse-nsflag", "NS/BH formation options: 0: original SSE; 1: Belczynski (2008); 2: Fryer (2012) rapid SN; 3: Fryer (2012) delayed SN; 4: Belczynski (2008); 5: no SN explosion"),
731  piflag(input_par_store, 1, "mobse-piflag", "PPSN condition (Spera et al. 2015)"),
732  //psflag(input_par_store, 1, "PPSN condition (Belczynski 2016): 0: no PPSN; 1: strong; (Leung 2019): 2: moderate; 3: weak"),
733  //kmech (input_par_store, 1, "Kick mechanism: 1: standard momentum-conserving; 2: convection-asymmetry-driven; 3: collapse-asymmerty-driven; 4: neutrino driven"),
734  //ecflag(input_par_store, 1, "if >0, ECS is switched on"),
735  pts1 (input_par_store, 0.05, "mobse-pts1", "time step of MS"),
736  pts2 (input_par_store, 0.01, "mobse-pts2", "time step of GB, CHeB, AGB, HeGB"),
737  pts3 (input_par_store, 0.02, "mobse-pts3", "time step of HG, HeMS"),
738  idum (input_par_store, 1234, "mobse-idum", "random number seed used by the kick routine"),
739  tscale(input_par_store, 1.0, "mobse-tscale", "Time scale factor from input data unit (IN) to Myr (time[Myr]=time[IN]*tscale)"),
740  rscale(input_par_store, 1.0, "mobse-rscale", "Radius scale factor from input data unit (IN) to Rsun (r[Rsun]=r[IN]*rscale)"),
741  mscale(input_par_store, 1.0, "mobse-msclae", "Mass scale factor from input data unit (IN) to Msun (m[Msun]=m[IN]*mscale)"),
742  vscale(input_par_store, 1.0, "mobse-vsclae", "Velocity scale factor from input data unit(IN) to km/s (v[km/s]=v[IN]*vscale)"),
743  z (input_par_store, 0.001, "mobse-metallicity", "Metallicity"),
744  print_flag(false) {}
745 #endif
746 
748 
754  int read(int argc, char *argv[], const int opt_used_pre=0) {
755  static int sse_flag=-1;
756  const struct option long_options[] = {
757  {neta.key, required_argument, &sse_flag, 0},
758  {bwind.key, required_argument, &sse_flag, 1},
759  {hewind.key, required_argument, &sse_flag, 2},
760  //{mxns.key, required_argument, &sse_flag, 3},
761  {alpha.key, required_argument, &sse_flag, 22},
762  {lambda.key, required_argument, &sse_flag, 23},
763  {beta.key, required_argument, &sse_flag, 24},
764  {xi.key, required_argument, &sse_flag, 25},
765  {bhwacc.key, required_argument, &sse_flag, 26},
766  {epsnov.key, required_argument, &sse_flag, 27},
767  {eddfac.key, required_argument, &sse_flag, 28},
768  {gamma.key, required_argument, &sse_flag, 29},
769 #if (defined BSEBBF) || (defined BSEEMP)
770  {sigma.key, required_argument, &sse_flag, 4},
771 #elif MOBSE
772  {sigma1.key, required_argument, &sse_flag, 4},
773  {sigma2.key, required_argument, &sse_flag, 5},
774 #endif
775  //{ifflag.key, required_argument, &sse_flag, 7},
776  {ceflag.key, required_argument, &sse_flag, 6},
777  {tflag.key, required_argument, &sse_flag, 7},
778  {wdflag.key, required_argument, &sse_flag, 8},
779  {bhflag.key, required_argument, &sse_flag, 9},
780  {nsflag.key, required_argument, &sse_flag, 10},
781 #if (defined BSEBBF) || (defined BSEEMP)
782  {psflag.key, required_argument, &sse_flag, 11},
783  {kmech.key, required_argument, &sse_flag, 12},
784  {ecflag.key, required_argument, &sse_flag, 13},
785 #elif MOBSE
786  {piflag.key, required_argument, &sse_flag, 11},
787 #endif
788  {pts1.key, required_argument, &sse_flag, 14},
789  {pts2.key, required_argument, &sse_flag, 15},
790  {pts3.key, required_argument, &sse_flag, 16},
791 #ifdef BSEEMP
792  {trackmode.key, required_argument, &sse_flag, 30},
793 #endif
794  {idum.key, required_argument, &sse_flag, 17},
795  {tscale.key, required_argument, &sse_flag, 18},
796  {rscale.key, required_argument, &sse_flag, 19},
797  {mscale.key, required_argument, &sse_flag, 20},
798  {vscale.key, required_argument, &sse_flag, 21},
799  {z.key, required_argument, 0, 'z'},
800  {"help", no_argument, 0, 'h'},
801  {0,0,0,0}
802  };
803 
804  int opt_used=opt_used_pre;
805  int copt;
806  int option_index;
807  std::string fname_par;
808  optind = 0;
809  while ((copt = getopt_long(argc, argv, "-z:p:h", long_options, &option_index)) != -1)
810  switch (copt) {
811  case 0:
812  switch (sse_flag) {
813  case 0:
814  neta.value = atof(optarg);
815  if(print_flag) neta.print(std::cout);
816  opt_used+=2;
817  break;
818  case 1:
819  bwind.value = atof(optarg);
820  if(print_flag) bwind.print(std::cout);
821  opt_used+=2;
822  break;
823  case 2:
824  hewind.value = atof(optarg);
825  if(print_flag) hewind.print(std::cout);
826  opt_used+=2;
827  break;
828  //case 3:
829  // mxns.value = atof(optarg);
830  // if(print_flag) mxns.print(std::cout);
831  // break;
832 #if (defined BSEBBF) || (defined BSEEMP)
833  case 4:
834  sigma.value = atof(optarg);
835  if(print_flag) sigma.print(std::cout);
836  opt_used+=2;
837  break;
838 #elif MOBSE
839  case 4:
840  sigma1.value = atof(optarg);
841  if(print_flag) sigma1.print(std::cout);
842  opt_used+=2;
843  break;
844  case 5:
845  sigma2.value = atof(optarg);
846  if(print_flag) sigma2.print(std::cout);
847  opt_used+=2;
848  break;
849 #endif
850  case 6:
851  ceflag.value = atof(optarg);
852  if(print_flag) ceflag.print(std::cout);
853  opt_used+=2;
854  break;
855  case 7:
856  tflag.value = atof(optarg);
857  if(print_flag) tflag.print(std::cout);
858  opt_used+=2;
859  break;
860  case 8:
861  wdflag.value = atof(optarg);
862  if(print_flag) wdflag.print(std::cout);
863  opt_used+=2;
864  break;
865  case 9:
866  bhflag.value = atof(optarg);
867  if(print_flag) bhflag.print(std::cout);
868  opt_used+=2;
869  break;
870  case 10:
871  nsflag.value = atof(optarg);
872  if(print_flag) nsflag.print(std::cout);
873  opt_used+=2;
874  break;
875 #if (defined BSEBBF) || (defined BSEEMP)
876  case 11:
877  psflag.value = atof(optarg);
878  if(print_flag) psflag.print(std::cout);
879  opt_used+=2;
880  break;
881  case 12:
882  kmech.value = atof(optarg);
883  if(print_flag) kmech.print(std::cout);
884  opt_used+=2;
885  break;
886  case 13:
887  ecflag.value = atof(optarg);
888  if(print_flag) ecflag.print(std::cout);
889  opt_used+=2;
890  break;
891 #elif MOBSE
892  case 11:
893  piflag.value = atof(optarg);
894  if(print_flag) piflag.print(std::cout);
895  opt_used+=2;
896  break;
897 #endif
898  case 14:
899  pts1.value = atof(optarg);
900  if(print_flag) pts1.print(std::cout);
901  opt_used+=2;
902  break;
903  case 15:
904  pts2.value = atof(optarg);
905  if(print_flag) pts2.print(std::cout);
906  opt_used+=2;
907  break;
908  case 16:
909  pts3.value = atof(optarg);
910  if(print_flag) pts3.print(std::cout);
911  opt_used+=2;
912  break;
913  case 17:
914  idum.value = atof(optarg);
915  if(print_flag) idum.print(std::cout);
916  opt_used+=2;
917  break;
918  case 18:
919  tscale.value = atof(optarg);
920  if(print_flag) tscale.print(std::cout);
921  opt_used+=2;
922  break;
923  case 19:
924  rscale.value = atof(optarg);
925  if(print_flag) rscale.print(std::cout);
926  opt_used+=2;
927  break;
928  case 20:
929  mscale.value = atof(optarg);
930  if(print_flag) mscale.print(std::cout);
931  opt_used+=2;
932  break;
933  case 21:
934  vscale.value = atof(optarg);
935  if(print_flag) vscale.print(std::cout);
936  opt_used+=2;
937  break;
938  case 22:
939  alpha.value = atof(optarg);
940  if(print_flag) alpha.print(std::cout);
941  opt_used+=2;
942  break;
943  case 23:
944  lambda.value = atof(optarg);
945  if(print_flag) lambda.print(std::cout);
946  opt_used+=2;
947  break;
948  case 24:
949  beta.value = atof(optarg);
950  if(print_flag) beta.print(std::cout);
951  opt_used+=2;
952  break;
953  case 25:
954  xi.value = atof(optarg);
955  if(print_flag) xi.print(std::cout);
956  opt_used+=2;
957  break;
958  case 26:
959  bhwacc.value = atof(optarg);
960  if(print_flag) bhwacc.print(std::cout);
961  opt_used+=2;
962  break;
963  case 27:
964  epsnov.value = atof(optarg);
965  if(print_flag) epsnov.print(std::cout);
966  opt_used+=2;
967  break;
968  case 28:
969  eddfac.value = atof(optarg);
970  if(print_flag) eddfac.print(std::cout);
971  opt_used+=2;
972  break;
973  case 29:
974  gamma.value = atof(optarg);
975  if(print_flag) gamma.print(std::cout);
976  opt_used+=2;
977  break;
978 #ifdef BSEEMP
979  case 30:
980  trackmode.value = atoi(optarg);
981  if(print_flag) trackmode.print(std::cout);
982  opt_used+=2;
983  break;
984 #endif
985  default:
986  break;
987  }
988  break;
989  case 'z':
990  z.value = atof(optarg);
991  if(print_flag) z.print(std::cout);
992  opt_used+=2;
993  break;
994  case 'p':
995  fname_par = optarg;
996  if(print_flag) {
997 #ifdef BSEBBF
998  std::string fbse_par = fname_par+".bse";
999 #elif MOBSE
1000  std::string fbse_par = fname_par+".mobse";
1001 #elif BSEEMP
1002  std::string fbse_par = fname_par+".bseEmp";
1003 #endif
1004  FILE* fpar_in;
1005  if( (fpar_in = fopen(fbse_par.c_str(),"r")) == NULL) {
1006  fprintf(stderr,"Error: Cannot open file %s.\n", fbse_par.c_str());
1007  abort();
1008  }
1009  input_par_store.readAscii(fpar_in);
1010  fclose(fpar_in);
1011  }
1012  opt_used+=2;
1013 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1014  input_par_store.mpi_broadcast();
1015  PS::Comm::barrier();
1016 #endif
1017  break;
1018  case 'h':
1019  if(print_flag){
1020 #ifdef BSEBBF
1021  std::cout<<"SSE/BSE options:"<<std::endl;
1022 #elif MOBSE
1023  std::cout<<"MOBSE options:"<<std::endl;
1024 #elif BSEEMP
1025  std::cout<<"BSEEMP options:"<<std::endl;
1026 #endif
1027  input_par_store.printHelp(std::cout, 2, 10, 23);
1028  }
1029  return -1;
1030  case '?':
1031  opt_used +=2;
1032  break;
1033  default:
1034  break;
1035  }
1036 
1037 #ifdef BSEBBF
1038  if(print_flag) std::cout<<"----- Finish reading input options of SSE/BSE -----\n";
1039 #elif MOBSE
1040  if(print_flag) std::cout<<"----- Finish reading input options of MOBSE -----\n";
1041 #elif BSEEMP
1042  if(print_flag) std::cout<<"----- Finish reading input options of BSEEMP -----\n";
1043 #endif
1044 
1045  return opt_used;
1046  }
1047 };
1048 
1050 
1053 class BSEManager{
1054 public:
1055  double z, zpars[20];
1056 #ifdef BSEEMP
1057  int trackmode;
1058 #endif
1059  double tscale;
1060  double rscale;
1061  double mscale;
1062  double vscale;
1063  const double year_to_day;
1064  const char* single_type[16];
1065  const char* binary_type[14];
1066 
1067  BSEManager(): z(0.0), zpars{0},
1068 #ifdef BSEEMP
1069  trackmode(0),
1070 #endif
1071  tscale(0.0), rscale(0.0), mscale(0.0), vscale(0.0), year_to_day(3.6525e8),
1072  single_type{"LMS", "MS", "HG", "GB", "CHeB", "FAGB", "SAGB", "HeMS", "HeHG", "HeGB", "HeWD", "COWD", "ONWD", "NS", "BH", "SN"},
1073  binary_type{"Unset", //0
1074  "Initial", //1
1075  "Type_change", //2
1076  "Start_Roche", //3
1077  "End_Roche", //4
1078  "Contact", //5
1079  "Start_Symbiotic", //6
1080  "End_Symbiotic", //7
1081  "Common_envelope", //8
1082  "Giant", //9
1083  "Coalescence", //10
1084  "Blue_straggler", //11
1085  "No_remain", //12
1086  "Disrupt" //13
1087  } {}
1088 
1089 
1090  bool checkParams() {
1091  assert(z>0.0);
1092 #ifdef BSEEMP
1093  assert(trackmode>0);
1094 #endif
1095  assert(tscale>0.0);
1096  assert(rscale>0.0);
1097  assert(mscale>0.0);
1098  assert(vscale>0.0);
1099  return true;
1100  }
1101 
1103  void dumpRandConstant(const char* _fname) {
1104  FILE* fin;
1105  if( (fin = fopen(_fname,"w")) == NULL) {
1106  fprintf(stderr,"Error: Cannot open file %s.\n", _fname);
1107  abort();
1108  }
1109  fprintf(fin, "%d %d %d ", value3_.idum, rand3_.idum2, rand3_.iy);
1110  for (int i=0; i<32; i++) fprintf(fin, "%d ", rand3_.ir[i]);
1111  fprintf(fin, "\n");
1112  fclose(fin);
1113  }
1114 
1116  void readRandConstant(const char* _fname) {
1117  FILE* fin;
1118  if( (fin = fopen(_fname,"r")) == NULL) {
1119  std::cerr<<_fname<<" not found.\n";
1120  }
1121  else {
1122  int rcount = fscanf(fin, "%d %d %d ", &value3_.idum, &rand3_.idum2, &rand3_.iy);
1123  for (int i=0; i<32; i++) rcount += fscanf(fin, "%d ", &rand3_.ir[i]);
1124  if(rcount<35) {
1125  std::cerr<<"Error: Data reading fails! requiring data number is 35, only obtain "<<rcount<<".\n";
1126  abort();
1127  }
1128  fclose(fin);
1129  }
1130  }
1131 
1132 #ifdef MOBSE
1133  static void printLogo(std::ostream & fout) {
1135  fout<<"---------------------------------------\n"
1136  <<" ╔╦╗╔═╗╔╗ ╔═╗╔═╗\n"
1137  <<" ║║║║ ║╠╩╗╚═╗║╣ \n"
1138  <<" ╩ ╩╚═╝╚═╝╚═╝╚═╝\n"
1139  <<"---------------------------------------"<<std::endl;
1140  fout<<" Online document: https://mobse-webpage.netlify.app/\n"
1141  <<std::endl;
1142  }
1143 #endif
1144 
1146  static void printReference(std::ostream & fout, const int offset=4) {
1147  for (int i=0; i<offset; i++) fout<<" ";
1148  fout<<"SSE: Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543\n";
1149  for (int i=0; i<offset; i++) fout<<" ";
1150  fout<<"BSE: Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897\n";
1151 #ifdef BSEBBF
1152  for (int i=0; i<offset; i++) fout<<" ";
1153  fout<<"Updated BSE: Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2020, A&A, 639, A41"
1154  <<std::endl;
1155 #elif BSEEMP
1156  for (int i=0; i<offset; i++) fout<<" ";
1157  fout<<"BSEEMP: Tanikawa A., Yoshida T., Kinugawa T., Takahashi K., Umeda H., 2020, MNRAS, 495, 4170\n";
1158  for (int i=0; i<offset; i++) fout<<" ";
1159  fout<<" Tanikawa A., Susa H., Yoshida T., Trani A.~A., Kinugawa T., 2021, ApJ, 910, 30\n";
1160 #elif MOBSE
1161  for (int i=0; i<offset; i++) fout<<" ";
1162  fout<<"MOBSE: Giacobbo N., Mapelli M. & Spera M., 2018, MNRAS, 474, 2959\n";
1163  fout<<"\t \t (Online document: https://mobse-webpage.netlify.app/)"
1164  <<std::endl;
1165 #endif
1166  }
1167 
1168  static std::string getSSEOutputFilenameSuffix() {
1169 #ifdef BSEBBF
1170  return std::string(".sse");
1171 #elif BSEEMP
1172  return std::string(".sseEmp");
1173 #elif MOBSE
1174  return std::string(".mosse");
1175 #endif
1176  }
1177 
1178  static std::string getBSEOutputFilenameSuffix() {
1179 #ifdef BSEBBF
1180  return std::string(".bse");
1181 #elif BSEEMP
1182  return std::string(".bseEmp");
1183 #elif MOBSE
1184  return std::string(".mobse");
1185 #endif
1186  }
1187 
1188  static std::string getSSEName() {
1189 #ifdef BSEBBF
1190  return std::string("SSE");
1191 #elif BSEEMP
1192  return std::string("SSEEMP");
1193 #elif MOBSE
1194  return std::string("MOSSE");
1195 #endif
1196  }
1197 
1198  static std::string getBSEName() {
1199 #ifdef BSEBBF
1200  return std::string("BSE");
1201 #elif BSEEMP
1202  return std::string("BSEEMP");
1203 #elif MOBSE
1204  return std::string("MOBSE");
1205 #endif
1206  }
1207 
1208  bool isMassTransfer(const int _binary_type) {
1209  return (_binary_type>=3&&_binary_type<=9);
1210  }
1211 
1213  //bool isKick(const int _binary_type) {
1214  // return (_binary_type==13);
1215  //}
1216 
1217  bool isMerger(const int _binary_type) {
1218  return (_binary_type>=10&&_binary_type<=12);
1219  }
1220 
1221  bool isDisrupt(const int _binary_type) {
1222  return (_binary_type==13);
1223  }
1224 
1226  void initial(const IOParamsBSE& _input, const bool _print_flag=false) {
1227  // common block
1228  value1_.neta = _input.neta.value;
1229  value1_.bwind = _input.bwind.value;
1230  value1_.hewind= _input.hewind.value;
1231 
1232  value2_.alpha = _input.alpha.value;
1233  value2_.lambda = _input.lambda.value;
1234 
1235 #if (defined BSEBBF) || (defined BSEEMP)
1236  value4_.sigma = _input.sigma.value;
1237 #elif MOBSE
1238  value4_.sigma1 = _input.sigma1.value;
1239  value4_.sigma2 = _input.sigma2.value;
1240 #endif
1241  value4_.mxns = 1.8;
1242 #if (defined BSEBBF) || (defined BSEEMP)
1243  if (_input.nsflag.value>0) value4_.mxns = 2.5;
1244 #elif MOBSE
1245  if (_input.nsflag.value>0) value4_.mxns = 3.0;
1246 #endif
1247  value4_.bhflag = _input.bhflag.value;
1248 
1249  value5_.beta = _input.beta.value;
1250  value5_.xi = _input.xi.value;
1251  value5_.bhwacc = _input.bhwacc.value;
1252  value5_.epsnov = _input.epsnov.value;
1253  value5_.eddfac = _input.eddfac.value;
1254  value5_.gamma = _input.gamma.value;
1255 
1256  flags_.ceflag = _input.ceflag.value;
1257  flags_.tflag = _input.tflag.value;
1258  //flags_.ifflag = _input.ifflag.value;
1259  flags_.wdflag = _input.wdflag.value;
1260  flags_.nsflag = _input.nsflag.value;
1261 #if (defined BSEBBF) || (defined BSEEMP)
1262  flags2_.psflag = _input.psflag.value;
1263  flags2_.kmech = _input.kmech.value;
1264  flags2_.ecflag = _input.ecflag.value;
1265 #elif MOBSE
1266  flags_.piflag = _input.piflag.value;
1267 #endif
1268 
1269  points_.pts1 = _input.pts1.value;
1270  points_.pts2 = _input.pts2.value;
1271  points_.pts3 = _input.pts3.value;
1272 
1273  tscale = _input.tscale.value;
1274  rscale = _input.rscale.value;
1275  mscale = _input.mscale.value;
1276  vscale = _input.vscale.value;
1277 
1278  // Set parameters which depend on the metallicity
1279  z = _input.z.value;
1280 #ifdef BSEEMP
1281  if (_print_flag&&(z>0.03))
1282  std::cerr<<"BSE warning! metallicity Z is not in (0.0, 0.03); given value:"<<z<<std::endl;
1283  trackmode = _input.trackmode.value;
1284  zcnsts_(&z, zpars, &trackmode);
1285 #else
1286  if (_print_flag&&(z<0.0001||z>0.03))
1287  std::cerr<<"BSE warning! metallicity Z is not in (0.0001, 0.03); given value:"<<z<<std::endl;
1288  zcnsts_(&z, zpars);
1289 #endif
1290  value3_.idum = (_input.idum.value>0)? -_input.idum.value: _input.idum.value;
1291 
1292 #ifdef PARTICLE_SIMULATOR_MPI_PARALLEL
1293  // add off set to random seed to avoid repeating random numbers
1294  value3_.idum += PS::Comm::getRank();
1295 #endif
1296 
1297  // collision matrix
1298  instar_();
1299 
1300  if (_print_flag) {
1301  printconst_();
1302  std::cout<<"z: "<<z<<" zpars: ";
1303  for (int i=0;i<20;i++) std::cout<<zpars[i]<<" ";
1304  std::cout<<std::endl;
1305 #ifdef BSEEMP
1306  std::cout<<"EMPTrack: "<<trackmode<<std::endl;
1307 #endif
1308  }
1309 
1310  }
1311 
1313  double getMass(StarParameter& _star) {
1314  return _star.mt/mscale;
1315  }
1316 
1318  double getMassLoss(StarParameterOut& _out) {
1319  return _out.dm/mscale;
1320  }
1321 
1323  double getMergerRadius(StarParameter& _star) {
1324  // use stellar radius as merger radius
1325  return _star.r/rscale;
1326  }
1327 
1329  double getStellarRadius(StarParameter& _star) {
1330  return _star.r/rscale;
1331  }
1332 
1334 
1336  double getSpeedOfLight() const {
1337  const double c = 2.99792458e5;
1338  return c/vscale;
1339  }
1340 
1342  double getTime(StarParameter& _star) {
1343  return _star.tphys/tscale;
1344  }
1345 
1347  double getDTMiss(StarParameterOut& _out) {
1348  return _out.dtmiss/tscale;
1349  }
1350 
1352  void printTypeChange(std::ostream& _fout, StarParameter& _star, StarParameterOut& _out, const int _width=4) {
1353  assert(_out.kw0>=0&&_out.kw0<16);
1354  assert(_star.kw>=0&&_star.kw<16);
1355  _fout<<" "<<std::setw(_width)<<single_type[_out.kw0]<<" -> "<<std::setw(_width)<<single_type[_star.kw];
1356  }
1357 
1359  void printBinaryEvent(std::ostream& _fout, const BinaryEvent& _bin_event) {
1360  int nmax = _bin_event.getEventNMax();
1361  for (int k=0; k<nmax; k++) {
1362  int type = _bin_event.getType(k);
1363  if(type>0) {
1364  _fout<<std::setw(16)<<binary_type[type];
1365  _bin_event.print(_fout, k);
1366  _fout<<std::endl;
1367  }
1368  else if(type<0) break;
1369  }
1370  }
1371 
1373  void printBinaryEventOne(std::ostream& _fout, const BinaryEvent& _bin_event, const int k) {
1374  int type = _bin_event.getType(k);
1375  assert(type>=0&&type<14);
1376  _fout<<std::setw(16)<<binary_type[type]<<" Init: ";
1377  if (k==0) _bin_event.print(_fout, _bin_event.getEventIndexInit());
1378  else _bin_event.print(_fout, k-1);
1379  _fout<<"\n"<<std::setw(16)<<" "<<" Final: ";
1380  _bin_event.print(_fout, k);
1381  }
1382 
1384  void printBinaryEventColumnOne(std::ostream& _fout, const BinaryEvent& _bin_event, const int k, const int _width=20, const bool print_type_name=true) {
1385  int type = _bin_event.getType(k);
1386  assert(type>=0&&type<14);
1387  if (print_type_name) _fout<<std::setw(16)<<binary_type[type];
1388  _fout<<std::setw(_width)<<type;
1389  if (k==0) _bin_event.printColumn(_fout, _bin_event.getEventIndexInit(), _width);
1390  else _bin_event.printColumn(_fout, k-1, _width);
1391  _bin_event.printColumn(_fout, k, _width);
1392  }
1393 
1395 
1399  double getVelocityChange(double* dv, StarParameterOut& _out) {
1400  for (int k=0; k<3; k++) dv[k] = _out.vkick[k]/vscale;
1401  return _out.vkick[3]/vscale;
1402  }
1403 
1405 
1412  int evolveStar(StarParameter& _star, StarParameterOut& _out, const double _dt, bool _unit_in_myr=false) {
1413  double tphysf = _dt*tscale + _star.tphys;
1414  if (_unit_in_myr) tphysf = _dt + _star.tphys;
1415  double dtp=tphysf*100.0+1000.0;
1416  _out.dm = _star.mt;
1417  _out.kw0 = _star.kw;
1418  int kw = _star.kw;
1419  evolv1_(&kw, &_star.m0, &_star.mt, &_star.r,
1420  &_star.lum, &_star.mc, &_star.rc, &_out.menv, &_out.renv,
1421  &_star.ospin, &_star.epoch,
1422  &_out.tm, &_star.tphys, &tphysf, &dtp, &z, zpars, _out.vkick);
1423  _star.kw = kw;
1424  _out.dm = _star.mt - _out.dm;
1425  _out.dtmiss = tphysf - _star.tphys;
1426 
1427  if (_star.kw<0) {
1428  //_star.kw = -_star.kw;
1429  return -1; // error
1430  }
1431  else if (_out.vkick[3]>0) return 2; // kick
1432  else if (_star.kw!=_out.kw0) return 1; // kw change
1433  else return 0;
1434  }
1435 
1437 
1450  int evolveBinary(StarParameter& _star1, StarParameter& _star2, StarParameterOut& _out1, StarParameterOut& _out2,
1451  double& _semi, double& _period, double& _ecc, BinaryEvent& _bse_event, const int& _binary_init_type, const double _dt_nb) {
1452  double tphys = std::max(_star1.tphys, _star2.tphys);
1453  double tphysf = _dt_nb*tscale + tphys;
1454  double dtp=tphysf*100.0+1000.0;
1455  double period_days = _period*tscale*year_to_day;
1456  double semi_rsun = _semi*rscale;
1457  // in case two component have different tphys, evolve to the same time first
1458  int event_flag = 0 ;
1459  _out1.dm = _star1.mt;
1460  _out2.dm = _star2.mt;
1461 
1462  // backup initial state
1463  _bse_event.recordEvent(_star1, _star2, semi_rsun, _ecc, _binary_init_type, _bse_event.getEventIndexInit());
1464 
1465  if (_star1.kw==15 || _star2.kw==15) {
1466  if (_star1.kw==15 && _star2.kw==15) {
1467  _star1.tphys = tphysf;
1468  _star2.tphys = tphysf;
1469  return 0;
1470  }
1471  if (_star1.kw==15) {
1472  double dt = tphysf - _star2.tphys;
1473  event_flag = evolveStar(_star2, _out2, dt, true);
1474  _star1.tphys = tphysf;
1475  }
1476  if (_star2.kw==15) {
1477  double dt = tphysf - _star1.tphys;
1478  event_flag = evolveStar(_star1, _out1, dt, true);
1479  _star2.tphys = tphysf;
1480  }
1481  if (event_flag<0) return event_flag;
1482  int event_index = 0;
1483  if (event_flag>0) _bse_event.recordEvent(_star1, _star2, semi_rsun, _ecc, 2, event_index++);
1484  _bse_event.setEventIndexEnd(event_index);
1485  return 0;
1486  }
1487 
1488  if (_star1.tphys<tphys) {
1489  double dt = tphys - _star1.tphys;
1490  event_flag = evolveStar(_star1, _out1, dt, true);
1491  }
1492  if (event_flag<0) return event_flag;
1493  if (_star2.tphys<tphys) {
1494  double dt = tphys - _star2.tphys;
1495  event_flag = evolveStar(_star2, _out2, dt, true);
1496  }
1497  if (event_flag<0) return event_flag;
1498 
1499  int kw[2];
1500  double m0[2],mt[2],r[2],lum[2],mc[2],rc[2],menv[2],renv[2],ospin[2],epoch[2],tm[2],vkick[8];
1501  for (int k =0; k<4; k++) {
1502  vkick[k] = _out1.vkick[k];
1503  vkick[k+4]= _out2.vkick[k];
1504  }
1505 
1506 
1507  kw[0] = _star1.kw;
1508  m0[0] = _star1.m0;
1509  mt[0] = _star1.mt;
1510  r[0] = _star1.r;
1511  mc[0] = _star1.mc;
1512  rc[0] = _star1.rc;
1513  ospin[0] = _star1.ospin;
1514  epoch[0] = _star1.epoch;
1515 
1516  kw[1] = _star2.kw;
1517  m0[1] = _star2.m0;
1518  mt[1] = _star2.mt;
1519  r[1] = _star2.r;
1520  mc[1] = _star2.mc;
1521  rc[1] = _star2.rc;
1522  ospin[1] = _star2.ospin;
1523  epoch[1] = _star2.epoch;
1524 
1525  evolv2_(kw, m0, mt, r, lum, mc, rc, menv, renv, ospin, epoch, tm, &tphys, &tphysf, &dtp, &z, zpars, &period_days, &_ecc, _bse_event.record[0], vkick);
1526  _period = period_days/year_to_day/tscale;
1527 
1528  _star1.kw = kw[0];
1529  _star1.m0 = m0[0];
1530  _star1.mt = mt[0];
1531  _star1.r = r[0];
1532  _star1.mc = mc[0];
1533  _star1.rc = rc[0];
1534  _star1.ospin = ospin[0];
1535  _star1.epoch = epoch[0];
1536  _star1.tphys = tphys;
1537  _star1.lum = lum[0];
1538 
1539  _star2.kw = kw[1];
1540  _star2.m0 = m0[1];
1541  _star2.mt = mt[1];
1542  _star2.r = r[1];
1543  _star2.mc = mc[1];
1544  _star2.rc = rc[1];
1545  _star2.ospin = ospin[1];
1546  _star2.epoch = epoch[1];
1547  _star2.tphys = tphys;
1548  _star2.lum = lum[1];
1549 
1550  _out1.menv = menv[0];
1551  _out1.renv = renv[0];
1552  _out1.tm = tm[0];
1553  _out1.dm = _star1.mt - _out1.dm;
1554  _out1.dtmiss = tphysf - _star1.tphys;
1555 
1556  _out2.menv = menv[1];
1557  _out2.renv = renv[1];
1558  _out2.tm = tm[1];
1559  _out2.dm = _star2.mt - _out2.dm;
1560  _out2.dtmiss = tphysf - _star2.tphys;
1561 
1562  for (int k=0; k<4; k++) _out1.vkick[k]=vkick[k];
1563  for (int k=0; k<4; k++) _out2.vkick[k]=vkick[k+4];
1564 
1565  if (kw[0]<0||kw[1]<0||(_star1.mt<0&&_star1.kw==15)||(_star2.mt<0&&_star2.kw==15)) {
1566  kw[0] = abs(kw[0]);
1567  kw[1] = abs(kw[1]);
1568  return -1; // error case
1569  }
1570  //else if (vkick[3]>0||vkick[7]>0) return 3; // kick
1571  //else if (isDisrupt(_binary_type)) return 4; // disrupt without kick
1572  //else if (isMerger(_binary_type)) return 5; // Merger
1573  //else if (isMassTransfer(_binary_type)) return 2; // orbit change
1574  //else if (_binary_type>0) return 1; // type change
1575  else return 0;
1576  }
1577 
1579  /* Use Roche Radius estimation
1580  @param[in] _star1: star parameter of first
1581  @param[in] _star2: star parameter of second
1582  @param[in] _semi: semi-major axis, [IN unit]
1583  @param[in] _ecc: eccentricity of binary, used for BSE
1584  \return true: Roche fill
1585  */
1586  bool isRocheFill(StarParameter& _star1, StarParameter& _star2, double& _semi, double& _ecc) {
1587  assert(_star2.mt>0);
1588  double q = _star1.mt/_star2.mt;
1589  double rl_over_semi1 = EstimateRocheRadiusOverSemi(q);
1590  q = 1.0/q;
1591  double rl_over_semi2 = EstimateRocheRadiusOverSemi(q);
1592  double rl_over_semi_max =std::max(rl_over_semi1,rl_over_semi2);
1593  double rad = _star1.r + _star2.r;
1594  double semi_rsun = _semi*rscale;
1595  double rcrit = std::max(rad,rl_over_semi_max*semi_rsun);
1596  double peri_rsun = semi_rsun*(1-_ecc);
1597  bool is_fill= (rcrit>0.5*peri_rsun);
1598  return is_fill;
1599  }
1600 
1601 
1603 
1612  void merge(StarParameter& _star1, StarParameter& _star2, StarParameterOut& _out1, StarParameterOut& _out2, double& _semi, double& _ecc) {
1613 
1614  double m0[2],mt[2],r[2],mc[2],rc[2],menv[2],renv[2],ospin[2],age[2],vkick[8];
1615  int kw[2];
1616 
1617  for (int k =0; k<4; k++) {
1618  vkick[k] = _out1.vkick[k];
1619  vkick[k+4]= _out2.vkick[k];
1620  }
1621 
1622  kw[0] = _star1.kw;
1623  m0[0] = _star1.m0;
1624  mt[0] = _star1.mt;
1625  r[0] = _star1.r;
1626  mc[0] = _star1.mc;
1627  rc[0] = _star1.rc;
1628  ospin[0] = _star1.ospin;
1629  age[0]= _star1.tphys-_star1.epoch;
1630 
1631  kw[1] = _star2.kw;
1632  m0[1] = _star2.m0;
1633  mt[1] = _star2.mt;
1634  r[1] = _star2.r;
1635  mc[1] = _star2.mc;
1636  rc[1] = _star2.rc;
1637  ospin[1] = _star2.ospin;
1638  age[1]= _star2.tphys-_star2.epoch;
1639 
1640  _out1.dm = _star1.mt;
1641  _out2.dm = _star2.mt;
1642 
1643  double semi_rsun = _semi*rscale;
1644 
1645  merge_(kw, m0, mt, r, mc, rc, menv, renv, ospin, age, &semi_rsun, &_ecc, vkick, zpars);
1646 
1647  _semi = semi_rsun/rscale;
1648 
1649  _star1.kw = kw[0];
1650  _star1.m0 = m0[0];
1651  _star1.mt = mt[0];
1652  _star1.r = r[0];
1653  _star1.mc = mc[0];
1654  _star1.rc = rc[0];
1655  _star1.ospin = ospin[0];
1656  _star1.epoch = _star1.tphys - age[0];
1657 
1658  _star2.kw = kw[1];
1659  _star2.m0 = m0[1];
1660  _star2.mt = mt[1];
1661  _star2.r = r[1];
1662  _star2.mc = mc[1];
1663  _star2.rc = rc[1];
1664  _star2.ospin = ospin[1];
1665  _star2.epoch = _star2.tphys - age[1];
1666 
1667  _out1.menv = menv[0];
1668  _out1.renv = renv[0];
1669  _out1.tm = 0;
1670  _out1.dm = _star1.mt - _out1.dm;
1671  _out1.dtmiss = 0;
1672 
1673  _out2.menv = menv[1];
1674  _out2.renv = renv[1];
1675  _out2.tm = 0;
1676  _out2.dm = _star2.mt - _out2.dm;
1677  _out2.dtmiss = 0;
1678 
1679  for (int k=0; k<4; k++) _out1.vkick[k]=vkick[k];
1680  for (int k=0; k<4; k++) _out2.vkick[k]=vkick[k+4];
1681 
1682  }
1683 
1685 
1689  double getTimeStepStar(StarParameter& _star) {
1690  if (_star.kw==15) return 1.0e30/tscale; // give very large value to avoid evolve
1691 
1692  // make a copy to avoid overwriting the origin values
1693  int kw = _star.kw;
1694  double m0 = _star.m0;
1695  double mt = _star.mt;
1696  double r = _star.r;
1697  double mc = _star.mc;
1698  double rc = _star.rc;
1699  double age = _star.tphys-_star.epoch;
1700  double dtm, dtr;
1701  // get next step
1702  trdot_(&kw, &m0, &mt, &r, &mc, &rc, &age, &dtm, &dtr, zpars);
1703 
1704  //assert(dtr>0.0);
1705  //assert(dtm>0.0);
1706 
1707  return std::min(dtr, dtm)/tscale;
1708  }
1709 
1711 
1719  double getTimeStepBinary(StarParameter& _star1, StarParameter& _star2,
1720  double& _semi, double& _ecc, int &_binary_type) {
1721 
1722  double dt = getTimeStepStar(_star1);
1723  dt = std::min(dt,getTimeStepStar(_star2));
1724 
1725  // mass transfer case
1726  if (_star1.kw>=10&&_star1.kw<15&&_star2.kw>=10&&_star2.kw<15&&_semi>0.0) {// GR effect
1727  double semi_rsun = _semi*rscale;
1728  dt = std::min(dt, EstimateGRTimescale(_star1, _star2, semi_rsun, _ecc));
1729  }
1730  else if (isMassTransfer(_binary_type)) {
1731  int kw[2];
1732  double m0[2],mt[2],r[2],mc[2],rc[2],age[2];
1733  kw[0] = _star1.kw;
1734  m0[0] = _star1.m0;
1735  mt[0] = _star1.mt;
1736  r[0] = _star1.r;
1737  mc[0] = _star1.mc;
1738  rc[0] = _star1.rc;
1739  age[0]= _star1.tphys-_star1.epoch;
1740 
1741  kw[1] = _star2.kw;
1742  m0[1] = _star2.m0;
1743  mt[1] = _star2.mt;
1744  r[1] = _star2.r;
1745  mc[1] = _star2.mc;
1746  rc[1] = _star2.rc;
1747  age[1]= _star2.tphys-_star2.epoch;
1748 
1749  double dtr;
1750  double semi_rsun = _semi*rscale;
1751  trflow_(kw,m0,mt,r,mc,rc,age,&dtr,&semi_rsun,&_ecc,zpars);
1752  dt = std::min(dt,dtr);
1753 
1754  }
1755 
1756  return dt/tscale;
1757  }
1758 
1760 
1772  bool isCallBSENeeded(StarParameter& _star1, StarParameter& _star2, double& _semi, double& _ecc, double& _dt) {
1773  if (_star1.kw>=10&&_star1.kw<15&&_star2.kw>=10&&_star2.kw<15) {
1774  // check GR effect
1775  double semi_rsun = _semi*rscale;
1776  double dt_gr = EstimateGRTimescale(_star1, _star2, semi_rsun, _ecc);
1777  if (dt_gr<_dt*tscale) return true;
1778  }
1779  else {
1780  // check Roche overflow
1781  if (isRocheFill(_star1, _star2, _semi, _ecc)) {
1782  int kw[2];
1783  double m0[2],mt[2],r[2],mc[2],rc[2],age[2];
1784  kw[0] = _star1.kw;
1785  m0[0] = _star1.m0;
1786  mt[0] = _star1.mt;
1787  r[0] = _star1.r;
1788  mc[0] = _star1.mc;
1789  rc[0] = _star1.rc;
1790  age[0]= _star1.tphys-_star1.epoch;
1791 
1792  kw[1] = _star2.kw;
1793  m0[1] = _star2.m0;
1794  mt[1] = _star2.mt;
1795  r[1] = _star2.r;
1796  mc[1] = _star2.mc;
1797  rc[1] = _star2.rc;
1798  age[1]= _star2.tphys-_star2.epoch;
1799 
1800  double dtr;
1801  double semi_rsun = _semi*rscale;
1802  trflow_(kw,m0,mt,r,mc,rc,age,&dtr,&semi_rsun,&_ecc,zpars);
1803  if (dtr<_dt*tscale) return true;
1804  }
1805  }
1806  return false;
1807  }
1808 };
StarParameter::mt
double mt
Initial stellar mass in solar units
Definition: bse_interface.h:260
StarParameter::lum
double lum
physical evolve time in Myr
Definition: bse_interface.h:267
BinaryEvent::printColumn
void printColumn(std::ostream &_fout, const int _index, const int _width=20) const
print data of class members using column style
Definition: bse_interface.h:606
StarParameterOut::dm
double dm
kick velocity for NS/BH formation
Definition: bse_interface.h:400
IOParamsBSE::nsflag
IOParams< long long int > nsflag
Definition: bse_interface.h:645
BinaryEvent::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: bse_interface.h:576
IOParams::value
Type value
Definition: io.hpp:43
BSEManager::mscale
double mscale
radius scaling factor from NB to Rsun
Definition: bse_interface.h:1061
BSEManager::printReference
static void printReference(std::ostream &fout, const int offset=4)
print reference to cite
Definition: bse_interface.h:1146
BSEManager::merge
void merge(StarParameter &_star1, StarParameter &_star2, StarParameterOut &_out1, StarParameterOut &_out2, double &_semi, double &_ecc)
merge two star using mix function, star 2 will becomes zero mass
Definition: bse_interface.h:1612
BSEManager::getBSEOutputFilenameSuffix
static std::string getBSEOutputFilenameSuffix()
Definition: bse_interface.h:1178
BSEManager::year_to_day
const double year_to_day
velocity scaling factor from NB to km/s
Definition: bse_interface.h:1063
BSEManager::rscale
double rscale
time scaling factor from NB to Myr (t[Myr]=t[NB]*tscale)
Definition: bse_interface.h:1060
IOParamsBSE
IO parameters manager for BSE based code.
Definition: bse_interface.h:619
IOParamsBSE::alpha
IOParams< double > alpha
Definition: bse_interface.h:625
StarParameterOut::vkick
double vkick[4]
Main sequence lifetime
Definition: bse_interface.h:399
IOParamsBSE::idum
IOParams< long long int > idum
Definition: bse_interface.h:659
galpy_pot_movie.fp
fp
Definition: galpy_pot_movie.py:131
BSEManager::evolveBinary
int evolveBinary(StarParameter &_star1, StarParameter &_star2, StarParameterOut &_out1, StarParameterOut &_out2, double &_semi, double &_period, double &_ecc, BinaryEvent &_bse_event, const int &_binary_init_type, const double _dt_nb)
call evolv2 for a binary
Definition: bse_interface.h:1450
BSEManager::printBinaryEventColumnOne
void printBinaryEventColumnOne(std::ostream &_fout, const BinaryEvent &_bin_event, const int k, const int _width=20, const bool print_type_name=true)
print binary event one in column
Definition: bse_interface.h:1384
IOParamsContainer::readAscii
void readAscii(FILE *_fin)
Definition: io.hpp:112
IOParamsBSE::input_par_store
IOParamsContainer input_par_store
Definition: bse_interface.h:621
BSEManager::BSEManager
BSEManager()
name of binary type return from BSE evolv2, notice if it is -1, it indicate the end of record
Definition: bse_interface.h:1067
BSEManager::isDisrupt
bool isDisrupt(const int _binary_type)
Definition: bse_interface.h:1221
BSEManager::printBinaryEvent
void printBinaryEvent(std::ostream &_fout, const BinaryEvent &_bin_event)
print binary event
Definition: bse_interface.h:1359
IOParamsContainer::printHelp
void printHelp(std::ostream &os, const int _offset_short_key=2, const int _offset_long_key=1, const int _width_key=23) const
Definition: io.hpp:198
BSEManager::isMassTransfer
bool isMassTransfer(const int _binary_type)
Definition: bse_interface.h:1208
galpy_pot_movie.dt
dt
Definition: galpy_pot_movie.py:134
BSEManager::getSSEName
static std::string getSSEName()
Definition: bse_interface.h:1188
IOParamsBSE::pts2
IOParams< double > pts2
Definition: bse_interface.h:654
StarParameter::printColumn
void printColumn(std::ostream &_fout, const int _width=20) const
print data of class members using column style
Definition: bse_interface.h:345
StarParameter::readAscii
void readAscii(FILE *fp)
read class data with ASCII format
Definition: bse_interface.h:299
IOParamsBSE::vscale
IOParams< double > vscale
Definition: bse_interface.h:663
BinaryEvent::print
void print(std::ostream &fout, const int index) const
for print in one line
Definition: bse_interface.h:549
StarParameterOut::print
void print(std::ostream &fout) const
for print in one line
Definition: bse_interface.h:439
BinaryEvent
BSE based code event recorder class.
Definition: bse_interface.h:493
StarParameterOut::StarParameterOut
StarParameterOut()
mass loss
Definition: bse_interface.h:402
IOParams::print
void print(std::ostream &os) const
Definition: io.hpp:54
BSEManager::getStellarRadius
double getStellarRadius(StarParameter &_star)
get stellar radius in NB unit
Definition: bse_interface.h:1329
BSEManager::getMass
double getMass(StarParameter &_star)
get current mass in NB unit
Definition: bse_interface.h:1313
BSEManager::getBSEName
static std::string getBSEName()
Definition: bse_interface.h:1198
IOParamsBSE::bwind
IOParams< double > bwind
Definition: bse_interface.h:623
BSEManager::getMassLoss
double getMassLoss(StarParameterOut &_out)
get mass loss in NB unit
Definition: bse_interface.h:1318
IOParamsBSE::eddfac
IOParams< double > eddfac
Definition: bse_interface.h:631
StarParameter::r
double r
Current mass in solar units (used for R)
Definition: bse_interface.h:261
StarParameterOut::printColumn
void printColumn(std::ostream &_fout, const int _width=20) const
print data of class members using column style
Definition: bse_interface.h:426
IOParams::key
const char * key
Definition: io.hpp:44
StarParameterOut::tm
double tm
radius of convective envelope
Definition: bse_interface.h:398
IOParamsBSE::beta
IOParams< double > beta
Definition: bse_interface.h:627
StarParameter::epoch
double epoch
spin of star
Definition: bse_interface.h:265
IOParamsBSE::epsnov
IOParams< double > epsnov
Definition: bse_interface.h:630
IOParamsBSE::bhflag
IOParams< long long int > bhflag
Definition: bse_interface.h:644
BSEManager::z
double z
Definition: bse_interface.h:1055
BSEManager::printBinaryEventOne
void printBinaryEventOne(std::ostream &_fout, const BinaryEvent &_bin_event, const int k)
print binary event one
Definition: bse_interface.h:1373
IOParamsBSE::ceflag
IOParams< long long int > ceflag
Definition: bse_interface.h:640
StarParameter::mc
double mc
Stellar radius in solar units
Definition: bse_interface.h:262
PIKG::min
T min(const Vector3< T > &v)
Definition: pikg_vector.hpp:150
BinaryEvent::getEventNMax
int getEventNMax() const
Maximum event number that can be recored in one call of evolv2.
Definition: bse_interface.h:534
BSEManager::getTime
double getTime(StarParameter &_star)
get evolved Time in NB unit
Definition: bse_interface.h:1342
BSEManager::checkParams
bool checkParams()
Definition: bse_interface.h:1090
BSEManager::getDTMiss
double getDTMiss(StarParameterOut &_out)
get the difference of required finishing time and actually evolved time in NB unit
Definition: bse_interface.h:1347
StarParameterOut
SSE/BSE based code star parameter for output.
Definition: bse_interface.h:393
BSEManager::readRandConstant
void readRandConstant(const char *_fname)
read BSE rand constant from file
Definition: bse_interface.h:1116
BSEManager::printTypeChange
void printTypeChange(std::ostream &_fout, StarParameter &_star, StarParameterOut &_out, const int _width=4)
print type change
Definition: bse_interface.h:1352
BSEManager::initial
void initial(const IOParamsBSE &_input, const bool _print_flag=false)
initial SSE/BSE based code global parameters
Definition: bse_interface.h:1226
BSEManager::getMergerRadius
double getMergerRadius(StarParameter &_star)
get merger radius in NB unit
Definition: bse_interface.h:1323
BSEManager::isRocheFill
bool isRocheFill(StarParameter &_star1, StarParameter &_star2, double &_semi, double &_ecc)
check Roche fill condition
Definition: bse_interface.h:1586
BSEManager::isMerger
bool isMerger(const int _binary_type)
notice kick priority is higher than others
Definition: bse_interface.h:1217
IOParamsBSE::pts1
IOParams< double > pts1
Definition: bse_interface.h:653
BSEManager::isCallBSENeeded
bool isCallBSENeeded(StarParameter &_star1, StarParameter &_star2, double &_semi, double &_ecc, double &_dt)
a simple check to determine whether call BSE is necessary by given a reference of time step
Definition: bse_interface.h:1772
IOParamsBSE::tflag
IOParams< long long int > tflag
Definition: bse_interface.h:641
StarParameter::initial
void initial(double _mass, int _kw=1, double _ospin=0.0, double _epoch=0.0)
Landmark luminosities
Definition: bse_interface.h:276
StarParameter::m0
double m0
stellar type
Definition: bse_interface.h:259
StarParameter::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: bse_interface.h:327
IOParamsBSE::rscale
IOParams< double > rscale
Definition: bse_interface.h:661
BSEManager::evolveStar
int evolveStar(StarParameter &_star, StarParameterOut &_out, const double _dt, bool _unit_in_myr=false)
call SSE evolv1 for single star
Definition: bse_interface.h:1412
IOParams< double >
BSEManager::tscale
double tscale
metallicity parameters
Definition: bse_interface.h:1059
StarParameter::kw
long long int kw
Definition: bse_interface.h:258
BSEManager::binary_type
const char * binary_type[14]
name of single type from SSE
Definition: bse_interface.h:1065
IOParamsBSE::bhwacc
IOParams< double > bhwacc
Definition: bse_interface.h:629
BSEManager
SSE/BSE interface manager.
Definition: bse_interface.h:1053
IOParamsBSE::neta
IOParams< double > neta
Definition: bse_interface.h:622
IOParamsBSE::gamma
IOParams< double > gamma
Definition: bse_interface.h:632
PIKG::max
T max(const Vector3< T > &v)
Definition: pikg_vector.hpp:143
IOParamsBSE::wdflag
IOParams< long long int > wdflag
Definition: bse_interface.h:643
BSEManager::getVelocityChange
double getVelocityChange(double *dv, StarParameterOut &_out)
get velocity change in NB unit
Definition: bse_interface.h:1399
BinaryEvent::getEventIndexInit
int getEventIndexInit() const
The index of initial event in record array (last one)
Definition: bse_interface.h:539
BinaryEvent::record
double record[20][9]
Definition: bse_interface.h:496
IOParamsBSE::hewind
IOParams< double > hewind
Definition: bse_interface.h:624
BSEManager::zpars
double zpars[20]
Definition: bse_interface.h:1055
StarParameter::writeAscii
void writeAscii(FILE *fp) const
write class data with ASCII format
Definition: bse_interface.h:291
IOParamsBSE::pts3
IOParams< double > pts3
Definition: bse_interface.h:655
BSEManager::vscale
double vscale
mass scaling factor from NB to Msun
Definition: bse_interface.h:1062
BinaryEvent::recordEvent
void recordEvent(const StarParameter &_p1, const StarParameter &_p2, const double _semi, const double _ecc, const int _type, const int _index)
set up the initial parameter of binary event based on the present status of a binary
Definition: bse_interface.h:501
IOParamsBSE::mscale
IOParams< double > mscale
Definition: bse_interface.h:662
StarParameter::tphys
double tphys
starting time of one evolution phase, age = tphys - epoch
Definition: bse_interface.h:266
BSEManager::getSSEOutputFilenameSuffix
static std::string getSSEOutputFilenameSuffix()
Definition: bse_interface.h:1168
IOParamsBSE::xi
IOParams< double > xi
Definition: bse_interface.h:628
BSEManager::getSpeedOfLight
double getSpeedOfLight() const
get speed of light in NB unit
Definition: bse_interface.h:1336
StarParameter::print
void print(std::ostream &fout) const
for print in one line
Definition: bse_interface.h:309
BSEManager::single_type
const char * single_type[16]
year to day
Definition: bse_interface.h:1064
StarParameterOut::dtmiss
double dtmiss
original type before evolution
Definition: bse_interface.h:395
IOParamsContainer
IO Params container.
Definition: io.hpp:82
BSEManager::getTimeStepBinary
double getTimeStepBinary(StarParameter &_star1, StarParameter &_star2, double &_semi, double &_ecc, int &_binary_type)
call BSE evolv2 for a binary
Definition: bse_interface.h:1719
StarParameterOut::kw0
long long int kw0
Definition: bse_interface.h:394
StarParameterOut::menv
double menv
required evolution time - actually evolved time
Definition: bse_interface.h:396
IOParamsBSE::tscale
IOParams< double > tscale
Definition: bse_interface.h:660
StarParameterOut::renv
double renv
mass of convective envelope
Definition: bse_interface.h:397
IOParamsBSE::read
int read(int argc, char *argv[], const int opt_used_pre=0)
reading parameters from GNU option API
Definition: bse_interface.h:754
BSEManager::dumpRandConstant
void dumpRandConstant(const char *_fname)
dump rand constant to file
Definition: bse_interface.h:1103
StarParameter::ospin
double ospin
core radius in solar units (output)
Definition: bse_interface.h:264
IOParamsBSE::z
IOParams< double > z
Definition: bse_interface.h:664
StarParameterOut::printColumnTitle
static void printColumnTitle(std::ostream &_fout, const int _width=20)
print titles of class members using column style
Definition: bse_interface.h:409
BSEManager::getTimeStepStar
double getTimeStepStar(StarParameter &_star)
get next time step to check in Myr
Definition: bse_interface.h:1689
StarParameter::printTitleWithMeaning
static int printTitleWithMeaning(std::ostream &_fout, const int _counter=0, const int _offset=0)
print column title with meaning (each line for one column)
Definition: bse_interface.h:364
StarParameter::rc
double rc
core mass in solar units
Definition: bse_interface.h:263
BinaryEvent::setEventIndexEnd
void setEventIndexEnd(const int index)
set binary type to -1 for the given event index to indicate the end of record
Definition: bse_interface.h:529
IOParamsBSE::lambda
IOParams< double > lambda
Definition: bse_interface.h:626
StarParameter
SSE/BSE based code star parameter for saving.
Definition: bse_interface.h:257
IOParamsBSE::print_flag
bool print_flag
Definition: bse_interface.h:666
BinaryEvent::getType
int getType(const int index) const
Get the binary type.
Definition: bse_interface.h:544