HfS, Hyperfine Structure Fitting Tool

HfS is a tool to fit the hyperfine structure of spectral lines, with multiple velocity components. The HfS_nh3 procedures included in HfS fit simultaneously the hyperfine structure of the NH$_3$ (J,K)= (1,1) and (2,2) transitions, and perform a standard analysis to derive $T_\mathrm{ex}$, NH$_3$ column density, $T_\mathrm{rot}$, and $T_\mathrm{k}$. HfS uses a Monte Carlo approach for fitting the line parameters. Especial attention is paid to the derivation of the parameter uncertainties. HfS includes procedures that make use of parallel computing for fitting spectra from a data cube.


INTRODUCTION
HfS (Hyperfine Structure) (Estalella 2016) is a tool to fit the hyperfine structure of spectral lines, with the possibility of fitting simultaneously multiple velocity components. The HfS nh3 procedures included in HfS fit simultaneously the hyperfine quadrupole and magnetic structure of the NH 3 (J, K) = (1, 1) and (2,2) inversion transitions. The assumptions made by HfS are that the beam filling factor f , the excitation temperature T ex , the hyperfine lines linewidth ∆V, and the central velocity V LSR are the same for all the hyperfine lines.
For HfS nh3 these assumptions hold for the hyperfine lines of both NH 3 inversion transitions, (1, 1) and (2,2). In addition, the results of the fit are used by HfS nh3 to derive physical parameters including the excitation temperature, NH 3 column density, rotational and kinetic temperature, with the assumption that the emitting region is homogeneous along the line of sight.
HfS is written in Fortran 77 and 90/95 (mainly for the dynamical storage of arrays). The graphic interface uses the PGplot Graphics Subroutine Library 2 , and some procedures have a multiprocessor version running under Open MPI 3 . The interactive procedures are menu driven, allowing to select options with the keyboard, and for some actions, with the mouse cursor and buttons.
One of the advantages of HfS when compared with other packages performing hyperfine fitting, like for instance, the widely used CLASS of GILDAS 4 is that it is well documented, easy to install, and with a simple interface. HfS can fit, in a single run, multiple velocity components to the spectra of a FITS data cube, and obtain the maps of the parameters fitted and derived from the fit, in a short computing time taking advantage of the multiple processors of current computers.
The structure of the paper is as follows: The fit parameters used by HfS are described in §2, the fitting strategy in §3, the calculation of the synthetic spectrum in §4, the line parame-ters derived from the fit parameters in §5, the NH 3 physical parameters derived by HfS nh3 in §6, the error estimation for the fit parameters in §7 and for the derived parameters in §8. A comparison with the results obtained with CLASS and the Rosolowsky et al. (2008) routine is presented in §9. The different procedures that compose HfS are described in §10. Finally, in several Appendices we describe the requisites and installation instructions (Appendices A to C), and examples of use of HfS, and of input and output files (Appendices D to G).

