{'fraction_non_alphanumeric': 0.03332109616256721, 'fraction_numerical': 0.018722871272161875, 'mean_word_length': 4.367798771299409, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling galactic spectra: I - A dynamical model for NGC3258

In this paper we present a method to analyse absorption line spectra of a galaxy designed to determine the stellar dynamics and the stellar populations by a direct fit to the spectra. This paper is the first one to report on the application of the method to data. The modelling results in the knowledge of distribution functions that are sums of basis functions. The practical implementation of the method is discussed and a new type of basis functions is introduced. With this method, a dynamical model for NGC 3258 is constructed. This galaxy can be successfully modelled with a potential containing 30% dark matter within 1r_e with a mass of 1.6x10^11 M_o. The total mass within 2r_e is estimated as 5x10^11 M_o, containing 63% dark matter. The model is isotropic in the centre, is radially anisotropic between 0.2 and 2 kpc (0.88 r_e) and becomes tangentially anisotropic further on. The photometry reveals the presence of a dust disk near the centre.

acknowledged that dynamical evidence for dark matter in elliptical galaxies requires the use of more information contained in the line-of-sight velocity distribution (hereafter LOSVD) than simply the mean velocity and velocity dispersion. The reason for this is the existence of a so called mass-anisotropy degeneracy. In a dynamical model it is possible to mimic the influence of a dark matter halo by tangential anisotropy. Therefore, also the anisotropy information in LOSVDs should be used. However, there is no uniformity in the parametrization of these LOSVDs and different algorithms may well have different biases . Moreover, it is not always straightforward to include these parameters in a fit, since not all of them are linearly dependent on the distribution function (hereafter DF) (Gerhard et al. (1998), De Bruyne et al. (2001). Another option is to use the complete LOSVD as dynamical input. This omits the need to approximate the LOSVD with a number of parameters, but in combination with the Schwarzschild orbit superposition method, this requires a careful smoothing both of the LOSVD and the dynamical model (Gebhardt et al. 2000). But there is still another way of dealing with the kinematic information of a galaxy, through direct modelling of the observed spectra.
In this paper we adopt a method to analyse spectra from elliptical galaxies that explores the dynamics of the galaxy and simultaneously offers a way to study the stellar populations in that galaxy. The method originates from the field of dynamical modelling and is outlined by De Rijcke & Dejonghe (1998), hereafter DD98. The main idea is that the different stellar populations that contribute to the integrated galaxy spectrum do not necessarily share the same kinematic characteristics. Hence, in terms of dynamical modelling, they should be associated with different distribution functions. A simultaneous determination of the best stellar mix and the best distribution functions associated with them de facto amounts to a dynamical model and a stellar population synthesis.
This paper is the first to present the application of the modelling method to data. The aim is to determine the gravitational potential of NGC 3258, by modelling the galaxy spectra using different stellar templates. The modelling strategy is outlined and commented, and a new family of dynamical components is presented. This is the first paper in a series of two, where a deeper interplay between dynamics and population synthesis is explored. In this first paper, the attention will be focused on the dynamical modelling aspect, while the second paper will go deeper into the stellar population analysis.
In the next section, the observations and data reduction steps are presented. The modelling strategy is outlined in section 3. Section 4 presents the dynamical components that were used. In section 5 the results of the modelling are presented. Section 6 contains a discussion of the results and the conclusions can be found in section 7.

OBSERVATIONS AND DATA REDUCTION
The E1 galaxy NGC 3258 is one of the brightest members of the Antlia Group. In its neighborhood, NGC 3257 can be found at 4.5 arcmin (about 21re) and NGC 3260 is at 2.6 arcmin (about 12re). We measured a systematic velocity of 2709 km/s, this puts NGC 3258 at a distance of 36.12 Mpc (taking H0 = 75km/s/Mpc). This means that a linear size of 1 kpc at the distance of NGC 3258 translates into an angular size of 5.7 ′′ .
NGC 3258 was observed with the ESO-NTT telescope in the nights of 27-28/2/2000. Spectra of the major axis covering the Ca II triplet around 8600Å were taken. During the observations in 2000, also imaging of NGC 3258 in the R band was performed. The detector used for this observation run was a Tektronix CCD with 2048 × 2047 pixels, 24µm×24µm in size and with a pixel scale of 0.27 ′′ /pix. Another observation run, in 1999, was used to obtain images in the B band.

Photometry
R band observations (in total 600 sec integration time) of NGC 3258 were obtained at the ESO-NTT using the red arm of EMMI during the spectroscopic run in februari 2000. An additional B band image (600 sec integration time in total) was obtained in february 1999 in the same configuration. The detector used for the B band images was a Tektronix CCD with 1024×1024 pixels, 24µm×24µm in size and with a pixel scale of 0.37 ′′ /pix.
The seeing (judged from a number of stars in the images) was 1.1 ′′ for the R band images and 1.7 ′′ for the B band images. The photometric zero point was established by comparison with published photometric data available on NED. Figure 1 shows an image of NGC 3258, isophotes are plotted in figure 2. In this figure, the central isophote seems to be a little displaced compared to the other isophotes. Both B and R images show a slightly off-centre nucleus, with a shift of about 0.45 ′′ (or 0.07 kpc) in the horizontal direction and a shift of about 0.1 ′′ (or 0.01 kpc) in the vertical direction. This is illustrated in the upper panel figure 5, where the shift of the centres of the isophotes is shown.
The photometric parameters were derived by fitting elliptic isophotes to the image. The intensity along the isophotes was expanded in a fourth order Fourier series, in order to reveal deviations   from pure elliptic isophotes, like in e.g. De . In figure 3, the main photometric parameters are shown: surface brightness, position angle, ellipticity and B-R profile as derived from the B and R surface brightness profiles in the left column. The coefficients indicating deviations from pure elliptic isophotes are shown in the right column. The values for these coefficients indicate that NGC 3258 is an almost perfect elliptical galaxy, with an ellipticity between 0.1 and 0.2.
A Sérsic profile (I(r) = I0 exp[−(r/r0) 1/n ], with I(r) the surface brightness profile, I0 the central surface brightness, r0 the scale length, an n the Sérsic parameter (Ciotti et al. 1995) was fitted to the surface brightness, the result of the fit is overplotted (in solid line) in the upper left panel of figure 3. The observed surface brightness profile can be well matched with a Sérsic profile with Sérsic parameter n = 5.47 ± 0.07. The effective radius is 12.95 ′′ or 2.26 kpc.
Values for B-R (derived from the photometric profiles) are shown in the lower left panel of figure 3. Within the central 5 ′′ , there seems to be a slight reddening of the nucleus. B-R changes from about 1 in the most central measurement to about 0.85 at 5 ′′ . Figure 4 shows the central part of the unsharp masked R image with a clearly visible dust feature. The estimated radius of the major axis of the disklike feature is about 1 ′′ -1.5 ′′ (0.17 kpc -0.26 kpc). Figure 5 shows a contour map of this region (lower panel) and a cut along a horizontal line through the centre (upper panel). This line does not coincide with any of the photometric axes. The effect of the dust lane on the apparent brightness is clearly visible in the upper panel of that figure. The brightest central point does not coincide with the centre of the isophotes for the main body of the galaxy. The same behaviour is seen in the B band images. The off-centre nucleus appearing in the B and R images is probably due to the obscuring effect of the dust lane.

Spectroscopy
For the spectra, grating #7 was used, having a dispersion of 0.66 A/pix. A slit width of 1.5 ′′ yielded a spectral resolution of 3.67Å FWHM, resulting in an instrumental dispersion of about 54 km/s in the region of the Ca II triplet. For NGC 3258, several exposures of 3600 sec were taken (in total 11 hours).
A number of standard stars (G dwarfs and K and M giants) were also observed, they are listed in table 1.
Standard reduction steps were applied to these spectra with ESO-MIDAS 1 . Dedicated exposures, taken during both runs, were used for bias correction, dark correction and flatfielding. Cosmic ray hits were removed using a top hat filter. After this correction, remaining cosmics were also removed by hand.
For the wavelength calibration, lamp spectra were taken just before or after each of the spectroscopic observations of the galaxy or the stellar template stars. For each row in the spectra a polynomial was constructed to transform the pixel scale into wavelength scale. The calibrated images with the Ca II triplet have a step of 0.3Å.
Airmass correction is applied, using the mean value of the airmass at the beginning and end of the exposure. The contribution of the sky to the spectra was estimated from an upper and lower region of the image, where there was no contribution of the light of the galaxy or template star. Several galaxy spectra were reduced separately, aligned and combined into one galaxy spectral image.

Kinematic parameters
In figure 6 kinematic parameters are presented. These parameters are derived from the spectra, using a χ 2 minimization technique in order to find the parameters of a truncated Gauss-Hermite series that gives the best approximation to the LOSVD. The spectra were spatially rebinned to satisfy a minimum S/N. Only data points from the side of the galaxy with mean streaming away from the observer are shown. A best fitting template mix (i.e. the one that gives the smallest χ 2 in the fitting) was determined for the central data point and this mix was then used throughout the whole galaxy. The spectra are analysed with a mix of a G5V dwarf and a K4III giant.
The kinematic parameters are determined following the method of van der Marel & Franx (1993) and implemented as in De , first fitting a Gaussian LOSVD to determine v and σ, and then keeping these lowest order parameters fixed to fit the other coefficients.
If v , σ, h3 and h4 are the lowest order parameters for the LOSVD, the true moments are approximated more accurately by for the first order moment and  for the second order moment, vp and σp are also called the 'corrected' values. These data are consistent with the kinematic parameters presented by Pellegrini et al. (1997) and Koprolin & Zeilinger (2000).

Other observations
NGC 3258 is a weak but extended radio source with a total flux of 18.3 mJy (VLA observations at 4.9 GHZ (6 cm) within 45 ′′ ) (Sadler et al. 1989), where the core (< 5 ′′ ) flux is only 7% of the total flux. A hot gaseous corona around NGC 3258 emits about 3 × 10 41 ergss −1 (Fabiano et al. 1992). There are no total mass estimates based on X-ray observations available. Bregman et al. (1998) studied this galaxy in the far-infrared, and derived a dust temperature of 24.5 K and a dust mass of 6 × 10 6 M⊙ from IRAS observations at 60 and 100µm. This is not unlike what is found in spiral galaxies. The mid-infrared properties of NGC 3258, based upon ISOCAM observations at 6.75, 9.63 and 15µm, indicate the presence of a hot T ≃ 260K) dust component (Ferrari et al. 2002). A mass of 48M⊙ of hot dust within a radius of 5.5 kpc was derived. This is at the lower end of the typical range of 10-400 M⊙ reported for early-type galaxies.

DYNAMICAL MODELLING STRATEGY
There is a large variety of ways to make a dynamical model and there also different strategies to use these models. For a recent review, see Dejonghe & De Bruyne (2003).
The data used for the modelling are the projected surface brightness and a set of spectra of the galaxy, taken at several positions along the major axis. The projected surface brightnes is used to obtain a graviational potential, this is described in more detail in 5.1.
An important issue for the modelling strategy applied in this paper, is the form of the data that are used to determine the DF. In this case, the fit is done directly to the spectra (after necessary data reduction steps). This means that the observed quantities that have to be reproduced by the model are pairs (wavelength, flux). The use of spectra as input for dynamical information is in contrast with more frequently used modelling techniques where a number of kinematic parameters are first derived from the spectra. The models are then fitted to these parameters that are moments of the distribution function. The modelling procedure used in this paper is described in detail in DD98. The input spectra are reproduced by a dynamical model that uses spectra of different template stars. An important characteristic of these models is that they provide a DF for each population that is used in the model. This is a function which gives the probability to find a star on an orbit characterized by the integrals of motion E (the energy), L (the total angular momentum) and Lz (the component of the angular momentum along the rotation axis of the system). This function holds all information on the dynamical state of the galaxy.
The resulting DF for each type of star is a sum of basic components and is constructed with a quadratic programming algorithm (Dejonghe 1989). There are two main requirements for the DF: (1) the difference between the observed quantities and the quantities calculated from the DF should be minimized (least squares minimization) and (2) the DF should be a positive function.
Modelling the spectra instead of kinematic parameters causes a considerable increase in the number of data points and hence an increase of the required computing time, but there are a number of advantages over using kinematic parameters.
• More traditional modelling involves a two step process (e.g. De Bruyne et al. (2001)): first determining kinematic parameters out of observed spectra and then calculating a model. Here this is turned into a one step process. Omitting the step of determining kinematic parameters and using the galactic spectra as input has the consequence that the kinematic information used as input for the modelling is independent of any parameterization. Also in the more traditional two step modelling process, there are some non-parametric methods to derive kinemtical profiles, see De Bruyne et al. (2003).
Of course, the dynamical profiles that are a result from the modelling depend on the assumptions concerning the gravitational potential and the dark matter halo, but this is the case for every dynamical modelling technique. In this paper, the gravitational potential is assumed to be spherical. This assumption is not a conceptual issue but is justified by the shape of the NGC 3258.
• Furthermore, since the basis distribution functions that are used in the fit are chosen freely by the algorithm from a large library, also the modelling procedure itself is non-parametric. The smoothness of the model is guaranteed by the use of smooth basis functions. Hence, considering also the previous remark, this method can be regarded as a non-parametric modelling method. Moreover, the model has considerably more freedom in determin-ing LOSVDs than can be achieved with a limited number of parameters as is the case when kinematic parameters are determined in a separate process.
• The modelling process based on quadratic programming, uses a χ 2 (see equation (10) in DD98) that has the statistical meaning of a goodness of fit indicator if the noise on the spectra is Poisson noise and the data points are weighted accordingly.
• The Hessian matrix (see equation (13) in DD98) involved in the quadratic programming can be used to calculate error bars on the distribution functions that have a statistical meaning and give an idea of the uncertainty on the model and on its moments. Tests with synthetic spectra (DD98) show that the error bars are a reliable measure for the spread between input model and fit, as caused by the adopted Poisson noise. Moreover, if the component library is sufficiently complete the systematic errors are negligible compared to the random errors.
• If the dynamical modelling uses different template stars, it can result in a population synthesis.
• The composition of the stellar mix is allowed to vary smoothly as a function of radius and is therefore directly connected to the internal dynamics. This minimizes stellar template mismatch, and allows in meantime to obtain a physical stellar population decomposition.
• Information at different wavelengths can be combined in an elegant way. This is an advantage of using kinematic information on the level of spectra. It is not clear how to interpret profiles of the same kinematic parameters taken from different parts of the spectrum that give different results (van der Marel 1994).

Templates
Five template stars that have discernible spectral properties were used: a G2V dwarf, a G5V dwarf, a K1III giant, a K4III giant and a M1III giant (see table 1 for their identification).
The dynamical models were created in an iterative process. First, models were calculated with a single template star, using a library of about 100 dynamical components. Afterward, the components that are selected for the single template models are combined into component libraries that contain a sample of dynamical components for different stellar templates. With these mixed libraries, new models were calculated.

DYNAMICAL COMPONENTS
Using a spherical potential in combination with rotating components for the DF, it is possible to obtain a flattened luminous density distribution, as is required to be able to fit the two dimensional surface brightness of this E1 galaxy well. We consider three integral models: F (E, L, Lz) for which σr, σ φ and σ θ are different. These models are linear combinations of three types of components.

Continuum terms
The continuum of the weighted sum of template spectra does not necessarily fit the continuum of the galaxy spectra. Therefore, low order polynomials are included in the component library to compensate for the small differences between both continua. Such a continuum component has the form with α and β real numbers, rc the cutoff radius (the component is set to zero for larger radii), jR(0) the central value of the surface brightness, jR(rc) the surface brightness at the cutoff radius, jR(rp) the surface brightness at the position of the line of sight, λ the wavelength at which the spectrum is evaluated, λ1 and λ2 bound the spectral region that is used in the fit. The value of α should be small so that these components only fit the continuum and not the absorption lines.

Fricke components
These components have the following augmented densitỹ with α and β positive integers, rc a scale length and E0 the energy level at which the component is truncated. All of these components are either isotropic (β = 0) or tangentially anisotropic (component parameter β > 0). However, the coefficients for these components can be negative, hence also radial anisotropic models can be constructed. These are spherical two integral (E, L) components.
They have equal velocity dispersions in the tangential directions, σ φ = σ θ , but a different radial dispersion in the case of nonisotropic components. A more detailed discussion of this family of components can be found in DD98.

Rotating components
These components add rotation to the model. Due to the ordered motion of the stars around the component's symmetry-axis, the tangential dispersion is usually lower than the radial dispersion. These rotating components are axisymmetric two integral (E, L z ′ ) components (with z ′ the symmetry axis of the component), for which σR = σ z ′ and the dispersion in the rotation direction σ φ is different.
In this section, a family of components is introduced that should enable one to reproduce both the rotation of the bulk of the galaxy and the kinematics of a small subsystem, rotating around an axis other than the main galaxy's rotation axis. The distribution function of the components of this family is given by with the energy per unit mass E = ψ(r) − v 2 M /2 − L 2 z ′ /(2̟ 2 ) and the z ′ -component of the angular momentum L z ′ = ̟v φ . The radial distance to the centre is measured by r, the distance to the centre measured in the equatorial plane is denoted by ̟. The velocity component in the meridional plane is vM = v 2 ̟ + v 2 z ′ , the velocity component perpendicular to the meridional plane is given by v φ . Here, (x ′ , y ′ , z ′ ) is a Cartesian reference frame, attached to the component, with the z ′ -axis the rotation axis of the component. The overall gravitational potential ψ(r) is taken to be spherical symmetric although the components themselves have an axisymmetric mass distribution. This is acceptable since a small kinematically distinct core probably does not significantly influence the global potential. Also, a slowly rotating elliptical galaxy can still be approximately round, justifying the use of a spherical potential. An axisymmetric potential generally allows E and Lz as integrals of motion, with Lz the component of the angular momentum parallel to the potentials symmetry axis. Therefore, these components can also be used in an axisymmetric potential if their rotation axis z ′ coincides with the symmetry axis z of the potential.
All the stars have the same sense of rotation, i.e. if irot = 1 the stars rotate counter-clockwise, if irot = −1 they rotate clockwise. The parameter α determines the distribution of the stars over orbits with different energies whereas β influences the way orbits with different angular momenta are populated. The energy cutoff E > E0 will prevent stars from reaching beyond a certain maximum radius. The angular momentum cutoff |L z ′ | > L0 on the other hand will impose a minimum radius that stars can reach.

The component's frame of reference
These components are to be employed in the modelling of slowly rotating, approximately spherical galaxies using a severally symmetric potential. Throughout this work, all dynamical quantities are expressed in spherical coordinates. These particular components however are axisymmetric and their characteristics will be most readily calculated in cylindrical coordinates. Therefore, we seek a transformation linking the spherical coordinates in which the properties of the galaxy as a whole are expressed to a cylindrical coordinate system, attached to each individual component of this family.
First, we introduce a cartesian reference frame (x, y, z) with its origin at the galaxy's centre. The components of a vector a in spherical coordinates (r, θ, ϕ) can be linked to its components in cartesian coordinates via the transformation with A = sin θ cos ϕ sin θ sin ϕ cos θ cos θ cos ϕ cos θ sin ϕ − sin θ − sin ϕ cos ϕ 0 Suppose the direction of the rotation axis of the system we are describing, whether it be a small subsystem or the galaxy's main body, is given by the angles (θ0, ϕ0) with respect to the cartesian reference frame. The cartesian coordinate system (x ′ , y ′ , z ′ ), fixed to the component, can be obtained by first rotating the (x, y, z)frame over an angle ϕ0 around the z-axis, yielding the auxiliary (x ′′ , y ′′ , z ′′ )-frame, followed by a rotation of the (x ′′ , y ′′ , z ′′ )frame over an angle θ0 around the y ′′ -axis : with B = cos θ0 cos ϕ0 − sin ϕ0 sin θ0 cos ϕ0 cos θ0 sin ϕ0 cos ϕ0 sin θ0 sin ϕ0 − sin θ0 0 cos θ0 .
To facilitate the calculations as much as possible, we will henceforth use cylindrical coordinates (̟, φ, z ′ ). These are connected to the (x ′ , y ′ , z ′ )-frame by the transformation One can now immediately make the connection between the spherical coordinates (r, θ, ϕ) in the galaxy's reference frame and the cylindrical coordinates (̟, φ, z ′ ) in the component's frame of reference using the matrix D = ABC :

The velocity moments
The axisymmetric velocity moments for a counter-clockwise rotating model (L z ′ 0, irot = 1) are defined as The moments of a clockwise rotating model (L z ′ 0, irot = −1) with the same parameters α and β are readily obtained sincẽ The distribution function is an even function of v̟ and v z ′ , therefore only the moments with even l and n values will be non-zero. One can define anisotropic moments that are connected to the true moments by the relatioñ µ + 2l,m,2n (ψ, ̟) = 1 π B(l + 1 2 , n + 1 2 )μ + 2(l+n),m (ψ, ̟).

Spatial density
The spatial density is given by The density is only non-zero inside the volume defined by ψ(r) − E0 > L 2 0 /(2̟ 2 ). Evidently, the boundary of this volume, where the density becomes exactly zero, is determined by For non-zero L0 and E0, this volume has approximately the shape of a torus. If L0 = 0, the density is non-zero inside the spherical region defined by r ψ −1 (E0). All models that share the same values for E0 and L0 fill the same volume, independent of α or β but the density distribution inside this volume does depend on α and β, see figure 7

Mean velocity For the mean rotation velocity, one finds
At the boundary of the component, the mean rotation velocity takes the value v φ = L0/̟ = 2(ψ(r) − E0). The parameter β has little impact on the rotation velocity profile, as shown in figure 8 Figure 8. The influence of the angular momentum cut-off on the kinematics of models with α = 6, β = 2, rmax = 6 kpc is demonstrated. The values of the inner radius are r min = 0 kpc (solid line), 0.5 kpc (dotted line), 1.0 kpc (dashed line) and 1.5 kpc (dashed-dotted line). In the left column, a number of intrinsic kinematic quantities are plotted. The radius R is measured in the equatorial plane, z ′ = 0. From top to bottom : the logarithm of the spatial luminosity density ρ(R), the radial (black) and the tangential velocity dispersion (gray), the rotation velocity v φ (R) and Binney's anisotropy parameter β(R). In the right column, the projected kinematics are presented. The radius rp lies along the long axis of the models. From top to bottom : the logarithm of the projected luminosity density ρp(rp), the projected velocity dispersion σp(rp), the projected streaming velocity vp(rp) and the Gauss-Hermite coefficients h 3 and h 4 (gray). and figure 9. The extreme values of the rotation velocity are determined solely by L0 and E0, making the rotation velocity a very sensitive function of these parameters.

Velocity dispersion
Since l and n are completely equivalent in equation (19), it is obvious that the velocity dispersions σ 2 ̟ (̟, z ′ ) and σ 2 z ′ (̟, z ′ ) are always equal. This is the case for any model with a distribution function that depends only on E and L z ′ . One finds that Figure 9. The influence of the parameter β on the kinematics of models with α = 6, r min = 0.5 kpc and rmax = 6 kpc is demonstrated. The values of β are 2 (full line), 4 (dotted line), 6 (dashed line) and 8 (dasheddotted line). In the left column, a number of intrinsic kinematic quantities are plotted. The radius R is measured in the equatorial plane, z ′ = 0. From top to bottom : the logarithm of the spatial density ρ(R), the radial (black) and the tangential velocity dispersion (gray), the rotation velocity v φ (R) and Binney's anisotropy parameter β(R). In the right column, the projected kinematics are presented. The radius rp lies along the long axis of the models. From top to bottom : the logarithm of the projected density ρp(rp), the projected velocity dispersion σp(rp), the projected streaming velocity vp(rp) and the Gauss-Hermite coefficients h 3 (black) and h 4 (gray).
Figures 8 and 9 also show that the anisotropy parameter β φ , with is almost independent of the components parameter β but very sensitive to the value of the angular momentum cutoff L0.

The spatial line-of-sight velocity distribution
For the calculation of the spatial LOSVD, the DF will be directly integrated over the two velocity components in the plane of the sky. The integration over the line of sight will have to be performed numerically since it depends on the shape of the potential, which is generally only known numerically. We decompose the velocity vector in a component parallel to the line of sight, vp, and two components in the plane of the sky, v1 and v2. The polar unit vector e φ can be written with respect to the vectors ( ez, e1, e2) with eg the unit vector along the line of sight and e1 and e2 unit vectors in the plane of the sky : with b 2 1 + b 2 2 + b 2 3 = 1. We now rewrite the distribution function as a function of vp, v1 and v2 : Here and in the following, the power on the angular momentum will be assumed to be an integer, denoted by n. LOSVDs of components with non-integer values for this parameter will have to be calculated numerically. The velocity components in the plane of the sky, v1 and v2, can be expressed in polar coordinates (v, ϑ) Taking this transformation into account and introducing the notation following expression can be obtained in which χ has been dropped, since the sign and the zero-point of ϑ are immaterial. The cosines b and b1 obey the relation b 2 + b 2 1 = 1, with b always positive. The distribution function has to be integrated over that part of the rectangle To relieve the notation throughout the remaining calculations, the following parameters will be used The calculations will proceed along different paths depending on whether P 0 or P 0. The area in the rectangle (31) over which the distribution function has to be integrated is further limited by the positivity constraint on the angular momentum, (32). This translates in the extra constraints The spatial LOSVD can then be calculated as with the integration interval defined by (31) and (34).

P
0 In this case, the distribution function has to be integrated over the grey shaded area in figure 10. Only if P < Ab will the second of the curves (34) intersect the rectangle (31). P Ab 0 The z ′ -component of the angular momentum is positive everywhere inside the rectangle (31). The integration over ϑ can be performed over the entire interval [0, π], and so π 0 (P + bv cos ϑ) n dϑ = (36) where [n/2] is the largest integer, less or equal to n/2. Now, the integration over v can be tackled. This integral reduces to a Betafunction and the spatial LOSVD can be written in the following elegant form : which is finite and well-behaved for all integer n and real α.

P Ab
In this case, cos ϑ0 0 and ϑ0 ∈ [ π 2 , π]. The second of the curves defined by (34) cuts out a piece of the rectangle (31) and hence the distribution function has to be integrated over the shaded area in figure 10 : The first auxiliary function h1 is simply a Beta-function The second function h2 is an incomplete Beta-function and can be rewritten in terms of a hypergeometric function : (41) Thus, we obtain the following expression for the spatial LOSVD φ(̟, z ′ , vp) = Γ(α + 1) 2 α+n/2 ̟ n A 2(α+1) P n n i=0 n i We first turn our attention to the second term of (42). Here, only negative odd powers of cos ϑ appear in the integral. The integration can be done analytically with the aid of Gradshteyn & Ryzhik (1965) and yields the result : As to the first term of (42), the even and odd powers of cos ϑ have to treated separately n i=0 n i Thus, we are lead to the following expression for the spatial LOSVD: In the limit P → Ab, one finds that ϑ0 → π and the complicated expression (45) reduces to the far easier formula (38), as it should. The rather daunting formula (45) will only be well-behaved for integer values of α since then the last summation (with summation variable k) terminates after a finite number of terms.
4.3.4.2 P < 0 In this case, cos ϑ0 = − P Ab > 0 and hence ϑ0 ∈ [0, π 2 ]. The distribution function now has to be integrated over the shaded area in the (v, ϑ)-plane presented in figure 11. Only if |P| < Ab will the first of the curves (34) intersect the rectangle (31). |P| Ab The integration area defined by the equations (34) and (31) vanishes, hence the spatial LOSVD is identically zero for all |P| Ab.
|P| < Ab The distribution function must be integrated over the grey shaded area in figure 11 : We start with the integration over the velocity in the plane of the sky, v, which can be performed analogously to (40) and (41) : Thus, we obtain an expression that is almost exactly like (42) φ(̟, z ′ , vp) = Γ(α + 1) 2 α+n/2 ̟ n A 2(α+1) P n × n i=0 n i The integrations over ϑ are handled analogously to (43) and (44) and we finally obtain This formula will only be used for integer values of α, as discussed in the case of (45). In the limit |P| → Ab, ϑ0 vanishes. The spatial LOSVD becomes identically zero and remains zero for all |P| Ab.

DYNAMICAL MODEL
The two most prominent absorption lines of the Ca II triplet were used as input for the modelling. The spectra were rebinned to reach sufficiently high S/N. The S/N was about 50 at radii less than 10 ′′ and from these spectra 182 points were included in the fit. For larger radii, the S/N was about 30 and only 91 points per spectrum were included in the fit.

Potential
The observed R-band surface brightness profile was deprojected into a spherical luminosity density. Assuming a mass-to-light ratio, this lead to a mass density that was used to calculate a potential by means of the Poisson equation. A model with constant M/L did not result in a satisfactory fit. Hence, a number of models with a dark matter halo were considered. For all these models, the dark halo component was added in the form of a spherical distribution of mass, where the amount of matter is increasing as a power of the radius. Different models have different fractions of total matter mass to light emitting mass at a chosen radius (here 2 kpc or 0.9re, this radius is only for scaling and has no physical meaning) and different total masses. In this way, a grid of models was created, and the smallest value for χ 2 was taken to be indicative for the best matching potential. The best results were obtained with a potential where the spatial density of the total amount of matter was 2.5 times the spatial density of the amount of luminous matter within 2 kpc. This scaling implies that the total mass within 1re is 1.6 × 10 11 M⊙, 70 % of which is luminous matter. The total mass at 2re is 5 × 10 11 M⊙, and 37 % of this is luminous matter (see also upper panel in figure  12). Figure 12 (lower panel) shows the total mass profile in solid line in comparison with the luminous matter mass profile in dashed line.
It was not necessary to include a central black hole in the model to obtain a good fit for the central spectra.

Dynamical model
The model that best fits the spectra around the two most prominent Ca II triplet lines, at 8542Å and 8662Å, is shown in figure  13. From these plots, no obvious shortcoming of the model can be indicated. Figure 14 shows some derived quantities of the model. The upper left panel shows the projected density on major and minor axis. As can be seen, the flattening of the galaxy is well reproduced. The right panels show calculated projected mean velocity and velocity dispersion profiles (solid lines), together with uncertainties on these quantities as derived from the model. The values and error estimates for the projected mean velocity vp and velocity dispersion σp, presented in 2.4, are plotted on top in dots. In figure 14 also the anisotropy parameters β φ = 1 − (σ 2 φ /σ 2 r ) (in solid line) and β θ = 1 − (σ 2 θ /σ 2 r ) (in dotted line) are displayed in the lower left panel. The parameter β φ , Binney's anisotropy parameter, is close to zero for the central parts, indicating an isotropic DF  there, becomes slightly posivite for intermediate radii and slightly negative for radii larger than 2 kpc, indicating that the model is tangentially anisotropic there.
A mix of G5V and K4III stars gives the best fit to the spectra in this wavelength range. The value for χ 2 obtained for the fit is 2267, the second best model has a χ 2 of 2383. The resulting template mix of this best model can be found figure 15. The relative densities are calculated from the spatial number densities of the model.
More models, with other template mixes, than shown here were calculated. The results mainly confirm the trends seen in figure 15. In the model with K4III and G2V stars e.g., the contribution of G2V stars to the density is larger than the G5V stars in the K4III-G5V panel, but the same radial behaviour occurs.

DISCUSSION
If a model is strongly tangential in the outer regions, most of the motion occurs along the line of sight and this may increase the projected velocity dispersion. Hence, the influence of dark matter may be underestimated. This is the strongest argument for including the fourth order shape parameter of the observed LOSVD as input parameter in a dynamical model, because it is believed to indicate whether a lot of orbits are observed near their tangent point or not. For the calculation of this dynamical model, we did not rely on the extracted LOSVD nor the derived kinematic parameters, but on the observed spectra.
The dynamical model is slightly tangentially anisotropic from roughly 1re on, as can be seen in figure 14.
Tangential anisotropy in the outer regions of the model is also seen in other models for e.g. NGC 3377 (Gebhardt et al. 2000),  NGC 1023  and also in the sample modeled by Kronawitter et al. (2000).
Estimating kinematic parameters through a dynamical model clearly reduces the amount of scatter in the kinematic profile, as is illustrated in figure 14 (right panels). This is easily understood since the model induces additional dynamical constraints on the parameters and as a result the values of the parameters taken at different radii are not independent. When the parameters are inferred directly from the observed spectra, there is no connection between the parameters at different radii. Figure 16 presents LOSVDs as derived from the observations and the model. The profiles calculated from the kinematic parameters presented in 2.4 are shown in dashed lines, the solid lines present the LOSVDs with error bars calculated from the model. It is clear that the LOSVDs calculated from the kinematic parameters all have an approximately Gaussian shape. This is of course due to the parametrization that is used to present them. The LOSVDs calculated from the model show a large variety of shapes. Some of these profiles cannot be reproduced by a Gauss-Hermite series truncated at fourth order (or even sixth order) and hence a perfect agreement with the LOSVDs calculated from the kinematic parameters cannot be expected.
The shapes of the LOSVDs calculated from the model can be better understood when looking at figure 17. In this figure the LOSVDs for the different stellar populations, together with the composed LOSVD (solid lines) are shown. The relative height of the LOSVDs for the two populations is in agreement with the information in figure 15. A full discussion of the role of the stellar populations will be deferred to a subsequent paper, where a fit will be presented that takes into account 2 wavelength regions.
Instead of comparing full profiles, one could argue in favor of comparing their most relevant parameters. Unfortunately, this is not a trivial exercise. The mathematical elegane of the Hermite expansion lies in the fact that it forms an orthogonal series, with coefficients that can be calculated by integrals that follow from the theory. However, this is almost never the way the Hermite coefficients are calculated in practice: they are generally determined with some kind of minimization after performing a convolution. This causes the loss of the orthogonality property. This becomes clear in an exercise where a stellar template spectrum convolved with a LOSVD with parameters −20km/s < v < 20km/s, σ = 300km/s , h3 = 0, h4 = 0, h5 = 0 and h6 = 0.1 is analysed in the same way as real galactic spectra and kinematic parameters up to h4 are determined. Since the template spectrum used to create this synthetic spectrum is equal to the one used to determine the parameters again, there is no influence of continuum subtraction or template mismatch on the results. The synthetic galaxy spectrum has a S/N of 130 in the centre and the light profile follows a de Vaucouleurs that resembles the light profile of NGC 3258. The result of this simulation is shown in figure 18. The values retrieved for v and h3 correspond to the original values. On the other hand, the values for σ and h4 are clearly higher than the values of the input LOSVDs.
This indicates that the kinematic parameters as derived from a Gauss-Hermite parametrization are in practice not independent, in the sense that power in the higher orders than considered in the expansion may contaminate lower order parameters. In the case of h4 this can give rise to wrong conclusions concerning the anisotropy in the stellar motion. Naively, one could conclude that it is better to go to higher order than h4 in the parametrization. Unfortunately, the signal-to-noise ratio in outer regions of a galaxy is hardly sufficient to derive reliable values for h4, so it does not make sense to try to determine higher order parameters. The useful information in LOSVDs derived from observations is thus very limited, due to the rather low signal-to-noise of the observations. Therefore, though it is a good idea to use a parametrization with a limited number of parameters to represent these LOSVD, there is, as a consequence, also only a limited number of shapes for LOSVDs possible, as ex- tracted from these observations. This on top of the fact that these parameters are not independent, which means that caution should be taken when interpreting them. On the other hand, LOSVDs calculated from dynamical models can have a much larger variety of shapes. It would be unfair to reject model LOSVDs because they do not match the ones that are derived from observations, simply because the information in LOSVDs derived from spectra is very limited, while the dynamical models offer much more freedom.

CONCLUSIONS
In this paper we use a method to analyse spectra from elliptical galaxies in order to retrieve dynamical information and to perform a stellar population synthesis at the same time. The method originates from the field of dynamical modelling and is presented by DD98. The main idea is that the different stellar populations that contribute to the integrated galaxy spectrum do not necessarily share the same kinematic characteristics. Hence, in terms of dynamical modelling, they should be given different distribution functions. This paper reports for the first time on the application of our modelling technique to data. The attention is focused on the presentation of the method and on establishing a mass estimate for the galaxy under study. A second paper will discuss the possibilities for stellar population studies and will go deeper into the distribution functions for separate stellar templates.
Modelling the spectra directly yields an estimate for the total amount of mass, and an analysis of the internal dynamics of the galaxy. This technique, using galaxy spectra as input data, also differs from more traditional modelling strategies that rely on kinematic parameters that are derived from the spectra in the sense that a two step process is turned into a one step process, thereby coming closer to the astrophysical issues.
Photometric and spectroscopic observations for NGC 3258 are presented. This galaxy has almost perfect elliptical isophotes, and has an effective radius of 12.95 ′′ (2.26 kpc). Images in B and R show a slightly off-centre nucleus. There seems to be a slight red-dening of the nucleus. The photometry also reveals the presence of a central dust disk (0.17 kpc -0.26 kpc). The R band photometry was used to derive a spherical potential for the luminous matter.
The photometry, together with galaxy spectra containing dynamical information up to 1.3re are combined into a dynamical model. The best dynamical model for NGC 3258 is the one with a dark matter halo that contributes 30% of the total mass at 1re and 63% of the total mass at 2re. From 1.45re outwards, the model is dominated by dark matter. The total mass within 1re is 1.6 × 10 11 M⊙ and (M/L)R ≈ 65, (M/L)R ≈ 30, whereas the total mass extrapolating to 2re is 5 × 10 11 M⊙ and (M/L)R ≈ 95,(M/L)B ≈ 43.5. In the region without dark matter (M/L)R ≈ 14 ((M/L)B ≈ 6.4). The modelling did not require to include a central black hole.
The anisotropy parameter (see figure 14) indicates that the model is isotropic in the centre, is radially anisotropic between 0.2 and 2 kpc (i.e. within 1re) and becomes slightly tangentially anisotropic further out.
A simulation has shown that for the analysis of LOSVDs derived from observations by means of a truncated Gauss-Hermite series, power in higher orders than considered in the expansion may contaminate lower order parameters. This is a worrysome result, since h4 profiles derived from observations tend to be interpreted in terms of anisotropy of the galaxy. If this parameter is contaminated by power in higher orders by some form of aliasing induced by the fitting process, this may lead to an erroneous interpretation or at least may make comparison of h4's obtained by different methods problematic. Moreover, due to the freedom in the dynamical model, one cannot expect that LOSVDs calculated from a dynamical model as presented in this paper can be approximated by Gauss-Hermite series that are truncated at a low order.

ACKNOWLEDGMENTS
VDB acknowledges financial support from FWO-Vlaanderen. WWZ acknowledges the support of the Austrian Science Fund (project P14783) and the support of the Bundesministerium für Bildung, Wissenschaft und Kultur.