Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/lwang-astro/mcluster_nb6pp
Browse files Browse the repository at this point in the history
  • Loading branch information
@lwang-astro
lwang-astro committed Mar 14, 2018
2 parents aa5a8e7 + db01ec3 commit 35f00c4
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 30 deletions.
6 changes: 5 additions & 1 deletion Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,14 @@ mcluster_sse: $(OBJECTS) $(LFLAGS)
$(CC) -c main.c -D SSE -lm
$(CC) $(OBJECTS) main.o -o mcluster_sse -lm $(CFLAGS)

mcluster_gpu: $(OBJECTS) $(LFLAGS) gpupot.gpu.o main.c
mcluster_ssegpu: $(OBJECTS) $(LFLAGS) gpupot.gpu.o main.c
$(CC) -c main.c -D SSE -D GPU -lm -I$(CUDA_PATH)/include
$(CC) $(OBJECTS) main.o gpupot.gpu.o -o mcluster_gpu -L$(CUDA_PATH)/lib64 -lcudart -lstdc++ -lm $(CFLAGS)

mcluster_gpu: gpupot.gpu.o main.c
$(CC) -c main.c -D GPU -I$(CUDA_PATH)/include
$(CC) main.o gpupot.gpu.o -o mcluster_gpu -L$(CUDA_PATH)/lib64 -lcudart -lstdc++ -lm

mcluster:
$(CC) -o mcluster main.c -lm