FIT PARAMETERS
The general HfS procedures fit simultaneously, for every velocity component of a transition with hyperfine structure, four independent parameters, • ∆V, hyperfine lines linewidth, assumed to be the same for all the hyperfine lines, • V LSR , main line central LSR velocity, • A * m ≡ A(1 − exp{−τ m }), main line peak intensity (for hyperfine lines wider than the hyperfine separation and the channel width), where A is the amplitude (see §5), • τ * m ≡ 1 − exp{−τ m }, where τ m is the optical depth of the main line.
The HfS nh3 procedures included in HfS fit simultaneously the hyperfine structure of a pair of spectra of the NH 3 inversion transitions (J, K) = (1, 1) and (2, 2) (see Table 1).
For each set of fit parameters, the optical depth of the (2,2) main line is obtained from the relation (1) TABLE 1 Relative velocity V hyp , and optical depth τ hyp with respect to that of the main line, τ m , of the hyperfine lines of the NH 3 (J, K) = (1, 1) and (2,2) inversion transitions (Mangum & Shirley 2015). (os: outer satellite, is: inner satellite, m: main) The use of A * m ≡ A(1 − exp{−τ m }) and τ * m ≡ 1 − exp{−τ m } as fit parameters instead of A or Aτ m , and τ m needs some justification. First, A * m is well determined even in the two limiting cases, τ m 1 (with A ill-determined) and τ m 1 (with Aτ m ill-determined). This is because A * m is a true physical magnitude, i.e. the peak intensity of the main line, provided that the the hyperfine lines that may compose the main line are blended. Second, τ * m is bounded, 0 < τ * m < 1, and the two bounds correspond to the optically thin and thick limiting cases, where the fit is degenerate, and no longer depends on τ m . In short, the two fit parameters are well determined in all cases, even in the optically thin and thick limits. The only drawback of the use of τ * m as fit parameter is for the case of extremely high opacity in the simultaneous fit of NH 3 (1, 1) and (2,2), when the intensity of the (2, 2) main and satellite lines is similar. In its present version, HfS nh3 is not able to deal accurately with τ 1m > 16, but this case is indeed very rare. In a survey of low-mass cores in Perseus (Rosolowsky et al. 2008), none of the sources has such a high value of optical depth, and in another survey of high-mass clumps in the Galactic Plane (Svoboda et al. 2016) only 1% of the sources show τ 1m > 16. Thus, HfS nh3 should work for the vast majority of cases.
The fitting procedure ends with a set of four or, for the HfS nh3 procedures, six values plus τ * 2m , of the fit param-eters for every velocity component, which minimize the fit residual χ 2 of the spectrum (see §3), and an estimation of the uncertainty of the fit parameters (see §7), σ(∆V), σ(V LSR ), σ(A * m ), σ(τ * m ), and, for the HfS nh3 procedures, σ(V LSR2 ), and σ(A * 2m ).
3. FITTING STRATEGY The fitting procedure is similar to that used in other fitting problems by Estalella et al. (2012) and Palau et al. (2014). HfS samples the space parameter of dimension m (four or six times the number of velocity components), defined by the parameters p 1 , . . . , p m , to find the minimum value of the fit residual χ 2 , where y obs i are the observed line intensities, for a total of N spectral channels, y mod i (p 1 , . . . , p m ) are the model line intensities, depending on m free parameters, and σ i are the errors of the observations.
The sampling strategy is based on that used in AGA (Asexual Genetic Algorithm) (Cantó, Curiel, Martínez-Gómez 2009). The procedure starts with a number of samples of the parameter space within the initial search range for each parameter (seeds). For each seed, the residual χ 2 is computed for a number of samples (descendants) within the search range centered on the seed value. The best samples, i.e. those with the lowest χ 2 , are kept as seeds for the next loop, for which the search ranges are decreased by a constant factor. The procedure is iterated, and is stopped after a given number of loops.
Several sampling methods of the m-dimensional parameter space are possible, i.e. regular grid, random, Halton or Sobol pseudo-random sequences (Halton 1964;Sobol 1967). HfS uses a Sobol pseudo-random sequence because it samples the m-dimensional parameter space more evenly than a purely random sequence, even for high values of m, and the convergence of the fitting procedure to the minimum of χ 2 is faster. 4. SYNTHETIC SPECTRUM Let us assume that the transition being fitted has n h hyperfine lines, and for each hyperfine j = 1, . . . n h , V j hyp is the velocity shift with respect to the main line, and τ j hyp is the ratio of the optical depth of the hyperfine and that of the main line, τ m . Let n c be the number of velocity components fitted, and for each velocity component i = 1, . . . n c , ∆V i , V i LSR , A i , and τ i m are the values of the fit or derived line parameters (see §5).
For each velocity channel k = 1, . . . N, with central velocity V k and channel width ∆V ch , the contribution of the velocity component i to the optical depth is calculated as where G is the integral over the velocity range of the channel, G is evaluated by means of the error function erf ( where If the difference between x + and x − is too low, within less than 1 part in 10 4 , Equation 5 can produce a large roundoff error, and G is evaluated as Finally, the intensity of the channel k of the synthetic spectrum is calculated as the sum of intensities of the different velocity components, 5. DERIVED LINE PARAMETERS From the values of the fit parameters, the following derived line parameters are calculated, which are necessary for calculating the synthetic spectra and the estimation of parameters with physical interest: where f is the beam filling factor, T ex the excitation temperature, T bg the background temperature, and is the Planck-corrected temperature. A is calculated from In the case of HfS nh3, the excitation temperature T ex and the filling factor f are assumed to be the same for the (1, 1) and (2,2) lines. Since the frequencies of both transitions are very close (see Table 2), these assumptions imply that the amplitude A is the same for both lines, within less than 1 in 10 4 .
Special care has to be taken when τ m τ * m 1, since the last expression involves the difference of 1 and a number near 1. In this case, a good approximation is the Taylor expansion In the case of HfS nh3, τ 1m and τ 2m are derived in the same way.
Aτ m , amplitude times the main line optical depth. -Calculated from Note that for τ * m 1, Aτ m A * m . In the case of HfS nh3, Aτ 1m and Aτ 2m are derived in the same way. 6. NH 3 DERIVED PHYSICAL PARAMETERS In addition to the derived line parameters, HfS nh3 derives physical parameters using the standard analysis of NH 3 (1, 1) and (2,2) observations, which assumes that the region observed is homogeneous along the line of sight. The physical parameters derived are the excitation temperature T ex , the NH 3 (1, 1) and (2,2) beam-averaged column densities f N(1, 1) and f N (2,2), the rotational temperature T rot , the NH 3 beam-averaged column density f N(NH 3 ), and the kinetic temperature T k .
Excitation temperature. -The excitation temperature is obtained from the amplitude A, where T ν = hν/k is the frequency in temperature units. The frequency of the (1, 1) transition is used, but since the frequency of the (2, 2) transition is very close (see Table 2), the result does not depend on which of the two frequencies is used. For values of T ex T bg > T ν , this expression simplifies to The value of the excitation temperature depends on the value assumed for the filling factor f . The usual assumption is that the filling factor f = 1. The value obtained with this assumption is a lower limit for the value of T ex . On the contrary, if we assume that f 1, T ex → ∞.
NH 3 (1, 1) and (2, 2) column densities. -The column density of the (J, K) = (1, 1) and (2, 2) levels (i.e. the sum of column densities of the two inversion levels of the corresponding (J, K) rotational level) can be given as (Anglada et al. 1995;Estalella & Anglada 1997) where A jk is the Einstein coefficient of the inversion transition of the rotational level (J, K), R m = τ tot /τ m is the ratio of total and main line optical depths of the inversion transition, and T ν = hν JK /k is the frequency of the inversion transition in temperature units (see Table 2). The last expression depends on the value of the filling factor f assumed to derive T ex . The explicit dependence on f is so that the beam-averaged column density can be expressed as The maximum value of f N(J, K) is obtained for f = 1 (the usual assumption to derive T ex ), while the minimum value is obtained for f 1. In the latter case (or for A T ex T bg > T ν ), the expression simplifies to  (2,2) inversion transition frequencies (Kukolich 1967), spontaneous emission Einstein coefficients (Osorio et al. 2009), ratios of total to main line optical depths (Mangum & Shirley 2015), and B(J, K) and C(J, K) coefficients of Anglada et al. (1995) recalculated with the improved values of the constants in this table. The values of the constants appearing in these equations are given in Table 2.
In practical units the two equations become (Anglada et al. 1995) and, for f The values of the constants B(J, K) and C(J, K) for the (1, 1) and (2,2) transitions are given in Table 2. The expression equivalent to these two equations, with an explicit dependence on f , is, Rotational temperature. -The rotational temperature is obtained from the ratio of (1, 1) and (2, 2) column densities, where E 11 , E 22 , and g 11 , g 22 , are the energies and degeneracies of the corresponding levels. In practical units the equation becomes (see Table 3),  (Poynter & Kakar 1975;Mangum & Shirley 2015). Note that the values of the energies reported by these authors are slightly different from those given in Ho & Townes (1983).
NH 3 column density. -The ammonia total column density is usually estimated assuming that the population ratios between pairs of levels are given by the same rotational temperature (CTEX approximation), and that only the metastable levels J = K, up to (J, K) = (3, 3), are populated. While this assumption is reasonable for low-mass dense cores with moderate temperatures, it is not appropriate for hot cores, where NH 3 inversion transitions (6,6) to (14,14) are detected (see for instance Goddi, Zhang, & Moscadelli 2015).
In addition, since the NH 3 (1, 1) and (2, 2) levels correspond to para-NH 3 , no information on the ortho-NH 3 is obtained from the (1, 1) and (2, 2) spectra. Thus, an ortho-to-para ratio has to be assumed. The usual assumption is to take an orthoto-para ratio of 1, although a value of ∼ 0.7 has been observed in some star-forming regions (see for instance the discussion in Faure et al. 2013). We will assume an ortho-to-para ratio of 1.
With the former assumptions, with the partition function Q given by In practical units (see Table 3 Kinetic temperature. -The kinetic temperature can be taken to be equal to the rotational temperature T rot , but a better estimation can be given taking into account collisional transitions to other levels. A 3-level approximation considering only the rotational levels (1, 1), (2,2), and (2, 1) where T 0 = (E 22 − E 11 )/k, and C(22 → 21) and C(22 → 11) are the collisional excitation and desexcitation rates between the corresponding levels. Calculations involving more rotational levels were performed by Danby et al. (1988), and more recently by Maret et al. (2009) with T 0 = 40.99 K, (E 21 − E 22 )/k = 16.26 K, and the numerical value 0.73 was determined to fit the results shown in Fig.  5 of Maret et al. (2009). Note that this relation implies that T rot is always below a value T rot = 40.99/ ln 1.73 = 74.8 K (see Fig.1). Given a value of T rot , the implicit equation must be solved to find T k . A possible iterative algorithm to solve the equation is 7. ERROR ESTIMATION OF THE FIT PARAMETERS 7.1. Parameter-space confidence region Let us assume that the best fit to the observed data is obtained for values p 0 k of the parameters, for which the residual χ 2 is minimum, where σ i is the error of y obs i . The uncertainty, σ(p k ), in the values derived for the parameters p k can be estimated as the projection of the confidence region of the m-dimensional space parameter for which χ 2 does not exceed the minimum value by an amount ∆(m, α), where α is the significance level (0 < α < 1). Following Avni (1976) and Wall & Jenkins (2003), the probability is that of a chi-square distribution with m degrees of freedom. Thus, ∆(m, α) is the increment of χ 2 such that if the observation is repeated a large number of times, a fraction α of times the values of the parameters fitted will be inside the confidence region, i.e. in the interval p k ± σ(p k ). The values of ∆(m, α) for significance levels equivalents to 1, 2, and 3 sigmas for a Gaussian error distribution, and different values of m are shown in Table 4. Assuming that the model fits well the observations, i.e. χ 2 min N − m, the condition can be written as This expression is useful since it can be given in terms of the weighted rms fit residual, σ, for which we obtain that the confidence region is given by the parameter values that increase the rms fit residual to This last expression can be used even when the errors of the observations are unknown.

7.2.
Modeling the fit residual χ 2 7.2.1. Quadratic approximation of the fit residual Let us assume that around its minimum value, the residual χ 2 can be approximated by a quadratic function, where x j = p j − p 0 j are the increment of the parameter values from their best-fit values. The confidence region of the m-dimensional parameter space will be the region inside the surface m i, j=1 In array form, the quadric equation can be expressed as with c = −∆(m, α), or, with the obvious definitions for arrays A, B, and X, In order to estimate the uncertainties of the fit parameters, σ(p k ), we have to calculate the projections of the quadric onto each axis k. If the fit residual is well behaved, we may expect that the quadric is an ellipsoid, i.e. its projections onto the plane defined by any pair of parameters is an ellipse, and the ellipsoid has finite projections onto any axis.

Projections of an ellipsoid
The projections of the ellipsoid onto each coordinate axis can be found as the intersections of the hyperplane perpendicular to the axis, tangent to the ellipsoid. The equation of the tangent hyperplane at a point (x 0 1 , . . . , x 0 m ) of the ellipsoid is given by (McConnell 2011) The equation of an hyperplane perpendicular to the k axis is and the coefficients of x j for j k have to be zero, Eqs. 45 form a system of m − 1 linear equations with m unknowns, the coordinates of the tangent point. We can consider that x 0 j ( j k) are the m − 1 unknowns of the system, which can be derived as a function of x 0 k . The system of equations can be written with arrays of size m, where A k is the array A with zeros in the k row and column, a ik = a ki = 0, and 1 in the diagonal term, a kk = 1, and Z k and B k have 1 and 0 respectively in the k row, From the system we can derive the solution giving x 0 j in terms of x 0 k , which can be expressed as where d kk = 1 and e kk = 0. The tangent point must fulfill the tangent hyperplane equation, Eq. 44. By substitution of Eq. 49 in the tangent hyperplane equation, and setting x k = x 0 k , we get a second degree The two solutions of the equation provide the two kcoordinates of the projections of the ellipsoid onto the k axis.

Case of a centered ellipsoid
Since we are assuming that we know that the minimum residual χ 2 is well determined, the quadratic function of Eq. 51 must have a minimum at the origin, and the linear terms in the variables vanish because the partial derivatives at the origin must be zero. This means that we can assume that the quadric is centered, and its equation becomes In this case, the equation of the tangent plane is simpler, since B = 0, C k = 0, and E k = 0. The second degree equation (Eq. 50) becomes and the projections are symmetric, ±x 0 k .

Coefficients of the ellipsoid
The centered ellipsoid depends on m(m+1)/2 symmetric a i j coefficients. Let us examine a sufficient number of constraints to derive these coefficients.
Diagonal coefficients a ii . -For each parameter i we can obtain constraints from the values of the residual χ 2 increment for different increments x i of the parameter p i , while keeping the rest of parameters constant, For each value of the increments of the parameter and the residual χ 2 increment, x n i , ∆ n i , we have a constraint on a ii , Although a single value is enough to determine a ii , at least two are recommended, above and below the best-fit value, i.e. with x i > 0 and x i < 0. In general, the best approximation for a ii (so that the sum of the squares of [∆ n i − a ii (x n i ) 2 ] is minimum) is given by where the sums are for all increments evaluated. Note that, provided that the ∆ n i are positive, the diagonal term a ii will always be positive. For a good characterization of the behavior of χ 2 , at least two values of ∆ n i should be close to ∆(m, α). The intersections of the ellipsoid with the coordinate axis i are ±(∆(m, α)/a ii ) 1/2 . For the case of statistically independent parameters, these intersections will coincide with the projections of the ellipsoid, since the coordinate axes are the principal axes of the ellipsoid and the cross-terms of the ellipsoid vanish. However, in general, there will be some dependence among the parameters, and the cross terms will not be zero.
Cross-coefficients a i j (i j). -The constraints to derive the cross terms a i j can be obtained from the value of the residual χ 2 for the simultaneous increment of the two parameters i and j, while keeping the rest of parameters constant, χ 2 (p 0 1 , . . . , p 0 i + x i , . . . , p 0 j + x j , . . . , p 0 m ) = χ 2 min + ∆ i j (56) For each pair of increments of the parameters and the residual χ 2 increment, x n i , x n j , ∆ n i j , we have a constraint on a i j , Taking into account Eq. 54, and defining δ n Although a single pair of values x i , x j , is enough to determine a i j , at least four are recommended, above and below the bestfit value for each parameter, i.e. with the four combinations of x i > 0, x i < 0, x j > 0, and x j < 0. In general, the best approximation for a i j (in the same sense as for the diagonal term) is given by where the sums are for all the pairs of increments evaluated.

Practical case
A practical implementation of the procedure is as follows. Let us call 1. For each parameter p i we estimate a positive and a negative increment, x + i and x − i , such that the increment of the residual χ 2 is close to ∆(m, α), The diagonal coefficient a ii (Eq. 54) is given by 2. For each pair of parameters p i , p j (i j) we calculate The cross-coefficient a i j (Eq. 59) is given by 3. For each parameter k we construct the arrays A k , Z k (Eq. 47), and solve the system of linear equations A k D k = Z k . Finally, we calculate the quadratic coefficient, m i=1 a ki d ik (Eq. 52), and the projection In some cases, the quadratic approximation of χ 2 is not good enough, and the determination of the projections can fail: for some parameter p k the array A k may have null determinant, or the quadratic coefficient of Eq. 52 may be negative. In these cases, a rough estimation of the uncertainty in p k can still be given as the intersection with the k axis, [∆(m, α)/a kk ] 1/2 . 8. ERROR ESTIMATION OF THE DERIVED PARAMETERS All the derived parameters depend on m fit parameters (four or six times the number of velocity component). Let us call the values of the fit parameters, and σ(p k ) their errors, found from the increase in the fit residual χ 2 (see §7). Let d be any of the parameters derived from the fit parameters, d = d(p 1 , . . . , p m ), for instance τ m or N(NH 3 ). For every fit parameter p 0 k (k = 1, . . . , m) we evaluate the values of the derived parameter when we increase and decrease the value of the k-th fit parameter by its error σ(p k ), . . , p 0 m ), k = 1, . . . , m (66) Assuming that the errors of the fit parameters are statistically independent, we can estimate the error σ(d) as 1.221 ± 0.194 −3.029 ± 0.061 6.2 ± 0.8 0.2 ± 0.5 7.1 ± 3.5 0.3 ± 0.9 HfS 1.240 ± 0.091 −3.030 ± 0.027 · · · · · · 7.1 ± 1.0 0.

COMPARISON OF HfS WITH OTHER ROUTINES
The fits obtained with HfS and HfS nh3 were compared with those obtained with other commonly used routines. In Tables 5 and 6 we show some examples of fits performed with HfS and CLASS for the NH 3 (1, 1) line, and the Rosolowsky routine (Rosolowsky et al. 2008) for the simultaneous fit of the NH 3 (1, 1) and (2,2) lines. The fits shown cover cases of low and high optical depth, and narrow and wide linewidths. As can be seen in the tables, there is in general agreement between the HfS results and the other routines.
Regarding the comparison with CLASS, in the case of linewidths lower than the channel width (spectra 3 and 4 in Table 5), CLASS gives a fixed value for the linewidth, higher than the channel width and much higher than the actual linewidth. As a consequence the values of Aτ m and τ m given by CLASS are scaled down by roughly the same factor, so that CLASS is consistent with HfS in the values of Aτ m ∆V and τ m ∆V. In general, the errors given by CLASS appear to be underestimated.
The comparison with the Rosolowsky routine was performed for the examples of fits given in Rosolowsky et al. (2008). The fitted parameters were taken from Table 3 of the electronic edition of Rosolowsky et al. (2008), and the raw spectral data were retrieved from the COMPLETE Web site. 5 For the HfS nh3 fits, the data were Hanning-smoothed with HFHW = 2 and fitted using Nksample = 400 (see below). For the parameters obtained specifically from the simultaneous fit of the NH 3 (1, 1) and (2,2) lines, there is a good agreement in the values obtained for T k , and the values of ammonia column density reported in Rosolowsky et al. (2008) lie between the two limiting cases given by HfS nh3, for filling factors f 1 and f = 1.
10. DESCRIPTION OF THE HfS PROCEDURES While the HfS nh3 procedures fit the NH 3 (1, 1) and (2,2) transitions simultaneously, the general HfS procedures fit a single transition, selected among those stored in the file hfs transitions.dat. This file has to be located in the working directory, or in the directory pointed at by the environment variable HFS DIR (see Appendix A). The first transition in the file is a single line at V LSR = 0, useful for fitting a single Gaussian line. Other transitions in the file include NH 3 (1, 1) and (2,2), NH 2 D(1 11 , 1 10 ), N 2 H + (1-0), CN, HCN, H 13 CN, C 2 H, C 2 D, C 17 O, and more transitions can be added easily to the file. The criterion for defining the "main component" of a transition is, in general, all the hyperfines with a velocity offset less than 0.001 km s −1 , and for NH 3 (1, 1) and (2,2), less than 0.6 km s −1 . Any transition not appearing in the file is assumed to be single.
An optional Hanning smoothing of the spectrum can be performed prior to fitting. You can select the Hanning filter half-width, HFHW. For HFHW = 0 no smoothing is performed. HFHW = 1 is the standard 3-point Hanning smoothing, and the resulting spectrum has half the initial number of channels. In general, the Hanning smoothing encompasses 2*HFHW+1 points, resulting in a final number of channels HFHW+1 times smaller.
The iterative process is controlled by two parameters, Nksample, the number of thousands of samples of the parameter space, and Final Range, the ratio of ranges of the last loop and initial search ranges. The number of loops is taken as nloop = Nksample 1/2 , and the number of seeds and descendants is taken as nseed = ndesc = (nloop × 1000) 1/2 , so that nloop × nseed × ndesc = Nksample × 1000. For each loop the ranges will be decreased a factor of f = Final Range 1/(nloop−1) .
The initial values for the fit are guessed from the intensity, position, and width of the data peak for the first component, and of the residual (data minus previous components) for the rest of components (up to a maximum of 9). The main line optical depth is set arbitrarily to 0.7 (1 − e −τ m = 0.5), or to an arbitrary low value (10 −6 ) for a single Gaussian fit.
Values for the initial search ranges are calculated by the procedures. For fitting a single Gaussian line the search range of τ * m is made 0 to keep it constant. Additional constraints for the search ranges of the parameters are A * m > 0; 0 < τ * m < 1; and ∆V > ∆V min , where ∆V min is the thermal linewidth for a kinetic temperature T k = T bg of a large molecule (for the general HfS procedures, 0.025 km s −1 for a mass of 200 m H ), or of NH 3 (0.086 km s −1 for the HfS nh3 procedures).
Note that the final value of a fit parameter can be outside the initial range for the parameter. For an initial range 2r, and a range decreasing factor per loop f , the final value of a 5 https://www.cfa.harvard.edu/COMPLETE/data_html_pages/ GBT_NH3.html parameter can can be up to roughly r/(1 − f ) apart from its initial value. For instance, for the default values Nksample = 200 (corresponding to nloop = 14), and Final Range = 0.05, we have f 0.8 and r/(1 − f ) 5r.
The different procedures that compose HfS are described in the following.
10.1. hfs fit, hfs nh3 These are interactive graphic procedures for fitting simultaneously multiple velocity components of an spectral line with hyperfine structure, or to a pair of NH 3 (1, 1) and (2,2), to spectra read from data files, and generating files with the synthetic spectra. See an example of a fitting run of hfs fit in Appendix D and of hfs nh3 in Appendix E.
You can set any number of velocity components, for which the procedures propose a first guess. Alternatively, you can use the cursor to add or delete components. The position of a new component is set at the cursor position, its intensity is the intensity at the cursor position, and the FWHM is estimated around the cursor position. The values of the initial search ranges can be changed. Any of the parameters can be kept constant by setting its range to 0.

INPUT
The spectra to fit are read from ASCII files with a pair of values (velocity, intensity) per line. Lines beginning with "!" or "#" are ignored. The file(s) can be given as argument(s) to hfs fit and hfs nh3: $ hfs_fit <source>.dat $ hfs_nh3 <source_11>.dat <source_22>.dat OUTPUT • hfs fit.log or hfs nh3.log, log file with the details of the fitting session.
• <source>.synt, or <source 11>.synt and <source 22>.synt, ASCII files with the synthesized spectra. The files have a header (lines beginning with "!") with the values of the parameters of the fitted spectra for each velocity component. The files are readable by GREG of the package GILDAS.The file is overwritten for every new fit. See an example in Appendix G.
• <source>.eps, plot showing the data, the components fitted, and the residual. See some examples in Fig. 2. The file is overwritten for every new fit.
10.2. hfs file. This is a batch procedure to fit the hyperfine structure of the same transition (and a single velocity component) for a set of data files, and generating a list of files with the synthetic spectra. The fit procedure is the same as that of hfs fit, but it is not interactive.

INPUT
The input is a parameter file that can be given as argument when running the procedure:

Iteration parameters: Nksample, Final Range.
3. and following lines: list of files, one file per line See an example of <file list>.par in Appendix F.

OUTPUT
• <file list>.log, log file with the details of the fitting process for all the files in <file list>.par.
• <file list>.out, ASCII file with the values and uncertainties of the parameters fitted for each file of the list. The first lines beginning with "!" are the header and give information about the transition and the column headers.
• <file list>.ps, PostScript file with plots of the data, the components fitted, and the residual for all the files in the list.
• <file #>.synt, an ASCII file for each file in <file list>.par, with the synthesized spectrum. The file has a header (lines beginning with "!") with the values of the parameters of the fitted spectrum.
10.3. hfs cube sp, hfs nh3 cube sp These are single processor batch procedures for fitting simultaneously multiple velocity components of spectra from 3-axes FITS data cubes. A subimage of the FITS data cube can be selected, and optional boxcar averaging of pixels and Hanning filtering of the spectra can be performed. The procedures are similar to hfs fit and hfs nh3, but they are not interactive. The maximum dimension of the image and the maximum number of channels of the data cube is arbitrary, i.e. the only limitation is the amount of memory available by the computer.

INPUT
The input is a parameter file that has to be given as argument when running the procedures: $ hfs_cube_sp <parameter>.par $ hfs_nh3_cube_sp <parameter>.par The parameter file for hfs cube sp has 10 lines (see an example in Appendix F): 1. Transition name, e.g. "C17O(1-0)".
2. FITS data cube file to read. The file must have 3 nondegenerate axes, in this order: x position, y position, velocity channel.
3. Rms of channels without emission, minimum SNR of spectra to be fitted. Components with a peak intensity below SNR times rms are not fitted.
5. Range of channels for each component: 2 × Ncomp values, with the first and last channel of the velocity range for each component. The channel ranges must be non overlapping. For the first component, 0 defaults to 1 (first), nchan (last).
The parameter file for hfs nh3 cube sp has 11 lines (see an example in Appendix F): 1. (1, 1) FITS data cube file to read. The file must have 3 non-degenerate axes, in this order: x position, y position, velocity channel.

2.
(1, 1) rms of channels without emission, minimum SNR of spectra to be fitted. Components with a peak intensity below SNR times rms are not fitted.
3. Number of velocity components to fit, Ncomp, between 1 and 9. 5. (2,2) FITS data cube file to read. The file must have the same geometry as that of the (1, 1) FITS file.
6. (2, 2) rms of channels without emission, minimum SNR of spectra to be fitted. Components with a peak intensity below SNR times rms are not fitted.
The 5 following lines are the same as for the case of hfs cube sp.

OUTPUT
• log/<parfile>.log, log file in folder log with the details of the fitting process for all the pixels of the subimage.
• <parfile> comp#.out, an ASCII file for each velocity component with the values of the parameters fitted and the line and physical parameters for each pixel of the subimage, and their uncertainty. The first lines beginning with "!" are the header and give information about the parameter file, FITS files, velocity component number, velocity range of the component, Hanning filtering applied, smoothing boxcar radius, and column headers. The files are used by hfs view, and are readable by GREG. See an example in Appendix G.
• ps/<parfile> <xoffset>.ps, PostScript files in folder ps, with plots of the data, the components fitted, and the residual for all pixels of the subimage with a given <xoffset> in arcsec.
• maps/<parfile> <parameter> comp#.fits, FITS files in folder maps, with maps, for each velocity component, of the parameters fitted and the line and physical parameters, and their uncertainty. Each FITS file has two planes, the first plane with the values of the parameter, and the second plane with the uncertainties. See in Fig. 3 an example of the the maps obtained from a 2-velocity-components fit of the C 18 O line for a FITS data cube.
10.4. hfs cube mp, hfs nh3 cube mp These are multiprocessor procedures that use Open MPI (see Appendix A) and run in parallel using a number of processors available in the machine, or in more than one host. The multiprocessor procedures are naked versions of the single-processor versions, hfs cube sp and hfs nh3 cube sp, without any graphic output. The instructions for running these procedures can be found in Appendix C.
10.5. hfs view, hfs nh3 view These are interactive graphic procedures for displaying spectra from FITS data cubes and the corresponding synthetic spectra fitted with hfs cube or hfs nh3 cube. A plot of the integrated intensity is shown, and you can select with the mouse the position for which the the spectrum (data and synthetic) is shown. You can select to show the integrated intensity for all channels, or for the channel ranges of each velocity component. The data and synthetic spectra at any pixel can be extracted in ASCII files.

INPUT
The input is the same parameter file used as input to hfs cube or hfs nh3 cube. From the information in this file, hfs view reads the corresponding FITS data cubes and the <parfile> comp#.out file for each velocity component (created by hfs cube or hfs nh3 cube), with the parameters of the fitted spectra.
10.6. hfs blanking This is an auxiliary procedure to flag the output FITS files of hfs cube and hfs nh3 cube, according to the parameter values, errors, or SNR. The FITS files have two planes, the first one with the values of the parameter at each pixel, and the second one with with the error of the value. The procedure reads the values and errors in the FITS files, and allows you to blank the pixels that fulfill a criterion on parameter error, value, or value/error, above or below a cutoff value.
INPUT <parfile> <parameter> comp#.fits, any FITS file created by hfs cube or hfs nh3 cube, for each velocity component and parameter.
OUTPUT <parfile> <parameter> comp# blank.fits, output 2axes FITS file, with the values of the parameter for nonflagged pixels, and NaN for flagged-out pixels. 10.7. hfs synt Auxiliary procedure to create a synthetic spectrum, with the option of adding Gaussian noise. You can select the transition, number of channels, spectral resolution, line parameters, and noise level.

OUTPUT
• hfs synt.synt, ASCII file with the synthetic spectrum generated. The file has a header (lines beginning with "!") with the values of the parameters used for generating the synthetic spectrum.
• hfs synt.eps, plot of the synthetic spectrum.
The author thanks Pau Estalella for suggesting the use of pseudo-random sequences, Ferran Sala for helpful discussions on the projections of a quadric, and Salvador Curiel for helping with the implementation of the multiprocessor procedures. Thanks also toÁlvaro Sánchez-Monge for reading the manuscript and, together with Aina Palau, Gemma Busquet and Carmen Juárez, for testing HfS, finding bugs, and suggesting improvements. This work has been partially supported by the Spanish MINECO grant AYA2014-57369-C3 (cofunded with FEDER funds) and MDM-2014-0369 of IC-CUB (Unidad de Excelencia 'María de Maeztu').