GalactICS - A Galaxy Model Building Package First do the following: 1. cd src 2. type 'make all' - This builds the programs. Note that some of the code is in fortran 77, using lines longer than 72 characters in some cases. The -e flag in the makefile allow for this for a Solaris f77 compiler. Other programs are written in C. Again, the linking between these routines works on solaris systems, but may need to be adjusted for other architectures. We have found that linking using f77 instead of ld will often automatically load the appropriate libraries. The graphics output by some of the programs (dbh, plotforce, diskdf, plothalo) use the PGPLOT library. It can be found at http://astro.caltech.edu/~tjp/pgplot/. Alternatively, remove all calls to routines with names starting with "PG", as well as the -lpgplot flag in the Makefile, and the programs should still run fine. 3. type 'make install' - This copies programs to the directory bin/ Now you're set. Test that the programs run: Go into directory Milky_Way/A Type "make galaxy" - this should build all the files and a small N-body galaxy - look at all of the in.* files to see the various parameters that go into building the models. And read below for more details. The distribution contains the input files for Milky_Way/A through D, the models that were used in Kuijken & Dubinski 1995. ----------- This package contains a set of programs and subroutines for building galaxy models including a disk, a bulge and halo. The details of the inner workings of the code are described in Kuijken and Dubinski 1995. There are 3 steps in building these models. A. Calculating the potential. B. Constructing a disk distribution function which will generate the given potential. C. And realizing each component with a self-consistent distribution of particle orbits. These actions are all performed by typing `make galaxy', which runs a succession of programs to end up with a set of N-body particle masses, positions and velocities representing your model. Descriptions of the individual steps: A. The Potential Program: dbh Sample input file: in.dbh Output: dbh.dat - contains tabulated values of the harmonic coefficients for the Legendre expansion of the density, potential and radial force at the specified radii for the entire model h.dat - same as above for halo only b.dat - same as above for bulge only mr.dat - gives mass and radial extent (or edge) of disk,bulge and halo Parameters in in.dbh: y #yes we want a halo (or no) -6.0 1.32 1 .1 0.8 #psi0, v0, q, (rc/rk)^2, ra y #yes we want a disk (or no) .867 1 5 .1 .5 #M_d, R_d, R_outer, z_d, dR_trunc y #yes we want a bulge (or no) 14.45 -3.7 .714 #rho_b, psi_cut, sig_b .01 5000 #delta_r, nr 10 #number of harmonics (even number) dbh.ps/ps # PGPLOT graphics device for the plots produced. The program asks for parameters describing each of the components. You have the option of including any combination of components (though I think models without a halo won't work). Halo Parameters: psi0 - central potential - the smaller (the more negative) this parameter the deeper the potential and the more extended the halo v0 - v0 = sqrt(2.0)*sigma0 where sigma0 is the central velocity dispersion. roughly the velocity where the halo rotation curve peaks q - an optional flattening parameter for the potential - generally 0.7 < q < 1.05 - q=1.0 will give a nearly spherical halo (rc/rk)^2 - a core smoothing parameter - ratio of the core radius to the derived King radius for halo only models set this to 1.0. For multicomponent models, this can be a smaller number 0.0 to 0.1. I've found that with this parameter=0.0 the program can crash. Ra - a scaling radius for the halo - The halo Ra radius is the radius at which the halo rotation curve, at its initial slope ignoring cutoffs and the other components, reaches v0. Disk Parameters: M_d - mass of the exponential disk ignoring cutoffs R_d - exponential scale length R_outer - outer radius where we begin to truncate the disk density z_d - disk scale height assuming a sech^2(z/zd) vertical density law dR_trunc - truncation width - the disk density smoothly drops to zero in the range R_outer < R < ~R_outer + 2*dR_trunc. Bulge Parameters: rho_b - bulge central density psi_cut - bulge cut-off potential psi0 < psi_cut < 0.0 - energy cut-off for the bulge sig_b - bulge central potential Potential parameters: dr - the width of the radial bins used to calculate the potential nr - number of radial bins - initially a guess since we don't know the radial extent of the system lmax - the largest value in the potential harmonic expansion - use lmax=2 to get a quick look at the mass profile and lmax=10 for the final calculation of the model Creating a galaxy model from these parameters is sort of a black art since the halo and bulge models are not parameterized in terms of their mass profiles but rather properties of their distribution functions. Changes in psi0, v0 etc. have weird but predictable effects on the mass profile. The halo is a flattened analogue of the King model so the concentration (R_tidal/R_core) is determined by the dimensionless central potential $\psi_0/\sigma_0^2$. The more negative the value the greater the concentration. The parameters $R_a$ and $v_0$, affect the scaling of the halo mass profile. The effect of different bulge parameters is more predictable. Decreasing the central velocity dispersion will create a more centrally concentrated bulge and decreasing the psi cut off will truncate the bulge and decrease its total mass. The disk is parameterized directly by its mass profile so its effect on the rotation curve is predictable ahead of time. Hit and miss seems to be a good strategy for finding a suitable profile. Generate a model to lmax=2 and then view the resulting rotation curve by typing 'make vr.dat'. This uses the program to generate the file vr.dat which tells you the contributions to the total rotation curve. Another useful file is 'mr.dat' which tells you the mass and radial extent of the disk bulge and halo. The program plotforce will also generate the rotation curves for you directly from the dbh.dat, b.dat and h.dat files. The potential is determined iteratively: starting from an initial guess at the potential, the density implied by the halo and bulge DFs is calculated, the disk density added, and the potential of that mass distribution is used as starting point for the next iteration. Initially only the monopole (l=0) components are calculated until the model converges, then one more harmonic is added per iteration up to the maximum requested, and once all harmonics are included the iterations are continued until the outer (tidal) radius of the halo is unchanged between iterations. At each iterations plots of the harmonic expansion coefficients are produced. If the tidal radius reported is "outside grid" for a large number of iterations, increase the number of radial bins or increase their size. Sometimes infinite tidal radii are also reported: this happens when the total mass of the model using the current guess for the potential is insufficient to generate a potential well as deep as requested. If this persists over many iterations, again increase the number or size of the radial bins. B. Disk distribution function Program: getfreqs Input files: dbh.dat h.dat b.dat Output: freqdbh.dat getfreqs tabulates various characteristic frequencies (omega, kappa etc.) in the equatorial plane for use by diskdf below. Program: diskdf Input files: freqdbh.dat dbh.dat in.diskdf Output files: cordbh.dat toomre.dat The program diskdf iteratively calculates the correction functions for the disk distribution function. These functions are multiplicative corrections to the surface density and vertical velocity dispersion which appear to leading order in the Shu (1969) distribution functions. See KD95 for details. It requires the sample parameters: .47 1.0 #central radial vel. dispersion, exponential scale length of sig_r^2 50 #number of radial intervals for correction functions 10 #number of iterations diskdf.ps/ps # PGPLOT device for plot of correction functions. It also outputs the Toomre Q as a function of radius in the file toomre.dat. C. Generating N-body realizations Programs: gendisk, genbulge, genhalo Input files: cordbh.dat dbh.dat gendisk parameters: 4000 #number of particles -1 #negative random integer seed 1 #1=yes we want to center 0=no we don't dbh.dat #multipole expansion data file genbulge parameters: 0.5 #streaming fraction 1000 #number of particles -1 #negative integer seed 1 #center the data 1=yes dbh.dat #harmonics file genhalo parameters: 0.5 #streaming fraction 6000 #number of particles -1 #negative integer seed for random number generator 1 #1=yes we want to center dbh.dat #multipole expansion data file The streaming fraction, f, sets the fraction of orbits with L_z > 0. The remaining fraction, 1-f, have L_z<0. With this parameter you can therefore vary the rotation of the bulge and the halo. f=0.5 refers to the non-rotating case. The N-body data are written to the stdout so the programs should be run as: gendisk < in.disk > disk genbulge < in.bulge > bulge genhalo < in.halo > halo Format is ascii with data arranged as: N_bodies time m_1 x_1 y_1 z_1 vx_1 vy_1 vz_1 m_2 x_2 y_2 z_2 vx_2 vy_2 vz_2 m_3 x_3 y_3 z_3 vx_3 vy_3 vz_3 . . . etc. There is a shell script 'mergerv' which can merge the disk, bulge and halo files into a single N-body file. The program tobinary turns the ascii files into a simple binary format, listing first the number of particles, then all their masses, then the time, and finally the x,y,z,vx,vy,vz coordinates for each particle.