Expand Down
111 changes: 83 additions & 28 deletions main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,17 +101,17 @@ int main (int argv, char **argc) {
//Binary parameters
int nbin = 0; //Number of primordial binary systems
double fbin = 0.0; //Primordial binary fraction, number of binary systems = 0.5*N*fbin, only used when nbin is set to 0
int pairing = 3; //Pairing of binary components; 0= random pairing, 1= ordered pairing for components with masses M>msort, 2= random but separate pairing for components with masses m>Msort; 3= Uniform distribution of mass ratio (0.1<q<1.0) for m>Msort and random pairing for remaining (Kiminki & Kobulnicky 2012; Sana et al., 2012; Kobulnicky et al. 2014; implemented by Long Wang)
int pairing = 3; //Pairing of binary components; 0= random pairing, 1= ordered pairing for components with masses M>msort, 2= random but separate pairing for components with masses m>Msort; 3= Use period distribution for M>msort from Sana et al. (2011,2012) and Oh et al. (2015).
double msort = 5.0; //Stars with masses > msort will be sorted and preferentially paired into binaries if pairing = 1
int adis = 3; //Semi-major axis distribution; 0= flat ranging from amin to amax, 1= based on Kroupa (1995) period distribution, 2= based on Duquennoy & Mayor (1991) period distribution, 3= based on Kroupa (1995) period distribution for M<Msort; based on Sana et al. (2012); Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015) period distribution for M>Msort (implemented by Long Wang)
int OBperiods = 1; //Use period distribution for massive binaries with M_primary > msort from Sana & Evans (2011) if OBperiods = 1
int adis = 3; //Semi-major axis distribution; 0= flat ranging from amin to amax, 1= based on Kroupa (1995) period distribution, 2= based on Duquennoy & Mayor (1991) period distribution, 3= based on Kroupa (1995) period distribution
int OBperiods = 2; //1: Use period distribution for massive binaries with M_primary > msort from Sana & Evans (2011); 2: Uniform distribution of mass ratio (0.1<q<1.0) for m>Msort and random pairing for remaining (Kiminki & Kobulnicky 2012; Sana et al., 2012; Kobulnicky et al. 2014; Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015) period distribution for M>Msort (implemented by Long Wang)
double amin = 0.0001; //Minimum semi-major axis for adis = 0 [pc]
double amax = 0.01; //Maximum semi-major axis for adis = 0 [pc]
#ifdef SSE
int eigen = 0; //Use Kroupa (1995) eigenevolution for pre-main sequence short-period binaries; =0 off, =1 on [use either eigenevolution or BSE; BSE recommended when using SSE]
int BSE = 1; //Apply binary star evolution using BSE (Hurley, Tout & Pols 2002) =0 off, =1 on [use either eigenevolution or BSE; BSE recommended when using SSE]
#else
int eigen = 1; //Use Kroupa (1995) eigenevolution for pre-main sequence short-period binaries; =0 off, =1 on [use either eigenevolution or BSE; BSE recommended when using SSE]
int eigen = 2; // 1. Use Kroupa (1995) eigenevolution for pre-main sequence short-period binaries; 2. Use modified eigenevolution from Belloni et al. (2017); 0. Off
int BSE = 0; //Apply binary star evolution using BSE (Hurley, Tout & Pols 2002) [needs special compiling and BSE]; =0 off, =1 on [use either eigenevolution or BSE; BSE recommended when using SSE]
#endif

Expand Down Expand Up @@ -500,7 +500,7 @@ int main (int argv, char **argc) {
if (pairing) {
if (pairing == 1) printf("\nApplying ordered pairing for stars with masses > %.1f Msun.\n",msort);
else if (pairing == 2) printf("\nApplying random pairing for stars with masses > %.1f Msun.\n",msort);
else if (pairing == 3) printf("\nApplying uniform mass ratio distribution for stars with masses > %.1f Msun.\n",msort);
else if (pairing == 3) printf("\nApplying Sana et al. (2012) mass atio distribution for stars with masses > %.1f Msun.\n",msort);
order(star, N, M, msort, pairing);
} else {
randomize(star, N);
Expand Down Expand Up @@ -3491,8 +3491,8 @@ int get_binaries(int nbin, double **star, double M, double rvir, int pairing, in

//Specify semi-major axis
if (((m1*M>=msort) || (m2*M>=msort)) && (OBperiods)) {
if (adis == 3) {
double lPmin = 0.15, alpha = 0.45; //Sana et al., (2012); Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015)
if (OBperiods==2) {
double lPmin = 0.15, alpha = 0.45; //Sana et al., (2012); Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015)
double ralpha = 1/alpha, eta = alpha/0.23;
lP = pow(pow(lPmin,alpha) + eta*drand48(),ralpha);
} else {
Expand Down Expand Up @@ -3562,7 +3562,7 @@ int get_binaries(int nbin, double **star, double M, double rvir, int pairing, in


if (!i && OBperiods && msort) {
if (adis == 3) {
if (OBperiods==2) {
printf("\nDeriving semi-major axis distribution for binaries with primary masses > %.3f Msun from Sana et al., (2012); Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015) period distribution.\n",msort);
} else {
printf("\nDeriving semi-major axis distribution for binaries with primary masses > %.3f Msun from Sana & Evans (2011) period distribution.\n",msort);
Expand All @@ -3571,20 +3571,36 @@ int get_binaries(int nbin, double **star, double M, double rvir, int pairing, in


//Specify eccentricity distribution
if (!i && OBperiods && msort) printf("\nApplying thermal eccentricity distribution for low-mass systems and Sana & Evans (2011) eccentricity distribution for high-mass systems.\n");
if (!i && OBperiods && msort) {
if (OBperiods==1) {
printf("\nApplying thermal eccentricity distribution for low-mass systems and Sana & Evans (2011) eccentricity distribution for high-mass systems.\n");
} else {
printf("\nApplying thermal eccentricity distribution for low-mass systems and Sana et al. (2012), Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015) eccentricity distribution for high-mass systems.\n");
}

}
else if (!i) printf("\nApplying thermal eccentricity distribution.\n");

if (((m1*M>=msort) || (m2*M>=msort)) && (OBperiods)) {
double k1, k2;
double elim = 0.8; //maximum eccentricity for high-mass binaries
double fcirc = 0.3; //fraction of circular orbits among high-mass binaries
k2 = (fcirc*exp(0.5*elim)-1.0)/(exp(0.5*elim)-1.0);
k1 = fcirc-k2;
ecc = drand48();
if (ecc < fcirc) ecc = 0.0;
else {
ecc = 2.0*log((ecc-k2)/k1);
}
if (OBperiods==1) {
double k1, k2;
double elim = 0.8; //maximum eccentricity for high-mass binaries
double fcirc = 0.3; //fraction of circular orbits among high-mass binaries
k2 = (fcirc*exp(0.5*elim)-1.0)/(exp(0.5*elim)-1.0);
k1 = fcirc-k2;
ecc = drand48();
if (ecc < fcirc) ecc = 0.0;
else {
ecc = 2.0*log((ecc-k2)/k1);
}
}
else if (OBperiods==2){
double Pmin=pow(10.0,0.15); // consistent with Sana et al. (2012)
double elimit = 1.0 - pow((P*365.25/Pmin),-2.0/3.0);
do{
ecc = pow(drand48(), 1.0/0.55); // f=0.55ecc^(-0.45)
}while(ecc>elimit);
}
} else {
ecc = sqrt(drand48()); // Thermal distribution f(e)=2e
}
Expand All @@ -3603,8 +3619,8 @@ int get_binaries(int nbin, double **star, double M, double rvir, int pairing, in
m2*=M;
abin*=rvir;
if (!(OBperiods && ((m1>=msort) || (m2>=msort)))) {
if (m1>=m2) eigenevolution(&m1, &m2, &ecc, &abin);
else eigenevolution(&m2, &m1, &ecc, &abin);
if (m1>=m2) eigenevolution(&m1, &m2, &ecc, &abin, eigen);
else eigenevolution(&m2, &m1, &ecc, &abin, eigen);
}
m1/=M;
m2/=M;
Expand Down Expand Up @@ -3888,6 +3904,30 @@ int order(double **star, int N, double M, double msort, int pairing){
// Find the second one based on uniform mass ratio
if (i<N-1) {
double mpair = (drand48()*0.9+0.1)*masses[i][0];
// second index
int k = -1;
// find closest one
int k1 = i+1;
while (k1<N&&mask[k1]==1) k1++;
if (k1<N) {
double mk = fabs(masses[k1][0]-mpair);
double dm = mk;
int k2 = k1;
do {
k1 = k2;
mk = dm;
k2++;
while (k2<N && mask[k2]==1) k2++;
if(k2<N) dm = fabs(masses[k2][0]-mpair);
else dm=mk;
} while(dm<mk);
k = k1;
//printf("mass_ratio(true,real): %f %f\n",mpair/masses[i][0], masses[k][0]/masses[i][0]);
if(dm/mpair>1e-2)
printf("Warning: dm too large: m1=%F, mp=%f, m2=%f, dm=%f, i=%d, k=%d\n",masses[i][0],mpair,masses[k][0],i,k);
}
//printf("mpair =%f, k=%d, i=%d, mass[k]=%f\n",mpair,k,i,masses[k][0]);
/*
int k = i+1;
while (masses[k][0] > mpair) {
k++;
Expand All @@ -3913,6 +3953,7 @@ int order(double **star, int N, double M, double msort, int pairing){
else k = k1;
}
}
*/
if (k>0) {
if(j>=N) {
perror("\n Error!: Pairing binary star exceed the maximum index!\n");
Expand Down Expand Up @@ -3945,6 +3986,7 @@ int order(double **star, int N, double M, double msort, int pairing){
// Store index for remaining stars
mmrand[ileft][0] = drand48();
mmrand[ileft][1] = masses[i][1];
mask[i] = 1;
ileft++;
}
}
Expand Down Expand Up @@ -4317,8 +4359,16 @@ double rtnewt (double ecc, double ma) {
exit(-1);
}

int eigenevolution(double *m1, double *m2, double *ecc, double *abin){
int eigenevolution(double *m1, double *m2, double *ecc, double *abin, int opt){
if(*m1<*m2) {
double mtmp = *m2;
*m2 = *m1;
*m1 = mtmp;
}
double alpha = 28.0;
double eta=0.25;
double m1crit = 1.0;
if(opt==2) alpha *= pow(*m1,eta);
double beta = 0.75;
double mtot,lper,lperi,ecci,mtoti,r0,rperi,qold,qnew;

Expand All @@ -4333,23 +4383,26 @@ int eigenevolution(double *m1, double *m2, double *ecc, double *abin){
/* Circularisation */
r0 = alpha*RSUN;
rperi = *abin*(1.0-*ecc);
alpha = -pow((r0/rperi),beta);
alpha = -pow((r0/rperi),beta); // alpha -> rho

if (*ecc > 0) {
*ecc = exp(alpha+log(*ecc));
}

/* pre-ms mass-transfer */
qold = *m1/ *m2;
if (qold > 1.0) qold = 1.0/qold;
//qold = *m1/ *m2;
//if (qold > 1.0) qold = 1.0/qold;
qold = *m2 / *m1;
alpha = -alpha;
if (opt==2&&*m1<m1crit) alpha *= 0.5;
if (alpha > 1.0) {
qnew = 1.0;
if(opt==2&&*m1<m1crit) qnew = 0.9 + 0.1*drand48();
else qnew = 1.0;
} else {
qnew = qold + (1.0-qold)*alpha;
}

*m1 = max(*m1,*m2);
//*m1 = max(*m1,*m2);
*m2 = qnew * *m1;

mtot = *m1 + *m2;
Expand Down Expand Up @@ -5132,7 +5185,9 @@ void help(double msort) {
printf(" -b <value> (binary fraction, specify either B or b) \n");
printf(" -p <0|1|2|3> (binary pairing, 0= random, 1= ordered for M>%.1f Msun,\n",msort);
printf(" 2= random but separate pairing for M>%.1f Msun)\n",msort);
printf(" 3= random but uniform distribution of mass ratio (0.1<q<1.0) for M>%.1f Msun)\n",msort);
printf(" 3= random but use period distribution from Sana et al., (2012);\n");
printf(" Oh, S., Kroupa, P., & Pflamm-Altenburg, J. (2015)\n");
printf(" for M>%.1f Msun)\n",msort);
printf(" -s <number> (seed for randomization; 0= randomize by timer) \n");
printf(" -t <0|1|2|3> (tidal field; 0= no tidal field, 1= near-field, \n");
printf(" 2= point-mass, 3= Milky-Way potential) \n");
Expand Down
2 changes: 1 addition & 1 deletion main.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,7 +192,7 @@ int segregate(double **star, int N, double S);
int energy_order(double **star, int N, int Nstars);
int randomize(double **star, int N);
double rtnewt (double ecc, double ma);
int eigenevolution(double *m1, double *m2, double *ecc, double *abin);
int eigenevolution(double *m1, double *m2, double *ecc, double *abin, int opt);
int radial_profile(double **star, int N, double rvir, double M,int create_radial_profile, int create_cumulative_profile, int code, int *NNBMAX, double *RS0, double *Rh2D, double *Rh3D, int NNBMAX_NBODY6);
int cmd(double **star, int l, double Rgal, double *abvmag, double *vmag, double *BV, double *Teff, double *dvmag, double *dBV);
int output0(char *output, int N, int NNBMAX, double RS0, double dtadj, double dtout, double tcrit, double rvir, double mmean, int tf, int regupdate, int etaupdate, int mloss, int bin, int esc, double M, double mlow, double mup, double MMAX, double epoch, double dtplot, double Z, int nbin, double Q, double *RG, double *VG, double rtide, int gpu, double **star, int sse, int seed, double extmass, double extrad, double extdecay, double extstart);
Expand Down

0 comments on commit 35f00c4

Please sign in to comment.