Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://www.stsci.edu/~marel/software/scalefree.html
Дата изменения: Wed Dec 9 07:56:04 1998 Дата индексирования: Tue Oct 2 12:58:09 2012 Кодировка: Поисковые слова: п п п п п п |
*********************************************************************** * Scale-free Modeling Software * * to predict the intrinsic and projected kinematical quantities of * * axisymmetric mass densities with anisotropic velocity distributions * * in spherical potentials. * *********************************************************************** *************************** PART I: GENERAL INFORMATION *************************** Authors ------- * Roeland P. van der Marel, 1994-1995 : development of code address : Space Telescope Science Institute Research Programs Office (RPO) 3700 San Martin Drive Baltimore, MD 21218 Tel : (+1) 410 338 4931 Fax : (+1) 410 338 4596 e-mail : marel@stsci.edu homepage : http://sol.stsci.edu/~marel/ * Jos H. J. de Bruijne 1994-1995 : testing and application of code address : Sterrewacht Leiden Postbus 9513 2300 RA Leiden The Netherlands Tel : (+31) 71 5275878 Fax : (+31) 71 5275819 e-mail : debruyne@strw.LeidenUniv.nl homepage : http://www.strw.leidenuniv.nl/~debruyne/ Summary of the method --------------------- We consider the case of a scale-free spheroidal mass density with an anisotropic velocity distribution in scale-free spherical potential. The assumption of a spherical potential has the advantage that all integrals of motion are known explicitly. Two families of phase-space distribution functions are considered. The `case I' distribution functions are anisotropic generalizations of the flattened f(E,L_z) model, which they include as a special case. The `case II' distribution functions generate flattened constant-anisotropy models. Free parameters control the radial power-law slopes of the mass density and potential, the flattening of the mass distribution, and the velocity dispersion anisotropy. The models can describe the outer parts of galaxies and the density cusp structure near a central black hole, but also provide general insight into the dynamical properties of flattened systems. Because of their simplicity they provide a useful complementary approach to the construction of flattened self-consistent three-integral models for elliptical galaxies. References ---------- The method and some applications are discussed in the following paper: `Scale-free dynamical models for galaxies: flattened densities in spherical potentials' de Bruijne J., van der Marel R.P., de Zeeuw P.T. MNRAS, 282, 909-925, 1996 which can be retrieved form R. van der Marel's homepage. Acknowledgments --------------- If you have found this software useful for your research, we would appreciate an acknowledgment to `use of the Scale-free Modeling Software developed by R.P. van der Marel and J.H.J. de Bruijne'. Bug Reports ----------- Please send bug reports and important comments or questions to marel@stsci.edu Caveats ------- This software comes without guarantee and on an `as is' basis. It is not being actively maintained or upgraded. All original testing of the code was done on a UNIX Sparc 2 workstation with the SunOS operating system in the period 1994-1995. In 1997 I ensured that the code would run on a Sun Ultra 170 workstation with the solaris operating system. However, little detailed testing of the code was performed on this operating system. It is not guaranteed that the software will work without problems on other platforms, with different operating systems and different compilers. However, the code is close to standard fortran, so if any, only minor revisions will be necessary. I will consider all requests for use of this software. However, I do not allow people who have received the software from me to distribute it to third parties. Anyone who uses this software is encouraged to send me an email with their address, so that I can send reports of bugs and revisions. Program structure ----------------- There are only a few files: README : This file. Contains instructions and a description of the software. scalefree.f : The fortran program that does all the calculations. the file also contains all necessary subroutines, some of which are from the Numerical Recipes book by Press et all. makefile : makefile that allows compilation of the program scalefree.f. examp/ : directory with example input files for the program scalefree.f. After installation (see below) there will be two additional files: scalefree.o : the object file generated from scalefree.f scalefree.e : the executable file generated from scalefree.f Installation ------------ To install the software issue the UNIX command prompt> source install.com This will construct the executable scalefree.e and the object file scalefree.o. Examples -------- The directory examp contains example input files. To run the program with the example input files, issue the UNIX commands: prompt> example.e < examp/scalefree.in1.eg prompt> example.e < examp/scalefree.in2.eg Known Bugs ---------- As discussed below, the calculations of the Gauss-Hermite moments of the projected velocity profiles (VPs) do not work well for certain parameter ranges. The intrinsic and projected velocity moments should always be highly accurate. ********************************* PART II: Discussion of the method ********************************* Construct galaxy models for flattened systems of test-particles in spherical potentials, as discussed in De Bruijne et al (1996). Let b be a reference length, rho_0 a reference mass density, and W a reference velocity. Introduce dimensionless units: r = r_true / b R = R_true / b z = z_true / b v = v_true / (SQRT(2)*W) L = L_true / (SQRT(2)*W*b) E = E_true / (W^2) Psi = Psi_true / (W^2) rho = rho_true / rho_0 f = f_true / [rho_0 (2 W^2)^{-3/2}] We consider two cases for the potential: KEPLERIAN POTENTIAL: Let M be the total mass of the galaxy, such that Psi_true = GM/r_true. Then choose W = SQRT(GM/b). In these units: Psi = 1 / r LOGARITHMIC POTENTIAL: Let Vc be the circular velocity, such that Psi_true = - Vc^2 ln (r_true/b). Then choose W = Vc. In these units: Psi = - ln r Let the mass density fall off with logarithmic slope gamma, and be stratified on spheroids with axial ratio q: rho = r^{-gamma} (sin(theta)^2 + [cos(theta)^2/q^2])^{-gamma/2} where theta is the polar angle such that: R = r sin(theta) and z = r cos(theta) Since the potential is spherical, the quantities E, L^2, L_z^2 are integrals of motion. Let L_max(E) be the maximum angular momentum that can be attained by a star at energy E = Psi - v^2. In the Keplerian potential: L_max(E)^2 = 1 / 4E In the Logarithmic potential: L_max(E)^2 = EXP(-2E -1) / 2 We consider even DFs that are separable functions or quasi-separable functions of E, zeta^2, eta^2: case I : f_e = g(E) zeta^{-2 beta} j(ecc^2 eta^2) case II : f_e = g(E) zeta^{-2 beta} h(ecc^2 eta^2 / zeta^2) where: zeta^2 = L^2/L_max(E)^2 eta^2 = L_z^2/L_max(E)^2 ecc^2 = 1 - q^2 and beta is a free parameter. With these ansatz's, the functions j and h are determined uniquely (see De Bruijne et al. 1996). To allow the construction of rotating models we adopt an odd part of the form: f_o = f_e * (2s-1) * (eta^2)^t The parameters s and t determine the amount of streaming in the model. The maximum streaming model has t=0 and s = 0 or 1. The total DF is f = f_e + f_o. All the intrinsic velocity moments for the models can be expressed as a power series in [ecc^2 SIN^2(theta)]. The program scalefree.f evaluates the velocity moments. These are then used to calculate the projected line-of-sight velocity moments on the sky as one-dimensional integrals. This can be done with essentially any given accuracy, either by Gauss-Legendre quadrature or Romberg integration (to be specified by the user). These parts of the program should work flawlessly for all allowed input parameter values. In addition to the calculation of the intrinsic and projected velocity moments, the program scalefree.f also reconstructs the projected line-of-sight VP shapes. This however is a rather more difficult problem than the calculation of the projected velocity moments themselves. The approach that we have used is to characterize each projected VP as a set of N values on an array. These values are then related to the (known) first N projected velocity moments through the solution of a so-called VanderMonde matrix. Once a VP is known on an array of velocity values, characteristic quantities such as the Gauss-Hermite moments (van der Marel & Franx 1993) can be evaluated in straightforward fashion. The difficulty lies in finding the solution of the VanderMonde matrix equation, which is an ill-conditioned problem. We have implemented two different approaches: (1) The first approach is to solve the VanderMonde matrix equation directly, without any regularization or other tampering with the equation itself. This yields solutions that are oscillatory on the scale of the velocity array spacing, due to the ill-conditioned nature of the problem (noise amplification). Although the VP itself is ridiculous, it turns out that quantities like the Gauss-Hermite moments (which are integrals over the VP) are often well and correctly determined. The question remains how many (N) projected velocity moments should be used to get the most accurate Gauss-Hermite moment determinations. The program chooses N by starting with a low value of N, and then solving the VanderMonde matrix equation repeatedly until the Gauss-Hermite moments have converged to within a preset tolerance. The reason that this approach works is because a Gauss-Hermite series bears a strong resemblance to a Fourier series. Gauss-Hermite moments h_i of higher order measure power in the VP on higher frequencies (Gerhard 1993). For high N, numerical errors cause the VanderMonde solution to become oscillatory on the scale of the velocity array spacing. This, however, does not influence the observables (gamma,V,sigma,h3,h4,....), because these only measure power on relatively low frequencies. (2) The second approach is to add a regularization term to the VanderMonde matrix equation, and then solve it in a least-squares sense (cf. Section 18.5 of Numerical Recipes). This works well, but leaves open the question how the regularization parameter must be chosen. The program gives two options: (a) Let the user give the regularization parameter by hand; (b) First start with the case of essentially no regularization. This yields an oscillatory solution with many local peaks. Then increase the amount of regularization in steps until the solution has become such that there are no more than 3 significant local maxima (a local maximum is considered significant if it exceeds the value of its neighbors on the grid by a preset tolerance eps times the absolute VP maximum). The accuracy of algorithm (1) was tested as described in Appendix A of de Bruijne et al (1996). To recap, the results are as follows. A program was written that calculates the VP for the spherical (q=1) case of the models by direct evaluation of the three-dimensional integral along the line of sight that defines the VP. The (gamma,V,sigma,h3,h4,....) calculated with the program were found to be accurate, except for the case of a logarithmic potential and large anisotropy (beta < -4.0 or beta > 0.8). This has two reasons: (i) the logarithmic potential has no escape velocity so that it is more difficult to represent the VP on a finite velocity array; and (ii) the VP becomes discontinuous in the limit beta --> -infinity (only circular orbits), and singular in the limit beta --> 1 (only radial orbits) (e.g., van der Marel & Franx 1993). A second test was provided by the f(E,L_z) case (beta=0), for which the three-dimensional VP integral could be calculated as in Qian et al. (1995). Again, the (gamma,V,sigma,h3,h4,....) calculated with the VanderMonde algorithm were generally found to be accurate, with the exception of rather flattened models q < 0.6 in a logarithmic potential. This is understood from the fact that the VPs of these models become contrived and double-peaked for large flattening (Dehnen & Gerhard 1994). In both test cases inaccurate results were accompanied by non-zero values of h_1 and h_2, which should be zero by definition. For the general case one may therefore take the values of h_1 and h_2 as an indicator of the numerical accuracy of the VanderMonde algorithm. From this it was found that the algorithm (1) fails only for the case of a logarithmic potential, large flattening and large anisotropy. If knowledge of the actual shape of the VP is required (rather than just a measurement of Gauss-Hermite moments), then algorithm (2) is the only possible solution. The algorithm (2) was implemented into the software after the de Bruijne et al. paper was written, and no detailed study of its accuracy was made. It is likely that the algorithm can be optimized further, in particular the choice of the optimal regularization parameter. Users may want work on this themselves. It appears useful for many application to compare the results of approach (1) to those of approach (2), as a test of the accuracy of the results. Note that algorithm (2) is likely to be inaccurate for the same set of parameter ranges for which algorithm (1) is inaccurate (see above). References ---------- de Bruijne J., van der Marel R.P., de Zeeuw P.T., 1996, MNRAS, 282, 909 Dehnen W., Gerhard O. E., 1994, MNRAS, 268, 1019 Gerhard O. E., 1993, MNRAS, 265, 213 Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical Recipes, Second Edition. Cambridge University Press, Cambridge Qian E. E., de Zeeuw P. T., van der Marel R. P., Hunter C., 1995, MNRAS, 274, 602 van der Marel R. P., Franx M., 1993, ApJ, 407, 525 *********************************************** PART III: Description of program in- and output *********************************************** The program works in straightforward fashion. Several questions are posed on the command line, and the user merely needs to let the program know what he wishes it to do. The command line information, combined with the discussion in this file, should be sufficient to figure out what needs to be done. After the user has entered the required model parameters, and has answered to the program how it must deal with the various numerical details, the user can let the program know what information he/she wants. If the user requests information on intrinsic properties, the program will return: rho v_ph v_r^2 v_th^2 v_ph^2 at a given angle theta in the meridional plane. If the user requests information on projected properties, the program will return for a given position angle from the major axis: rho_p v_p v^2_p v^3_p v^4_p as well as the VP characteristics: true (gam,V,sig) Gauss-fit (gam,V,sig) h0, h1, h2, h3, ... and the actual VanderMonmde matrix VP solution VP(v) (in a table). If additional properties are required (such asor h_10), the fortran code can be easily modified to return those. Users are advised to run the program with the input files in examp/ to get a better idea of the workings of the program. The model in these input files is a spherical model with density proportional to r^{-4} in a Kepler potential. The results for this case can be compared to those in Table 1 of van der Marel & Franx (1993), who calculated the VP analytically as an integral over the line of sight. The input files in examp/ are for this model, but for the two different approaches of solving the vanderMonde matrix equation to get the VP. The Gauss-Hermite moments obtained with both approaches agree well, and with the results in van der Marel & Franx. Roeland van der Marel.