The bulge and halo are a King and lowered Evans distribution resp, whereas the disk a 3D generalization of Shu’s planar disk to include a 3rd integral of motion to account for the vertical structure. The construction occurs in 3 steps. First the potential is calculated, after which the disk distribution function is constructed which generates the given potential. Finally, each component is realized with a self-consistent distribution of particle orbits.

Although one has the option of including any combination of components, leaving out the halo probably won’t work.

Two new versions have been published since the current 1995 K&D version:
Widrow & Dubinski (2005) and Widrow, Pym and Dubinski (2008). These are not
(yet) supported in the current **mkkd95** program.

**out=**- output snapshot (a rundirectory $out.tmpdir is also created). No default.
**ndisk=**- Number of
particles in disk. Set them to 0 to leave out this component. Notice these
are the first set (of three) of particles in the snapshot. [8000]

**nbulge=**- Number of particles in bulge, the second set of particles in the output snapshot. [4000]
**nhalo=**- Number of particles in halo. Notice these are the thirds and last set of particles in the snapshot. [6000]
**psi0=**- in.dbh: Psi0 (HALO) central potential - the smaller (the more negative) this parameter the deeper the potential and the more extended the halo [-4.6]
**v0=**- in.dbh: v0 = sqrt(2.0)*sigma0 where sigma0 is the central velocity dispersion. roughly the velocity where the halo rotation curve peaks [1.42]
**q=**- in.dbh: q an optional flattening parameter for the potential - generally 0.7 < q < 1.05 - q=1.0 will give a nearly spherical halo [1]
**rck2=**- in.dbh: (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. [0.1]
**ra=**- in.dbh: 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. [0.8]

**md=**- in.dbh: M_d - mass of the exponential disk ignoring cutoffs [0.867]
**rd=**- in.dbh: R_d - exponential scale length [1]
**router=**- in.dbh: R_outer - outer radius where we begin to truncate the disk density [5]
**zd=**- in.dbh: z_d - disk scale height assuming a sech^2(z/zd) vertical density law [0.1]
**drtrunc=**- in.dbh: dR_trunc - truncation width - the disk density smoothly drops to zero in the range R_outer < R < ~R_outer + 2*dR_trunc. [0.5]
**rhob=**- in.dbh: rho_b - bulge central density [14.45]
**psicut=**- in.dbh:
psi_cut - bulge cut-off potential psi0 < psi_cut < 0.0 - energy cut-off
for the bulge

[-2.3] **sigb=**- in.dbh: sig_b - bulge central potential [0.714]
**dr=**- in.dbh: delta_r the width of the radial bins used to calculate the potential [0.01]
**nr=**- in.dbh: nr - number of radial bins - initially a guess since we don’t know the radial extent of the system [2400]
**lmax=**- in.dbh: 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 number of harmonics [10]

**sigvr0=**- in.diskdf: central radial velocity dispersion [0.47]
**sigr0=**- in.diskdf: scalelength of sig_r^2 [1.0]
**ncorr=**- in.diskdf: number of intervals for correction function [50]
**niter=**- in.diskdf: number of iterations [10]
**fstreamb=**- in.bulge: the streaming fraction of stars with L_z > 0. This controls the rotation of the bulge and the halo, with f=0.5 the non-rotating case. [0.75]
**fstreamh=**- in.halo: see above, for the halo. [0.5]
**iseed=**- in.- Random Seed - KD95 convention. Notice this seed is shared by all 3 components. [-1]
**zerocm=**- in.- Center the snapshot? Notice again this is shared by all three components. [t]
**bin=**- directory in which KD95 binaries live.
Unless your binaries have been added to your search PATH, you need to
specify the directory where these live. [.]

**model=A|B|C|D**- Base K&D95 model. The default model parameters above are already for model-A, but with this keyword another one of the 4 published models can be quickly loaded. After that specific parameters can still be overriden by setting their value explicitly. By default model is not set via this method.
**nmodel=**- Number of models computed. Useful if you want to stack a large number of models to increase signal to noise. Default: 1.
**cleanup=t|f**- Cleanup the temporary run directory after usage. The contents of this directory
is probably only useful for debugging, it contains the output of all the
individual programs and temporary files before final conversion to the
*out*file. Default: t

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 using the *vcirc*
program, in which the contributions to the total rotation curve are tabulated.
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.

**dbh** calculates the potential. From in.dbh is computes dbh.dat, h.dat,
b.dat and mr.dat.

**getfreqs** tabulates various characteristic frequencies (omega,
kappa etc.) in the equatorial plane for use by **diskdf**. Input files are dbh.dat
h.dat b.dat, and it generates freqdbh.dat

**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. The keywords sigrv0, sigr0, ncorr, niter
are used for this. It also outputs the Toomre Q as a function of radius
in the file toomre.dat.

**gendisk, genbulge, genhalo** assembled the respective
components using an input file.

**mergerv** is a small shell script that merges
the 3 ascii files.

% mkkd95 A0.dat % snapmstat A0.dat sort=f 0 0:7999 = 8000 Mass= 0.000108822 TotMas= 0.87058 CumMas= 0.87058 1 8000:11999 = 4000 Mass= 0.00010631 TotMas= 0.425242 CumMas= 1.29582 2 12000:17999 = 6000 Mass= 0.000819365 TotMas= 4.91619 CumMas= 6.21201 % snapplot A0.dat color=’i<8000?2.0/16.0:(i<12000?3.0/16.0:4.0/16.0)’ yvar=z % snapxyz A0.dat - color=’i<8000?1:(i<12000?2:4)’ | xyzview - maxpoint=18000 nfast=18000 scale=8 fullscreen=t % mkkd95 B0.dat model=B % snapplot B0.dat color=’i<1000?2.0/16.0:(i<2000?3.0/16.0:4.0/16.0)’ yvar=z % mkkd95 C0.dat model=C % snapplot C0.dat color=’i<4000?2.0/16.0:(i<6000?3.0/16.0:4.0/16.0)’ yvar=z % mkkd95 D0.dat model=D % snapplot D0.dat color=’i<1000?2.0/16.0:(i<2000?3.0/16.0:4.0/16.0)’ yvar=z

http://www.astro.rug.nl/~kuijken/galactics.html original GalactICS distribution 1995MNRAS.277.1341K - Kuijken, K.; Dubinski, J. Nearly Self-Consistent Disc / Bulge / Halo Models for Galaxies Widrow & Dubinski 2005 (version 2) Widrow, Pym and Dubinski 2008 (version 3) makegalaxy https://bitbucket.org/lutorm/makegalaxy

NEMO/src/nbody/init/mkkd95.c NEMO/usr/kuijken/GalactICS-exp/ $out.tmdir/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 $out.tmdir/h.dat same as above, but only for halo $out.tmdir/b.dat same as above, but only for bulge $out.tmdir/mr.dat mass and radial extent (or edge) of disk, bulge and halo

06-Mar-04V1.0 Created, using kd95’s README filePJT 11-mar-04V1.2 added nmodel= and warns about using model=PJT 23-mar-04V1.4 use logfile in tmpdir, added cleanup=, some key reorderPJT 19-jul-06V1.5 fix POSIX problems and document order of particles better 27-jul-06fixed documentation on vcirc usage, disabled pgplot for ia64PJT 5-aug-06merged two previous (CVS) docsPJT