diff --git a/Makefile b/Makefile index 2ba153a..dd72d2b 100644 --- a/Makefile +++ b/Makefile @@ -1,43 +1,39 @@ -IGNORE: +# Compiler FC = gfortran -O2 #-ffast-math #FC = gfortran-mp-4.4 -O2 #-ffast-math #use of -ffast-math may cause problems on some machines -LFLAGS = const_bse.h zdata.h +INCPATH=include -CC = gcc -O2 -fopenmp -Wall #-ffast-math +CC = gcc -O2 -fopenmp -Wall #-ffast-math #use "-D NOOMP" and remove "-fopenmp" for compilation without OpenMP -#CFLAGS = -L/opt/local/lib/gcc44/ -lgfortran CFLAGS = -L/usr/lib/ -lgfortran -CUFLAGS= -O3 -D WITH_CUDA5 -CUDA_PATH = /usr/local/cuda +CUFLAGS= -O3 -D WITH_CUDA5 -L$(CUDA_PATH)/lib64 -lcudart -lstdc++ -lm $(CFLAGS) +CUDA_PATH = $(shell which nvcc | rev | cut -d'/' -f3- | rev) SDK_PATH=$(CUDA_PATH)/samples -.f.o: - $(FC) -c $< +SOURCES := $(shell find . -type f -name '*.f') +INCLUDES := $(shell find . -type f -name '*.h') +OBJECTS := $(SOURCES:.f=.o) -SOURCE = \ -deltat.f evolv1.f hrdiag.f kick.f mlwind.f mrenv.f \ -ran3.f star.f zcnsts.f zfuncs.f\ -comenv.f corerd.f dgcore.f evolv2.f gntage.f \ -instar.f mix.f rl.f +%.o:%.f + $(FC) -c $^ -o $@ -OBJECTS = $(SOURCE:.f=.o) - -mcluster_sse: $(OBJECTS) $(LFLAGS) +mcluster_sse: $(OBJECTS) $(CC) -c main.c -D SSE -lm - $(CC) $(OBJECTS) main.o -o mcluster_sse -lm $(CFLAGS) + $(CC) $(OBJECTS) main.o -o mcluster_sse -lm $(CFLAGS) -mcluster_gpu: $(OBJECTS) $(LFLAGS) gpupot.gpu.o main.c +mcluster_gpu: $(OBJECTS) $(INCPATH)/gpupot.gpu.o $(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) + $(CC) $(OBJECTS) main.o gpupot.gpu.o -o mcluster_gpu $(CUFLAGS) -mcluster: +mcluster: + @echo $(OBJECTS) $(CC) -o mcluster main.c -lm -gpupot.gpu.o: gpupot.gpu.cu cuda_pointer.h - nvcc -c $(CUFLAGS) -Xcompiler "-fPIC -O3 -Wall" -I$(SDK_PATH)/common/inc -I. gpupot.gpu.cu +$(INCPATH)/gpupot.gpu.o: $(INCPATH)/gpupot.gpu.cu $(INCPATH)/cuda_pointer.h + nvcc -c $(CUFLAGS) -Xcompiler "-fPIC -O3 -Wall" -I$(SDK_PATH)/common/inc -I$(INCPATH) $(INCPATH)/gpupot.gpu.cu clean: - rm -f *.o mcluster_sse mcluster mcluster_gpu + rm -f $(INCPATH)/*.o *.o mcluster_sse mcluster mcluster_gpu diff --git a/README b/README deleted file mode 100644 index f35f827..0000000 --- a/README +++ /dev/null @@ -1,388 +0,0 @@ -************************************************************** -McLuster - a tool to make a star cluster - -Developed by Andreas H.W. Kuepper at AIfA Bonn, Germany, -in collaboration with Pavel Kroupa & Holger Baumgardt. - -25 March 2013 - -Detailed information on the code and the parameters can -be found in the following paper: - -Kuepper A.H.W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300 -"Mass segregation and fractal substructure in young massive clusters: -(I) the McLuster code and method calibration" - -This reference should be used for any publications that -utilize McLuster. - -For comments, questions and bug reports contact: -ahwkuepper@gmail.com -************************************************************** - -COMPILING THE CODE: - -For compiling the basic McLuster version use the MAKEFILE by -typing - -> make mcluster - -or simply type - -> cc -o mcluster main.c -O2 -fopenmp -lm - -This should work on most computers on which a recent -C-compiler is installed. - -For the SSE-version of McLuster you have to use the MAKEFILE. -Try compiling by typing - -> make mcluster_sse - -You may have to change the C- or FORTRAN-compiler in the -MAKEFILE by editing the variable CC or FC, respectively. -Furthermore, you may have to give the correct directory of -your gfortran installation to the C-compiler. This can be -done by changing the path in the CFLAGS variable - --L/usr/lib/ - -to the correct path of your gfortran installation. - -Note that some warnings like the following may appear: - -ld warning: alignment lost in merging tentative definition - -which can be ignored. - -************************************************************** - -HELP: - -Get a quick help by typing - -> mcluster -h - -for the basic version, or - -> mcluster_sse -h - -respectively for the SSE version. - -************************************************************** - -INTRODUCTION: - -McLuster can be run from the command line by passing arguments -to the code which specify the desired cluster, e.g., - -> mcluster -N1000 -R2.0 -P0 - -would give a cluster with a Plummer profile and a half-mass -radius of 2 pc, consisting of 1000 stars. For parameters -which are not specified default values are taken, which can be -changed within the code. Go to the top of the 'main' routine -for a complete list of available parameters with descriptions. - - -By typing - -> mcluster -h - -you get the allowed arguments which can be passed to McLuster -which will be described in detail in the following: - - - -N (number of stars) - -The number of stars may vary from 3 to ~10^6. Remember though -that some procedures within the code require ~N^2 computational -steps which may take incredibly long. Furthermore, processes -like mass segregation and fractalization may temporarily -need a lot of memory/or may temporarily generate more stars than -the required number. The parameter NMAX within the main -routine may have to be set to a higher value then. - - - -M (mass of cluster; specify either N or M) - -If you specify the total mass of your cluster then McLuster -will generate stars until the total mass reaches M. - - - -P <0|1|2|3|-1> (density profile; 0= Plummer, 1= King (1966), - 2= Subr et al. (2007) mass-segregated, - 3= 2-dimensional EFF/Nuker template, - -1= no density gradient) - -The density profile is used to generate the positions & -velocities of the stars. The final cluster is then scaled to -exactly match the half-mass radius which you specified. In case -you chose the EFF/Nuker profile McLuster, will integrate the density -distribution and determine the expected half-mass radius of -that specific distribution. The final cluster is then scaled -to this expectation value. - -For the generation of the Subr et al. profile the routine -PLUMIX is used, whose most recent version was incorporated -in McLuster. - -Note that some choices for the profile necessitate further -parameters to be specified. - - - -W <1-12> (W0 parameter for King model) - -In case you choose the King profile you"ll have to specify its -concentration which is done by fixing its W0 value. - - - -R (half-mass radius [pc], ignored for P = 3) - -The half-mass radius of the cluster has to be specified for all -density profiles except for the EFF/Nuker profile. For the latter it -gets determined automatically according to your choice of the -scale radius and the cut-off radius. All radii are set in pc! -By setting the half-mass radius to -1, McLuster will use the -cluster mass-half-mass radius relation from Marks & Kroupa (2012) -to estimate an initial half-mass radius according to your choice -of mass. - - - -r (scale radius of EFF/Nuker template [pc]) - -This specifies the scale radius, see Elson, Fall & Freeman (1987) -and or Lauer et al. (1995) for Nuker template. - - - -c (cut-off radius of EFF/Nuker template [pc]) - -Since the EFF and Nuker profiles extend infinitely they have to be -cut off at some radius. - - - -g (power-law slope(s) of EFF/Nuker template; use - once for EFF template; use three times for Nuker - template (outer slope, inner slope, transition) - -With this option you can set the power-law slopes of the EFF and Nuker -profiles. Use once for EFF to set the outer slope. Use three times for -outer slope, inner slope and the transition parameter. Note that the EFF -profile is a special case of the Nuker profile with the inner slope set -to 0.0 and the transition parameter set to 2.0. - - - -S <0.0-1.0> (degree of mass segregation; 0.0= no segregation) - -Mass segregation can be added to any choice of profile. In case of -the Subr et al. profile the S value should lie between 0.0 and 0.5 -in order to produce reasonable clusters (see Subr, Kroupa & Baumgardt -2007 for further details). All other density profiles allow values -from 0.0 to 1.0, where 1.0 corresponds to total segregation, i.e. -the most massive star sits most deeply within the cluster potential, -followed by the second most massive star, etc. The segregation of -stars is done following Baumgardt, de Marchi & Kroupa (2008). - - - -D <1.6-3.0> (fractal dimension; 3.0= no fractalization) - -Any density profile can be fractalized. The fractal dimension -specifies the degree of fractalization. Values of 2.6 are -reasonable, and values lower than 1.6 should be omitted. The -fractalization is realized similar as described in Whitworth & -Goodwin (2004). If you choose P = -1 you can get a fractal -distribution of stars without explicit density gradient. - - - -T (tcrit in N-body units) - -This specifies the time after which the N-body calculations -should be stopped. For details see Sverre Aarseth's NbodyX -handbook. - - - -Q (virial ratio) - -The virial ratio can be any positive value. 0.5 means virial -equlibrium, values >0.5 mean initial expansion, and values -<0.5 result in collapsing clusters. - - - -C <0|1|3|4> (code; 0= Nbody6, 1= Nbody4, - 3= table of stars, 4= Nbody7) - -If you want to use the output of McLuster as input for your -N-body computations you should specify the corresponding code. -You can also generate an ascii table of stars instead. - - - -A (dtadj in N-body units) - -Sets the adjustment time-step for Sverre's NbodyX codes. - - - -O (deltat in N-body units) - -Sets the output interval of NbodyX. - - - -G <0|1> (GPU usage; 0= no GPU, 1= use GPU) - -Specify whether you want to use Nbody6 with a GPU. - - - -o (output name of cluster model) - -Change the name of the output file(s). - - - -f <0|1|2|3|4> (IMF; 0= no IMF, 1= Kroupa (2001), - 2= user defined, 3= equal to "=1" but with optimal - sampling, 4= L3 IMF from Maschberger 2012) - -Set f to 0 for a cluster consisting of single-mass stars. For f -set to 1 the canonical two-part power-law Kroupa IMF is used from -0.08 Msun to 100 Msun. Those default values can be changed within -the main routine of the code or with the -m option. The first -m -will specify the lower mass limit, the second one will set the -upper mass limit. -If f is set to 2 you can specify the limiting masses and slopes -between those limits for up to 10 power-law slopes and 11 mass limits. -If you need more, change the parameters MAX_AN and MAX_MN within -main.h accordingly. -If f is set to 4 you can also use the -m option to set upper and -lower mass limits, and you can also use the -a option to specify -the power law slopes and the mu parameter of the L3 IMF. First -a -will set alpha, the second will set beta and the third is for mu. - - - -a (IMF slope for user defined IMF, may be used - multiple times, from low mass to high mass; or - for parameters alpha, beta and mu of L3 IMF) - - -m (IMF mass limits for the IMF; for user defined IMF - it may be used multiple times; always set limits - from low mass to high mass [Msun]) - -Specify the mass limits and the power-law slopes of the IMF as follows - -> mcluster -f2 -m 0.08 -a -1.3 -m 0.5 -a -2.3 -m 150.0 - -This would yield a Kroupa IMF ranging from 0.08 Msun to 150 Msun. - -> mcluster -f4 -m 0.01 -m 100.0 -a 2.3 -a 1.4 -a 0.2 - -This would give a L3 IMF from 0.01 to 100 Msun with alpha = 2.3, -beta = 1.4 and mu = 0.2. - - - -B (number of binary systems) - - -b (binary fraction, specify either B or b) - -For the binary content you can specify the fraction of stars -in binaries, 0.0 <= b <= 1.0, or the number of binary systems, -0 <= B <= N/2. Note that the total mass of the cluster may exceed -the mass you specified initially when using high binary fractions -and eigenevolution. This is due to the form of the eigenevolution -algorithm. Also, do not use eigenevolution when using BSE! - - - -p <0|1|2> (binary pairing, 0= random, 1= ordered for M>5Msun, - 2= random but separate pairing for M>5Msun) - -If p is 0 then the components of each binary are randomly drawn -from the stars of the cluster. If p is 1 then first all stars above -5 Msun are put in binaries, where the most massive star is put -together with the second most massive star, the third with the -forth, and so on until the desired binary fraction/number of -binaries is reached, or until no more stars above 5 Msun are -available. In this case the following stars are again drawn -randomly. The limit between random and ordered sampling can be -set within the main routine by changin the variable 'msort'. In -addition, binaries with primary masses above msort can be chosen -to follow a period distribution as given in Sana & Evans (2011) by -setting the parameter 'OBperiods' in the main routine to 1 (recommended). -If p is set to 2 then stars above 5 Msun are paired randomly among each -other and the rest is also paired randomly. - - -s (seed for randomization; 0= randomize by timer) - -If s is set to 0 then each call of McLuster will generate an -independent realization of the specified cluster parameters (random -number generator is seeded with the local time). Any other integer -number will result in exactly the same realization (except for rounding -differences). - - - -t <1|2|3> (tidal field; 1= near-field, 2= point-mass, - 3= Milky-Way potential) - -The tidal field can be specified according to the NbodyX manual. -By setting the cluster's galactocentric radius and its orbital -velocity, McLuster calculates the theoretical tidal radius and -cuts off the Plummer profile at this radius (if P is set to 0). - - - -e (epoch for stellar evolution [Myr]) - -This parameter can be used to set the cluster age but only with -mcluster_sse. McLuster will use the SSE routines by Hurley et al. -to evolve the cluster stars from their initial masses to the -masses at the specified age. It also calculates the according -luminosites, stellar radii, etc. If you are interested in a -cluster of age zero but also need those stellar parameters, then -use mcluster_sse and set e to a small value such as 0.1. Within -the main.c file you can also activate BSE for binary evolution of -primordial binaries. Therefore, set BSE = 1. - - - -Z (metallicity [0.0001-0.03, 0.02 = solar]) - -Use Z to set the cluster metallicity. Alternatively, you can set -the parameter FeH within the main routine to your desired [Fe/H] -value such that McLuster will calculate the according Z value. - - - -X (galactocentric radius vector, use 3x, [pc]) - -Use X three times to specify the position of the cluster within -the galaxy, e.g., - -> mcluster -X 8500 -X 0 -X 0 - -to set the cluster to a position of x = 8.5 kpc, y = z = 0 kpc. - - - -V (cluster velocity vector, use 3x, [km/s]) - -With V you can specify the cluster's orbital velocity, e.g., - -> mcluster -V 0 -V 220 -V 0 - -for a cluster orbital velocity of 220 km/s within the x-y- -plane. - - - -u <0|1> (output units; 0= Nbody, 1= astrophysical) - -When you generate a table of stars you may be interested in -astrophysical units rather than Nbody units. - - - -h (display this help) - - -? (display this help) - -Get help with those options. - - -Further parameters can be specified within the code. Note -that you have to compile the code again after any change. -Therefore, first type - -> make clean - -ignore any error message, and compile the code again. - -************************************************************** diff --git a/README.md b/README.md new file mode 100644 index 0000000..6fed2fa --- /dev/null +++ b/README.md @@ -0,0 +1,326 @@ +## McLuster - a tool to make a star cluster + + * Developer: Andreas H.W. Kuepper at [AIfA Bonn, Germany](https://astro.uni-bonn.de/en) + * Collaborators: [Pavel Kroupa](https://astro.uni-bonn.de/~pavel/) & + [Holger Baumgardt](https://people.smp.uq.edu.au/HolgerBaumgardt/) + * Date: 25 March 2013 + +Detailed information on the code and the parameters can +be found in the following paper: + +> Kuepper A.H.W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300 +> [Mass segregation and fractal substructure in young massive clusters: +> (I) the McLuster code and method calibration](https://arxiv.org/abs/1107.2395) + +This reference should be used for any publications that +utilize McLuster. + +For comments, questions and bug reports contact: +ahwkuepper@gmail.com + +### Compiling the code + +For compiling the basic McLuster version use the MAKEFILE by typing +`make mcluster` or simply type: +```bash +cc -o mcluster main.c -O2 -fopenmp -lm +``` + +This should work on most computers on which a recent +C-compiler is installed. + +For the SSE-version of McLuster you have to use the MAKEFILE. +Try compiling by typing `make mcluster_sse`. + +You may have to change the C or FORTRAN compiler in the +MAKEFILE by editing the variable CC or FC, respectively. +Furthermore, you may have to give the correct directory of +your gfortran installation to the C compiler. This can be +done by changing the path in the CFLAGS variable `-L/usr/lib/` +to the correct path of your gfortran installation. + +Note that some warnings like the following may appear: +ld warning: alignment lost in merging tentative definition +which can be ignored. + +### Display available commands + +Get a quick help by typing +```bash +mcluster -h +``` +for the basic version, or +```bash +mcluster_sse -h +``` +respectively for the SSE version. + +### How to use it + +McLuster can be run from the command line by passing arguments +to the code which specify the desired cluster, e.g., +```bash +mcluster -N1000 -R2.0 -P0 +``` +would give a cluster with a Plummer profile and a half-mass +radius of 2 pc, consisting of 1000 stars. For parameters +which are not specified default values are taken, which can be +changed within the code. Go to the top of the 'main' routine +for a complete list of available parameters with descriptions. + +By typing +```bash +mcluster -h +``` +you get the allowed arguments which can be passed to McLuster +which will be described in detail in the following: + +* *-N (number of stars)* + + The number of stars may vary from 3 to ~10^6. Remember though + that some procedures within the code require ~N^2 computational + steps which may take incredibly long. Furthermore, processes + like mass segregation and fractalization may temporarily + need a lot of memory/or may temporarily generate more stars than + the required number. The parameter NMAX within the main + routine may have to be set to a higher value then. + +* *-M (mass of cluster; specify either N or M)* + + If you specify the total mass of your cluster then McLuster + will generate stars until the total mass reaches M. + +* *-P <0|1|2|3|-1> (density profile; 0= Plummer, 1= King (1966), 2= Subr et al. (2007) mass-segregated, 3= 2-dimensional EFF/Nuker template,-1= no density gradient)* + + The density profile is used to generate the positions & + velocities of the stars. The final cluster is then scaled to + exactly match the half-mass radius which you specified. In case + you chose the EFF/Nuker profile McLuster, will integrate the density + distribution and determine the expected half-mass radius of + that specific distribution. The final cluster is then scaled + to this expectation value. + + For the generation of the Subr et al. profile the routine + PLUMIX is used, whose most recent version was incorporated + in McLuster. + + Note that some choices for the profile necessitate further + parameters to be specified. + +* *-W <1-12> (W0 parameter for King model)* + + In case you choose the King profile you"ll have to specify its + concentration which is done by fixing its W0 value. + +* *-R (half-mass radius [pc], ignored for P = 3)* + + The half-mass radius of the cluster has to be specified for all + density profiles except for the EFF/Nuker profile. For the latter it + gets determined automatically according to your choice of the + scale radius and the cut-off radius. All radii are set in pc! + By setting the half-mass radius to -1, McLuster will use the + cluster mass-half-mass radius relation from Marks & Kroupa (2012) + to estimate an initial half-mass radius according to your choice + of mass. + +* *-r (scale radius of EFF/Nuker template [pc])* + + This specifies the scale radius, see Elson, Fall & Freeman (1987) + and or Lauer et al. (1995) for Nuker template. + +* *-c (cut-off radius of EFF/Nuker template [pc])* + + Since the EFF and Nuker profiles extend infinitely they have to be + cut off at some radius. + +* *-g (power-law slope(s) of EFF/Nuker template; use once for EFF template; use three times for Nuker template (outer slope, inner slope, transition)* + + With this option you can set the power-law slopes of the EFF and Nuker + profiles. Use once for EFF to set the outer slope. Use three times for + outer slope, inner slope and the transition parameter. Note that the EFF + profile is a special case of the Nuker profile with the inner slope set + to 0.0 and the transition parameter set to 2.0. + +* *-S <0.0-1.0> (degree of mass segregation; 0.0= no segregation)* + + Mass segregation can be added to any choice of profile. In case of + the Subr et al. profile the S value should lie between 0.0 and 0.5 + in order to produce reasonable clusters (see Subr, Kroupa & Baumgardt + 2007 for further details). All other density profiles allow values + from 0.0 to 1.0, where 1.0 corresponds to total segregation, i.e. + the most massive star sits most deeply within the cluster potential, + followed by the second most massive star, etc. The segregation of + stars is done following Baumgardt, de Marchi & Kroupa (2008). + +* *-D <1.6-3.0> (fractal dimension; 3.0= no fractalization)* + + Any density profile can be fractalized. The fractal dimension + specifies the degree of fractalization. Values of 2.6 are + reasonable, and values lower than 1.6 should be omitted. The + fractalization is realized similar as described in Whitworth & + Goodwin (2004). If you choose P = -1 you can get a fractal + distribution of stars without explicit density gradient. + +* *-T (tcrit in N-body units)* + + This specifies the time after which the N-body calculations + should be stopped. For details see Sverre Aarseth's NbodyX + handbook. + +* *-Q (virial ratio)* + + The virial ratio can be any positive value. 0.5 means virial + equlibrium, values >0.5 mean initial expansion, and values + <0.5 result in collapsing clusters. + +* *-C <0|1|3|4> (code; 0= Nbody6, 1= Nbody4, 3= table of stars, 4= Nbody7)* + + If you want to use the output of McLuster as input for your + N-body computations you should specify the corresponding code. + You can also generate an ascii table of stars instead. + +* *-A (dtadj in N-body units)* + + Sets the adjustment time-step for Sverre's NbodyX codes. + +* *-O (deltat in N-body units)* + + Sets the output interval of NbodyX. + +* *-G <0|1> (GPU usage; 0= no GPU, 1= use GPU)* + + Specify whether you want to use Nbody6 with a GPU. + +* *-o (output name of cluster model)* + + Change the name of the output file(s). + +* *-f <0|1|2|3|4> (IMF; 0= no IMF, 1= Kroupa (2001), 2= user defined, 3= equal to "=1" but with optimal sampling, 4= L3 IMF from Maschberger 2012)* + + Set f to 0 for a cluster consisting of single-mass stars. For f + set to 1 the canonical two-part power-law Kroupa IMF is used from + 0.08 Msun to 100 Msun. Those default values can be changed within + the main routine of the code or with the -m option. The first -m + will specify the lower mass limit, the second one will set the + upper mass limit. + + If f is set to 2 you can specify the limiting masses and slopes + between those limits for up to 10 power-law slopes and 11 mass limits. + If you need more, change the parameters MAX_AN and MAX_MN within + main.h accordingly. + + If f is set to 4 you can also use the -m option to set upper and + lower mass limits, and you can also use the -a option to specify + the power law slopes and the mu parameter of the L3 IMF. First -a + will set alpha, the second will set beta and the third is for mu. + +* *-a (IMF slope for user defined IMF, may be used multiple times, from low mass to high mass; or for parameters alpha, beta and mu of L3 IMF)* +* *-m (IMF mass limits for the IMF; for user defined IMF it may be used multiple times; always set limits from low mass to high mass [Msun])* + + Specify the mass limits and the power-law slopes of the IMF as follows + + ```bash + mcluster -f2 -m 0.08 -a -1.3 -m 0.5 -a -2.3 -m 150.0 + ``` + This would yield a Kroupa IMF ranging from 0.08 Msun to 150 Msun. + ```bash + mcluster -f4 -m 0.01 -m 100.0 -a 2.3 -a 1.4 -a 0.2 + ``` + This would give a L3 IMF from 0.01 to 100 Msun with alpha = 2.3, + beta = 1.4 and mu = 0.2. + +* *-B (number of binary systems)* +* *-b (binary fraction, specify either B or b)* + + For the binary content you can specify the fraction of stars + in binaries, *0.0 <= b <= 1.0*, or the number of binary systems, + *0 <= B <= N/2*. Note that the total mass of the cluster may exceed + the mass you specified initially when using high binary fractions + and eigenevolution. This is due to the form of the eigenevolution + algorithm. Also, do not use eigenevolution when using BSE! + +* *-p <0|1|2> (binary pairing, 0= random, 1= ordered for M>5Msun, 2= random but separate pairing for M>5Msun)* + + If p is 0 then the components of each binary are randomly drawn + from the stars of the cluster. If p is 1 then first all stars above + 5 Msun are put in binaries, where the most massive star is put + together with the second most massive star, the third with the + forth, and so on until the desired binary fraction/number of + binaries is reached, or until no more stars above 5 Msun are + available. In this case the following stars are again drawn + randomly. The limit between random and ordered sampling can be + set within the main routine by changin the variable 'msort'. In + addition, binaries with primary masses above msort can be chosen + to follow a period distribution as given in Sana & Evans (2011) by + setting the parameter 'OBperiods' in the main routine to 1 (recommended). + If p is set to 2 then stars above 5 Msun are paired randomly among each + other and the rest is also paired randomly. + +* *-s (seed for randomization; 0= randomize by timer)* + + If s is set to 0 then each call of McLuster will generate an + independent realization of the specified cluster parameters (random + number generator is seeded with the local time). Any other integer + number will result in exactly the same realization (except for rounding + differences). + +* *-t <1|2|3> (tidal field; 1= near-field, 2= point-mass, 3= Milky-Way potential)* + + The tidal field can be specified according to the NbodyX manual. + By setting the cluster's galactocentric radius and its orbital + velocity, McLuster calculates the theoretical tidal radius and + cuts off the Plummer profile at this radius (if P is set to 0). + +* *-e (epoch for stellar evolution [Myr])* + + This parameter can be used to set the cluster age but only with + mcluster_sse. McLuster will use the SSE routines by Hurley et al. + to evolve the cluster stars from their initial masses to the + masses at the specified age. It also calculates the according + luminosites, stellar radii, etc. If you are interested in a + cluster of age zero but also need those stellar parameters, then + use mcluster_sse and set e to a small value such as 0.1. Within + the main.c file you can also activate BSE for binary evolution of + primordial binaries. Therefore, set BSE = 1. + +* *-Z (metallicity [0.0001-0.03, 0.02 = solar])* + + Use Z to set the cluster metallicity. Alternatively, you can set + the parameter FeH within the main routine to your desired [Fe/H] + value such that McLuster will calculate the according Z value. + +* *-X (galactocentric radius vector, use 3x, [pc])* + + Use X three times to specify the position of the cluster within + the galaxy, e.g., + ```bash + mcluster -X 8500 -X 0 -X 0 + ``` + to set the cluster to a position of x = 8.5 kpc, y = z = 0 kpc. + +* *-V (cluster velocity vector, use 3x, [km/s])* + + With V you can specify the cluster's orbital velocity, e.g., + ```bash + mcluster -V 0 -V 220 -V 0 + ``` + for a cluster orbital velocity of 220 km/s within the x-y- + plane. + +* *-u <0|1> (output units; 0= Nbody, 1= astrophysical)* + + When you generate a table of stars you may be interested in + astrophysical units rather than Nbody units. + +* *-h (display this help)* +* *-? (display this help)* + + Get help with those options. + +Further parameters can be specified within the code. Note +that you have to compile the code again after any change. +Therefore, first type +```bash +make clean +``` +ignore any error message, and compile the code again. diff --git a/McLusterManual.pdf b/doc/McLusterManual.pdf similarity index 100% rename from McLusterManual.pdf rename to doc/McLusterManual.pdf diff --git a/McLusterManual.tex b/doc/McLusterManual.tex similarity index 100% rename from McLusterManual.tex rename to doc/McLusterManual.tex diff --git a/mn2e.cls b/doc/mn2e.cls similarity index 100% rename from mn2e.cls rename to doc/mn2e.cls diff --git a/comenv.f b/include/comenv.f similarity index 100% rename from comenv.f rename to include/comenv.f diff --git a/const_bse.h b/include/const_bse.h similarity index 100% rename from const_bse.h rename to include/const_bse.h diff --git a/corerd.f b/include/corerd.f similarity index 100% rename from corerd.f rename to include/corerd.f diff --git a/cuda_pointer.h b/include/cuda_pointer.h similarity index 100% rename from cuda_pointer.h rename to include/cuda_pointer.h diff --git a/deltat.f b/include/deltat.f similarity index 100% rename from deltat.f rename to include/deltat.f diff --git a/dgcore.f b/include/dgcore.f similarity index 100% rename from dgcore.f rename to include/dgcore.f diff --git a/evolv1.f b/include/evolv1.f similarity index 100% rename from evolv1.f rename to include/evolv1.f diff --git a/evolv2.f b/include/evolv2.f similarity index 100% rename from evolv2.f rename to include/evolv2.f diff --git a/gntage.f b/include/gntage.f similarity index 100% rename from gntage.f rename to include/gntage.f diff --git a/gpupot.gpu.cu b/include/gpupot.gpu.cu similarity index 100% rename from gpupot.gpu.cu rename to include/gpupot.gpu.cu diff --git a/hrdiag.f b/include/hrdiag.f similarity index 100% rename from hrdiag.f rename to include/hrdiag.f diff --git a/instar.f b/include/instar.f similarity index 100% rename from instar.f rename to include/instar.f diff --git a/kick.f b/include/kick.f similarity index 100% rename from kick.f rename to include/kick.f diff --git a/main.h b/include/main.h similarity index 100% rename from main.h rename to include/main.h diff --git a/mix.f b/include/mix.f similarity index 100% rename from mix.f rename to include/mix.f diff --git a/mlwind.f b/include/mlwind.f similarity index 100% rename from mlwind.f rename to include/mlwind.f diff --git a/mrenv.f b/include/mrenv.f similarity index 100% rename from mrenv.f rename to include/mrenv.f diff --git a/ran3.f b/include/ran3.f similarity index 100% rename from ran3.f rename to include/ran3.f diff --git a/rl.f b/include/rl.f similarity index 100% rename from rl.f rename to include/rl.f diff --git a/star.f b/include/star.f similarity index 100% rename from star.f rename to include/star.f diff --git a/zcnsts.f b/include/zcnsts.f similarity index 100% rename from zcnsts.f rename to include/zcnsts.f diff --git a/zdata.h b/include/zdata.h similarity index 100% rename from zdata.h rename to include/zdata.h diff --git a/zfuncs.f b/include/zfuncs.f similarity index 100% rename from zfuncs.f rename to include/zfuncs.f diff --git a/main.c b/main.c index 4b1e1f7..01d0390 100644 --- a/main.c +++ b/main.c @@ -36,7 +36,6 @@ * star[][14] = Zstar, metallicity of star * ***************************************************************************/ - #include #include #include @@ -44,5105 +43,6232 @@ #include #include #include + #ifdef GPU #include #include #endif -#include "main.h" -//use OpenMP if not specified otherwise -#ifndef NOOMP -#include -#endif +// Use OpenMP if not specified otherwise +#ifndef NOOMP +#include +#endif + +#include "include/main.h" + + +int main (int argv, char **argc) { + + /******************* + * Input variables * + *******************/ + + // Basic physical parameters + + // Number of stars, Mcl will be set to 0 if specified! + int N = 0; + // Total mass of the cluster, only used when N is set to 0, necessary for + // usage of maximum stellar mass relation of Weidner & Kroupa (2006) + double Mcl = 1000.0; + // Density profile; + // =0 Plummer model, + // =1 King model (based on king0.f by D.C. Heggie), + // =2 mass segregated Subr model (Subr et al. 2007), + // =3 2-dimensional EFF template (Elson, Fall & Freeman 1987) or + // Nuker template (Lauer et al. 1995), + // =-1 no density gradient + int profile = 0; + // King's W0 paramter [0.3-12.0] + double W0 = 5.0; + // Fraction of mass segregation for profile; =0.0 unsegregated, + // =1.0 completely segregated (take maximally S=0.99 for profile=2) + double S = 0.0; + // Fractal dimension; =3.0 unfractal, =2.6 2/8 fractal, =2.0 4/8 fractal, + // =1.6 6/8 fractal, (2^D children per parent, following Goodwin & + // Whitworth 2004) + double D = 3.0; + // Initial virial ratio; =0.5 virial equilibrium, <0.5 collapsing, + // >0.5 expanding + double Q = 0.5; + // Half-mass radius [pc], ignored if profile = 3, set =-1 for using Marks + // & Kroupa (2012) Mcl-Rh relation + double Rh = 0.8; + // Power-law slopes of EFF/Nuker templates (outer slope, inner slope, + // transition); set gamma[1] = 0.0 and gamma[2] = 2.0 for EFF (profile = 3) + double gamma[3] = {2.0, 0.0, 2.0}; + // Scale radius of EFF/Nuker template (profile = 3) [pc] + double a = 1.0; + // Cut-off radius for EFF/Nuker template (profile = 3) [pc] + double Rmax = 100.0; + // Simulation time [N-body units (Myr in Nbody6 custom)] + double tcrit = 100.0; + // Tidal field: =0 no tidal field, + // =1 Near-field approximation, + // =2 point-mass galaxy, + // =3 Allen & Santillan (1991) MW potential + // (or Sverre's version of it) + int tf = 3; + // Initial Galactic coordinates of the cluster [pc] + double RG[3] = {8500.0,0.0,0.0}; + // Initial velocity of the cluster [km/s] + double VG[3] = {0.0,220.0,0.0}; + + //Mass function parameters + + // 0 = single mass stars; + // 1 = use Kroupa (2001) mass function; + // 2 = use multi power law (based on mufu.c by L.Subr) + int mfunc = 1; + // Stellar mass in case of single-mass cluster + double single_mass = 1.0; + // Lower mass limit for mfunc = 1 & mfunc = 4 + double mlow = 0.08; + // Upper mass limit for mfunc = 1 & mfunc = 4 + double mup = 100.0; + // alpha slopes for mfunc = 2 + double alpha[MAX_AN] = {-1.35, -2.35, -2.7, 0.0, 0.0}; + // mass limits for mfunc = 2 + double mlim[MAX_MN] = {0.08, 0.5, 4.0, 100.0, 0.0, 0.0}; + // alpha slope for mfunc = 4 (L3 mass function, Maschberger 2012) + double alpha_L3 = 2.3; + // beta slope for mfunc = 4 + double beta_L3 = 1.4; + // mu parameter for mfunc = 4 + double mu_L3 = 0.2; + // Usage of Weidner & Kroupa (2006) relation for most massive star; + // =0 off, =1 on + int weidner = 0; + //Stellar evolution; 0 = off, 3 = Eggleton, Tout & Hurley [KZ19] + int mloss = 3; + // Use random kick velocity and present-day escape velocity to determine + // retention of compact remnants & evolved binary components + // (only for SSE/BSE version); =0 off, =1 on + int remnant = 1; + // Age of the cluster, i.e. star burst has been ... Myr before + // [e.g. 1000.0, default = 0.0] [needs special compiling and SSE by + // Hurley, Pols & Tout] + double epoch = 0.0; + // Metallicity [0.0001-0.03, 0.02 = solar] + double Z = 0.02; + // Metallicity [Fe/H], only used when Z is set to 0 + double FeH = -1.41; + // Usage of Prantzos (2007) relation for the life-times of stars. + // Set upper mass limit to Lifetime(mup) >= epoch + int prantzos = 0; + + // Binary parameters + + // Number of primordial binary systems + int nbin = 0; + // Primordial binary fraction, number of binary systems = 0.5*N*fbin, + // only used when nbin is set to 0 + double fbin = 0.0; + // 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.1Msort and + // random pairing for remaining (Kiminki & Kobulnicky 2012; + // Sana et al., 2012; Kobulnicky et al. 2014; implemented by Long Wang) + int pairing = 3; + // Stars with masses > msort will be sorted and preferentially paired into + // binaries if pairing = 1 + double msort = 5.0; + // 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 MMsort (implemented by Long Wang) + int adis = 3; + // Use period distribution for massive binaries with M_primary > msort + // from Sana & Evans (2011) if OBperiods = 1 + int OBperiods = 1; + // Minimum semi-major axis for adis = 0 [pc] + double amin = 0.0001; + // Maximum semi-major axis for adis = 0 [pc] + double amax = 0.01; +#ifdef SSE + // 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 = 0; + // Apply binary star evolution using BSE (Hurley, Tout & Pols 2002) + // =0 off, + // =1 on [use either eigenevolution or BSE; + // BSE recommended when using SSE] + int BSE = 1; +#else + // 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 = 1; + // 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] + int BSE = 0; +#endif + + // Gas parameters (only used for Nbody6 input) + + // external Plummer (gas) sphere mass [Msun] + double extmass = 0.0; + // external Plummer (gas) sphere scale factor [pc] + double extrad = 0.0; + // decay time for gas expulsion (set 0 for no decay) [Myr] + double extdecay = 0.0; + // delay time for start of gas expulsion [Myr] + double extstart = 0.0; + + // Code parameters + + // Nbody version: + // =0 Nbody6, + // =1 Nbody4, + // =2 Nbody6 custom, + // =3 only create output list of stars, + // =4 Nbody7 (not yet fully functional), + // =5 Nbody6++GPU + int code = 3; + // Number seed for random number generator; + // =0 for randomization by local time + unsigned int seed = 0; + // Name of output files + char *output = "test"; + // DTADJ [N-body units (Myr in Nbody6 custom)], energy-check time step + double dtadj = 1.0; + // DELTAT [N-body units (Myr in Nbody6 custom)], + // output interval, must be multiple of DTADJ + double dtout = 1.0; + // DTPLOT [N-body units (Myr in Nbody6 custom)], + // output of HRdiagnostics, should be multiple of DTOUT, + // set to zero if output not desired + double dtplot = 1.0; + // Use of GPU, + // 0= off, + // 1= on + int gpu = 0; + // Update of regularization parameters during computation; + // 0 = off, + // 0 > on + int regupdate = 0; + // Update of ETAI & ETAR during computation; + // 0 = off, + // 0 > on + int etaupdate = 0; + // Removal of escapers; + // 0 = no removal, + // 1 = regular removal at 2*R_tide; + // 2 = removal and output in ESC + int esc = 2; + // Units of McLuster output; + // 0= Nbody-Units, + // 1= astrophysical units + int units = 1; + + // McLuster internal parameters + + // Make cluster half-mass radius exactly match the expected/desired value; + // =0 off, + // =1 on (recommended) + int match = 1; + // Force spherical symmetry for fractal clusters; + // =0 off, + // =1 on (recommended) + int symmetry = 1; + // Make energy check at end of McLuster; + // =0 off, + // =1 on + int check = 0; + // Creates a radial density profile and prints it to the screen; + // =0 off, + // =1 on + int create_radial_profile = 1; + // Creates a radial cumulative profile and prints it to the screen; + // =0 off, + // =1 on + int create_cumulative_profile = 1; + // Distance of cluster from sun for artificial CMD with + // observational errors [pc] + double Rgal = 10000.0; + // Solar metallicity + double Zsun = 0.02; + // Maximum number of stars & orbits allowed in McLuster + int NMAX = 1500000; + // Maximum number of neighbours allowed in NBODY6 + int NNBMAX_NBODY6 = 500; + // Maximum stellar mass allowed in McLuster [Msun] + double upper_IMF_limit = 150.0; + // Counter for number of alpha slopes for mfunc = 2 + int an = 0; + // Counter for number of mass limits for mfunc = 1, 2 & 4 + int mn = 0; + // Counter for components of galactocentric radius vector + int xn = 0; + // Counter for components of cluster velocity vector + int vn = 0; + // Counter for external potential input parameters + int xx = 0; + // Counter for EFF/Nuker profile parameters + int gn = 0; + // Input array for external potential parameters + double extgas[4]; + + + // SSE internal parameters (see Hurley, Pols & Tout 2000) + + // Reimers mass-loss coefficent (neta*4x10^-13; 0.5 normally) + value1_.neta = 0.5; + // Binary enhanced mass loss parameter (inactive for single) + value1_.bwind = 0.0; + // Helium star mass loss factor (1.0 normally) + value1_.hewind = 1.0; + // Maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1) + value1_.mxns = 3.0; + // Time-step parameter in evolution phase: MS (0.05) + points_.pts1 = 0.05; + // Time-step parameter in evolution phase: GB, CHeB, AGB, HeGB (0.01) + points_.pts2 = 0.01; + // Time-step parameter in evolution phase: HG, HeMS (0.02) + points_.pts3 = 0.02; + // Kick velocities + value4_.sigma = 190.0; + // bhflag > 0 allows velocity kick at BH formation + value4_.bhflag = 1; + + // BSE internal parameters (see Hurley, Pols & Tout 2002) + + // ceflag > 0 activates spin-energy correction in common-envelope (0) + // #defunct# + flags_.ceflag = 0; + // tflag > 0 activates tidal circularisation (1) + flags_.tflag = 1; + // ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0) + flags_.ifflag = 0; + // nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ,572,407(1) + flags_.nsflag = 1; + // wdflag > 0 uses modified-Mestel cooling for WDs (0) + flags_.wdflag = 1; + + // beta is wind velocity factor: proportional to vwind**2 (1/8) + value5_.beta = 1.0/8.0; + // xi is the wind accretion efficiency factor (1.0) + value5_.xi = 1.0; + // acc2 is the Bondi-Hoyle wind accretion factor (3/2) + value5_.acc2 = 3.0/2.0; + // epsnov is the fraction of accreted matter retained in nova + // eruption (0.001) + value5_.epsnov = 0.001; + // eddfac is Eddington limit factor for mass transfer (1.0) + value5_.eddfac = 1.0; + // gamma is the angular momentum factor for mass lost during Roche (-1.0) + value5_.gamma = -1.0; + // alpha1 is the common-envelope efficiency parameter (1.0) + value2_.alpha1 = 1.0; + // lambda is the binding energy factor for common envelope evolution (0.5) + value2_.lambda = 0.5; + + + // Command line input + int option; + while ((option = getopt(argv, argc, "N:M:P:W:R:r:c:g:S:D:T:Q:C:A:O:G:o:f:a:m:B:b:p:s:t:e:Z:X:x:V:u:h:?")) != -1) switch (option) + { + case 'N': N = atoi(optarg); Mcl = 0.0; break; + case 'M': Mcl = atof(optarg); N = 0; break; + case 'P': profile = atoi(optarg); break; + case 'W': W0 = atof(optarg); break; + case 'R': Rh = atof(optarg); break; + case 'r': a = atof(optarg); break; + case 'c': Rmax = atof(optarg); break; + case 'g': + if (gn < 3) { gamma[gn] = atof(optarg); gn++; break; } + case 'S': S = atof(optarg); break; + case 'D': D = atof(optarg); break; + case 'T': tcrit = atof(optarg); break; + case 'Q': Q = atof(optarg); break; + case 'C': code = atoi(optarg); break; + case 'A': dtadj = atof(optarg); break; + case 'O': dtout = atof(optarg); break; + case 'G': gpu = atoi(optarg); break; + case 'o': output = optarg; break; + case 'f': mfunc = atoi(optarg); break; + case 'a' : + if (an < MAX_AN) { + alpha[an] = atof(optarg); + if (an == 0) alpha_L3 = atof(optarg); + if (an == 1) beta_L3 = atof(optarg); + if (an == 2) mu_L3 = atof(optarg); + an++; + break; + } else { printf("\nError: Number of alphas exceeded maximum limit of %d\n", MAX_AN); return 1; } + case 'm' : + if (mn < MAX_MN) { + mlim[mn] = atof(optarg); + if (mn == 0) mlow = atof(optarg); + if (mn == MAX_MN-1) mup = atof(optarg); + mn++; + break; + } else { printf("\nError: Number of mass params exceded maximum limit of %d\n", MAX_MN); return 1; } + case 'B': nbin = atoi(optarg); break; + case 'b': fbin = atof(optarg); break; + case 'p': pairing = atoi(optarg); break; + case 's': seed = atoi(optarg); break; + case 't': tf = atoi(optarg); break; + case 'e': epoch = atof(optarg); break; + case 'Z': Z = atof(optarg); break; + case 'X' : + if (xn < 3) { RG[xn] = atof(optarg); xn++; break; } + case 'V' : + if (vn < 3) { VG[vn] = atof(optarg); vn++; break; } + case 'x' : + if (xx < 4) { extgas[xx] = atof(optarg); xx++; break; } + case 'u': units = atoi(optarg); break; + case ':': + case 'h': help(msort); return 1; + case '?': help(msort); return 1; + }; + + if (mn-1 > 0) mup = mlim[mn-1]; + + // print summary of input parameters to .info file + info(output, N, Mcl, profile, W0, S, D, Q, Rh, gamma, a, Rmax, tcrit, tf, + RG, VG, mfunc, single_mass, mlow, mup, alpha, mlim, alpha_L3, beta_L3, + mu_L3, weidner, mloss, remnant, epoch, Z, prantzos, nbin, fbin, + pairing, msort, adis, amin, amax, eigen, BSE, extmass, extrad, + extdecay, extstart, code, seed, dtadj, dtout, dtplot, gpu, regupdate, + etaupdate, esc, units, match, symmetry, OBperiods); + + /********* + * Start * + *********/ + + printf("\n-----START---- \n"); + +#ifdef NOOMP + clock_t t1, t2; + // start stop-watch + t1 = clock(); +#else + double t1, t2; +#pragma omp parallel + { + if (omp_get_thread_num() == 0) + if (omp_get_num_threads() > 1) + printf("\n\nUsing OpenMP with %d threads\n", + omp_get_num_threads()); + + t1 = omp_get_wtime(); //start stop-watch + } +#endif + + if (seed) { + // initialize random number generator by seed + srand48(seed); + } + else { + // initialize random number generator by local time + seed = (unsigned) time(NULL); + seed %= 100000; + srand48(seed); + } + + printf ("\n\nRandom seed = %i\n\n\n", seed); + if (seed) { + // idum is the random number seed used in the kick routine. + value3_.idum = seed; + } + else { + value3_.idum = 10000000.0*drand48(); + } + + int i,j; + // Total mass [M_sun] + double M; + // Mean stellar mass [M_sun] + double mmean; + // Maximum neighbour number (Nbody6 only) + int NNBMAX; + // Initial radius of neighbour sphere [pc], Nbody6 only + double RS0; + // Tidal radius [pc] + double rtide; + // Angular velocity of cluster around the galaxy + double omega; + // Virial radius [pc] + double rvir; + // For CoM correction + double cmr[7]; + // King-, Plummer radius + double rking, rplummer; + // most massive star + double MMAX; + // time-scale factor + double tscale; + // kinetic energy + double ekin = 0.0; + // potential energy + double epot = 0.0; + // velocity dispersion + double sigma = 0.0; + // KZ(22) parameter (binaries yes/no) + int bin; + // (evolved stellar population yes/no) + int sse; + // mass function parameters for mfunc = 2 + double submass[MAX_AN], subcount[MAX_AN], norm[MAX_AN], N_tmp, M_tmp; + // actual 2D/3D half-mass radius of the model + double Rh2D, Rh3D; + + // set half-mass radius temporarily to scale radius for computation of + // escape velocity + if (profile == 3) + Rh = a; + + if ((Mcl) && (N)) { + printf("\nWARNING:" + "specify either Mcl (-M) or N (-N)!\n\n"); + exit (1); + } else if ((!Mcl) && (mfunc == 3)) { + printf("\nWARNING:" + "specify Mcl (-M) when using optimal sampling (-f 3)!\n\n"); + exit (1); + } + + if (xx > 0) { + if ((xx == 4) && (extgas[2])) { + extmass = extgas[0]; + extrad = extgas[1]; + extdecay = extgas[2]; + extstart = extgas[3]; + } else { + printf("\nWARNING: Insufficient or incorrect parameters specified" + " for external Plummer potential!\n\n"); + exit (1); + } + } + + + /*********************** + * Generate star array * + ***********************/ + + int columns = 15; + double **star; + star = (double **)calloc(NMAX,sizeof(double *)); + for (j=0;j epoch\n"); + while (Lifetime(MMAX) < sqrt(pow(epoch,2))) { + MMAX -= 0.01; + } + } + + + /******************* + * Generate masses * + *******************/ + + printf("\n\n-----GENERATE MASSES----- \n"); + + // This loop has to be called multiple times with different N (or Mcl), + // Z and epoch in order to create multiple stellar populations make sure + // that epoch and Z are stored together with the stars, then concatenate + // all arrays for now, all Z and epochs are the same, + // i.e. a single stellar population + if (mfunc == 1) { + printf("\nMaximum stellar mass set to: %.2f\n", MMAX); + generate_m1(&N, star, mlow, mup, &M, &mmean, MMAX, Mcl, epoch, Z, Rh, + remnant); + } else if (mfunc == 2) { + + if (mn) { + for (i = mn+1; i < MAX_MN; i++) mlim[i] = 0.0; + } else { + for (i=0; i= mn) + an = mn - 1; + + mn = an + 1; + + if (!mn){ + printf("\nError: at least one mass limit has to be specified\n"); + return 1; + } else if (mn == 1) { + single_mass = mlim[0]; + printf("\nSetting stellar masses to %g solar mass\n",single_mass); + + if (!N) + N = Mcl/single_mass; + + for (j=0;j= 0; i--) { + norm[i] = norm[i+1] * pow(mlim[i+1], alpha[i+1] - alpha[i]); + subcount[i] = norm[i] * subint(mlim[i], mlim[i+1], + alpha[i] + 1.); + N_tmp += subcount[i]; + submass[i] = norm[i] * subint(mlim[i], mlim[i+1], + alpha[i] + 2.); + M_tmp += submass[i]; + } + + generate_m2(an, mlim, alpha, Mcl, M_tmp, subcount, &N, &mmean, &M, + star, MMAX, epoch, Z, Rh, remnant); + } + } else if (mfunc == 3) { + printf("\nMaximum stellar mass set to: %.2f\n",MMAX); + generate_m3(&N, star, mlow, mup, &M, &mmean, MMAX, Mcl); + randomize(star, N); + } else if (mfunc == 4) { + printf("\nMaximum stellar mass set to: %.2f\n", MMAX); + printf("\nUsing L3 IMF (Maschberger 2012)\n"); + generate_m4(&N, star, alpha_L3, beta_L3, mu_L3, mlow, mup, &M, &mmean, + MMAX, Mcl, epoch, Z, Rh, remnant); + } else { + printf("\nSetting stellar masses to %.1f solar mass\n", single_mass); + if (!N) N = Mcl/single_mass; + for (j=0;j %.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); + } + + order(star, N, M, msort, pairing); + } else { + randomize(star, N); + printf("\nApplying random pairing.\n"); + } + N -= nbin; + for (j=0;j mhighest) + mhighest = star[i][0]; + } + mmeancom /= N; + // number of necessary pos & vel pairs for Baumgardt et al. (2008) + // mass segregation routine + int Nseg = ceil(N*mmeancom/mlowest); + int Nunseg = N; + + double *Mcum; + Mcum = (double *)calloc(N,sizeof(double)); + + // sort masses when mass segregation parameter > 0 + if ((S) && !(profile == 2)) { + printf("\nApplying mass segregation with S = %f\n",S); + segregate(star, N, S); + + //calculate cumulative mass function Mcum + for (i=0;i0.5) + p[3] = gamma[0]; + //inner power-law slope (0.0 for EFF template) + p[4] = gamma[1]; + //transition parameter (2.0 for EFF template) + p[5] = gamma[2]; + + generate_profile(N, star, Rmax, M, p, &Rh, D, symmetry); + printf("\nRh = %.1f pc\n", Rh); + } else if (profile == -1) { + printf("\nGenerating fractal distribution with parameters: " + "N = %i\t Rh = %.3f\t D = %.2f\n", N, Rh, D); + fractalize(D, N, star, 0, symmetry); + rvir = Rh; + } else { + printf("\nGenerating Plummer model with parameters: " + "N = %i\t Rh = %.3f\t D = %.2f\n", N, Rh, D); + rvir = Rh/0.772764; + rplummer = Rh/1.305; + generate_plummer(N, star, rtide, rvir, D, symmetry); + } + + // Apply Baumgardt et al. (2008) mass segregation + if (!(profile == 2) && (S)) { + double **m_temp; + m_temp = (double **)calloc(Nunseg,sizeof(double *)); + + for (j=0;jtphys no data will be stored + double dtp = epoch+1; + // metallicity + double z = Z; + // metallicity parameters + double zpars[20]; + // kick velocity for compact remnants + double vkick; + // escape velocity of cluster + double vesc; + // number of ejected compact remnants + int lostremnants = 0; + // mass of ejected compact remnants + double lostremnantsmass = 0.0; + + if (Mcl && remnant && (Rh != -1)) { + vesc = sqrt(2.0*G*Mcl/Rh); + printf("Escape velocity of cluster = %.4f km/s\n", vesc); + } else if ((!remnant) || (Rh == -1)) { + vesc = 1.0E10; + printf("Keeping all compact remnants\n"); + } else { + vesc = sqrt(2.0*0.4**N/Rh); + printf("Estimated escape velocity of cluster assuming mean stellar " + "mass of 0.4 Msun = %.4f km/s\n", vesc); + } + + for (i=0; i<20; i++) + zpars[i] = 0; + + // get metallicity parameters + zcnsts_(&z,zpars); + + printf("\nSetting up stellar population with Z = %.4f.\n",Z); + + if (epoch) + printf("\nEvolving stellar population for %.1f Myr.\n",epoch); + + // set up mass function parameters + ty = 2; + alpha1 = 1.3; + alpha2 = 2.3; + + c1 = 1.0-alpha1; + c2 = 1.0-alpha2; + + k1 = 2.0/c1*(pow(0.5,c1)-pow(mlow,c1)); + if (mlow>0.5) { + k1 = 0; + k2 = 1.0/c2*(pow(mup,c2)-pow(mlow,c2)); + } else + k2 = k1 + 1.0/c2*(pow(mup,c2)-pow(0.5,c2)); + if (mup<0.5) { + k1 = 2.0/c1*(pow(mup,c1)-pow(mlow,c1)); + k2 = k1; + } + + // determine theoretical mean mass from mass function + c1 = 2.0-alpha1; + c2 = 2.0-alpha2; + + if (mlow != mup) { + int mlow_c1 = pow(mlow, c1); + int mlow_c2 = pow(mlow, c2); + int mup_c1 = pow(mup, c1); + int mup_c2 = pow(mup, c2); + + if (mlow>0.5) { + mth = (1.0/c2*(mup_c2-mlow_c2))/k2; + } else if (mup<0.5) { + mth = (2.0/c1*(mup_c1-mlow_c1))/k2; + } else + mth = (2.0/c1*(pow(0.5,c1)-mlow_c1)+1.0/c2*(mup_c2-pow(0.5,c2)))/k2; + } else { + mth = mlow; + } + + if (!*N) { + *N = max(floor((Mcl-MMAX)/mth), 1); + if (!epoch) printf("Estimated number of necessary stars: %i\n", *N); + *N = 1; + } + + c1 = 1.0-alpha1; + c2 = 1.0-alpha2; + *mmean = 0.0; + *M = 0.0; + double mostmassive = 0.0; + + for (i=0; i<*N; i++) { + do{ + do { + xx = drand48(); + if (xx MMAX); + + // evolve star for deltat = epoch with SSE + // (Hurley, Pols & Tout 2002) + tphys = 0.0; + kw = 1; + // initial mass + mass = star[i][0]; + // actual mass + mt = mass; + ospin = 0.0; + tphysf = epoch; + epoch1 = 0.0; + vkick = 0.0; + //printf("MASS %.2f", mass); + star[i][7] = mass; + + evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, + &epoch1, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick); + // printf("-> %.2f\n", mass); + // if (vkick) printf("KICK: %.5f\n", vkick); + lostremnants++; + lostremnantsmass += mt; + } while (vkick > vesc); + + lostremnants--; + lostremnantsmass -= mt; + + star[i][0] = mt; + star[i][8] = kw; + star[i][9] = epoch1; + star[i][10] = ospin; + star[i][11] = r; + star[i][12] = lum; + + if (star[i][0] > mostmassive) + mostmassive = star[i][0]; + + *M += star[i][0]; + + if ((i==*N-1) && (*Mtphys no data will be stored + double dtp = epoch+1; + // metallicity + double z = Z; + // metallicity parameters + double zpars[20]; + // kick velocity for compact remnants + double vkick; + // escape velocity of cluster + double vesc; + // number of ejected compact remnants + int lostremnants = 0; + // mass of ejected compact remnants + double lostremnantsmass = 0.0; + + if (Mcl && remnant && (Rh != -1)) { + vesc = sqrt(2.0*G*Mcl/Rh); + printf("Escape velocity of cluster = %.4f km/s\n", vesc); + } else if ((!remnant) || (Rh == -1)) { + vesc = 1.0E10; + printf("Keeping all compact remnants\n"); + } else { + vesc = sqrt(2.0*0.4**N/Rh); + printf("Estimated escape velocity of cluster assuming mean stellar " + "mass of 0.4 Msun = %.4f km/s\n", vesc); + } + + for (i=0; i<20; i++) + zpars[i] = 0; + + // get metallicity parameters + zcnsts_(&z,zpars); + + printf("\nSetting up stellar population with Z = %.4f.\n",Z); + + if (epoch) + printf("\nEvolving stellar population for %.1f Myr.\n", epoch); + + double tmp, ml, mup; + double mostmassive = 0.0; + *mmean = 0.0; + *M = 0.0; + + if (!*N) + *N = 1; + + for (i = 0; i < an; i++) { + printf("# <%.2f , %.2f> .. %.2f\n", mlim[i], mlim[i+1], alpha[i]); + } + + for (i = 1; i < an; i++) + subcount[i] += subcount[i-1]; + + for (i = 0; i < *N; i++) { + do { + do { + tmp = drand48() * subcount[an-1]; + for (j = 0; (j < an) && (subcount[j] < tmp); j++); + if (alpha[j] != -1.) { + ml = pow(mlim[j], 1. + alpha[j]); + mup = pow(mlim[j+1], 1. + alpha[j]); + } else { + ml = log(mlim[j]); + mup = log(mlim[j+1]); + } + tmp = ml + drand48() * (mup - ml); + + if (alpha[j] != -1.) + star[i][0] = pow(tmp, 1. / (1. + alpha[j])); + else + star[i][0] = exp(tmp); + + } while (star[i][0] > MMAX); + //printf("%8.4f\n", star[i][0]); + + // evolve star for deltat = epoch with SSE + // (Hurley, Pols & Tout 2002) + tphys = 0.0; + kw = 1; + // initial mass + mass = star[i][0]; + // actual mass + mt = mass; + ospin = 0.0; + vkick = 0.0; + tphysf = epoch; + epoch1 = 0.0; + //printf("MASS %.2f", mass); + star[i][7] = mass; + if (epoch) { + evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, + &ospin, &epoch1, &tms, &tphys, &tphysf, &dtp, &z, zpars, + &vkick); + } + //printf("-> %.2f\n", mass); + //printf("KICK: %.5f\n", vkick); + lostremnants++; + lostremnantsmass += mt; + } while (vkick > vesc); + + lostremnants--; + lostremnantsmass -= mt; + + star[i][0] = mt; + star[i][8] = kw; + star[i][9] = epoch1; + star[i][10] = ospin; + star[i][11] = r; + star[i][12] = lum; + + if (star[i][0] > mostmassive) + mostmassive = star[i][0]; + + *M += star[i][0]; + + if ((i==*N-1) && (*M mb) { + m = pow(pow(star[i-1][0], 2.0-a2) - star[i-1][0]/k2*(2.0-a2), + 1.0/(2.0-a2)); + + if (m < mb) { + m = pow(pow(mb, 2.0-a1) - (2.0-a1)/pow(mb, a1-a2)* + (star[i-1][0]/k2 + pow(mb, 2.0-a2)/(2.0-a2) - 1.0/(2.0-a2)* + (pow(star[i-1][0], 2.0-a2))), 1.0/(2.0-a1)); + } + } else { + m = pow((a1-2.0)/k1*star[i-1][0] + + pow(star[i-1][0], 2.0-a1), 1.0/(2.0-a1)); + } + + if (mtphys no data will be stored + double dtp = epoch+1; + // metallicity + double z = Z; + // metallicity parameters + double zpars[20]; + // kick velocity for compact remnants + double vkick; + // escape velocity of cluster + double vesc; + // number of ejected compact remnants + int lostremnants = 0; + // mass of ejected compact remnants + double lostremnantsmass = 0.0; + double mostmassive = 0.0; + + double G_l; + double G_u; + double u; + double x; + int i; + double mth; + + double gamma; + double mpeak; + // variables for the mean + double B_l; + double B_u; + double a; + double b; + double beta_log; + int ifault; + + printf("\nL3 Parameters:\n"); + printf("alpha = %1.2f beta = %1.2f mu = %1.2f m_low = %3.3f m_up = " + "%3.3f\n",alpha,beta,mu,mlow,mup); + + mpeak = mu*pow((beta-1.0),1.0/(alpha-1.0)); + gamma = alpha + beta*(1.0-alpha); + + printf("'Peak' mass = %3.3f Effective low mass exponent gamma = " + "%1.2f\n\n",mpeak,gamma); + + // normalisation + G_l = pow(1.0 + pow(mlow/mu, 1.0-alpha) , 1.0-beta); + G_u = pow(1.0 + pow(mup/mu, 1.0-alpha) , 1.0-beta); + + + // Calculate the mean mass + a = (2.0-alpha)/(1.0-alpha); + b = beta - a; + + // logarithm of the beta function necessary to calculate + // the incomplete beta fn + beta_log = alogam( a, &ifault ) + alogam( b, &ifault ) - + alogam( a + b, &ifault ); + + x = pow(mlow/mu,1.0-alpha); + x = x/(1.0+x); + B_l = betain(x, a, b, beta_log, &ifault)*exp(beta_log); + + x = pow(mup/mu,1.0-alpha); + x = x/(1.0+x); + B_u = betain(x, a, b, beta_log, &ifault)*exp(beta_log); + mth = mu*(1.0-beta)*(B_u - B_l)/(G_u - G_l) ; + + if (Mcl && remnant && (Rh != -1)) { + vesc = sqrt(2.0*G*Mcl/Rh); + printf("Escape velocity of cluster = %.4f km/s\n", vesc); + } else if ((!remnant) || (Rh == -1)) { + vesc = 1.0E10; + printf("Keeping all compact remnants\n"); + } else { + vesc = sqrt(2.0*0.4**N/Rh); + printf("Estimated escape velocity of cluster assuming mean stellar " + "mass of 0.4 Msun = %.4f km/s\n", vesc); + } + for (i=0; i<20; i++) + zpars[i] = 0; + + zcnsts_(&z,zpars); //get metallicity parameters + + printf("\nSetting up stellar population with Z = %.4f.\n",Z); + + if (epoch) printf("\nEvolving stellar population for %.1f Myr.\n",epoch); + + if (!*N) { + *N = max(floor((Mcl-mup)/mth), 1); + if (!epoch) + printf("Estimated number of necessary stars: %i\n", *N); + *N = 1; + } + + // generate random masses + for (i=0; i<*N; i++) { + do { + u = drand48(); + x = u*(G_u-G_l)+G_l; + x = pow(x,1/(1-beta)) - 1.0; + star[i][0] = mu*pow(x,1.0/(1.0-alpha)); + + tphys = 0.0; + kw = 1; + // initial mass + mass = star[i][0]; + // actual mass + mt = mass; + ospin = 0.0; + tphysf = epoch; + epoch1 = 0.0; + vkick = 0.0; + //printf("MASS %.2f", mass); + star[i][7] = mass; + evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, + &epoch1, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick); + //printf("-> %.2f\n", mass); + //if (vkick) printf("KICK: %.5f\n", vkick); + lostremnants++; + lostremnantsmass += mt; + } while (vkick > vesc); + + lostremnants--; + lostremnantsmass -= mt; + + star[i][0] = mt; + star[i][8] = kw; + star[i][9] = epoch1; + star[i][10] = ospin; + star[i][11] = r; + star[i][12] = lum; + + if (star[i][0] > mostmassive) + mostmassive = star[i][0]; + + *M += star[i][0]; + + if ((i==*N-1) && (*M=3.0) { + for (i=0;ircut); + + // velocities + do { + a[4] = drand48(); + a[5] = drand48(); + a[6] = pow(a[4],2)*pow(1.0 - pow(a[4],2),3.5); + } while (0.1*a[5]>a[6]); + + a[8] = a[4]*sqrt(2.0)/pow(1.0 + ri*ri,0.25); + a[6] = drand48(); + a[7] = drand48(); + + star[i][6] = (1.0 - 2.0*a[6])*a[8]; + star[i][4] = sqrt(a[8]*a[8] - pow(star[i][6],2))*cos(TWOPI*a[7]); + star[i][5] = sqrt(a[8]*a[8] - pow(star[i][6],2))*sin(TWOPI*a[7]); + } + } else { + double xcut; + xcut = 0.000; + while (1.0/sqrt(pow(xcut,-2.0/3.0) - 1.0)<=rcut) xcut+=0.00001; + + fractalize(D, N, star, 1, symmetry); + for (i=0;i12.0) { + printf("W0 too large\n"); + return 0; + } else if (W0 < 0.2) { + printf("W0 too small\n"); + return 0; + } + + /*************************** + * INTERPOLATE KING VALUES * + ***************************/ + + // step size for interpolation + h = pow(W0-rmin,2)/(1.0*M-1.0); + // central density + den = densty(W0); + + //x is King's W, y[0] is z**2, y[1] is 2*z*dz/dW, where z is scaled radius + //so x(y[0]) is energy as function of radius^2 and x(y[1]) is the derivative + + xstart = 0.000001; + ystart0 = 2.0*xstart/3.0; + ystart1 = -2.0/3.0; + + x1 = W0 - xstart; + x2 = 0.0; + + // integrate Poisson's eqn + printf("Integrating Poisson's equation\n"); + + odeint(ystart0, ystart1, x1, x2, den, &kount, xp, yp, M, KMAX); + + printf("No of integration steps = %i\n",kount); + + // interpolate yking + for (k=0;k xp[0]) { + yking[k][0] = (W0-x[k])*yp[0][0]/(W0 - xp[0]); + yking[k][1] = yp[0][1]; + printf("+"); + } else { + i = 0; + + do { + if ((x[k]-xp[i])*(x[k]-xp[i+1]) <= 0.0) { + yking[k][0] = yp[i][0] + (yp[i+1][0]-yp[i][0])* + (x[k]-xp[i])/(xp[i+1]-xp[i]); + yking[k][1] = yp[i][1] + (yp[i+1][1]-yp[i][1])* + (x[k]-xp[i])/(xp[i+1]-xp[i]); + goto jump; + } else { + i++; + } + } while (i= kount) { + yking[k][0] = yp[kount][0]; + yking[k][1] = yp[kount][1]; + } + + // get mass as a function of radius + mass[k] = -2.0*pow(yking[k][0],1.5)/yking[k][1]; + + + //calculate total potential energy + if (k == 0) { + pot = -0.5*pow(mass[k],2)/sqrt(yking[k][0]); + } else { + pot -= 0.5*(pow(mass[k],2) - pow(mass[k-1],2))/(0.5 + *(sqrt(yking[k][0]) + sqrt(yking[k-1][0]))); + } + } + + /****************************** + * DETERMINE BASIC QUANTITIES * + ******************************/ + + // Total mass + totmas = -2.0*pow(yking[M-2][0],1.5)/yking[M-2][1]; + + // Half-mass radius + k=0; + + do { + k++; + } while (mass[k] < 0.5*totmas); + + zh = yking[k][0] - (yking[k][0]-yking[k-1][0])* + (mass[k]-0.5*totmas)/(mass[k]-mass[k-1]); + + *rh = sqrt(zh); + + // Virial radius and King radius + *rvir = -pow(totmas,2)/(2.0*pot); + *rking = sqrt(yking[M-2][0]); + + //Central velocity dispersion**2 (3-dimensional) + ve = sqrt(2.0*W0); + cg1 = (-pow(ve,3)*exp(-pow(ve,2)/2.0) - 3.0*ve*exp(-pow(ve,2)/2.0) + + 3.0/2.0*sqrt(2.0*PI)*erf(ve*sqrt(2.0)/2.0) - pow(ve,5)* + exp(-pow(ve,2)/2.0)/5.0)/(-ve*exp(-pow(ve,2)/2.0) + + sqrt(2.0*PI)*erf(ve*sqrt(2.0)/2.0)/2.0 - pow(ve,3)* + exp(-pow(ve,2)/2.0)/3.0); + + printf("\nTheoretical values:\n"); + printf("Total mass (King units) = %g\n", totmas); + printf("Viriral radius (King units) = %g\n", *rvir); + printf("Half-mass radius (King units) = %g\t(Nbody units) = %g\n", + *rh, *rh/ *rvir); + printf("Edge radius (King units) = %g\t(Nbody units) = %g\n", + *rking, *rking/ *rvir); + printf("Concentration = %g\n", log10(*rking)); + printf("Core radius (King units) = %g\t(Nbody units) = %g\n", + 1.0, 1.0/ *rvir); + printf("3d velocity dispersion**2: %g (central)\t %g (global)\n", + cg1, -pot/totmas); + + /*************************** + * GENERATE STAR POSITIONS * + ***************************/ + + printf("\nGenerating Stars:\n"); + + if (D>=3.0) { + for (i=0;iM) printf("WARNING: failing iteration\n"); + } while (mass[j] <= fmass); + + r2 = yking[j-1][0] + (fmass-mass[j-1])*(yking[j][0] - + yking[j-1][0])/(mass[j]-mass[j-1]); + } + + r = sqrt(r2); + costh = 2.0*drand48()-1.0; + sinth = sqrt(1.0-pow(costh,2)); + phi = 2.0*PI*drand48(); + xstar = r*sinth*cos(phi); + ystar = r*sinth*sin(phi); + zstar = r*costh; + + /**************** + * CHOOSE SPEED * + ****************/ + + if (r < *rking) { + if (r < sqrt(yking[0][0])) { + w1 = x[0]; + w = W0 - (r2/yking[0][0])*(W0 - w1); + } else { + j = 0; + do { + j++; + } while (r > sqrt(yking[j][0])); + wj1 = x[j-1]; + wj = x[j]; + w = wj1 + (r2-yking[j-1][0])*(wj-wj1)/(yking[j][0]-yking[j-1][0]); + } + } else { + printf("radius too big\n"); + } + + vmax = sqrt(2.0*w); + do { + speed = vmax*drand48(); + fstar = pow(speed,2)*(exp(-0.5*pow(speed,2))-exp(-1.0*w)); + } while (fstar < 2.0*drand48()/exp(1.0)); + + costh = 2.0*drand48()-1.0; + phi = 2.0*PI*drand48(); + sinth = sqrt(1.0-pow(costh,2)); + ustar = speed*sinth*cos(phi); + vstar = speed*sinth*sin(phi); + wstar = speed*costh; + mstar = star[i][0]; + + //printf("i: %i\tr=%g\tm=%.5f\tx=%.5f y=%.5f z=%.5f\tvx=%.5f " + // "vy=%.5f vz=%.5f\n",i,sqrt(r2),mstar,xstar,ystar,zstar,ustar, + // vstar,wstar); + + coord[i][0] = xstar; + coord[i][1] = ystar; + coord[i][2] = zstar; + coord[i][3] = ustar; + coord[i][4] = vstar; + coord[i][5] = wstar; + } + + for (i=0;iM) printf("WARNING: failing iteration\n"); + } while (mass[j] <= fmass); + + r2 = yking[j-1][0] + (fmass-mass[j-1])*(yking[j][0] - + yking[j-1][0])/(mass[j]-mass[j-1]); + } + + r = sqrt(r2); + + /**************** + * CHOOSE SPEED * + ****************/ + + if (r < *rking) { + if (r < sqrt(yking[0][0])) { + w1 = x[0]; + w = W0 - (r2/yking[0][0])*(W0 - w1); + } else { + j = 0; + do { + j++; + } while (r > sqrt(yking[j][0])); + wj1 = x[j-1]; + wj = x[j]; + w = wj1 + (r2-yking[j-1][0])*(wj-wj1)/(yking[j][0]- + yking[j-1][0]); + } + } else { + printf("radius too big\n"); + } + + vmax = sqrt(2.0*w); + do { + speed = vmax*drand48(); + fstar = pow(speed,2)*(exp(-0.5*pow(speed,2))-exp(-1.0*w)); + } while (fstar < 2.0*drand48()/exp(1.0)); + + for (k=1;k<4;k++) + star[i][k] *= r/ri; + + for (k=4;k<7;k++) + star[i][k] *= speed; + } + + } + + for (j=0;j= 0.0) { + h = sqrt(pow(H1,2)); + } else { + h = - sqrt(pow(H1,2)); + } //step size + + // z^2 + y[0] = ystart0; + // 2*z*dz/dW where z is scaled radius + y[1] = ystart1; + + xsav = x-DXSAV*2.0; + + // limit integration to MAXSTP steps + for (i=0;i sqrt(pow(DXSAV,2))) { + + if (*kount < KMAX-1) { + xp[*kount] = x; + + for (j=0;j<2;j++) { + yp[*kount][j] = y[j]; + } + + *kount = *kount + 1; + xsav = x; + } + } //store x, y1 and y2 if the difference in x is smaller as the + // desired output step size DXSAV + + if (((x+h-x2)*(x+h-x1)) > 0.0) h = x2-x; + + // do a Runge-Kutta step + rkqc(y,dydx,&x,&h,den,yscal, &hdid, &hnext, TOL); + + if ((x-x2)*(x2-x1) >= 0.0) { + ystart0 = y[0]; + ystart1 = y[1]; + + xp[*kount] = x; + for (j=0;j<2;j++) { + yp[*kount][j] = y[j]; + } + return 0; + *kount = *kount +1; + } + + if (sqrt(pow(hnext,2)) < HMIN) { + printf("Stepsize smaller than minimum.\n"); + return 0; + } + + h = hnext; + } + + return 0; +} + +int derivs(double x, double *y, double *dydx, double den){ + + double rhox; + + if (x >= 0.0) { + rhox =-sqrt(x)*(x+1.5)+0.75*sqrt(PI)*exp(x)*erf(sqrt(x)); + } else { + rhox = 0.0; + } + + dydx[0]= y[1]; + dydx[1] = 0.25*pow(y[1],2)*(6.0+9.0*y[1]*rhox/den)/y[0]; + + return 0; +} + +int rkqc(double *y,double *dydx, double *x, double *h, double den, + double *yscal, double *hdid, double *hnext, double TOL) { + + double safety = 0.9; + double fcor = 0.0666666667; + double errcon = 6.0E-4; + double pgrow = -0.20; + double pshrnk = -0.25; + double xsav; + int i; + double ysav[2], dysav[2], ytemp[2]; + double errmax; + double hh; + + xsav = *x; + + for (i=0;i<2;i++) { + ysav[i] = y[i]; + dysav[i] = dydx[i]; + } + + do { + hh = 0.5**h; + rk4(xsav, ysav, dysav, hh, ytemp, den); + + *x = xsav + hh; + // find derivative + derivs(*x,ytemp,dydx,den); + rk4(*x, ytemp, dydx, hh, y, den); + + *x = xsav + *h; + if (*x == xsav) { + printf("ERROR: Stepsize not significant in RKQC.\n"); + return 0; + } + rk4(xsav, ysav, dysav, *h, ytemp, den); + + errmax = 0.0; + for (i=0;i<2;i++) { + ytemp[i] = y[i] - ytemp[i]; + errmax = max(errmax, sqrt(pow(ytemp[i]/yscal[i],2))); + } + errmax /= TOL; + // if integration error is too large, decrease h + if (errmax > 1.0) + *h = safety**h*(pow(errmax,pshrnk)); + + } while (errmax > 1.0); + + + *hdid = *h; + if (errmax > errcon) { + // increase step size for next step + *hnext = safety**h*(pow(errmax,pgrow)); + } else { + // if integration error is very small increase step size significantly + // for next step + *hnext = 4.0**h; + } + + for (i=0;i<2;i++) { + y[i] += ytemp[i]*fcor; + } + + return 0; + +} + +int rk4(double x, double *y, double *dydx, double h, double *yout, + double den) { + + double hh, h6, xh; + double yt[2], dyt[2], dym[2]; + int i; + + hh=h*0.5; + h6=h/6.0; + xh=x+hh; + + for (i=0;i<2;i++) { + yt[i] = y[i] + hh*dydx[i]; + } + + // find derivative + derivs(xh,yt,dyt,den); + for (i=0;i<2;i++) { + yt[i] = y[i] + hh*dyt[i]; + } + + // find derivative + derivs(xh,yt,dym,den); + for (i=0;i<2;i++) { + yt[i] = y[i] + h*dym[i]; + dym[i] += dyt[i]; + } + + // find derivative + derivs(x+h,yt,dyt,den); + for (i=0;i<2;i++) { + yout[i] = y[i] + h6*(dydx[i]+dyt[i]+2.0*dym[i]); + } + + return 0; +} + +int generate_subr(int N, double S, double **star, double rtide, double rvir) { + + long i, j; + double r1, r2, tmp, rcut, rtemp; + double U_tot, U_tot1, K_tot; + double UT_sub, M_sub; + double maxim, beta; + int N_star; + double M_tot; + star1 = NULL; + star2 = NULL; + + if (S<0.5) { + printf("Setting cut-off radius of Subr model to the approximate " + "tidal radius\n"); + + // cut-off radius for Plummer sphere = tidal radius + rcut = rtide/rvir; + } else { + printf("Attention! Be aware that for about S>0.5 the Subr mass " + "segregation may generate a virially hot system\n" + "Setting cut-off radius of Subr model to twice the approximate " + "tidal radius\n"); + // cut-off radius for Plummer sphere = tidal radius + rcut = 2.0*rtide/rvir; + } + + M_tot = 0.; + N_star = N; + star1 = (struct t_star1 *)realloc(star1, N_star * sizeof(struct t_star1)); + star2 = (struct t_star2 *)realloc(star2, N_star * sizeof(struct t_star2)); + + for (i=0; i 0.) + star1[i].Ui = star1[i].mass * pow(M_sub, -S); + else + star1[i].Ui = star1[i].mass; + + for (j = 0; j < i; j++) { + tmp = star1[i].Ui * star1[j].Ui; + star2[i].UT_add += tmp; + star1[i].U_tmp += tmp; + star1[j].U_tmp += tmp; + }; + U_tot1 += star2[i].UT_add; + }; + + UT_sub = U_tot = K_tot = 0.; + + for (i = 0; i < N_star; i++) { + + if ((i/1000)*1000 == i) + printf("Generating orbit #%li\n", i); + + star2[i].UT_add *= 0.5 / U_tot1; + star2[i].UT = star1[i].U_tmp * 0.5 / U_tot1; + UT_sub += star2[i].UT_add; + + do { + position(i, U_tot, UT_sub, S); + rtemp = sqrt(pow(star1[i].r[0],2)+pow(star1[i].r[1],2)+ + pow(star1[i].r[2],2)); + } while (rtemp>rcut); + U_tot += star1[i].mass * star2[i].U_sub; + }; + + for (i = 0; i < N_star; i++) { + maxim = 0.25 * star1[i].mass / star2[i].UT; + find_beta(&maxim, &beta); + if (beta > 0.) + { + do + { + r1 = maxim * drand48(); + r2 = drand48(); + r2 *= r2; + tmp = 1. - r2; + } + while (r1 > r2 * exp(beta * log(tmp))); + } + else + { + do + { + r1 = maxim * drand48(); + r2 = sin(M_PI_2 * drand48()); + r2 *= r2; + tmp = sqrt(1. - r2); + } + while (r1 > r2 * exp((2.* beta + 1.) * log(tmp))); + }; + + star2[i].E = r2 * (star2[i].U + star2[i].U_sub); + + tmp = sqrt(2. * star2[i].E); + isorand(star2[i].v); + star2[i].v[0] *= tmp; + star2[i].v[1] *= tmp; + star2[i].v[2] *= tmp; + K_tot += star1[i].mass * star2[i].E; + }; + + for (i = 0; i < N_star; i++){ + //printf("%.12f %18.14f %18.14f %18.14f %16.12f %16.12f " + // "%16.12f\n",star1[i].mass, star1[i].r[0], star1[i].r[1], + // star1[i].r[2], star2[i].v[0], star2[i].v[1], star2[i].v[2]); + star[i][0] = star1[i].mass; + star[i][1] = star1[i].r[0]; + star[i][2] = star1[i].r[1]; + star[i][3] = star1[i].r[2]; + star[i][4] = star2[i].v[0]; + star[i][5] = star2[i].v[1]; + star[i][6] = star2[i].v[2]; + star[i][7] = star1[i].m0; + star[i][8] = star1[i].kw; + star[i][9] = star1[i].epoch; + star[i][10] = star1[i].spin; + star[i][11] = star1[i].rstar; + star[i][12] = star1[i].lum; + star[i][13] = star1[i].epochstar; + star[i][14] = star1[i].zstar; + + }; + + return 0; +} + +void quick(int start, int stop) { + int i, j; + double temp, median; + + median = 0.5 * (star1[start].mass + star1[stop].mass); + + i = start; + j = stop; + while (i <= j) + { + while ((i <= stop) && (star1[i].mass > median)) + i++; + + while ((j >= start) && (star1[j].mass < median)) + j--; + + if (j > i) + { + SWAP(i,j); + i++; + j--; + } + else if (i == j) + { + i++; + j--; + }; + }; + + if (start + 1 < j) + quick(start, j); + else if ((start < j) && (star1[start].mass < star1[j].mass)) + SWAP(start, j); + + if (i + 1 < stop) + quick(i, stop); + else if ((i < stop) && (star1[i].mass < star1[stop].mass)) + SWAP(i, stop) +} + +void position(int id, double U_sub, double UT_sub, double S){ + + long i; + double rx, ry, rz, r1; + double rfac; + + rfac = pow(star2[id].M_sub, 2. * S) / (1. - S); + + star2[id].ntry = 0; + + do + { + star2[id].U_sub = 0.; + star2[id].ntry++; + + r1 = drand48(); + r1 = exp(-2. * log(r1) / 3.); + r1 = RSCALE * sqrt(1. / (r1 - 1.)); + + r1 *= rfac; + + isorand(star1[id].r); + star1[id].r[0] *= r1; + star1[id].r[1] *= r1; + star1[id].r[2] *= r1; + star2[id].rad = r1; + + for (i = 0; i < id; i++) + { + rx = star1[i].r[0] - star1[id].r[0]; + ry = star1[i].r[1] - star1[id].r[1]; + rz = star1[i].r[2] - star1[id].r[2]; + r1 = sqrt(rx*rx + ry*ry + rz*rz); + star1[i].U_tmp = star1[id].mass / r1; + star2[id].U_sub += star1[i].mass / r1; + }; + r1 = U_sub + star1[id].mass * star2[id].U_sub; + } + while (fabs(r1 - UT_sub) > 0.1 * (r1 + UT_sub) / sqrt(1. + id)); + + for (i = 0; i < id; i++) star2[i].U += star1[i].U_tmp; +} + +void find_beta(double *avg, double *beta){ + + int i; + double a1, a2, I1, I2, J1, J2; + double alo, ah, Ilo, Ih, Jlo, Jh; + double tmpavg; + + if (*avg > 0.75) + { + *beta = -0.5; + *avg = 1.; + return; + }; + + Ilo = I1 = 3. * M_PI / 16.; + Ih = I2 = 0.2; + Jlo = J1 = M_PI / 4; + Jh = J2 = 1. / 3.; + alo = a1 = -0.5; + ah = a2 = 0.; + + tmpavg = I2 / J2; + i = 0; + while (tmpavg > *avg) + { + if (!(i & 1)) + { + a1 += 1.; + I1 /= 1. + 5. / (2. * a1); + J1 /= 1. + 3. / (2. * a1); + Ilo = I2; Ih = I1; + Jlo = J2; Jh = J1; + alo = a2; ah = a1; + } + else + { + a2 += 1.; + I2 /= 1. + 5. / (2. * a2); + J2 /= 1. + 3. / (2. * a2); + Ilo = I1; Ih = I2; + Jlo = J1; Jh = J2; + alo = a1; ah = a2; + }; + tmpavg = Ih / Jh; + i++; + }; + + *beta = alo + (ah - alo) * (*avg - Ilo / Jlo) / (tmpavg - Ilo / Jlo); + + if (*beta > 0.) + *avg = exp(*beta * log(*beta / (1. + *beta))) / (1. + *beta); + else + *avg = 2. * exp((2.* *beta + 1.) * log((2.* *beta + 1.) / (2.* *beta + + 3.))) / (2.* *beta + 3.); +} + +void isorand(double *r) { + + double rad; + + do + { + r[0] = 2. * drand48() - 1.; + r[1] = 2. * drand48() - 1.; + r[2] = 2. * drand48() - 1.; + rad = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]; + } + while ((rad > 1.) || (rad < 0.0001)); + + rad = sqrt(rad); + r[0] /= rad; + r[1] /= rad; + r[2] /= rad; +} + +// smallest to largest +int cmpmy(double *x1, double *x2) { + if(*x1<*x2) return -1; + return 1; +} + +// largest to smallest +int cmpmy_reverse(double *x1, double *x2) { + if(*x1>*x2) return -1; + return 1; +} + +double generate_profile (int N, double **star, double Rmax, double Mtot, + double *p, double *Rh, double D, int symmetry) { + + double Mnorm; + double r, Rmin = 1E-05; + // 51 steps is more than sufficient (in most cases) + int steps = 51, i, h, j; + double r_array[steps], rho_array[steps], sigma_array[steps], + M_array[steps], X_array[steps]; + double random[7], x,y,z, ri, vx, vy, vz; + double r_norm, v_norm; + + Mnorm = M(Rmax, p); + p[1] = Mtot/Mnorm; + + double stepsize; + stepsize = (log10(Rmax)-log10(Rmin))/(steps-1); + + printf("\nIntegrating profile...\n"); + for (i=0;i=3.0) { + + for (i = 0; i< N; i++){ + + if ((i/1000)*1000 == i) + printf("Generating orbit #%i\n", i); + + do { + random[0] = drand48(); + + j = 0; + while (X_array[j]Rmax); + + double sigma; + + j = 0; + while (r_array[j] arg2) errt = arg1; + else errt = arg2; + if (errt <= *err) { + *err=errt; + ans=a[j][i]; + } + } + if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*(*err)) break; + } + + return ans; +} + +double midexp(double (*func)(double, double*), double aa, double bb, int n, + double *p) { + + double x,tnm,sum,del,ddel,a,b; + static double s; + int it,j; + + b=exp(-aa); + a=0.0; + if (n == 1) { + x = 0.5*(a+b); + return (s=(b-a)*func(-log(x),p)/x); + } else { + for(it=1,j=1;j Elson, Fall & Freeman (1987) + // for p[4] = 0 & p[5] = 2 + return p[1]*pow(2.0,0.5*(p[3]-p[4]))*pow(x/p[2],-p[4])* + pow(1.0+pow(x/p[2],p[5]),-(p[3]-p[4])/p[5]); +} + +double rho(double r, double *p) { + int j = 0; + double s,st,ost,os,stemp; + if (r < p[2]) { + ost = os = -1.0e30; + for (j=1;j<=JMAX;j++) { + if (p[4]) st=midinf(kernel,r,p[2],j,p); + else st=midpnt(kernel,r,p[2],j,p); + s=(4.0*st-ost)/3.0; + if (fabs(s-os) < EPS*fabs(os)) break; + if (s == 0.0 && os == 0.0 && j > 6) break; + os=s; + ost=st; + } + stemp = s; + ost = os = -1.0e30; + for (j=1;j<=JMAX;j++) { + st=midinf(kernel,p[2],BIGNUMBER,j,p); + s=(4.0*st-ost)/3.0; + if (fabs(s-os) < EPS*fabs(os)) break; + if (s == 0.0 && os == 0.0 && j > 6) break; + os=s; + ost=st; + } + s += stemp; + } else { + ost = os = -1.0e30; + for (j=1;j<=JMAX;j++) { + st=midinf(kernel,r,BIGNUMBER,j,p); + s=(4.0*st-ost)/3.0; + if (fabs(s-os) < EPS*fabs(os)) break; + if (s == 0.0 && os == 0.0 && j > 6) break; + os=s; + ost=st; + } + } + + // no negative density! + return max(s * -1.0/PI, 0.0); +} + +double rho_kernel (double x, double *p) { + double s; + s = 4.0*PI*x*x*rho(x,p); + return s; +} + +double M(double r, double *p){ + int j = 0; + double s,st,ost,os; + + ost = os = -1.0e30; + for (j=1;j<=JMAX;j++) { + st=midsql(rho_kernel,0.0,r,j,p); + s=(4.0*st-ost)/3.0; + if (fabs(s-os) < EPS*fabs(os)) break; + if (s == 0.0 && os == 0.0 && j > 6) break; + os=s; + ost=st; + } + return s; + +} + +double sigma_kernel (double x, double *p) { + double s; + s = 1.0*rho(x,p)*M(x,p)/(x*x); + return s; +} + +double sigma(double r, double *p) { + int j = 0; + double s,st,ost,os,t; + + ost = os = -1.0e30; + for (j=1;j<=JMAX;j++) { + st=midexp(sigma_kernel,r,BIGNUMBER,j,p); + s=(4.0*st-ost)/3.0; + if (fabs(s-os) < EPS*fabs(os)) break; + if (s == 0.0 && os == 0.0 && j > 6) break; + os=s; + ost=st; + } + + if (rho(r,p)) t = sqrt(G*s/rho(r,p)); + else t = 0.0; + + return t; + +} + +double get_gauss(void){ + double random[2],p,q; + do { + random[0] = 2.0*drand48()-1.0; + random[1] = 2.0*drand48()-1.0; + q = random[0]*random[0]+random[1]*random[1]; + } while (q>1.0); + + p = sqrt(-2.0*log(q)/q); + return random[0]*p; + +} + +double fractalize(double D, int N, double **star, int radial, int symmetry) { + int i, j, h, Nparent, Nparentlow; + int Ntot = 128.0*pow(8,ceil(log(N)/log(8))); + int Ntotorg = Ntot; + double l = 2.0; + double prob = pow(2.0, D-3.0); + double scatter; + if (radial) scatter = 0.01; + else scatter = 0.1; + double vx, vy, vz; + double vscale; + int subi; + double morescatter = 0.1; //0.1 looks good + + printf("\nFractalizing initial conditions...\n\n"); + + double **star_temp; //temporary array for fractalized structure + star_temp = (double **)calloc(Ntot,sizeof(double *)); + for (j=0;j 1.0)); + for (h=1;h<7;h++) star[i][h] = star_temp[j][h]; + star_temp[j][0] = 0.0; + } + + double r, r_norm, vnorm = 0.0, start[4]; + if (radial) { + for (i=0;i1.0); + for (h=1;h<4;h++) star[i][h] = start[h]; + } + } + for (j=0;j=msort) || (m2*M>=msort)) && (OBperiods)) { + if (adis == 3) { + // Sana et al., (2012); Oh, S., Kroupa, P., & + // Pflamm-Altenburg, J. (2015) + double lPmin = 0.15, alpha = 0.45; + double ralpha = 1/alpha, eta = alpha/0.23; + lP = pow(pow(lPmin,alpha) + eta*drand48(),ralpha); + } else { + // derive from Sana & Evans (2011) period distribution for + // massive binaries + + // parameters of Sana & Evans (2011) period distribution in + // days (eq. 5.1) + double lPmin = 0.3, lPmax = 3.5, lPbreak = 1.0; + // fraction of binaries with periods below Pbreak + double Fbreak = 0.5; + double xperiod = drand48(); + if (xperiod <= Fbreak) + lP = xperiod *(lPbreak-lPmin)/Fbreak + lPmin; + else + lP = (xperiod - Fbreak)*(lPmax-lPbreak)/(1.0-Fbreak) + + lPbreak; + } + + // days + P = pow(10.0,lP); + // yr + P /= 365.25; + + abin = pow((m1+m2)*M*P*P,(1.0/3.0));//AU + abin /= 206264.806;//pc + abin /= rvir;//Nbody units + } else if (adis == 0) { + // flat semi-major axis distribution + if (!i) + printf("\nApplying flat semi-major axis distribution with " + "amin = %g and amax = %g.\n", amin, amax); + if (!i) + amin /= rvir; + if (!i) + amax /= rvir; + abin = amin+drand48()*(amax-amin); + } else if (adis == 1 || adis == 3) { + // derive from Kroupa (1995) period distribution + if (!i) { + printf("\nDeriving semi-major axis distribution from " + "Kroupa (1995) period distribution.\n"); + } + // parameters of Kroupa (1995) period distribution + double Pmin = 10, delta = 45, eta = 2.5; + do { + lP = log10(Pmin) + sqrt(delta*(exp(2.0*drand48()/eta)-1)); + } while (lP > 8.43); + + P = pow(10.0,lP);//days + P /= 365.25;//yr + + abin = pow((m1+m2)*M*P*P,(1.0/3.0));//AU + abin /= 206264.806;//pc + abin /= rvir;//Nbody units + } else { + // derive from Duquennoy & Mayor (1991) period distribution + + // mean of Duquennoy & Mayor (1991) period distribution + // [log days] + lPmean = 4.8; + // full width half maximum of Duquennoy & Mayor (1991) + // period distribution [log days] + lPsigma = 2.3; + if (!i) { + printf("\nDeriving semi-major axis distribution from " + "Duquennoy & Mayor (1991) period distribution.\n"); + } + do { + //generate two random numbers in the interval [-1,1] + u1 = 2.0*drand48()-1.0; + u2 = 2.0*drand48()-1.0; + + //combine the two random numbers + q = u1*u1 + u2*u2; + } while (q > 1.0); + + p = sqrt(-2.0*log(q)/q); + x1 = u1*p; + x2 = u2*p; + lP1 = lPsigma*x1 + lPmean; + lP2 = lPsigma*x2 + lPmean; + + P = pow(10,lP1);//days + P /= 365.25;//yr + + abin = pow((m1+m2)*M*P*P,(1.0/3.0));//AU + abin /= 206264.806;//pc + abin /= rvir;//Nbody units + } + + + if (!i && OBperiods && msort) { + if (adis == 3) { + 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); + } + } + + // 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"); + } + else if (!i) { + printf("\nApplying thermal eccentricity distribution.\n"); + } + + if (((m1*M>=msort) || (m2*M>=msort)) && (OBperiods)) { + double k1, k2; + // maximum eccentricity for high-mass binaries + double elim = 0.8; + // fraction of circular orbits among high-mass binaries + double fcirc = 0.3; + 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 { + // Thermal distribution f(e)=2e + ecc = sqrt(drand48()); + } + //ecc = 0.0; // all circular + + //Apply Kroupa (1995) eigenevolution + eccold = ecc; + abinold = abin; + m1old = m1; + m2old = m2; + + if (eigen) { + if (!i) { + printf("\nApplying Kroupa (1995) eigenevolution for " + "short-period binaries\n"); + } + // temporary re-scaling + m1*=M; + m2*=M; + abin*=rvir; + + if (!(OBperiods && ((m1>=msort) || (m2>=msort)))) { + if (m1>=m2) + eigenevolution(&m1, &m2, &ecc, &abin); + else + eigenevolution(&m2, &m1, &ecc, &abin); + } + + m1/=M; + m2/=M; + abin/=rvir; + } + + // Apply Binary Star Evolution (Hurley, Tout & Pols 2002) + if (BSE) { + // stellar type + int kw[2] = {1, 1}; + // initial mass + double mass[2]; + // actual mass + double mt[2]; + // radius + double r[2] = {0.0, 0.0}; + // luminosity + double lum[2] = {0.0, 0.0}; + // core mass + double mc[2] = {0.0, 0.0}; + // core radius + double rc[2] = {0.0, 0.0}; + // envelope mass + double menv[2] = {0.0, 0.0}; + // envelope radius + double renv[2] = {0.0, 0.0}; + // spin + double ospin[2] = {0.0, 0.0}; + // time spent in current evolutionary state + double epoch1[2] = {0.0, 0.0}; + // main-sequence lifetime + double tms[2] = {0.0, 0.0}; + // initial age + double tphys = 0.0; + // final age + double tphysf = epoch; + // data store value, if dtp>tphys no data will be stored + double dtp = epoch+1; + // metallicity + double z = Z; + + // initial mass of primary + mass[0] = m1*M; + // initial mass of secondary + mass[1] = m2*M; + // actual mass + mt[0] = mass[0]; + // actual mass + mt[1] = mass[1]; + vkick[0] = 0.0; + vkick[1] = 0.0; + + abin *= rvir; //pc + abin *= 206264.806; //AU + P = sqrt(pow(abin, 3.0)/(m1+m2));//yr + P *= 365.25;//days + + evolv2_(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin,epoch1,tms, + &tphys,&tphysf,&dtp,&z,zpars,&P,&eccold, vkick); + + P /= 365.25;//yr + + abin = pow((m1+m2)*P*P,(1.0/3.0));//AU + abin /= 206264.806;//pc + abin /= rvir;//Nbody units + + // Nbody units + m1 = mt[0]/M; + // Nbody units + m2 = mt[1]/M; + + star[2*i][8] = kw[0]; + star[2*i+1][8] = kw[1]; + star[2*i][9] = epoch1[0]; + star[2*i+1][9] = epoch1[1]; + star[2*i][10] = ospin[0]; + star[2*i+1][10] = ospin[1]; + star[2*i][11] = r[0]; + star[2*i+1][11] = r[1]; + star[2*i][12] = lum[0]; + star[2*i+1][12] = lum[1]; + } + + //pos & vel in binary frame + ea = rtnewt(ecc, drand48()); + rop[0] = abin*(cos(ea) - ecc); + rop[1] = abin*sqrt(1.0-ecc*ecc)*sin(ea); + + mm = sqrt((m1+m2)/pow(abin,3)); + eadot = mm/(1.0 - ecc*cos(ea)); + vop[0] = -abin*sin(ea)*eadot; + vop[1] = abin*sqrt(1.0-ecc*ecc)*cos(ea)*eadot; + + //Convert to cluster frame + cosi = 2.0*drand48()-1.0; + inc = acos(cosi); + node = 2.0*PI*drand48(); + peri = 2.0*PI*drand48(); + + pmat[0][0] = cos(peri)*cos(node) - sin(peri)*sin(node)*cosi; + pmat[1][0] = cos(peri)*sin(node) + sin(peri)*cos(node)*cosi; + pmat[2][0] = sin(peri)*sin(inc); + pmat[0][1] = -sin(peri)*cos(node) - cos(peri)*sin(node)*cosi; + pmat[1][1] = -sin(peri)*sin(node) + cos(peri)*cos(node)*cosi; + pmat[2][1] = cos(peri)*sin(inc); + + for (j=0;j<3;j++) { + rrel[j] = pmat[j][0]*rop[0] + pmat[j][1]*rop[1]; + vrel[j] = pmat[j][0]*vop[0] + pmat[j][1]*vop[1]; + } + + for (j=0;j<3;j++) { + // Star2 pos + star[2*i+1][j+1] = star[2*i][j+1] + m1/(m1+m2)*rrel[j]; + // Star2 vel + star[2*i+1][j+4] = star[2*i][j+4] + m1/(m1+m2)*vrel[j]; + // Star1 pos + star[2*i][j+1] -= m2/(m1+m2)*rrel[j]; + // Star1 vel + star[2*i][j+4] -= m2/(m1+m2)*vrel[j]; + } + + star[2*i][0] = m1; + star[2*i+1][0] = m2; + + } while ((sqrt(vkick[0]) > vesc) && (sqrt(vkick[1]) > vesc)); + + //printf("%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf" + // "\t%lf\n", m1/(m1+m2)*rrel[0], m1/(m1+m2)*rrel[1], + // m1/(m1+m2)*rrel[2], m1/(m1+m2)*vrel[0], m1/(m1+m2)*vrel[1], + // m1/(m1+m2)*vrel[2], -m2/(m1+m2)*rrel[0], -m2/(m1+m2)*rrel[1], + // -m2/(m1+m2)*rrel[2], -m2/(m1+m2)*vrel[0], -m2/(m1+m2)*vrel[1], + // -m2/(m1+m2)*vrel[2]); + + //printf("%i\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t" + // "%.10lf\t%.10lf\t%.10lf\n",i,m1*M, m2*M, P, abin*rvir, ecc, + // m1old*M, m2old*M, abinold*rvir, eccold); + + //printf("%i\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t" + // %.10lf\t%.10lf\n",i,star[2*i][4], star[2*i+1][4], star[2*i][1], + // star[2*i+1][1], star[2*i][1]-star[2*i+1][1], m2old*M, + // abinold*rvir, eccold); + + if (!star[2*i][0] || !star[2*i+1][0]) { + // one binary less because one component went supernova + nbin--; + + // temporarily save surviving companion + double startemp[15]; + if (!star[2*i][0]) { + for (k=0;k<15;k++) + startemp[k] = star[2*i+1][k]; + } else { + for (k=0;k<15;k++) + startemp[k] = star[2*i][k]; + } + + for (j=2*i+2;j<*N;j++) { + for (k=0;k<15;k++) { + //move all remaining stars two positions up in the array + star[j-2][k] = star[j][k]; + } + } + + // reduce total number of stars by one, i.e. remove massless + // supernova remnant from computations + *N = *N-1; + + for (k=0;k<15;k++) { + //make surviving companion the last particle in array + star[*N-1][k] = startemp[k]; + } + + //reduce binary counter + i--; + } + } + + return 0; +} + +// largest up +void shellsort(double **array, int N, int k) { + int i,j,l,n; + N = N-1; + double swap[k]; + //guess distance n + for (n = 1; n <= N/9; n = 3*n+1); + for (; n > 0; n /= 3) { + for (i = n; i <= N; i++) { + for (l=0; l= n) && (array[j-n][0] < swap[0])); j -= n) { + for (l=0; l 0; n /= 3) { + for (i = n; i <= N; i++) { + for (l=0; l= n) && (array[j-n][0] > swap[0])); j -= n) { + for (l=0; l 0; n /= 3) { + for (i = n; i <= N; i++) { + swap = array[i]; + for (j = i; ((j >= n) && (array[j-n] < swap)); j -= n) { + array[j] = array[j-n]; + } + array[j] = swap; + } + } +} + +//smallest up +void shellsort_reverse_1d(double *array, int N) { + int i,j,n; + N = N-1; + double swap; + //guess distance n + for (n = 1; n <= N/9; n = 3*n+1); + for (; n > 0; n /= 3) { + for (i = n; i <= N; i++) { + swap = array[i]; + for (j = i; ((j >= n) && (array[j-n] > swap)); j -= n) { + array[j] = array[j-n]; + } + array[j] = swap; + } + } +} + +int order(double **star, int N, double M, double msort, int pairing){ + + int i,j; + int Nhighmass; + int columns = 15; + double **star_temp; + star_temp = (double **)calloc(N,sizeof(double *)); + for (j=0;j= msort) { + // First star in binarie + int m_tmp = masses[i][1]; + star_temp[j][0] = star[m_tmp][0]; + star_temp[j][1] = star[m_tmp][1]; + star_temp[j][2] = star[m_tmp][2]; + star_temp[j][3] = star[m_tmp][3]; + star_temp[j][4] = star[m_tmp][4]; + star_temp[j][5] = star[m_tmp][5]; + star_temp[j][6] = star[m_tmp][6]; + star_temp[j][7] = star[m_tmp][7]; + star_temp[j][8] = star[m_tmp][8]; + star_temp[j][9] = star[m_tmp][9]; + star_temp[j][10] = star[m_tmp][10]; + star_temp[j][11] = star[m_tmp][11]; + star_temp[j][12] = star[m_tmp][12]; + star_temp[j][13] = star[m_tmp][13]; + star_temp[j][14] = star[m_tmp][14]; + mask[i] = 1; + j++; + // Find the second one based on uniform mass ratio + if (i mpair) { + k++; + if (k>=N) { + k--; + break; + } + } + int k1 = k-1; + int k2 = k; + if (k1 == i) { + while(k1i && mask[k1]==1) + k1--; + + while(k20) { + if(j>=N) { + perror("\n Error!: Pairing binary star exceed " + "the maximum index!\n"); + exit(0); + } + int k_tmp = masses[k][1]; + // Second star in binarie + star_temp[j][0] = star[k_tmp][0]; + star_temp[j][1] = star[k_tmp][1]; + star_temp[j][2] = star[k_tmp][2]; + star_temp[j][3] = star[k_tmp][3]; + star_temp[j][4] = star[k_tmp][4]; + star_temp[j][5] = star[k_tmp][5]; + star_temp[j][6] = star[k_tmp][6]; + star_temp[j][7] = star[k_tmp][7]; + star_temp[j][8] = star[k_tmp][8]; + star_temp[j][9] = star[k_tmp][9]; + star_temp[j][10] = star[k_tmp][10]; + star_temp[j][11] = star[k_tmp][11]; + star_temp[j][12] = star[k_tmp][12]; + star_temp[j][13] = star[k_tmp][13]; + star_temp[j][14] = star[k_tmp][14]; + mask[k] = 1; + j++; + } else { + printf("Mass %g did not find good pair, expected pair " + "mass: %g\n",masses[i][0],mpair); + } + } + } + else { + // Store index for remaining stars + mmrand[ileft][0] = drand48(); + mmrand[ileft][1] = masses[i][1]; + ileft++; + } + } + } + + // Error check, ileft+j should be N + if(ileft+j!=N) { + fprintf(stderr,"\n Error! M>Msort binaries + MMsort: %d, MMsort: %d",j/2); + } + + // Randomize indexs for random pairing + shellsort(mmrand,ileft,2); + + // Random pairing remaining binaries + for(i=0;i= msort) { + star_temp[j][0] = star[(int) masses[i][1]][0]; + star_temp[j][1] = star[(int) masses[i][1]][1]; + star_temp[j][2] = star[(int) masses[i][1]][2]; + star_temp[j][3] = star[(int) masses[i][1]][3]; + star_temp[j][4] = star[(int) masses[i][1]][4]; + star_temp[j][5] = star[(int) masses[i][1]][5]; + star_temp[j][6] = star[(int) masses[i][1]][6]; + star_temp[j][7] = star[(int) masses[i][1]][7]; + star_temp[j][8] = star[(int) masses[i][1]][8]; + star_temp[j][9] = star[(int) masses[i][1]][9]; + star_temp[j][10] = star[(int) masses[i][1]][10]; + star_temp[j][11] = star[(int) masses[i][1]][11]; + star_temp[j][12] = star[(int) masses[i][1]][12]; + star_temp[j][13] = star[(int) masses[i][1]][13]; + star_temp[j][14] = star[(int) masses[i][1]][14]; + j++; + } + } + if (msort) + Nhighmass = j; + + for (i=0;i0.5 expanding - double Rh = 0.8; //Half-mass radius [pc], ignored if profile = 3, set =-1 for using Marks & Kroupa (2012) Mcl-Rh relation - double gamma[3] = {2.0, 0.0, 2.0}; //Power-law slopes of EFF/Nuker templates (outer slope, inner slope, transition); set gamma[1] = 0.0 and gamma[2] = 2.0 for EFF (profile = 3) - double a = 1.0; //Scale radius of EFF/Nuker template (profile = 3) [pc] - double Rmax = 100.0; //Cut-off radius for EFF/Nuker template (profile = 3) [pc] - double tcrit = 100.0; //Simulation time [N-body units (Myr in Nbody6 custom)] - int tf = 3; //Tidal field: =0 no tidal field, =1 Near-field approximation, =2 point-mass galaxy, =3 Allen & Santillan (1991) MW potential (or Sverre's version of it) - double RG[3] = {8500.0,0.0,0.0}; //Initial Galactic coordinates of the cluster [pc] - double VG[3] = {0.0,220.0,0.0}; //Initial velocity of the cluster [km/s] - - //Mass function parameters - int mfunc = 1; //0 = single mass stars; 1 = use Kroupa (2001) mass function; 2 = use multi power law (based on mufu.c by L.Subr) - double single_mass = 1.0; //Stellar mass in case of single-mass cluster - double mlow = 0.08; //Lower mass limit for mfunc = 1 & mfunc = 4 - double mup = 100.0; //Upper mass limit for mfunc = 1 & mfunc = 4 - double alpha[MAX_AN] = {-1.35, -2.35, -2.7, 0.0, 0.0}; //alpha slopes for mfunc = 2 - double mlim[MAX_MN] = {0.08, 0.5, 4.0, 100.0, 0.0, 0.0}; //mass limits for mfunc = 2 - double alpha_L3 = 2.3; //alpha slope for mfunc = 4 (L3 mass function, Maschberger 2012) - double beta_L3 = 1.4; //beta slope for mfunc = 4 - double mu_L3 = 0.2; //mu parameter for mfunc = 4 - int weidner = 0; //Usage of Weidner & Kroupa (2006) relation for most massive star; =0 off, =1 on - int mloss = 3; //Stellar evolution; 0 = off, 3 = Eggleton, Tout & Hurley [KZ19] - int remnant = 1; //Use random kick velocity and present-day escape velocity to determine retention of compact remnants & evolved binary components (only for SSE/BSE version); =0 off, =1 on - double epoch = 0.0; //Age of the cluster, i.e. star burst has been ... Myr before [e.g. 1000.0, default = 0.0] [needs special compiling and SSE by Hurley, Pols & Tout] - double Z = 0.02; //Metallicity [0.0001-0.03, 0.02 = solar] - double FeH = -1.41; //Metallicity [Fe/H], only used when Z is set to 0 - int prantzos = 0; //Usage of Prantzos (2007) relation for the life-times of stars. Set upper mass limit to Lifetime(mup) >= epoch - - //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.1Msort and random pairing for remaining (Kiminki & Kobulnicky 2012; Sana et al., 2012; Kobulnicky et al. 2014; implemented by Long Wang) - 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 MMsort (implemented by Long Wang) - int OBperiods = 1; //Use period distribution for massive binaries with M_primary > msort from Sana & Evans (2011) if OBperiods = 1 - 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 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 - - //Gas parameters (only used for Nbody6 input) - double extmass = 0.0; //external Plummer (gas) sphere mass [Msun] - double extrad = 0.0; //external Plummer (gas) sphere scale factor [pc] - double extdecay = 0.0; //decay time for gas expulsion (set 0 for no decay) [Myr] - double extstart = 0.0; //delay time for start of gas expulsion [Myr] - - //Code parameters - int code = 3; //Nbody version: =0 Nbody6, =1 Nbody4, =2 Nbody6 custom, =3 only create output list of stars, =4 Nbody7 (not yet fully functional), =5 Nbody6++GPU - unsigned int seed = 0; //Number seed for random number generator; =0 for randomization by local time - char *output = "test"; //Name of output files - double dtadj = 1.0; //DTADJ [N-body units (Myr in Nbody6 custom)], energy-check time step - double dtout = 1.0; //DELTAT [N-body units (Myr in Nbody6 custom)], output interval, must be multiple of DTADJ - double dtplot = 1.0; //DTPLOT [N-body units (Myr in Nbody6 custom)], output of HRdiagnostics, should be multiple of DTOUT, set to zero if output not desired - int gpu = 0; //Use of GPU, 0= off, 1= on - int regupdate = 0; //Update of regularization parameters during computation; 0 = off, 0 > on - int etaupdate = 0; //Update of ETAI & ETAR during computation; 0 = off, 0 > on - int esc = 2; //Removal of escapers; 0 = no removal, 1 = regular removal at 2*R_tide; 2 = removal and output in ESC - int units = 1; //Units of McLuster output; 0= Nbody-Units, 1= astrophysical units - - //McLuster internal parameters - int match = 1; //Make cluster half-mass radius exactly match the expected/desired value; =0 off, =1 on (recommended) - int symmetry = 1; //Force spherical symmetry for fractal clusters; =0 off, =1 on (recommended) - int check = 0; //Make energy check at end of McLuster; =0 off, =1 on - int create_radial_profile = 1; //Creates a radial density profile and prints it to the screen; =0 off, =1 on - int create_cumulative_profile = 1; //Creates a radial cumulative profile and prints it to the screen; =0 off, =1 on - double Rgal = 10000.0; //Distance of cluster from sun for artificial CMD with observational errors [pc] - double Zsun = 0.02; //Solar metallicity - int NMAX = 1500000; //Maximum number of stars & orbits allowed in McLuster - int NNBMAX_NBODY6 = 500; //Maximum number of neighbours allowed in NBODY6 - double upper_IMF_limit = 150.0; //Maximum stellar mass allowed in McLuster [Msun] - int an = 0; //Counter for number of alpha slopes for mfunc = 2 - int mn = 0; //Counter for number of mass limits for mfunc = 1, 2 & 4 - int xn = 0; //Counter for components of galactocentric radius vector - int vn = 0; //Counter for components of cluster velocity vector - int xx = 0; //Counter for external potential input parameters - int gn = 0; //Counter for EFF/Nuker profile parameters - double extgas[4]; //Input array for external potential parameters - - - //SSE internal parameters (see Hurley, Pols & Tout 2000) - value1_.neta = 0.5; //Reimers mass-loss coefficent (neta*4x10^-13; 0.5 normally) - value1_.bwind = 0.0; //Binary enhanced mass loss parameter (inactive for single) - value1_.hewind = 1.0; //Helium star mass loss factor (1.0 normally) - value1_.mxns = 3.0; //Maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1) - points_.pts1 = 0.05; //Time-step parameter in evolution phase: MS (0.05) - points_.pts2 = 0.01; //Time-step parameter in evolution phase: GB, CHeB, AGB, HeGB (0.01) - points_.pts3 = 0.02; //Time-step parameter in evolution phase: HG, HeMS (0.02) - value4_.sigma = 190.0; //Kick velocities - value4_.bhflag = 1; //bhflag > 0 allows velocity kick at BH formation - - //BSE internal parameters (see Hurley, Pols & Tout 2002) - flags_.ceflag = 0; //ceflag > 0 activates spin-energy correction in common-envelope (0) #defunct# - flags_.tflag = 1; //tflag > 0 activates tidal circularisation (1) - flags_.ifflag = 0; //ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0) - flags_.nsflag = 1; //nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1) - flags_.wdflag = 1; //wdflag > 0 uses modified-Mestel cooling for WDs (0) - - value5_.beta = 1.0/8.0; //beta is wind velocity factor: proportional to vwind**2 (1/8) - value5_.xi = 1.0; //xi is the wind accretion efficiency factor (1.0) - value5_.acc2 = 3.0/2.0; //acc2 is the Bondi-Hoyle wind accretion factor (3/2) - value5_.epsnov = 0.001; //epsnov is the fraction of accreted matter retained in nova eruption (0.001) - value5_.eddfac = 1.0; //eddfac is Eddington limit factor for mass transfer (1.0) - value5_.gamma = -1.0; //gamma is the angular momentum factor for mass lost during Roche (-1.0) - - value2_.alpha1 = 1.0; //alpha1 is the common-envelope efficiency parameter (1.0) - value2_.lambda = 0.5; //lambda is the binding energy factor for common envelope evolution (0.5) - - - //Command line input - int option; - while ((option = getopt(argv, argc, "N:M:P:W:R:r:c:g:S:D:T:Q:C:A:O:G:o:f:a:m:B:b:p:s:t:e:Z:X:x:V:u:h:?")) != -1) switch (option) - { - case 'N': N = atoi(optarg); Mcl = 0.0; break; - case 'M': Mcl = atof(optarg); N = 0; break; - case 'P': profile = atoi(optarg); break; - case 'W': W0 = atof(optarg); break; - case 'R': Rh = atof(optarg); break; - case 'r': a = atof(optarg); break; - case 'c': Rmax = atof(optarg); break; - case 'g': - if (gn < 3) { gamma[gn] = atof(optarg); gn++; break; } - case 'S': S = atof(optarg); break; - case 'D': D = atof(optarg); break; - case 'T': tcrit = atof(optarg); break; - case 'Q': Q = atof(optarg); break; - case 'C': code = atoi(optarg); break; - case 'A': dtadj = atof(optarg); break; - case 'O': dtout = atof(optarg); break; - case 'G': gpu = atoi(optarg); break; - case 'o': output = optarg; break; - case 'f': mfunc = atoi(optarg); break; - case 'a' : - if (an < MAX_AN) { - alpha[an] = atof(optarg); - if (an == 0) alpha_L3 = atof(optarg); - if (an == 1) beta_L3 = atof(optarg); - if (an == 2) mu_L3 = atof(optarg); - an++; - break; - } else { printf("\nError: Number of alphas exceeded maximum limit of %d\n", MAX_AN); return 1; } - case 'm' : - if (mn < MAX_MN) { - mlim[mn] = atof(optarg); - if (mn == 0) mlow = atof(optarg); - if (mn == MAX_MN-1) mup = atof(optarg); - mn++; - break; - } else { printf("\nError: Number of mass params exceded maximum limit of %d\n", MAX_MN); return 1; } - case 'B': nbin = atoi(optarg); break; - case 'b': fbin = atof(optarg); break; - case 'p': pairing = atoi(optarg); break; - case 's': seed = atoi(optarg); break; - case 't': tf = atoi(optarg); break; - case 'e': epoch = atof(optarg); break; - case 'Z': Z = atof(optarg); break; - case 'X' : - if (xn < 3) { RG[xn] = atof(optarg); xn++; break; } - case 'V' : - if (vn < 3) { VG[vn] = atof(optarg); vn++; break; } - case 'x' : - if (xx < 4) { extgas[xx] = atof(optarg); xx++; break; } - case 'u': units = atoi(optarg); break; - case ':': - case 'h': help(msort); return 1; - case '?': help(msort); return 1; - }; + for (i=0;i 0) mup = mlim[mn-1]; + do { + l++; + if (star_temp[l][0]) { + j++; + } + } while (l 1) printf("\n\nUsing OpenMP with %d threads\n", omp_get_num_threads()); + for (j=0;j 0) { - if ((xx == 4) && (extgas[2])) { - extmass = extgas[0]; - extrad = extgas[1]; - extdecay = extgas[2]; - extstart = extgas[3]; - } else { - printf("\nWARNING: Insufficient or incorrect parameters specified for external Plummer potential!\n\n"); - exit (1); - } - } - - - /*********************** - * Generate star array * - ***********************/ - - int columns = 15; - double **star; - star = (double **)calloc(NMAX,sizeof(double *)); - for (j=0;j epoch\n"); - while (Lifetime(MMAX) < sqrt(pow(epoch,2))) { - MMAX -= 0.01; - } - } - - - /******************* - * Generate masses * - *******************/ - - printf("\n\n-----GENERATE MASSES----- \n"); - - //This loop has to be called multiple times with different N (or Mcl), Z and epoch in order to create multiple stellar populations - //make sure that epoch and Z are stored together with the stars, then concatenate all arrays - //for now, all Z and epochs are the same, i.e. a single stellar population - if (mfunc == 1) { - printf("\nMaximum stellar mass set to: %.2f\n", MMAX); - generate_m1(&N, star, mlow, mup, &M, &mmean, MMAX, Mcl, epoch, Z, Rh, remnant); - } else if (mfunc == 2) { - if (mn) { - for (i = mn+1; i < MAX_MN; i++) mlim[i] = 0.0; - } else { - for (i=0; i= mn) an = mn - 1; - mn = an + 1; - if (!mn){ - printf("\nError: at least one mass limit has to be specified\n"); - return 1; - } else if (mn == 1) { - single_mass = mlim[0]; - printf("\nSetting stellar masses to %g solar mass\n",single_mass); - if (!N) N = Mcl/single_mass; - for (j=0;j= 0; i--) { - norm[i] = norm[i+1] * pow(mlim[i+1], alpha[i+1] - alpha[i]); - subcount[i] = norm[i] * subint(mlim[i], mlim[i+1], alpha[i] + 1.); - N_tmp += subcount[i]; - submass[i] = norm[i] * subint(mlim[i], mlim[i+1], alpha[i] + 2.); - M_tmp += submass[i]; - } - generate_m2(an, mlim, alpha, Mcl, M_tmp, subcount, &N, &mmean, &M, star, MMAX, epoch, Z, Rh, remnant); - } - } else if (mfunc == 3) { - printf("\nMaximum stellar mass set to: %.2f\n",MMAX); - generate_m3(&N, star, mlow, mup, &M, &mmean, MMAX, Mcl); - randomize(star, N); - } else if (mfunc == 4) { - printf("\nMaximum stellar mass set to: %.2f\n", MMAX); - printf("\nUsing L3 IMF (Maschberger 2012)\n"); - generate_m4(&N, star, alpha_L3, beta_L3, mu_L3, mlow, mup, &M, &mmean, MMAX, Mcl, epoch, Z, Rh, remnant); - } else { - printf("\nSetting stellar masses to %.1f solar mass\n", single_mass); - if (!N) N = Mcl/single_mass; - for (j=0;j %.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); - order(star, N, M, msort, pairing); - } else { - randomize(star, N); - printf("\nApplying random pairing.\n"); - } - N -= nbin; - for (j=0;j mhighest) mhighest = star[i][0]; - } - mmeancom /= N; - int Nseg = ceil(N*mmeancom/mlowest); //number of necessary pos & vel pairs for Baumgardt et al. (2008) mass segregation routine - int Nunseg = N; - - double *Mcum; - Mcum = (double *)calloc(N,sizeof(double)); - - if ((S) && !(profile == 2)) {//sort masses when mass segregation parameter > 0 - printf("\nApplying mass segregation with S = %f\n",S); - segregate(star, N, S); - for (i=0;i0.5) - p[4] = gamma[1];//inner power-law slope (0.0 for EFF template) - p[5] = gamma[2];//transition parameter (2.0 for EFF template) - generate_profile(N, star, Rmax, M, p, &Rh, D, symmetry); - printf("\nRh = %.1f pc\n", Rh); - } else if (profile == -1) { - printf("\nGenerating fractal distribution with parameters: N = %i\t Rh = %.3f\t D = %.2f\n", N, Rh, D); - fractalize(D, N, star, 0, symmetry); - rvir = Rh; - } else { - printf("\nGenerating Plummer model with parameters: N = %i\t Rh = %.3f\t D = %.2f\n", N, Rh, D); - rvir = Rh/0.772764; - rplummer = Rh/1.305; - generate_plummer(N, star, rtide, rvir, D, symmetry); - } - - - //Apply Baumgardt et al. (2008) mass segregation - if (!(profile == 2) && (S)) { - double **m_temp; - m_temp = (double **)calloc(Nunseg,sizeof(double *)); - for (j=0;jtphys no data will be stored - double z = Z; //metallicity - double zpars[20]; //metallicity parameters - double vkick; //kick velocity for compact remnants - double vesc; //escape velocity of cluster - int lostremnants = 0; //number of ejected compact remnants - double lostremnantsmass = 0.0; //mass of ejected compact remnants - if (Mcl && remnant && (Rh != -1)) { - vesc = sqrt(2.0*G*Mcl/Rh); - printf("Escape velocity of cluster = %.4f km/s\n", vesc); - } else if ((!remnant) || (Rh == -1)) { - vesc = 1.0E10; - printf("Keeping all compact remnants\n"); - } else { - vesc = sqrt(2.0*0.4**N/Rh); - printf("Estimated escape velocity of cluster assuming mean stellar mass of 0.4 Msun = %.4f km/s\n", vesc); - } - for (i=0; i<20; i++) zpars[i] = 0; - zcnsts_(&z,zpars); //get metallicity parameters - - printf("\nSetting up stellar population with Z = %.4f.\n",Z); - - if (epoch) printf("\nEvolving stellar population for %.1f Myr.\n",epoch); - - //set up mass function parameters - ty = 2; - alpha1 = 1.3; - alpha2 = 2.3; - - c1 = 1.0-alpha1; - c2 = 1.0-alpha2; - - k1 = 2.0/c1*(pow(0.5,c1)-pow(mlow,c1)); - if (mlow>0.5) { - k1 = 0; - k2 = 1.0/c2*(pow(mup,c2)-pow(mlow,c2)); - } else - k2 = k1 + 1.0/c2*(pow(mup,c2)-pow(0.5,c2)); - if (mup<0.5) { - k1 = 2.0/c1*(pow(mup,c1)-pow(mlow,c1)); - k2 = k1; - } - - - //determine theoretical mean mass from mass function - c1 = 2.0-alpha1; - c2 = 2.0-alpha2; - - if (mlow != mup) { - if (mlow>0.5) { - mth = (1.0/c2*(pow(mup,c2)-pow(mlow,c2)))/k2; - } else if (mup<0.5) { - mth = (2.0/c1*(pow(mup,c1)-pow(mlow,c1)))/k2; - } else - mth = (2.0/c1*(pow(0.5,c1)-pow(mlow,c1))+1.0/c2*(pow(mup,c2)-pow(0.5,c2)))/k2; - } else { - mth = mlow; - } - - if (!*N) { - *N = max(floor((Mcl-MMAX)/mth), 1); - if (!epoch) printf("Estimated number of necessary stars: %i\n", *N); - *N = 1; - } - - - c1 = 1.0-alpha1; - c2 = 1.0-alpha2; - *mmean = 0.0; - *M = 0.0; - double mostmassive = 0.0; - - for (i=0; i<*N; i++) { - do{ - do { - xx = drand48(); - if (xx MMAX); - - //evolve star for deltat = epoch with SSE (Hurley, Pols & Tout 2002) - tphys = 0.0; - kw = 1; - mass = star[i][0]; //initial mass - mt = mass; //actual mass - ospin = 0.0; - tphysf = epoch; - epoch1 = 0.0; - vkick = 0.0; - //printf("MASS %.2f", mass); - star[i][7] = mass; - evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, &epoch1, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick); - //printf("-> %.2f\n", mass); - //if (vkick) printf("KICK: %.5f\n", vkick); - lostremnants++; - lostremnantsmass += mt; - } while (vkick > vesc); - lostremnants--; - lostremnantsmass -= mt; - star[i][0] = mt; - star[i][8] = kw; - star[i][9] = epoch1; - star[i][10] = ospin; - star[i][11] = r; - star[i][12] = lum; - if (star[i][0] > mostmassive) mostmassive = star[i][0]; - *M += star[i][0]; - if ((i==*N-1) && (*Mtphys no data will be stored - double z = Z; //metallicity - double zpars[20]; //metallicity parameters - double vkick; //kick velocity for compact remnants - double vesc; //escape velocity of cluster - int lostremnants = 0; //number of ejected compact remnants - double lostremnantsmass = 0.0; //mass of ejected compact remnants - if (Mcl && remnant && (Rh != -1)) { - vesc = sqrt(2.0*G*Mcl/Rh); - printf("Escape velocity of cluster = %.4f km/s\n", vesc); - } else if ((!remnant) || (Rh == -1)) { - vesc = 1.0E10; - printf("Keeping all compact remnants\n"); - } else { - vesc = sqrt(2.0*0.4**N/Rh); - printf("Estimated escape velocity of cluster assuming mean stellar mass of 0.4 Msun = %.4f km/s\n", vesc); - } - - for (i=0; i<20; i++) zpars[i] = 0; - zcnsts_(&z,zpars); //get metallicity parameters - - printf("\nSetting up stellar population with Z = %.4f.\n",Z); - - if (epoch) printf("\nEvolving stellar population for %.1f Myr.\n",epoch); - - double tmp, ml, mup; - double mostmassive = 0.0; - *mmean = 0.0; - *M = 0.0; - if (!*N) *N = 1; - - for (i = 0; i < an; i++) { - printf("# <%.2f , %.2f> .. %.2f\n", mlim[i], mlim[i+1], alpha[i]); - } - - for (i = 1; i < an; i++) subcount[i] += subcount[i-1]; - - for (i = 0; i < *N; i++) { - do { - do { - tmp = drand48() * subcount[an-1]; - for (j = 0; (j < an) && (subcount[j] < tmp); j++); - if (alpha[j] != -1.) { - ml = pow(mlim[j], 1. + alpha[j]); - mup = pow(mlim[j+1], 1. + alpha[j]); - } else { - ml = log(mlim[j]); - mup = log(mlim[j+1]); - } - tmp = ml + drand48() * (mup - ml); - if (alpha[j] != -1.) star[i][0] = pow(tmp, 1. / (1. + alpha[j])); - else star[i][0] = exp(tmp); - } while (star[i][0] > MMAX); - //printf("%8.4f\n", star[i][0]); - - //evolve star for deltat = epoch with SSE (Hurley, Pols & Tout 2002) - tphys = 0.0; - kw = 1; - mass = star[i][0]; //initial mass - mt = mass; //actual mass - ospin = 0.0; - vkick = 0.0; - tphysf = epoch; - epoch1 = 0.0; - //printf("MASS %.2f", mass); - star[i][7] = mass; - if (epoch) evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, &epoch1, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick); - //printf("-> %.2f\n", mass); - //printf("KICK: %.5f\n", vkick); - lostremnants++; - lostremnantsmass += mt; - } while (vkick > vesc); - lostremnants--; - lostremnantsmass -= mt; - star[i][0] = mt; - star[i][8] = kw; - star[i][9] = epoch1; - star[i][10] = ospin; - star[i][11] = r; - star[i][12] = lum; - - if (star[i][0] > mostmassive) mostmassive = star[i][0]; - *M += star[i][0]; - if ((i==*N-1) && (*M mb) { - m = pow(pow(star[i-1][0], 2.0-a2) - star[i-1][0]/k2*(2.0-a2), 1.0/(2.0-a2)); - if (m < mb) { - m = pow(pow(mb, 2.0-a1) - (2.0-a1)/pow(mb, a1-a2)*(star[i-1][0]/k2 + pow(mb, 2.0-a2)/(2.0-a2) - 1.0/(2.0-a2)*(pow(star[i-1][0], 2.0-a2))), 1.0/(2.0-a1)); - } - } else { - m = pow((a1-2.0)/k1*star[i-1][0] + pow(star[i-1][0], 2.0-a1), 1.0/(2.0-a1)); - } - if (mtphys no data will be stored - double z = Z; //metallicity - double zpars[20]; //metallicity parameters - double vkick; //kick velocity for compact remnants - double vesc; //escape velocity of cluster - int lostremnants = 0; //number of ejected compact remnants - double lostremnantsmass = 0.0; //mass of ejected compact remnants - double mostmassive = 0.0; - - - double G_l; - double G_u; - double u; - double x; - int i; - double mth; - - double gamma; - double mpeak; - // variables for the mean - double B_l; - double B_u; - double a; - double b; - double beta_log; - int ifault; - - printf("\nL3 Parameters:\n"); - printf("alpha = %1.2f beta = %1.2f mu = %1.2f m_low = %3.3f m_up = %3.3f\n",alpha,beta,mu,mlow,mup); - - mpeak = mu*pow((beta-1.0),1.0/(alpha-1.0)); - gamma = alpha + beta*(1.0-alpha); - - printf("'Peak' mass = %3.3f Effective low mass exponent gamma = %1.2f\n\n",mpeak,gamma); - - // normalisation - G_l = pow(1.0 + pow(mlow/mu, 1.0-alpha) , 1.0-beta); - G_u = pow(1.0 + pow(mup/mu, 1.0-alpha) , 1.0-beta); - - - // Calculate the mean mass - a = (2.0-alpha)/(1.0-alpha); - b = beta - a; - - // logarithm of the beta function necessary to calculate the incomplete beta fn - beta_log = alogam( a, &ifault ) + alogam( b, &ifault ) - alogam( a + b, &ifault ); - - x = pow(mlow/mu,1.0-alpha); - x = x/(1.0+x); - B_l = betain(x, a, b, beta_log, &ifault)*exp(beta_log); - - x = pow(mup/mu,1.0-alpha); - x = x/(1.0+x); - B_u = betain(x, a, b, beta_log, &ifault)*exp(beta_log); - mth = mu*(1.0-beta)*(B_u - B_l)/(G_u - G_l) ; - - if (Mcl && remnant && (Rh != -1)) { - vesc = sqrt(2.0*G*Mcl/Rh); - printf("Escape velocity of cluster = %.4f km/s\n", vesc); - } else if ((!remnant) || (Rh == -1)) { - vesc = 1.0E10; - printf("Keeping all compact remnants\n"); - } else { - vesc = sqrt(2.0*0.4**N/Rh); - printf("Estimated escape velocity of cluster assuming mean stellar mass of 0.4 Msun = %.4f km/s\n", vesc); - } - for (i=0; i<20; i++) zpars[i] = 0; - zcnsts_(&z,zpars); //get metallicity parameters - - printf("\nSetting up stellar population with Z = %.4f.\n",Z); - - if (epoch) printf("\nEvolving stellar population for %.1f Myr.\n",epoch); - - if (!*N) { - *N = max(floor((Mcl-mup)/mth), 1); - if (!epoch) printf("Estimated number of necessary stars: %i\n", *N); - *N = 1; - } - - // generate random masses - for (i=0; i<*N; i++) { - do { - u = drand48(); - x = u*(G_u-G_l)+G_l; - x = pow(x,1/(1-beta)) - 1.0; - star[i][0] = mu*pow(x,1.0/(1.0-alpha)); - - tphys = 0.0; - kw = 1; - mass = star[i][0]; //initial mass - mt = mass; //actual mass - ospin = 0.0; - tphysf = epoch; - epoch1 = 0.0; - vkick = 0.0; - //printf("MASS %.2f", mass); - star[i][7] = mass; - evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, &epoch1, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick); - //printf("-> %.2f\n", mass); - //if (vkick) printf("KICK: %.5f\n", vkick); - lostremnants++; - lostremnantsmass += mt; - } while (vkick > vesc); - lostremnants--; - lostremnantsmass -= mt; - star[i][0] = mt; - star[i][8] = kw; - star[i][9] = epoch1; - star[i][10] = ospin; - star[i][11] = r; - star[i][12] = lum; - if (star[i][0] > mostmassive) mostmassive = star[i][0]; - - *M += star[i][0]; - - if ((i==*N-1) && (*M=3.0) { - - for (i=0;ircut); //reject particles beyond tidal radius - - //velocities - do { - a[4] = drand48(); - a[5] = drand48(); - a[6] = pow(a[4],2)*pow(1.0 - pow(a[4],2),3.5); - } while (0.1*a[5]>a[6]); - - a[8] = a[4]*sqrt(2.0)/pow(1.0 + ri*ri,0.25); - a[6] = drand48(); - a[7] = drand48(); - - star[i][6] = (1.0 - 2.0*a[6])*a[8]; - star[i][4] = sqrt(a[8]*a[8] - pow(star[i][6],2))*cos(TWOPI*a[7]); - star[i][5] = sqrt(a[8]*a[8] - pow(star[i][6],2))*sin(TWOPI*a[7]); - } - } else { - double xcut; - xcut = 0.000; - while (1.0/sqrt(pow(xcut,-2.0/3.0) - 1.0)<=rcut) xcut+=0.00001; - - fractalize(D, N, star, 1, symmetry); - for (i=0;i12.0) { - printf("W0 too large\n"); - return 0; - } else if (W0 < 0.2) { - printf("W0 too small\n"); - return 0; - } - - - - - /*************************** - * INTERPOLATE KING VALUES * - ***************************/ - - h = pow(W0-rmin,2)/(1.0*M-1.0); //step size for interpolation - den = densty(W0); //central density - - //x is King's W, y[0] is z**2, y[1] is 2*z*dz/dW, where z is scaled radius - //so x(y[0]) is energy as function of radius^2 and x(y[1]) is the derivative - - xstart = 0.000001; - ystart0 = 2.0*xstart/3.0; - ystart1 = -2.0/3.0; - - x1 = W0 - xstart; - x2 = 0.0; - - //integrate Poisson's eqn - printf("Integrating Poisson's equation\n"); - - odeint(ystart0, ystart1, x1, x2, den, &kount, xp, yp, M, KMAX); - - printf("No of integration steps = %i\n",kount); - - - - //interpolate yking - for (k=0;k xp[0]) { - yking[k][0] = (W0-x[k])*yp[0][0]/(W0 - xp[0]); - yking[k][1] = yp[0][1]; - printf("+"); - } else { - - i = 0; - - do { - if ((x[k]-xp[i])*(x[k]-xp[i+1]) <= 0.0) { - yking[k][0] = yp[i][0] + (yp[i+1][0]-yp[i][0])*(x[k]-xp[i])/(xp[i+1]-xp[i]); - yking[k][1] = yp[i][1] + (yp[i+1][1]-yp[i][1])*(x[k]-xp[i])/(xp[i+1]-xp[i]); - goto jump; - } else { - i++; - } - } while (i= kount) { - yking[k][0] = yp[kount][0]; - yking[k][1] = yp[kount][1]; - } - - - //get mass as a function of radius - mass[k] = -2.0*pow(yking[k][0],1.5)/yking[k][1]; - - - //calculate total potential energy - if (k == 0) { - pot = -0.5*pow(mass[k],2)/sqrt(yking[k][0]); - } else { - pot -= 0.5*(pow(mass[k],2) - pow(mass[k-1],2))/(0.5*(sqrt(yking[k][0]) + sqrt(yking[k-1][0]))); - } - } - - - - - /****************************** - * DETERMINE BASIC QUANTITIES * - ******************************/ - - //Total mass - totmas = -2.0*pow(yking[M-2][0],1.5)/yking[M-2][1]; - - //Half-mass radius - k=0; - - do { - k++; - } while (mass[k] < 0.5*totmas); - - zh = yking[k][0] - (yking[k][0]-yking[k-1][0])*(mass[k]-0.5*totmas)/(mass[k]-mass[k-1]); - *rh = sqrt(zh); - - //Virial radius and King radius - *rvir = -pow(totmas,2)/(2.0*pot); - *rking = sqrt(yking[M-2][0]); - - //Central velocity dispersion**2 (3-dimensional) - ve = sqrt(2.0*W0); - cg1 = (-pow(ve,3)*exp(-pow(ve,2)/2.0) - 3.0*ve*exp(-pow(ve,2)/2.0) + 3.0/2.0*sqrt(2.0*PI)*erf(ve*sqrt(2.0)/2.0) - pow(ve,5)*exp(-pow(ve,2)/2.0)/5.0)/(-ve*exp(-pow(ve,2)/2.0) + sqrt(2.0*PI)*erf(ve*sqrt(2.0)/2.0)/2.0 - pow(ve,3)*exp(-pow(ve,2)/2.0)/3.0); - - printf("\nTheoretical values:\n"); - printf("Total mass (King units) = %g\n", totmas); - printf("Viriral radius (King units) = %g\n", *rvir); - printf("Half-mass radius (King units) = %g\t(Nbody units) = %g\n", *rh, *rh/ *rvir); - printf("Edge radius (King units) = %g\t(Nbody units) = %g\n", *rking, *rking/ *rvir); - printf("Concentration = %g\n", log10(*rking)); - printf("Core radius (King units) = %g\t(Nbody units) = %g\n", 1.0, 1.0/ *rvir); - printf("3d velocity dispersion**2: %g (central)\t %g (global)\n", cg1, -pot/totmas); - - - - /*************************** - * GENERATE STAR POSITIONS * - ***************************/ - - printf("\nGenerating Stars:\n"); - - if (D>=3.0) { - for (i=0;iM) printf("WARNING: failing iteration\n"); - } while (mass[j] <= fmass); - - r2 = yking[j-1][0] + (fmass-mass[j-1])*(yking[j][0] - yking[j-1][0])/(mass[j]-mass[j-1]); - } - - - r = sqrt(r2); - costh = 2.0*drand48()-1.0; - sinth = sqrt(1.0-pow(costh,2)); - phi = 2.0*PI*drand48(); - xstar = r*sinth*cos(phi); - ystar = r*sinth*sin(phi); - zstar = r*costh; - - - /**************** - * CHOOSE SPEED * - ****************/ - - if (r < *rking) { - if (r < sqrt(yking[0][0])) { - w1 = x[0]; - w = W0 - (r2/yking[0][0])*(W0 - w1); - } else { - j = 0; - do { - j++; - } while (r > sqrt(yking[j][0])); - wj1 = x[j-1]; - wj = x[j]; - w = wj1 + (r2-yking[j-1][0])*(wj-wj1)/(yking[j][0]-yking[j-1][0]); - } - } else { - printf("radius too big\n"); - } - - vmax = sqrt(2.0*w); - do { - speed = vmax*drand48(); - fstar = pow(speed,2)*(exp(-0.5*pow(speed,2))-exp(-1.0*w)); - } while (fstar < 2.0*drand48()/exp(1.0)); - - costh = 2.0*drand48()-1.0; - phi = 2.0*PI*drand48(); - sinth = sqrt(1.0-pow(costh,2)); - ustar = speed*sinth*cos(phi); - vstar = speed*sinth*sin(phi); - wstar = speed*costh; - mstar = star[i][0]; - - //printf("i: %i\tr=%g\tm=%.5f\tx=%.5f y=%.5f z=%.5f\tvx=%.5f vy=%.5f vz=%.5f\n",i,sqrt(r2),mstar,xstar,ystar,zstar,ustar,vstar,wstar); - - - coord[i][0] = xstar; - coord[i][1] = ystar; - coord[i][2] = zstar; - coord[i][3] = ustar; - coord[i][4] = vstar; - coord[i][5] = wstar; - - } - - for (i=0;iM) printf("WARNING: failing iteration\n"); - } while (mass[j] <= fmass); - - r2 = yking[j-1][0] + (fmass-mass[j-1])*(yking[j][0] - yking[j-1][0])/(mass[j]-mass[j-1]); - } - - r = sqrt(r2); - - - /**************** - * CHOOSE SPEED * - ****************/ - - if (r < *rking) { - if (r < sqrt(yking[0][0])) { - w1 = x[0]; - w = W0 - (r2/yking[0][0])*(W0 - w1); - } else { - j = 0; - do { - j++; - } while (r > sqrt(yking[j][0])); - wj1 = x[j-1]; - wj = x[j]; - w = wj1 + (r2-yking[j-1][0])*(wj-wj1)/(yking[j][0]-yking[j-1][0]); - } - } else { - printf("radius too big\n"); - } - - - vmax = sqrt(2.0*w); - do { - speed = vmax*drand48(); - fstar = pow(speed,2)*(exp(-0.5*pow(speed,2))-exp(-1.0*w)); - } while (fstar < 2.0*drand48()/exp(1.0)); - - for (k=1;k<4;k++) star[i][k] *= r/ri; - for (k=4;k<7;k++) star[i][k] *= speed; - } - - } - - - for (j=0;j= 0.0) { - h = sqrt(pow(H1,2)); - } else { - h = - sqrt(pow(H1,2)); - } //step size - - y[0] = ystart0; //z^2 - y[1] = ystart1; //2*z*dz/dW where z is scaled radius - - xsav = x-DXSAV*2.0; - - for (i=0;i sqrt(pow(DXSAV,2))) { - if (*kount < KMAX-1) { - xp[*kount] = x; - for (j=0;j<2;j++) { - yp[*kount][j] = y[j]; - } - *kount = *kount + 1; - xsav = x; - } - } //store x, y1 and y2 if the difference in x is smaller as the desired output step size DXSAV - - if (((x+h-x2)*(x+h-x1)) > 0.0) h = x2-x; - - rkqc(y,dydx,&x,&h,den,yscal, &hdid, &hnext, TOL); //do a Runge-Kutta step - - if ((x-x2)*(x2-x1) >= 0.0) { - ystart0 = y[0]; - ystart1 = y[1]; - - xp[*kount] = x; - for (j=0;j<2;j++) { - yp[*kount][j] = y[j]; - } - return 0; - *kount = *kount +1; - } - - if (sqrt(pow(hnext,2)) < HMIN) { - printf("Stepsize smaller than minimum.\n"); - return 0; - } - - h = hnext; - } - - return 0; -} + //copying back to original array + for (i=0;i= 0.0) { - rhox =-sqrt(x)*(x+1.5)+0.75*sqrt(PI)*exp(x)*erf(sqrt(x)); - } else { - rhox = 0.0; - } - - dydx[0]= y[1]; - dydx[1] = 0.25*pow(y[1],2)*(6.0+9.0*y[1]*rhox/den)/y[0]; - - return 0; -} + for (j=0;j 1.0) *h = safety**h*(pow(errmax,pshrnk)); //if integration error is too large, decrease h - - } while (errmax > 1.0); - - - *hdid = *h; - if (errmax > errcon) { - *hnext = safety**h*(pow(errmax,pgrow));//increase step size for next step - } else { - *hnext = 4.0**h;//if integration error is very small increase step size significantly for next step - } - - for (i=0;i<2;i++) { - y[i] += ytemp[i]*fcor; - } - - return 0; - + return 0; } -int rk4(double x, double *y, double *dydx, double h, double *yout, double den){ - double hh, h6, xh; - double yt[2], dyt[2], dym[2]; - int i; - - hh=h*0.5; - h6=h/6.0; - xh=x+hh; - for (i=0;i<2;i++) { - yt[i] = y[i] + hh*dydx[i]; - } - - derivs(xh,yt,dyt,den); //find derivative - for (i=0;i<2;i++) { - yt[i] = y[i] + hh*dyt[i]; - } - - derivs(xh,yt,dym,den); //find derivative - for (i=0;i<2;i++) { - yt[i] = y[i] + h*dym[i]; - dym[i] += dyt[i]; - } - - derivs(x+h,yt,dyt,den); //find derivative - for (i=0;i<2;i++) { - yout[i] = y[i] + h6*(dydx[i]+dyt[i]+2.0*dym[i]); - } - - return 0; +double rtnewt (double ecc, double ma) { + + double x1,x2,xacc,rtnewt,f,df,dx; + int j,jmax; + + x1 = 0; + x2 = 2*PI; + xacc = 1E-6; + jmax = 20; + ma = 2*PI*ma; + + rtnewt=.5*(x1+x2); + for (j=1;j<=jmax;j++) { + f = ma - rtnewt + ecc*sin(rtnewt); + df = -1 + ecc*cos(rtnewt); + dx=f/df; + rtnewt=rtnewt-dx; + if ((x1-rtnewt)*(rtnewt-x2)<0) + printf("jumped out of brackets\n"); + if(abs(dx)0.5 the Subr mass segregation may generate a virially hot system\nSetting cut-off radius of Subr model to twice the approximate tidal radius\n"); - rcut = 2.0*rtide/rvir; //cut-off radius for Plummer sphere = tidal radius - } - M_tot = 0.; - N_star = N; - star1 = (struct t_star1 *)realloc(star1, N_star * sizeof(struct t_star1)); - star2 = (struct t_star2 *)realloc(star2, N_star * sizeof(struct t_star2)); - - for (i=0; i 0.) star1[i].Ui = star1[i].mass * pow(M_sub, -S); - else star1[i].Ui = star1[i].mass; - for (j = 0; j < i; j++) { - tmp = star1[i].Ui * star1[j].Ui; - star2[i].UT_add += tmp; - star1[i].U_tmp += tmp; - star1[j].U_tmp += tmp; - }; - U_tot1 += star2[i].UT_add; - }; - - UT_sub = U_tot = K_tot = 0.; - - for (i = 0; i < N_star; i++) { - if ((i/1000)*1000 == i) printf("Generating orbit #%li\n", i); - star2[i].UT_add *= 0.5 / U_tot1; - star2[i].UT = star1[i].U_tmp * 0.5 / U_tot1; - UT_sub += star2[i].UT_add; - do { - position(i, U_tot, UT_sub, S); - rtemp = sqrt(pow(star1[i].r[0],2)+pow(star1[i].r[1],2)+pow(star1[i].r[2],2)); - } while (rtemp>rcut); - U_tot += star1[i].mass * star2[i].U_sub; - }; - - for (i = 0; i < N_star; i++) { - maxim = 0.25 * star1[i].mass / star2[i].UT; - find_beta(&maxim, &beta); - if (beta > 0.) - { - do - { - r1 = maxim * drand48(); - r2 = drand48(); - r2 *= r2; - tmp = 1. - r2; - } - while (r1 > r2 * exp(beta * log(tmp))); - } - else - { - do - { - r1 = maxim * drand48(); - r2 = sin(M_PI_2 * drand48()); - r2 *= r2; - tmp = sqrt(1. - r2); - } - while (r1 > r2 * exp((2.* beta + 1.) * log(tmp))); - }; - - star2[i].E = r2 * (star2[i].U + star2[i].U_sub); - - tmp = sqrt(2. * star2[i].E); - isorand(star2[i].v); - star2[i].v[0] *= tmp; - star2[i].v[1] *= tmp; - star2[i].v[2] *= tmp; - K_tot += star1[i].mass * star2[i].E; - }; - - for (i = 0; i < N_star; i++){ - //printf("%.12f %18.14f %18.14f %18.14f %16.12f %16.12f %16.12f\n",star1[i].mass, star1[i].r[0], star1[i].r[1], star1[i].r[2], star2[i].v[0], star2[i].v[1], star2[i].v[2]); - star[i][0] = star1[i].mass; - star[i][1] = star1[i].r[0]; - star[i][2] = star1[i].r[1]; - star[i][3] = star1[i].r[2]; - star[i][4] = star2[i].v[0]; - star[i][5] = star2[i].v[1]; - star[i][6] = star2[i].v[2]; - star[i][7] = star1[i].m0; - star[i][8] = star1[i].kw; - star[i][9] = star1[i].epoch; - star[i][10] = star1[i].spin; - star[i][11] = star1[i].rstar; - star[i][12] = star1[i].lum; - star[i][13] = star1[i].epochstar; - star[i][14] = star1[i].zstar; - - }; - - return 0; -} +int eigenevolution(double *m1, double *m2, double *ecc, double *abin){ + double alpha = 28.0; + double beta = 0.75; + double mtot,lper,lperi,ecci,mtoti,r0,rperi,qold,qnew; -void quick(int start, int stop) { - int i, j; - double temp, median; - - median = 0.5 * (star1[start].mass + star1[stop].mass); - - i = start; - j = stop; - while (i <= j) - { - while ((i <= stop) && (star1[i].mass > median)) i++; - while ((j >= start) && (star1[j].mass < median)) j--; - if (j > i) - { - SWAP(i,j); - i++; - j--; - } - else if (i == j) - { - i++; - j--; - }; - }; - - if (start + 1 < j) quick(start, j); - else if ((start < j) && (star1[start].mass < star1[j].mass)) SWAP(start, j); - - if (i + 1 < stop) quick(i, stop); - else if ((i < stop) && (star1[i].mass < star1[stop].mass)) SWAP(i, stop) -} + *abin *= PARSEC; + mtot = *m1+*m2; -void position(int id, double U_sub, double UT_sub, double S){ - long i; - double rx, ry, rz, r1; - double rfac; - - rfac = pow(star2[id].M_sub, 2. * S) / (1. - S); - - star2[id].ntry = 0; - - do - { - star2[id].U_sub = 0.; - star2[id].ntry++; - - r1 = drand48(); - r1 = exp(-2. * log(r1) / 3.); - r1 = RSCALE * sqrt(1. / (r1 - 1.)); - - r1 *= rfac; - - isorand(star1[id].r); - star1[id].r[0] *= r1; - star1[id].r[1] *= r1; - star1[id].r[2] *= r1; - star2[id].rad = r1; - - for (i = 0; i < id; i++) - { - rx = star1[i].r[0] - star1[id].r[0]; - ry = star1[i].r[1] - star1[id].r[1]; - rz = star1[i].r[2] - star1[id].r[2]; - r1 = sqrt(rx*rx + ry*ry + rz*rz); - star1[i].U_tmp = star1[id].mass / r1; - star2[id].U_sub += star1[i].mass / r1; - }; - r1 = U_sub + star1[id].mass * star2[id].U_sub; - } - while (fabs(r1 - UT_sub) > 0.1 * (r1 + UT_sub) / sqrt(1. + id)); - - for (i = 0; i < id; i++) star2[i].U += star1[i].U_tmp; -} + lperi = sqrt(pow(*abin,3)/mtot*4.0*PI*PI/GBIN); -void find_beta(double *avg, double *beta){ - int i; - double a1, a2, I1, I2, J1, J2; - double alo, ah, Ilo, Ih, Jlo, Jh; - double tmpavg; - - if (*avg > 0.75) - { - *beta = -0.5; - *avg = 1.; - return; - }; - - Ilo = I1 = 3. * M_PI / 16.; - Ih = I2 = 0.2; - Jlo = J1 = M_PI / 4; - Jh = J2 = 1. / 3.; - alo = a1 = -0.5; - ah = a2 = 0.; - - tmpavg = I2 / J2; - i = 0; - while (tmpavg > *avg) - { - if (!(i & 1)) - { - a1 += 1.; - I1 /= 1. + 5. / (2. * a1); - J1 /= 1. + 3. / (2. * a1); - Ilo = I2; Ih = I1; - Jlo = J2; Jh = J1; - alo = a2; ah = a1; - } - else - { - a2 += 1.; - I2 /= 1. + 5. / (2. * a2); - J2 /= 1. + 3. / (2. * a2); - Ilo = I1; Ih = I2; - Jlo = J1; Jh = J2; - alo = a1; ah = a2; - }; - tmpavg = Ih / Jh; - i++; - }; - - *beta = alo + (ah - alo) * (*avg - Ilo / Jlo) / (tmpavg - Ilo / Jlo); - - if (*beta > 0.) *avg = exp(*beta * log(*beta / (1. + *beta))) / (1. + *beta); - else *avg = 2. * exp((2.* *beta + 1.) * log((2.* *beta + 1.) / (2.* *beta + 3.))) / (2.* *beta + 3.); -} + ecci = *ecc; + mtoti = mtot; -void isorand(double *r) { - double rad; - - do - { - r[0] = 2. * drand48() - 1.; - r[1] = 2. * drand48() - 1.; - r[2] = 2. * drand48() - 1.; - rad = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]; - } - while ((rad > 1.) || (rad < 0.0001)); - - rad = sqrt(rad); - r[0] /= rad; - r[1] /= rad; - r[2] /= rad; -} + /* Circularisation */ + r0 = alpha*RSUN; + rperi = *abin*(1.0-*ecc); + alpha = -pow((r0/rperi),beta); -int cmpmy(double *x1, double *x2) { //smallest to largest - if(*x1<*x2) return -1; - return 1; -} + if (*ecc > 0) { + *ecc = exp(alpha+log(*ecc)); + } -int cmpmy_reverse(double *x1, double *x2) { //largest to smallest - if(*x1>*x2) return -1; - return 1; -} + /* pre-ms mass-transfer */ + qold = *m1/ *m2; + if (qold > 1.0) qold = 1.0/qold; + alpha = -alpha; + if (alpha > 1.0) { + qnew = 1.0; + } else { + qnew = qold + (1.0-qold)*alpha; + } -double generate_profile (int N, double **star, double Rmax, double Mtot, double *p, double *Rh, double D, int symmetry) { - - double Mnorm; - double r, Rmin = 1E-05; - int steps = 51, i, h, j; //51 steps is more than sufficient (in most cases) - double r_array[steps], rho_array[steps], sigma_array[steps], M_array[steps], X_array[steps]; - double random[7], x,y,z, ri, vx, vy, vz; - double r_norm, v_norm; - - Mnorm = M(Rmax, p); - p[1] = Mtot/Mnorm; - - double stepsize; - stepsize = (log10(Rmax)-log10(Rmin))/(steps-1); - - printf("\nIntegrating profile...\n"); - for (i=0;i=3.0) { - - for (i = 0; i< N; i++){ - - if ((i/1000)*1000 == i) printf("Generating orbit #%i\n", i); - - do { - random[0] = drand48(); - - j = 0; - while (X_array[j]Rmax); //reject particles beyond Rmax - - double sigma; - - j = 0; - while (r_array[j] arg2) errt = arg1; - else errt = arg2; - if (errt <= *err) { - *err=errt; - ans=a[j][i]; - } - } - if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*(*err)) break; - } - - return ans; -} + mtot = *m1 + *m2; -double midexp(double (*func)(double, double*), double aa, double bb, int n, double *p) { - double x,tnm,sum,del,ddel,a,b; - static double s; - int it,j; - - b=exp(-aa); - a=0.0; - if (n == 1) { - x = 0.5*(a+b); - return (s=(b-a)*func(-log(x),p)/x); - } else { - for(it=1,j=1;j=0.5*M) && !(*Rh3D)) *Rh3D = rarray[i][0]; + i++; + } else { + j++; + } + } -double rhoR (double x, double *p) { //user-defined 2d-density profile - // return p[1]*pow(1.0+x/p[2]*x/p[2],-0.5*p[3]); //Elson, Fall & Freeman (1987) - return p[1]*pow(2.0,0.5*(p[3]-p[4]))*pow(x/p[2],-p[4])*pow(1.0+pow(x/p[2],p[5]),-(p[3]-p[4])/p[5]); //Nuker (Lauer et al. 1995) -> Elson, Fall & Freeman (1987) for p[4] = 0 & p[5] = 2 - -} + j = 0; i = 0; Mtemp = 0.0; + while ((j < noofradialbins) && (i < N)) { + if (rarray2D[i][0] < rprofile[j][0]) {//2D binning in astrophysical units + rprofile[j][6] += 1.0; + rprofile[j][7] += rarray2D[i][1]; + rprofile[j][8] += rarray2D[i][2]; + Mtemp += rarray2D[i][1]; + if ((Mtemp>=0.5*M) && !(*Rh2D)) *Rh2D = rarray2D[i][0]; + i++; + } else { + j++; + } + } -double rho(double r, double *p) { - int j = 0; - double s,st,ost,os,stemp; - if (r < p[2]) { - ost = os = -1.0e30; - for (j=1;j<=JMAX;j++) { - if (p[4]) st=midinf(kernel,r,p[2],j,p); - else st=midpnt(kernel,r,p[2],j,p); - s=(4.0*st-ost)/3.0; - if (fabs(s-os) < EPS*fabs(os)) break; - if (s == 0.0 && os == 0.0 && j > 6) break; - os=s; - ost=st; - } - stemp = s; - ost = os = -1.0e30; - for (j=1;j<=JMAX;j++) { - st=midinf(kernel,p[2],BIGNUMBER,j,p); - s=(4.0*st-ost)/3.0; - if (fabs(s-os) < EPS*fabs(os)) break; - if (s == 0.0 && os == 0.0 && j > 6) break; - os=s; - ost=st; - } - s += stemp; - } else { - ost = os = -1.0e30; - for (j=1;j<=JMAX;j++) { - st=midinf(kernel,r,BIGNUMBER,j,p); - s=(4.0*st-ost)/3.0; - if (fabs(s-os) < EPS*fabs(os)) break; - if (s == 0.0 && os == 0.0 && j > 6) break; - os=s; - ost=st; - } - } - - return max(s * -1.0/PI, 0.0); //no negative density! - -} + if (create_radial_profile) { + printf("\nRadial density profile:\n\n# R [pc] N [1/pc^3] " + "M [Msun/pc^3] L [Lsun/pc^3] N [1/pc^2] M [Msun/pc^2] " + "L [Lsun/pc^2] \n"); + for (j=0;j=0.5*nmax)) nhalf = rarray[j][0]; + if (!(mhalf) && (mtemp>=0.5*mmax)) mhalf = rarray[j][0]; + if (!(lhalf) && (ltemp>=0.5*lmax)) lhalf = rarray[j][0]; + ntemp2D += 1.0; + mtemp2D += rarray2D[j][1]; + ltemp2D += rarray2D[j][2]; + if (!(nhalf2D) && (ntemp2D>=0.5*nmax)) nhalf2D = rarray2D[j][0]; + if (!(mhalf2D) && (mtemp2D>=0.5*mmax)) mhalf2D = rarray2D[j][0]; + if (!(lhalf2D) && (ltemp2D>=0.5*lmax)) lhalf2D = rarray2D[j][0]; + } + //print radii + printf("\nHalf-number radius = %.4f pc (%.4f pc, 2D)\nHalf-mass " + "radius = %.4f pc (%.4f pc, 2D)\nHalf-light radius = %.4f pc " + "(%.4f pc, 2D)\n",nhalf,nhalf2D,mhalf,mhalf2D,lhalf,lhalf2D); -double rho_kernel (double x, double *p) { - double s; - s = 4.0*PI*x*x*rho(x,p); - return s; -} + } -double M(double r, double *p){ - int j = 0; - double s,st,ost,os; - - ost = os = -1.0e30; - for (j=1;j<=JMAX;j++) { - st=midsql(rho_kernel,0.0,r,j,p); - s=(4.0*st-ost)/3.0; - if (fabs(s-os) < EPS*fabs(os)) break; - if (s == 0.0 && os == 0.0 && j > 6) break; - os=s; - ost=st; - } - return s; - -} + //estimate NNBMAX and RS0 (Nbody6 only) + if ((code == 0) || (code == 2) || code == 4 || (code == 5)) { + *NNBMAX = 2.0*sqrt(N); + if (*NNBMAX < 30) *NNBMAX = 30; + if (N<=*NNBMAX) *NNBMAX = 0.5*N; + if (*NNBMAX > NNBMAX_NBODY6) *NNBMAX = NNBMAX_NBODY6; + *RS0 = rarray[*NNBMAX][0]; + printf("\nEstimating appropriate NNBMAX = %i and RS0 = %f [pc]\n", + *NNBMAX,*RS0); + } -double sigma_kernel (double x, double *p) { - double s; - s = 1.0*rho(x,p)*M(x,p)/(x*x); - return s; -} + for (j=0;j 6) break; - os=s; - ost=st; - } - - if (rho(r,p)) t = sqrt(G*s/rho(r,p)); - else t = 0.0; - - return t; - + return 0; } -double get_gauss(void){ - double random[2],p,q; - do { - random[0] = 2.0*drand48()-1.0; - random[1] = 2.0*drand48()-1.0; - q = random[0]*random[0]+random[1]*random[1]; - } while (q>1.0); - - p = sqrt(-2.0*log(q)/q); - return random[0]*p; - -} +int cmd(double **star, int l, double Rgal, double *abvmag, double *vmag, + double *BV, double *Teff, double *dvmag, double *dBV) { + + double lTeff, BC, kb; + double bvc[8], bcc[8]; + double dbmag; + double BCsun, abvmagsun; + + // Stefan-Boltzmann constant in Lsun Rsun^-2 K^-4 + kb = 5.6704E-08*0.5*1.3914E9*0.5*1.3914E9/3.846E26; + + bvc[0] = -654597.405559323; + bvc[1] = 1099118.61158915; + bvc[2] = -789665.995692672; + bvc[3] = 314714.220932623; + bvc[4] = -75148.4728506455; + bvc[5] = 10751.803394526; + bvc[6] = -853.487897283685; + bvc[7] = 28.9988730655392; + + bcc[0] = -4222907.80590972; + bcc[1] = 7209333.13326442; + bcc[2] = -5267167.04593882; + bcc[3] = 2134724.55938336; + bcc[4] = -518317.954642773; + bcc[5] = 75392.2372207101; + bcc[6] = -6082.7301194776; + bcc[7] = 209.990478646363; + + BCsun = 0.11; //sun's bolometric correction + abvmagsun = 4.83; //sun's absolute V magnitude + + if (star[l][11] && (star[l][8]<14)) { + *Teff = pow(star[l][12]/(4.0*PI*star[l][11]*star[l][11]*kb),0.25); + if ((*Teff>3000.0) && (*Teff<55000.0)) { + + lTeff = log10(*Teff); + + *BV = bvc[0] + bvc[1]*lTeff + bvc[2]*pow(lTeff,2) + + bvc[3]*pow(lTeff,3) + bvc[4]*pow(lTeff,4) + + bvc[5]*pow(lTeff,5) + bvc[6]*pow(lTeff,6) + + bvc[7]*pow(lTeff,7); + + BC = bcc[0] + bcc[1]*lTeff + bcc[2]*pow(lTeff,2) + + bcc[3]*pow(lTeff,3) + bcc[4]*pow(lTeff,4) + + bcc[5]*pow(lTeff,5) + bcc[6]*pow(lTeff,6) + + bcc[7]*pow(lTeff,7); + + if (star[l][12]) + *abvmag = -2.5*log10(star[l][12])-BC+BCsun+abvmagsun; + + *vmag = *abvmag + 5.0*log10(Rgal) - 5.0; + + + double rand1, rand2, prand; + + do { + rand1 = 2.0*drand48()-1.0; + rand2 = 2.0*drand48()-1.0; + } while (rand1*rand1+rand2*rand2 > 1.0); + + prand = sqrt(-2.0*log(rand1*rand1+rand2*rand2)/(rand1*rand1+ + rand2*rand2)); + *dvmag = rand1*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, + 0.4*(*vmag-25.0)),2)); + dbmag = rand2*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, + 0.4*(*vmag-25.0)),2)); + *dBV = *dvmag + dbmag; + + } else { + *vmag = 9999.9; + *abvmag = 9999.9; + *BV = 9999.9; + *dvmag = 0.0; + *dBV = 0.0; + } + } else { + *Teff = 0.0; + *vmag = 9999.9; + *abvmag = 9999.9; + *BV = 9999.9; + *dvmag = 0.0; + *dBV = 0.0; + } -double fractalize(double D, int N, double **star, int radial, int symmetry) { - int i, j, h, Nparent, Nparentlow; - int Ntot = 128.0*pow(8,ceil(log(N)/log(8))); - int Ntotorg = Ntot; - double l = 2.0; - double prob = pow(2.0, D-3.0); - double scatter; - if (radial) scatter = 0.01; - else scatter = 0.1; - double vx, vy, vz; - double vscale; - int subi; - double morescatter = 0.1; //0.1 looks good - - printf("\nFractalizing initial conditions...\n\n"); - - double **star_temp; //temporary array for fractalized structure - star_temp = (double **)calloc(Ntot,sizeof(double *)); - for (j=0;j 1.0)); - for (h=1;h<7;h++) star[i][h] = star_temp[j][h]; - star_temp[j][0] = 0.0; - } - - double r, r_norm, vnorm = 0.0, start[4]; - if (radial) { - for (i=0;i1.0); - for (h=1;h<4;h++) star[i][h] = start[h]; - } - } - for (j=0;j=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) - double ralpha = 1/alpha, eta = alpha/0.23; - lP = pow(pow(lPmin,alpha) + eta*drand48(),ralpha); - } else { - //derive from Sana & Evans (2011) period distribution for massive binaries - double lPmin = 0.3, lPmax = 3.5, lPbreak = 1.0; //parameters of Sana & Evans (2011) period distribution in days (eq. 5.1) - double Fbreak = 0.5; //fraction of binaries with periods below Pbreak - double xperiod = drand48(); - if (xperiod <= Fbreak) - lP = xperiod *(lPbreak-lPmin)/Fbreak + lPmin; - else - lP = (xperiod - Fbreak)*(lPmax-lPbreak)/(1.0-Fbreak) + lPbreak; - } - - P = pow(10.0,lP);//days - P /= 365.25;//yr - - abin = pow((m1+m2)*M*P*P,(1.0/3.0));//AU - abin /= 206264.806;//pc - abin /= rvir;//Nbody units - } else if (adis == 0) { - //flat semi-major axis distribution - if (!i) printf("\nApplying flat semi-major axis distribution with amin = %g and amax = %g.\n", amin, amax); - if (!i) amin /= rvir; - if (!i) amax /= rvir; - abin = amin+drand48()*(amax-amin); - } else if (adis == 1 || adis == 3) { - //derive from Kroupa (1995) period distribution - if (!i) printf("\nDeriving semi-major axis distribution from Kroupa (1995) period distribution.\n"); - double Pmin = 10, delta = 45, eta = 2.5; //parameters of Kroupa (1995) period distribution - do { - lP = log10(Pmin) + sqrt(delta*(exp(2.0*drand48()/eta) - 1.0)); - } while (lP > 8.43); - - P = pow(10.0,lP);//days - P /= 365.25;//yr - - abin = pow((m1+m2)*M*P*P,(1.0/3.0));//AU - abin /= 206264.806;//pc - abin /= rvir;//Nbody units - } else { - //derive from Duquennoy & Mayor (1991) period distribution - lPmean = 4.8; //mean of Duquennoy & Mayor (1991) period distribution [log days] - lPsigma = 2.3; //full width half maximum of Duquennoy & Mayor (1991) period distribution [log days] - if (!i) printf("\nDeriving semi-major axis distribution from Duquennoy & Mayor (1991) period distribution.\n"); - do { - //generate two random numbers in the interval [-1,1] - u1 = 2.0*drand48()-1.0; - u2 = 2.0*drand48()-1.0; - - //combine the two random numbers - q = u1*u1 + u2*u2; - } while (q > 1.0); - - p = sqrt(-2.0*log(q)/q); - x1 = u1*p; - x2 = u2*p; - lP1 = lPsigma*x1 + lPmean; - lP2 = lPsigma*x2 + lPmean; - - P = pow(10,lP1);//days - P /= 365.25;//yr - - abin = pow((m1+m2)*M*P*P,(1.0/3.0));//AU - abin /= 206264.806;//pc - abin /= rvir;//Nbody units - } - +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) { + + //Open output files + char PARfile[50], NBODYfile[50], SSEfile[50]; + FILE *PAR, *NBODY, *SSE12; + sprintf(PARfile, "%s.input",output); + PAR = fopen(PARfile,"w"); + sprintf(NBODYfile, "%s.fort.10",output); + NBODY = fopen(NBODYfile,"w"); + + int hrplot = 0; + if (dtplot) hrplot = 1; + if (sse) { + sprintf(SSEfile, "%s.fort.12",output); + SSE12 = fopen(SSEfile,"w"); + hrplot = 2; + } - if (!i && OBperiods && msort) { - if (adis == 3) { - 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); - } - } - - - //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"); - 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); - } - } else { - ecc = sqrt(drand48()); // Thermal distribution f(e)=2e - } - //ecc = 0.0; // all circular - - - - //Apply Kroupa (1995) eigenevolution - eccold = ecc; - abinold = abin; - m1old = m1; - m2old = m2; - if (eigen) { - if (!i) printf("\nApplying Kroupa (1995) eigenevolution for short-period binaries\n"); - m1*=M;//temporary re-scaling - m2*=M; - abin*=rvir; - if (!(OBperiods && ((m1>=msort) || (m2>=msort)))) { - if (m1>=m2) eigenevolution(&m1, &m2, &ecc, &abin); - else eigenevolution(&m2, &m1, &ecc, &abin); - } - m1/=M; - m2/=M; - abin/=rvir; - } - - - //Apply Binary Star Evolution (Hurley, Tout & Pols 2002) - if (BSE) { - int kw[2] = {1, 1}; //stellar type - double mass[2]; //initial mass - double mt[2]; //actual mass - double r[2] = {0.0, 0.0}; //radius - double lum[2] = {0.0, 0.0}; //luminosity - double mc[2] = {0.0, 0.0}; //core mass - double rc[2] = {0.0, 0.0}; //core radius - double menv[2] = {0.0, 0.0}; //envelope mass - double renv[2] = {0.0, 0.0}; //envelope radius - double ospin[2] = {0.0, 0.0}; //spin - double epoch1[2] = {0.0, 0.0}; //time spent in current evolutionary state - double tms[2] = {0.0, 0.0}; //main-sequence lifetime - double tphys = 0.0; //initial age - double tphysf = epoch;//final age - double dtp = epoch+1; //data store value, if dtp>tphys no data will be stored - double z = Z; //metallicity - - mass[0] = m1*M; //initial mass of primary - mass[1] = m2*M; //initial mass of secondary - mt[0] = mass[0]; //actual mass - mt[1] = mass[1]; //actual mass - vkick[0] = 0.0; - vkick[1] = 0.0; - - abin *= rvir; //pc - abin *= 206264.806; //AU - P = sqrt(pow(abin, 3.0)/(m1+m2));//yr - P *= 365.25;//days - - evolv2_(kw,mass,mt,r,lum,mc,rc,menv,renv,ospin,epoch1,tms,&tphys,&tphysf,&dtp,&z,zpars,&P,&eccold, vkick); - - P /= 365.25;//yr - - abin = pow((m1+m2)*P*P,(1.0/3.0));//AU - abin /= 206264.806;//pc - abin /= rvir;//Nbody units - - m1 = mt[0]/M; //Nbody units - m2 = mt[1]/M; //Nbody units - - star[2*i][8] = kw[0]; - star[2*i+1][8] = kw[1]; - star[2*i][9] = epoch1[0]; - star[2*i+1][9] = epoch1[1]; - star[2*i][10] = ospin[0]; - star[2*i+1][10] = ospin[1]; - star[2*i][11] = r[0]; - star[2*i+1][11] = r[1]; - star[2*i][12] = lum[0]; - star[2*i+1][12] = lum[1]; - } - - //pos & vel in binary frame - ea = rtnewt(ecc, drand48()); - rop[0] = abin*(cos(ea) - ecc); - rop[1] = abin*sqrt(1.0-ecc*ecc)*sin(ea); - - mm = sqrt((m1+m2)/pow(abin,3)); - eadot = mm/(1.0 - ecc*cos(ea)); - vop[0] = -abin*sin(ea)*eadot; - vop[1] = abin*sqrt(1.0-ecc*ecc)*cos(ea)*eadot; - - //Convert to cluster frame - cosi = 2.0*drand48()-1.0; - inc = acos(cosi); - node = 2.0*PI*drand48(); - peri = 2.0*PI*drand48(); - - pmat[0][0] = cos(peri)*cos(node) - sin(peri)*sin(node)*cosi; - pmat[1][0] = cos(peri)*sin(node) + sin(peri)*cos(node)*cosi; - pmat[2][0] = sin(peri)*sin(inc); - pmat[0][1] = -sin(peri)*cos(node) - cos(peri)*sin(node)*cosi; - pmat[1][1] = -sin(peri)*sin(node) + cos(peri)*cos(node)*cosi; - pmat[2][1] = cos(peri)*sin(inc); - - for (j=0;j<3;j++) { - rrel[j] = pmat[j][0]*rop[0] + pmat[j][1]*rop[1]; - vrel[j] = pmat[j][0]*vop[0] + pmat[j][1]*vop[1]; - } - - for (j=0;j<3;j++) { - star[2*i+1][j+1] = star[2*i][j+1] + m1/(m1+m2)*rrel[j]; //Star2 pos - star[2*i+1][j+4] = star[2*i][j+4] + m1/(m1+m2)*vrel[j]; //Star2 vel - star[2*i][j+1] -= m2/(m1+m2)*rrel[j]; //Star1 pos - star[2*i][j+4] -= m2/(m1+m2)*vrel[j]; //Star1 vel - } - - star[2*i][0] = m1; - star[2*i+1][0] = m2; - - } while ((sqrt(vkick[0]) > vesc) && (sqrt(vkick[1]) > vesc)); - - //printf("%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", m1/(m1+m2)*rrel[0], m1/(m1+m2)*rrel[1], m1/(m1+m2)*rrel[2], m1/(m1+m2)*vrel[0], m1/(m1+m2)*vrel[1], m1/(m1+m2)*vrel[2], -m2/(m1+m2)*rrel[0], -m2/(m1+m2)*rrel[1], -m2/(m1+m2)*rrel[2], -m2/(m1+m2)*vrel[0], -m2/(m1+m2)*vrel[1], -m2/(m1+m2)*vrel[2]); - //printf("%i\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\n",i,m1*M, m2*M, P, abin*rvir, ecc, m1old*M, m2old*M, abinold*rvir, eccold); - //printf("%i\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\n",i,star[2*i][4], star[2*i+1][4], star[2*i][1], star[2*i+1][1], star[2*i][1]-star[2*i+1][1], m2old*M, abinold*rvir, eccold); - - if (!star[2*i][0] || !star[2*i+1][0]) { - nbin--; //one binary less because one component went supernova - - double startemp[15]; //temporarily save surviving companion - if (!star[2*i][0]) { - for (k=0;k<15;k++) startemp[k] = star[2*i+1][k]; - } else { - for (k=0;k<15;k++) startemp[k] = star[2*i][k]; - } - - for (j=2*i+2;j<*N;j++) { - for (k=0;k<15;k++) star[j-2][k] = star[j][k]; //move all remaining stars two positions up in the array - } - - *N = *N-1; //reduce total number of stars by one, i.e. remove massless supernova remnant from computations - - for (k=0;k<15;k++) star[*N-1][k] = startemp[k]; //make surviving companion the last particle in array - - i--; //reduce binary counter - } - - } - - return 0; + //write to .PAR file + fprintf(PAR,"1 5000000.0 0\n"); + fprintf(PAR,"%i 1 10 %i %i 1\n",N,seed,NNBMAX); + fprintf(PAR,"0.02 0.02 %.8f %.8f %.8f %.8f 1.0E-03 %.8f %.8f\n", + RS0,dtadj,dtout,tcrit,rvir,mmean); + fprintf(PAR,"2 2 1 0 1 0 2 0 0 2\n"); + fprintf(PAR,"0 %i 0 %i 2 %i %i 0 %i 3\n", + hrplot,tf,regupdate,etaupdate,mloss); + fprintf(PAR,"0 %i %i 0 1 2 0 1 0 1\n",bin,esc); + fprintf(PAR,"0 0 0 2 1 0 0 2 0 3\n"); + fprintf(PAR,"0 0 0 0 0 0 0 0 0 0\n"); + fprintf(PAR,"1.0E-5 1.0E-4 0.2 1.0 1.0E-06 0.001\n"); + fprintf(PAR,"2.350000 %.8f %.8f %i 0 %.8f %.8f %.8f\n", + MMAX,mlow,nbin,Z,epoch,dtplot); + fprintf(PAR,"%.2f 0.0 0.0 0.00000 0.125\n",Q); + if (tf == 2) { + fprintf(PAR,"%.8e %.8f\n", + M1pointmass,sqrt(RG[0]*RG[0]+RG[1]*RG[1]+RG[2]*RG[2])/1000.0); + } else if (tf == 3) { + //old version: + //fprintf(PAR,"%.6e %.6e %.6f %.6f %.6f %.6f %.6f %.6f %.6f %.6f " + // "%.6f %.6f\n",M1allen,M2allen,a2allen,b2allen,VCIRC,RCIRC, + // RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); + + //new version including bulge potential: + fprintf(PAR,"%.6e %.6e %.6f %.6f %.6f %.6f %.6e %.6f %.6f\n", + M1allen,M2allen,a2allen,b2allen,VCIRC,RCIRC, GMB, AR, GAM); + fprintf(PAR,"%.6f %.6f %.6f %.6f %.6f %.6f\n", + RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); + } + if (tf > 2) { + fprintf(PAR,"%.6f %.6f %.6f %.6f\n", extmass,extrad,extdecay,extstart); + } -} - -void shellsort(double **array, int N, int k) {//largest up - int i,j,l,n; - N = N-1; - double swap[k]; - //guess distance n - for (n = 1; n <= N/9; n = 3*n+1); - for (; n > 0; n /= 3) { - for (i = n; i <= N; i++) { - for (l=0; l= n) && (array[j-n][0] < swap[0])); j -= n) { - for (l=0; l 0; n /= 3) { - for (i = n; i <= N; i++) { - for (l=0; l= n) && (array[j-n][0] > swap[0])); j -= n) { - for (l=0; l 0; n /= 3) { - for (i = n; i <= N; i++) { - swap = array[i]; - for (j = i; ((j >= n) && (array[j-n] < swap)); j -= n) { - array[j] = array[j-n]; - } - array[j] = swap; - } - } -} -void shellsort_reverse_1d(double *array, int N) {//smallest up - int i,j,n; - N = N-1; - double swap; - //guess distance n - for (n = 1; n <= N/9; n = 3*n+1); - for (; n > 0; n /= 3) { - for (i = n; i <= N; i++) { - swap = array[i]; - for (j = i; ((j >= n) && (array[j-n] > swap)); j -= n) { - array[j] = array[j-n]; - } - array[j] = swap; - } - } -} + //write to .fort.10 file + int j; + for (j=0;j= msort) { - // First star in binarie - star_temp[j][0] = star[(int) masses[i][1]][0]; - star_temp[j][1] = star[(int) masses[i][1]][1]; - star_temp[j][2] = star[(int) masses[i][1]][2]; - star_temp[j][3] = star[(int) masses[i][1]][3]; - star_temp[j][4] = star[(int) masses[i][1]][4]; - star_temp[j][5] = star[(int) masses[i][1]][5]; - star_temp[j][6] = star[(int) masses[i][1]][6]; - star_temp[j][7] = star[(int) masses[i][1]][7]; - star_temp[j][8] = star[(int) masses[i][1]][8]; - star_temp[j][9] = star[(int) masses[i][1]][9]; - star_temp[j][10] = star[(int) masses[i][1]][10]; - star_temp[j][11] = star[(int) masses[i][1]][11]; - star_temp[j][12] = star[(int) masses[i][1]][12]; - star_temp[j][13] = star[(int) masses[i][1]][13]; - star_temp[j][14] = star[(int) masses[i][1]][14]; - mask[i] = 1; - j++; - // Find the second one based on uniform mass ratio - if (i mpair) { - k++; - if (k>=N) { - k--; - break; - } - } - int k1 = k-1; - int k2 = k; - if(k1==i) { - while(k1i && mask[k1]==1) k1--; - while(k20) { - if(j>=N) { - perror("\n Error!: Pairing binary star exceed the maximum index!\n"); - exit(0); - } - // Second star in binarie - star_temp[j][0] = star[(int) masses[k][1]][0]; - star_temp[j][1] = star[(int) masses[k][1]][1]; - star_temp[j][2] = star[(int) masses[k][1]][2]; - star_temp[j][3] = star[(int) masses[k][1]][3]; - star_temp[j][4] = star[(int) masses[k][1]][4]; - star_temp[j][5] = star[(int) masses[k][1]][5]; - star_temp[j][6] = star[(int) masses[k][1]][6]; - star_temp[j][7] = star[(int) masses[k][1]][7]; - star_temp[j][8] = star[(int) masses[k][1]][8]; - star_temp[j][9] = star[(int) masses[k][1]][9]; - star_temp[j][10] = star[(int) masses[k][1]][10]; - star_temp[j][11] = star[(int) masses[k][1]][11]; - star_temp[j][12] = star[(int) masses[k][1]][12]; - star_temp[j][13] = star[(int) masses[k][1]][13]; - star_temp[j][14] = star[(int) masses[k][1]][14]; - mask[k] = 1; - j++; - } else { - printf("Mass %g did not find good pair, expected pair mass: %g\n",masses[i][0],mpair); - } - } - } - else { - // Store index for remaining stars - mmrand[ileft][0] = drand48(); - mmrand[ileft][1] = masses[i][1]; - ileft++; - } + //write to .fort.12 file + if (sse) { + for (j=0;jMsort binaries + MMsort: %d, MMsort: %d",j/2); - } + } - // Randomize indexs for random pairing - shellsort(mmrand,ileft,2); + fclose(PAR); + fclose(NBODY); + if (bin == 5) { + fclose(SSE12); + printf("\nData written to %s, %s and %s\n", + PARfile, NBODYfile, SSEfile); + } else { + printf("\nData written to %s and %s\n", PARfile, NBODYfile); + } - // Random pairing remaining binaries - for(i=0;i= msort) {//copying to front of temporary array if massive - star_temp[j][0] = star[(int) masses[i][1]][0]; - star_temp[j][1] = star[(int) masses[i][1]][1]; - star_temp[j][2] = star[(int) masses[i][1]][2]; - star_temp[j][3] = star[(int) masses[i][1]][3]; - star_temp[j][4] = star[(int) masses[i][1]][4]; - star_temp[j][5] = star[(int) masses[i][1]][5]; - star_temp[j][6] = star[(int) masses[i][1]][6]; - star_temp[j][7] = star[(int) masses[i][1]][7]; - star_temp[j][8] = star[(int) masses[i][1]][8]; - star_temp[j][9] = star[(int) masses[i][1]][9]; - star_temp[j][10] = star[(int) masses[i][1]][10]; - star_temp[j][11] = star[(int) masses[i][1]][11]; - star_temp[j][12] = star[(int) masses[i][1]][12]; - star_temp[j][13] = star[(int) masses[i][1]][13]; - star_temp[j][14] = star[(int) masses[i][1]][14]; - j++; - } - } - if (msort) Nhighmass = j; - - for (i=0;i 0) { - *ecc = exp(alpha+log(*ecc)); - } - - /* pre-ms mass-transfer */ - qold = *m1/ *m2; - if (qold > 1.0) qold = 1.0/qold; - alpha = -alpha; - if (alpha > 1.0) { - qnew = 1.0; - } else { - qnew = qold + (1.0-qold)*alpha; - } - - *m1 = max(*m1,*m2); - *m2 = qnew * *m1; - - mtot = *m1 + *m2; - - lper = lperi*pow((1.0-ecci)/(1.0-*ecc),1.5); - lper = lper*sqrt(mtoti/mtot); - - *abin = pow(mtot*lper*lper*GBIN/4.0/PI/PI,1.0/3.0); - *abin /= PARSEC; - - return 0; } -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 i, j; - *Rh2D = 0.0; - *Rh3D = 0.0; - double Mtemp; - double **rarray; - rarray = (double **)calloc(N,sizeof(double *)); - for (j=0;j=0.5*M) && !(*Rh3D)) *Rh3D = rarray[i][0]; - i++; - } else { - j++; - } - } - - j = 0; i = 0; Mtemp = 0.0; - while ((j < noofradialbins) && (i < N)) { - if (rarray2D[i][0] < rprofile[j][0]) {//2D binning in astrophysical units - rprofile[j][6] += 1.0; - rprofile[j][7] += rarray2D[i][1]; - rprofile[j][8] += rarray2D[i][2]; - Mtemp += rarray2D[i][1]; - if ((Mtemp>=0.5*M) && !(*Rh2D)) *Rh2D = rarray2D[i][0]; - i++; - } else { - j++; - } - } - - if (create_radial_profile) { - printf("\nRadial density profile:\n\n# R [pc] N [1/pc^3] M [Msun/pc^3] L [Lsun/pc^3] N [1/pc^2] M [Msun/pc^2] L [Lsun/pc^2] \n"); - for (j=0;j=0.5*nmax)) nhalf = rarray[j][0]; - if (!(mhalf) && (mtemp>=0.5*mmax)) mhalf = rarray[j][0]; - if (!(lhalf) && (ltemp>=0.5*lmax)) lhalf = rarray[j][0]; - ntemp2D += 1.0; - mtemp2D += rarray2D[j][1]; - ltemp2D += rarray2D[j][2]; - if (!(nhalf2D) && (ntemp2D>=0.5*nmax)) nhalf2D = rarray2D[j][0]; - if (!(mhalf2D) && (mtemp2D>=0.5*mmax)) mhalf2D = rarray2D[j][0]; - if (!(lhalf2D) && (ltemp2D>=0.5*lmax)) lhalf2D = rarray2D[j][0]; - } - printf("\nHalf-number radius = %.4f pc (%.4f pc, 2D)\nHalf-mass radius = %.4f pc (%.4f pc, 2D)\nHalf-light radius = %.4f pc (%.4f pc, 2D)\n",nhalf,nhalf2D,mhalf,mhalf2D,lhalf,lhalf2D); //print radii - - } - - - //estimate NNBMAX and RS0 (Nbody6 only) - if ((code == 0) || (code == 2) || code == 4 || (code == 5)) { - *NNBMAX = 2.0*sqrt(N); - if (*NNBMAX < 30) *NNBMAX = 30; - if (N<=*NNBMAX) *NNBMAX = 0.5*N; - if (*NNBMAX > NNBMAX_NBODY6) *NNBMAX = NNBMAX_NBODY6; - *RS0 = rarray[*NNBMAX][0]; - printf("\nEstimating appropriate NNBMAX = %i and RS0 = %f [pc]\n",*NNBMAX,*RS0); - } - - for (j=0;j3000.0) && (*Teff<55000.0)) { - - lTeff = log10(*Teff); - - *BV = bvc[0] + bvc[1]*lTeff + bvc[2]*pow(lTeff,2) + bvc[3]*pow(lTeff,3) + bvc[4]*pow(lTeff,4) + bvc[5]*pow(lTeff,5) + bvc[6]*pow(lTeff,6) + bvc[7]*pow(lTeff,7); - - BC = bcc[0] + bcc[1]*lTeff + bcc[2]*pow(lTeff,2) + bcc[3]*pow(lTeff,3) + bcc[4]*pow(lTeff,4) + bcc[5]*pow(lTeff,5) + bcc[6]*pow(lTeff,6) + bcc[7]*pow(lTeff,7); - - if (star[l][12]) *abvmag = -2.5*log10(star[l][12])-BC+BCsun+abvmagsun; - - *vmag = *abvmag + 5.0*log10(Rgal) - 5.0; - - - double rand1, rand2, prand; - - do { - rand1 = 2.0*drand48()-1.0; - rand2 = 2.0*drand48()-1.0; - } while (rand1*rand1+rand2*rand2 > 1.0); - - prand = sqrt(-2.0*log(rand1*rand1+rand2*rand2)/(rand1*rand1+rand2*rand2)); - *dvmag = rand1*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2)); - dbmag = rand2*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2)); - *dBV = *dvmag + dbmag; - - } else { - *vmag = 9999.9; - *abvmag = 9999.9; - *BV = 9999.9; - *dvmag = 0.0; - *dBV = 0.0; - } - } else { - *Teff = 0.0; - *vmag = 9999.9; - *abvmag = 9999.9; - *BV = 9999.9; - *dvmag = 0.0; - *dBV = 0.0; - } - - return 0; -} - -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){ - - //Open output files - char PARfile[50], NBODYfile[50], SSEfile[50]; - FILE *PAR, *NBODY, *SSE12; - sprintf(PARfile, "%s.input",output); - PAR = fopen(PARfile,"w"); - sprintf(NBODYfile, "%s.fort.10",output); - NBODY = fopen(NBODYfile,"w"); - - int hrplot = 0; - if (dtplot) hrplot = 1; - if (sse) { - sprintf(SSEfile, "%s.fort.12",output); - SSE12 = fopen(SSEfile,"w"); - hrplot = 2; - } - - //write to .PAR file - fprintf(PAR,"1 5000000.0 0\n"); - fprintf(PAR,"%i 1 10 %i %i 1\n",N,seed,NNBMAX); - fprintf(PAR,"0.02 0.02 %.8f %.8f %.8f %.8f 1.0E-03 %.8f %.8f\n",RS0,dtadj,dtout,tcrit,rvir,mmean); - fprintf(PAR,"2 2 1 0 1 0 2 0 0 2\n"); - fprintf(PAR,"0 %i 0 %i 2 %i %i 0 %i 3\n",hrplot,tf,regupdate,etaupdate,mloss); - fprintf(PAR,"0 %i %i 0 1 2 0 1 0 1\n",bin,esc); - fprintf(PAR,"0 0 0 2 1 0 0 2 0 3\n"); - fprintf(PAR,"0 0 0 0 0 0 0 0 0 0\n"); - fprintf(PAR,"1.0E-5 1.0E-4 0.2 1.0 1.0E-06 0.001\n"); - fprintf(PAR,"2.350000 %.8f %.8f %i 0 %.8f %.8f %.8f\n",MMAX,mlow,nbin,Z,epoch,dtplot); - fprintf(PAR,"%.2f 0.0 0.0 0.00000 0.125\n",Q); - if (tf == 2) { - fprintf(PAR,"%.8e %.8f\n",M1pointmass,sqrt(RG[0]*RG[0]+RG[1]*RG[1]+RG[2]*RG[2])/1000.0); - } else if (tf == 3) { - //old version: - //fprintf(PAR,"%.6e %.6e %.6f %.6f %.6f %.6f %.6f %.6f %.6f %.6f %.6f %.6f\n",M1allen,M2allen,a2allen,b2allen,VCIRC,RCIRC,RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); - - //new version including bulge potential: - fprintf(PAR,"%.6e %.6e %.6f %.6f %.6f %.6f %.6e %.6f %.6f\n",M1allen,M2allen,a2allen,b2allen,VCIRC,RCIRC, GMB, AR, GAM); - fprintf(PAR,"%.6f %.6f %.6f %.6f %.6f %.6f\n", RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); - } - if (tf > 2) fprintf(PAR,"%.6f %.6f %.6f %.6f\n",extmass,extrad,extdecay,extstart); - - - - - //write to .fort.10 file - int j; - for (j=0;j 2) { + fprintf(PAR,"%.2f %.2f %.2f %.2f\n", extmass,extrad,extdecay,extstart); + } + + //write to .NBODY file + int j; + for (j=0;j 2) fprintf(PAR,"%.2f %.2f %.2f %.2f\n",extmass,extrad,extdecay,extstart); - - - - //write to .NBODY file - int j; - for (j=0;j 2) fprintf(PAR,"%.6f %.6f %.6f %.6f\n",extmass,extrad,extdecay,extstart); - fprintf(PAR,"10000.0 2 0\n"); - - //write to .fort.10 file - int j; - for (j=0;j 2) + fprintf(PAR,"%.6f %.6f %.6f %.6f\n",extmass,extrad,extdecay,extstart); + + fprintf(PAR,"10000.0 2 0\n"); + + //write to .fort.10 file + int j; + for (j=0;j0?2:0)); - fprintf(PAR,"0 %i 0 %i 2 %i %i 0 %i 6\n",hrplot,tf,regupdate,etaupdate,mloss); - fprintf(PAR,"0 6 %i 0 1 2 2 0 0 1\n", esc); - fprintf(PAR,"1 0 3 2 1 0 0 2 0 0\n"); - fprintf(PAR,"0 0 0 0 0 1 2 0 0 0\n"); - fprintf(PAR,"1.0E-5 1.0E-4 0.2 1.0 1.0E-06 0.001 0.125\n"); - fprintf(PAR,"2.350000 %.8f %.8f %i 0 %.8f %.8f %.8f\n",MMAX,mlow,nbin,Z,epoch,dtplot); - fprintf(PAR,"%.2f 0.0 0.0 0.00000\n",Q); -// if (tf == 1) { -// rtide = pow(1.0*M/(3.0*M1pointmass),1.0/3.0)*sqrt(RG[0]*RG[0]+RG[1]*RG[1]+RG[2]*RG[2]); -// fprintf(PAR,"%i %.8f\n",0,sqrt(1.0/(3.0*pow(rtide,3)))); - if (tf == 2) { - fprintf(PAR,"%.8e %.8f\n",M1pointmass,sqrt(RG[0]*RG[0]+RG[1]*RG[1]+RG[2]*RG[2])/1000.0); - } else if (tf == 3) { - //fprintf(PAR,"%.6e %.6f %.6e %.6f %.6f %.6e %.6f %.6f %.6f %.6f %.6f %.6f %.6f\n",M1allen,b1allen,M2allen,a2allen,b2allen,M3allen,a3allen,RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); - //new version including bulge potential: - fprintf(PAR,"%.6e %.6f %.6e %.6f %.6f %.6e %.6f %.6f %.6f %.6f\n",M1allen,b1allen,M2allen,a2allen,b2allen,M3allen,a3allen, GMB, AR, GAM); - fprintf(PAR,"%.6f %.6f %.6f %.6f %.6f %.6f\n", RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); - } - if (tf > 2) fprintf(PAR,"%.2f %.2f %.2f %.2f\n",extmass,extrad,extdecay,extstart); - - - - //write to .NBODY file - int j; - for (j=0;j0?2:0)); + fprintf(PAR,"0 %i 0 %i 2 %i %i 0 %i 6\n", + hrplot,tf,regupdate,etaupdate,mloss); + fprintf(PAR,"0 6 %i 0 1 2 2 0 0 1\n", esc); + fprintf(PAR,"1 0 3 2 1 0 0 2 0 0\n"); + fprintf(PAR,"0 0 0 0 0 1 2 0 0 0\n"); + fprintf(PAR,"1.0E-5 1.0E-4 0.2 1.0 1.0E-06 0.001 0.125\n"); + fprintf(PAR,"2.350000 %.8f %.8f %i 0 %.8f %.8f %.8f\n", + MMAX,mlow,nbin,Z,epoch,dtplot); + fprintf(PAR,"%.2f 0.0 0.0 0.00000\n",Q); +// if (tf == 1) { +// rtide = pow(1.0*M/(3.0*M1pointmass),1.0/3.0)* +// sqrt(RG[0]*RG[0]+RG[1]*RG[1]+RG[2]*RG[2]); +// fprintf(PAR,"%i %.8f\n",0,sqrt(1.0/(3.0*pow(rtide,3)))); + if (tf == 2) { + fprintf(PAR,"%.8e %.8f\n",M1pointmass, + sqrt(RG[0]*RG[0]+RG[1]*RG[1]+RG[2]*RG[2])/1000.0); + } else if (tf == 3) { + //fprintf(PAR,"%.6e %.6f %.6e %.6f %.6f %.6e %.6f %.6f %.6f %.6f " + // "%.6f %.6f %.6f\n",M1allen,b1allen,M2allen,a2allen,b2allen,M3allen, + // a3allen,RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); + //new version including bulge potential: + fprintf(PAR,"%.6e %.6f %.6e %.6f %.6f %.6e %.6f %.6f %.6f %.6f\n", + M1allen,b1allen,M2allen,a2allen,b2allen,M3allen,a3allen, + GMB, AR, GAM); + fprintf(PAR,"%.6f %.6f %.6f %.6f %.6f %.6f\n", + RG[0]/1000.0,RG[1]/1000.0,RG[2]/1000.0,VG[0],VG[1],VG[2]); + } + if (tf > 2) + fprintf(PAR,"%.2f %.2f %.2f %.2f\n",extmass,extrad,extdecay,extstart); + + //write to .NBODY file + int j; + for (j=0;j (number of stars) \n"); - printf(" -M (mass of cluster; specify either N or M) \n"); - printf(" -P <0|1|2|3|-1> (density profile; 0= Plummer, 1= King (1966), \n"); - printf(" 2= Subr et al. (2007) mass-segregated, \n"); - printf(" 3= 2-dimensional EFF/Nuker template, \n"); - printf(" -1= no density gradient) \n"); - printf(" -W <1-12> (W0 parameter for King model) \n"); - printf(" -R (half-mass radius [pc], ignored for P = 3) \n"); - printf(" -r (scale radius of EFF/Nuker template [pc]) \n"); - printf(" -c (cut-off radius of EFF/Nuker template [pc]) \n"); - printf(" -g (power-law slope(s) of EFF/Nuker template; use \n"); - printf(" once for EFF template; use three times for Nuker \n"); - printf(" template (outer slope, inner slope, transition) \n"); - printf(" -S <0.0-1.0> (degree of mass segregation; 0.0= no segregation)\n"); - printf(" -D <1.6-3.0> (fractal dimension; 3.0= no fractality) \n"); - printf(" -T (tcrit in N-body units) \n"); - printf(" -Q (virial ratio) \n"); - printf(" -C <0|1|3|5> (code; 0= Nbody6, 1= Nbody4, 3= table of stars, 5= Nbody6++) \n"); - printf(" -A (dtadj in N-body units) \n"); - printf(" -O (deltat in N-body units) \n"); - printf(" -G <0|1> (GPU usage; 0= no GPU, 1= use GPU) \n"); - printf(" -o (output name of cluster model) \n"); - printf(" -f <0|1|2|3|4> (IMF; 0= no IMF, 1= Kroupa (2001), \n"); - printf(" 2= user defined, 3= Kroupa (2001) with optimal sampling,\n"); - printf(" 4= L3 IMF (Maschberger 2012)) \n"); - printf(" -a (IMF slope; for user defined IMF, may be used \n"); - printf(" multiple times, from low mass to high mass; \n"); - printf(" for L3 IMF use three times for alpha, beta and mu)\n"); - printf(" -m (IMF mass limits, has to be used multiple times \n"); - printf(" (at least twice), from low mass to high mass [Msun])\n"); - printf(" -B (number of binary systems) \n"); - printf(" -b (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%.1f Msun)\n",msort); - printf(" -s (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"); - printf(" -e (epoch for stellar evolution [Myr]) \n"); - printf(" -Z (metallicity [0.0001-0.03, 0.02 = solar]) \n"); - printf(" -X (galactocentric radius vector, use 3x, [pc]) \n"); - printf(" -V (cluster velocity vector, use 3x, [km/s]) \n"); - printf(" -x (specify external (gas) Plummer potential, use 4x, \n"); - printf(" 1. gas mass [Msun], 2. Plummer radius [pc] \n"); - printf(" 3. decay time for gas expulsion [Myr], 4. delay \n"); - printf(" time for start of gas expulsion [Myr]) \n"); - printf(" -u <0|1> (output units; 0= Nbody, 1= astrophysical) \n"); - printf(" -h (display this help) \n"); - printf(" -? (display this help) \n"); - printf(" \n"); - printf(" Examples: mcluster -N 1000 -R 0.8 -P 1 -W 3.0 -f 1 -B 100 -o test1 \n"); - printf(" mcluster -f 2 -m 0.08 -a -1.35 -m 0.5 -a -2.7 -m 100.0 \n"); - printf(" mcluster -t 3 -X 8500 -X 0 -X 0 -V 0 -V 220 -V 0 \n"); - printf(" mcluster -D 1.6 -Q 0.4 -P -1 \n"); - printf(" \n"); - -} - + printf("\n Usage: mcluster -[N|M|P|W|R|r|c|g|S|D|T|Q|C|A|O|G|o|f|a|m|B|b|p|s|t|e|Z|X|V|x|u|h|?]\n"); + printf(" \n"); + printf(" -N (number of stars) \n"); + printf(" -M (mass of cluster; specify either N or M) \n"); + printf(" -P <0|1|2|3|-1> (density profile; 0= Plummer, 1= King (1966), \n"); + printf(" 2= Subr et al. (2007) mass-segregated, \n"); + printf(" 3= 2-dimensional EFF/Nuker template, \n"); + printf(" -1= no density gradient) \n"); + printf(" -W <1-12> (W0 parameter for King model) \n"); + printf(" -R (half-mass radius [pc], ignored for P = 3) \n"); + printf(" -r (scale radius of EFF/Nuker template [pc]) \n"); + printf(" -c (cut-off radius of EFF/Nuker template [pc]) \n"); + printf(" -g (power-law slope(s) of EFF/Nuker template; use \n"); + printf(" once for EFF template; use three times for Nuker \n"); + printf(" template (outer slope, inner slope, transition) \n"); + printf(" -S <0.0-1.0> (degree of mass segregation; 0.0= no segregation)\n"); + printf(" -D <1.6-3.0> (fractal dimension; 3.0= no fractality) \n"); + printf(" -T (tcrit in N-body units) \n"); + printf(" -Q (virial ratio) \n"); + printf(" -C <0|1|3|5> (code; 0= Nbody6, 1= Nbody4, 3= table of stars, 5= Nbody6++) \n"); + printf(" -A (dtadj in N-body units) \n"); + printf(" -O (deltat in N-body units) \n"); + printf(" -G <0|1> (GPU usage; 0= no GPU, 1= use GPU) \n"); + printf(" -o (output name of cluster model) \n"); + printf(" -f <0|1|2|3|4> (IMF; 0= no IMF, 1= Kroupa (2001), \n"); + printf(" 2= user defined, 3= Kroupa (2001) with optimal sampling,\n"); + printf(" 4= L3 IMF (Maschberger 2012)) \n"); + printf(" -a (IMF slope; for user defined IMF, may be used \n"); + printf(" multiple times, from low mass to high mass; \n"); + printf(" for L3 IMF use three times for alpha, beta and mu)\n"); + printf(" -m (IMF mass limits, has to be used multiple times \n"); + printf(" (at least twice), from low mass to high mass [Msun])\n"); + printf(" -B (number of binary systems) \n"); + printf(" -b (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%.1f Msun)\n",msort); + printf(" -s (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"); + printf(" -e (epoch for stellar evolution [Myr]) \n"); + printf(" -Z (metallicity [0.0001-0.03, 0.02 = solar]) \n"); + printf(" -X (galactocentric radius vector, use 3x, [pc]) \n"); + printf(" -V (cluster velocity vector, use 3x, [km/s]) \n"); + printf(" -x (specify external (gas) Plummer potential, use 4x, \n"); + printf(" 1. gas mass [Msun], 2. Plummer radius [pc] \n"); + printf(" 3. decay time for gas expulsion [Myr], 4. delay \n"); + printf(" time for start of gas expulsion [Myr]) \n"); + printf(" -u <0|1> (output units; 0= Nbody, 1= astrophysical) \n"); + printf(" -h (display this help) \n"); + printf(" -? (display this help) \n"); + printf(" \n"); + printf(" Examples: mcluster -N 1000 -R 0.8 -P 1 -W 3.0 -f 1 -B 100 -o test1 \n"); + printf(" mcluster -f 2 -m 0.08 -a -1.35 -m 0.5 -a -2.7 -m 100.0 \n"); + printf(" mcluster -t 3 -X 8500 -X 0 -X 0 -V 0 -V 220 -V 0 \n"); + printf(" mcluster -D 1.6 -Q 0.4 -P -1 \n"); + printf(" \n"); +}