document: the goal of the development of the code was to have a simple and efficient tool for the computation of adiabatic oscillation frequencies and eigenfunctions for general stellar models , emphasizing also the accuracy of the results . 
   not surprisingly , given the long development period , the simplicity is now less evident . 
   however , the code offers considerable flexibility in the choice of integration method as well as ability to determine all frequencies of a given model , in a given range of degree and frequency . 
   the choice of variables describing the equilibrium model and oscillations was to a large extent inspired by @xcite . 
   as discussed in section  [ sec : eqmodel ] the equilibrium model is defined in terms of a minimal set of dimensionless variables , as well as by mass and radius of the model . 
   fairly extensive documentation of the code , on which the present paper in part is based , is provided with the distribution packagejcd / adipack.n ] . 
   @xcite provided an extensive review of adiabatic stellar oscillations , emphasizing applications to helioseismology , and discussed many aspects and tests of the aarhus package , whereas @xcite carried out careful tests and comparisons of results on polytropic models ; this includes extensive tables of frequencies which can be used for comparison with other codes . 
   the equilibrium model is defined in terms of the following dimensionless variables : @xmath0 here @xmath1 is distance to the centre , @xmath2 is the mass interior to @xmath1 , @xmath3 is the photospheric radius of the model and @xmath4 is its mass ; also , @xmath5 is the gravitational constant , @xmath6 is pressure , @xmath7 is density , and @xmath8 , the derivative being at constant specific entropy . 
   in addition , the model file defines @xmath4 and @xmath3 , as well as central pressure and density , in dimensional units , and scaled second derivatives of @xmath6 and @xmath7 at the centre ( required from the expansions in the central boundary condition ) ; finally , for models with vanishing surface pressure , assuming a polytropic relation between @xmath6 and @xmath7 in the near - surface region , the polytropic index is specified . 
   the following relations between the variables defined here and more `` physical '' variables are often useful : @xmath9 we may also express the characteristic frequencies for adiabatic oscillations in terms of these variables . thus if @xmath10 is the buoyancy frequency , @xmath11 is the lamb frequency at degree @xmath12 and @xmath13 is the acoustical cut - off frequency for an isothermal atmosphere , we have @xmath14 where @xmath15 is the adiabatic sound speed , and @xmath16 is the pressure scale height , @xmath17 being the gravitational acceleration . 
   finally it may be noted that the squared sound speed is given by @xmath18 these equations also define the dimensionless characteristic frequencies @xmath19 , @xmath20 and @xmath21 as well as the dimensionless sound speed @xmath22 , which are often useful .      as is well known the displacement vector of nonradial ( spheroidal ) modes 
   can be written in terms of polar coordinates @xmath23 as @xmath24 \exp ( - { { \rm i}}\omega t ) \right\ } \ ; . \nonumber\end{aligned}\ ] ] here @xmath25 is a spherical harmonic of degree @xmath12 and azimuthal order @xmath2 , @xmath26 being co - latitude and @xmath27 longitude ; @xmath28 is an associated legendre function , and @xmath29 is a suitable normalization constant . 
   also , @xmath30 , @xmath31 , and @xmath32 are unit vectors in the @xmath1 , @xmath26 , and @xmath27 directions . finally , @xmath33 is time and @xmath34 is the angular frequency of the mode . 
   similarly , e.g. , the eulerian perturbation to pressure may be written @xmath35 \ ; . \label{eq : e2.2}\ ] ] as the oscillations are adiabatic ( and only conservative boundary conditions are considered ) @xmath34 is real , and the amplitude functions @xmath36 , 
   @xmath37 , @xmath38 , etc . can be chosen to be real . 
   the equations of adiabatic stellar oscillations , in the nonradial case , are expressed in terms of the following variables : , @xmath39 results from the earlier use of an unconventional sign convention for @xmath40 ; now , as usual , @xmath40 is defined such that the perturbed poisson equation has the form @xmath41 , where @xmath42 is the eulerian density perturbation . ] 
   @xmath43 here @xmath40 is the perturbation to the gravitational potential . also , we introduce the dimensionless frequency @xmath44 by @xmath45 corresponding to eqs [ eq : buoy ]  [ eq : cutoff ] . these quantities satisfy the following equations : @xmath46   y_1 + ( a - 1 )   y_2 + \eta a   y_3 \ ; , \\ \label{eq : ea.3 } x   { { { \rm d}}y_3   \over { { \rm d}}x } & = & y_3 + y_4 \ ; , \\ \label{eq : ea.4 } x   { { { \rm d}}y_4   \over { { \rm d}}x } & = & - a u   y_1    - u { v_g    \over \eta } y_2 \\ & & + [ l ( l + 1 ) + u(a - 2 ) +   u v_g ]   y_3 + 2(1 - u )   y_4 \ ; . \nonumber\end{aligned}\ ] ] here @xmath47 , and the notation is otherwise as defined in eq . 
   [ eq : fivea ] . in the @xcite approximation , where the perturbation to the gravitational potential is neglected , 
   the terms in @xmath48 are neglected in eqs  [ eq : ea.1 ] and [ eq : ea.2 ] and eqs  [ eq : ea.3 ] and [ eq : ea.4 ] are not used . 
   the dependent variables @xmath49 in the nonradial case have been chosen in such a way that for @xmath50 they all vary as @xmath51 for @xmath52 . for large @xmath12 a considerable ( and fundamentally unnecessary ) computational effort 
   would be needed to represent this variation sufficiently accurately with , e.g. , a finite difference technique , if these variables were to be used in the numerical integration . 
   instead i introduce a new set of dependent variables by @xmath53 these variables are then @xmath54 in @xmath55 near the centre . 
   they are used in the region where the variation in the @xmath56 is dominated by the @xmath51 behaviour , for @xmath57 , say , where @xmath58 is determined on the basis of the asymptotic properties of the solution . this transformation permits calculating modes of arbitrarily high degree in a complete model . for radial oscillations only @xmath59 and @xmath60 
   are used , where @xmath59 is defined as above , and @xmath61 here the equations become @xmath62 y_1 + a   y_2 \ ; . \label{eq : ea.6}\end{aligned}\ ] ]    the equations are solved on the interval @xmath63 $ ] in @xmath55 . here , 
   in the most common case involving a complete stellar model @xmath64 , where @xmath65 is a suitably small number such that the series expansion around @xmath66 is sufficiently accurate ; however , the code can also deal with envelope models with arbitrary @xmath67 , typically imposing @xmath68 at the bottom of the envelope . 
   the outermost point is defined by @xmath69 where @xmath70 is the surface radius , including the atmosphere ; thus , typically , @xmath71 . 
   the centre of the star , @xmath72 , is obviously a singular point of the equations . 
   as discussed , e.g. , by @xcite boundary conditions at this point are obtained from a series expansion , in the present code to second significant order . in the general case 
   this defines two conditions at the innermost non - zero point in the model . for radial oscillations , or nonradial oscillations in the cowling approximation 
   , one condition is obtained . 
   the surface in a realistic model is typically defined at a suitable point in the stellar atmosphere , with non - zero pressure and density . here 
   the simple condition of vanishing lagrangian pressure perturbation is implemented and sometimes used . 
   however , more commonly a condition between pressure perturbation and displacement is established by matching continuously to the solution in an isothermal atmosphere extending continuously from the uppermost point in the model . 
   a very similar condition was presented by @xcite . in addition , in the full nonradial case a condition is obtained from the continuous match of @xmath40 and its derivative to the vacuum solution outside the star .    in full polytropic models , or other models with vanishing surface pressure , the surface is also a singular point . in this case 
   a boundary condition at the outermost non - singular point is obtained from a series expansion , assuming a near - surface polytropic behaviour ( see * ? ? ? * for details ) . 
   the code also has the option of considering truncated ( e.g. , envelope ) models although at the moment only in the cowling approximation or for radial oscillations . in this case 
   the innermost boundary condition is typically the vanishing of the radial displacement @xmath73 although other options are available . 
   the numerical problem can be formulated generally as that of solving @xmath74 with the boundary conditions @xmath75 @xmath76 here the order @xmath77 of the system is 4 for the full nonradial case , and 2 for radial oscillations or nonradial oscillations in the cowling approximation . 
   this system only allows non - trivial solutions for selected values of @xmath78 which is thus an eigenvalue of the problem . 
   the programme permits solving these equations with two basically different techniques , each with some variants . 
   the first is a shooting method , where solutions satisfying the boundary conditions are integrated separately from the inner and outer boundary , and the eigenvalue is found by matching these solutions at a suitable inner fitting point @xmath79 . 
   the second technique is to solve the equations together with a normalization condition and all boundary conditions using a relaxation technique ; the eigenvalue is then found by requiring continuity of one of the eigenfunctions at an interior matching point . 
   for simplicity i do not distinguish between @xmath80 and @xmath56 ( cf . section  [ sec : eq ] ) in this section . 
   it is implicitly understood that the dependent variable ( which is denoted @xmath56 ) is @xmath80 for @xmath57 and @xmath56 for @xmath81 . 
   the numerical treatment of the transition between @xmath80 and @xmath56 has required a little care in the coding . 
   it is convenient here to distinguish between @xmath77 = 2 and @xmath77 = 4 . for @xmath77 = 2 the differential eqs  [ eq : e3.1 ] 
   have a unique ( apart from normalization ) solution @xmath82 satisfying the inner boundary conditions [ eq : e3.2 ] , and a unique solution @xmath83 satisfying the outer boundary conditions [ eq : e3.3 ] . 
   these are obtained by numerical integration of the equations . 
   the final solution can then be represented as @xmath84 . 
   the eigenvalue is obtained by requiring that the solutions agree at a suitable matching point @xmath85 , say . 
   thus @xmath86 these equations clearly have a non - trivial solution @xmath87 only when their determinant vanishes , i.e. , when @xmath88 equation [ eq : e3.5 ] is therefore the eigenvalue equation . 
   for @xmath77 = 4 there are two linearly independent solutions satisfying the inner boundary conditions , and two linearly independent solutions satisfying the outer boundary conditions . 
   the former set @xmath89 is chosen by setting @xmath90 and the latter set @xmath91 is chosen by setting @xmath92 the inner and outer boundary conditions are such that , given @xmath59 and @xmath48 , @xmath60 and @xmath39 may be calculated from them ; thus eqs  [ eq : e3.6 ] and [ eq : e3.7 ] completely specify the solutions , which are obtained by integrating from the inner or outer boundary . 
   the final solution can then be represented as @xmath93 at the fitting point @xmath79 continuity of the solution requires that @xmath94 this set of equations only has a non - trivial solution if @xmath95 where , e.g. , @xmath96 . 
   thus eq . 
   [ eq : e3.9 ] is the eigenvalue equation in this case . 
   clearly @xmath97 as defined in either eq . 
   [ eq : e3.5 ] or eq . 
   [ eq : e3.9 ] is a smooth function of @xmath78 , and the eigenfrequencies are found as the zeros of this function . 
   this is done in the programme using a standard secant technique . 
   however , the programme also has the option for scanning through a given interval in @xmath78 to look for changes of sign of @xmath97 , possibly iterating for the eigenfrequency at each change of sign . 
   thus it is possible to search a given region of the spectrum completely automatically . 
   the programme allows the use of two different techniques for solving the differential equations . 
   one is the standard second - order centred difference technique , where the differential equations are replaced by the difference equations @xmath98 
   , \quad i = 1 , \ldots , i \ ; . 
   \label{eq : e3.11}\ ] ] here i have introduced a mesh @xmath99 in @xmath55 , where @xmath10 is the total number of mesh points ; @xmath100 , and @xmath101 . 
   these equations allow the solution at @xmath102 to be determined from the solution at @xmath103 . 
   the second technique was proposed by @xcite ; here on each mesh interval @xmath104 we consider the equations @xmath105 with constant coefficients , where @xmath106 . 
   these equations may be solved analytically on the mesh intervals , and the complete solution is obtained by continuous matching at the mesh points . 
   this technique clearly permits the computation of solutions varying arbitrarily rapidly , i.e. , the calculation of modes of arbitrarily high order . on the other hand solving eqs  [ eq : e3.12 ] involves finding the eigenvalues and eigenvectors of the coefficient matrix , and therefore becomes very complex and time consuming for higher - order systems . 
   thus in practice it has only been implemented for systems of order 2 , i.e. , radial oscillations or nonradial oscillations in the cowling approximation . 
   if one of the boundary conditions is dropped , the difference equations , with the remaining boundary condition and a normalization condition , constitute a set of linear equations for the @xmath107 which can be solved for any value of @xmath44 ; this set may be solved efficiently by forward elimination and backsubstitution ( e.g. , * ? ? ? 
   * ) , with a procedure very similar to the so - called henyey technique ( e.g. , * ? ? 
   * see also christensen - dalsgaard 2007 ) used in stellar modelling . 
   the eigenvalue is then found by requiring that the remaining boundary condition , which effectively takes the role of @xmath108 , be satisfied . 
   however , as both boundaries , at least in a complete model , are either singular or very nearly singular , the removal of one of the boundary conditions tends to produce solutions that are somewhat ill - behaved , in particular for modes of high degree . 
   this in turn is reflected in the behaviour of @xmath97 as a function of @xmath44 . 
   this problem is avoided in a variant of the relaxation technique where the difference equations are solved separately for @xmath109 and @xmath110 , by introducing a double point @xmath111 in the mesh . the solution is furthermore required to satisfy the boundary conditions [ eq : e3.2 ] and [ eq : e3.3 ] , a suitable normalization condition ( e.g. @xmath112 ) , and continuity of all but one of the variables at @xmath85 , e.g. , @xmath113 ( when @xmath77 = 2 clearly only the first continuity condition is used ) we then set @xmath114 and the eigenvalues are found as the zeros of @xmath97 , regarded as a function of @xmath78 . with this definition 
   , @xmath97 may have singularities with discontinuous sign changes that are not associated with an eigenvalue , and hence a little care is required in the search for eigenvalues . 
   however , close to an eigenvalue @xmath97 is generally well - behaved , and the secant iteration may be used without problems .    as implemented here 
   the shooting technique is considerably faster than the relaxation technique , and so should be used whenever possible ( notice that both techniques may use the difference eqs  [ eq : e3.11 ] and so they are numerically equivalent , in regions of the spectrum where they both work ) . for _ second - order systems _ the shooting technique can probably always be used ; the integrations of the inner and outer solutions should cause no problems , and the matching determinant in eq . 
   [ eq : e3.5 ] is well - behaved . for _ fourth - order systems _ 
   , however , this needs not be the case . for modes where the perturbation to the gravitational potential has little effect on the solution , the two solutions @xmath115 and @xmath116 , and similarly the two solutions @xmath117 and @xmath118 , are almost linearly dependent , and so the matching determinant nearly vanishes for any value of @xmath78 . 
   this is therefore the situation where the relaxation technique may be used with advantage . 
   this applies , in particular , to the calculation of modes of moderate and high degree which are essential to helioseismology .      to make full use of the increasingly accurate observed frequencies the computed frequencies should clearly at the very least match the observational accuracy , for a given model . only in this way 
   do the frequencies provide a faithful representation of the properties of the model , in comparisons with the observations 
   . however , since the numerical errors in the computed frequencies are typically highly systematic , they may affect the asteroseismic inferences even if they are smaller than the random errors in the observations , and hence more stringent requirements should be imposed on the computations . 
   also , the fact that solar - like oscillations , and several other types of asteroseismically interesting modes , tend to be of high radial order complicates reaching the required precision . 
   the numerical techniques discussed so far are generally of second order . 
   this yields insufficient precision in the evaluation of the eigenfrequencies , unless a very dense mesh is used in the computation ( see also * ? ? ? 
   the code may apply two techniques to improve the precision . 
   one technique ( cf . 
   * ) uses the fact that the frequency approximately satisfies a variational principle @xcite . 
   the variational expression may formally be written as @xmath119 where @xmath120 and @xmath121 are integrals over the equilibrium model depending on the eigenfunction , here represented by @xmath122 . 
   the variational property implies that any error @xmath123 in @xmath122 induces an error in @xmath124 that is @xmath125 . 
   thus by substituting the computed eigenfunction into the variational expression a more precise determination of @xmath78 should result . 
   this has indeed been confirmed @xcite . 
   the second technique uses explicitly that the difference scheme [ eq : e3.11 ] , which is used by one version of the shooting technique , and the relaxation technique , is of second order . consequently the truncation errors in the eigenfrequency and eigenfunction scale as @xmath126 . 
   if @xmath127 and @xmath128 are the eigenfrequencies obtained from solutions with @xmath129 and @xmath10 meshpoints , the leading - order error term therefore cancels in @xmath130 \ ; . \label{eq : e3.18}\ ] ] this procedure , known as _ richardson extrapolation _ , was used by @xcite . 
   it provides an estimate of the eigenfrequency that is substantially more accurate than @xmath131 , although of course at some added computational expense . 
   indeed , since the error in the representation [ eq : e3.11 ] depends only on even powers of @xmath132 , the leading term of the error in @xmath133 is @xmath134 .    even with these techniques 
   the precision of the computed frequencies may be inadequate if the mesh used in stellar - evolution calculations is used also for the computation of the oscillations . 
   the number of meshpoints is typically relatively modest and the distribution may not reflect the requirement to resolve properly the eigenfunctions of the modes . 
   @xcite discussed techniques to redistribute the mesh in a way that takes into account the asymptotic behaviour of the eigenfunctions ; a code to do so , based on four - point lagrangian interpolation , is included in the adipls distribution package . 
   on the other hand , for computing low - order modes ( as are typically relevant for , say , @xmath135 scuti or @xmath136 cephei stars ) , the original mesh of the evolution calculation may be adequate . 
   it is difficult to provide general recommendations concerning the required number of points or the need for redistribution , since this depends strongly on the types of modes and the properties of the stellar model . 
   it is recommended to carry out experiments varying the number and distribution of points to obtain estimates of the intrinsic precision of the computation ( e.g. , * ? ? ? 
   * ; * ? ? ? 
   in the latter case , considering simple polytropic models , it was found that 4801 points yielded a relative precision substantially better than @xmath137 for high - order p modes , when richardson extrapolation was used .    in the discussion of the frequency calculation it is important to distinguish between _ precision _ and _ accuracy _ , 
   the latter obviously referring to the extent to which the computed frequencies represent what might be considered the ` true ' frequencies of the model . 
   in particular , the manipulations required to derive eq .  [ eq : varprinc ] and to demonstrate its variational property depend on the equation of hydrostatic support being satisfied . 
   if this is not the case , as might well happen in an insufficiently careful stellar model calculation , the value determined from the variational principle may be quite precise , in the sense of numerically stable , but still unacceptably far from the correct value . 
   indeed , a comparison between @xmath138 and @xmath133 provides some measure of the reliability of the computed frequencies ( e.g. * ? ? ? 
   the programme finds the order of the mode according to the definition developed by @xcite and @xcite , based on earlier work by @xcite . 
   specifically , the order is defined by @xmath139 here the sum is over the zeros @xmath140 in @xmath59 ( excluding the centre ) , and @xmath141 is the sign function , @xmath142 if @xmath143 and @xmath144 if @xmath145 . for a complete model that includes the centre @xmath146 for radial oscillations and @xmath147 for nonradial oscillations . 
   thus the lowest - order radial oscillation has order @xmath148 . 
   although this is contrary to the commonly used convention of assigning order 0 to the fundamental radial oscillation , the convention used here is in fact the more reasonable , in the sense that it ensures that @xmath149 is invariant under a continuous variation of @xmath12 from 0 to 1 . with this definition @xmath150 for p modes , @xmath151 for f modes , and @xmath152 for g modes , at least in simple models . 
   it has been found that this procedure has serious problems for dipolar modes in centrally condensed models ( e.g. , * ? ? ? 
   * ; * ? ? ? 
   * ; * ? ? ? 
   the eigenfunctions @xmath59 are shifted such that nodes disappear or otherwise provide spurious results when eq . 
   [ eq : e4.1 ] is used to determine the mode order . a procedure that does not suffer from this difficulty has recently been developed by @xcite 
   ; i discuss it further in section  [ sec : develop ] .    a powerful measure of the characteristics of a mode is provided by the _ normalized inertia_. the code calculates this as @xmath153   \rho   r^2   { { \rm d}}r \over    m [ \xi_r ( r_{\rm phot } ) ^2 + l(l+1 ) { \xi_{\rm h}}(r_{\rm phot } ) ^2 ] } \nonumber \\    & = & { \int_{x_1}^{{x_{\rm s } } } \left [ y_1 ^ 2    + y_2 ^ 2   / l ( l + 1 ) \right ]   q u { { \rm d}}x / x    \over 4 \pi [ y_1 ( x_{\rm phot } ) ^2 + y_2 ( x_{\rm phot } ) ^2/l(l+1 ) ] } \ ; .\end{aligned}\ ] ] ( for radial modes the terms in @xmath60 are not included . ) here @xmath154 and @xmath155 are the distance of the innermost mesh point from the centre and the surface radius , respectively , and @xmath156 is the fractional photospheric radius . 
   the normalization at the photosphere is to some extent arbitrary , of course , but reflects the fact that many radial - velocity observations use lines formed relatively deep in the atmosphere . 
   a more common definition of the inertia is @xmath157 where @xmath158 is the so - called _ 
   mode mass_.    the code has the option to output the eigenfunctions , in the form of @xmath159 . 
   in addition ( or instead ) the displacement eigenfunctions can be output in a form indicating the region where the mode predominantly resides , in an energetical sense , as @xmath160 ( for radial modes only @xmath161 is found ) . 
   these are defined in such a way that @xmath162   { { \rm d}}x / x      \over 4 \pi [ y_1 ( x_{\rm phot } ) ^2 + y_2 ( x_{\rm phot } ) ^2/l(l+1 ) ] } \ ; . \label{eq : e4.4}\ ] ] the form provided by the @xmath163 is also convenient , e.g. , for computing rotational splittings @xmath164 ( e.g. , * ? ? ? * ) , where @xmath165 is the frequency of a mode of radial order @xmath149 , degree @xmath12 and azimuthal order @xmath2 . for slow rotation 
   the splittings are obtained from first - order perturbation analysis as @xmath166 characterized by _ kernels _ 
   @xmath167 , where in general the angular velocity @xmath168 depends on both @xmath1 and @xmath26 . 
   the code has built in the option to compute kernels for first - order rotational splitting in the special case where @xmath168 depends only on @xmath1 . 
   several revisions of the code have been implemented in preliminary form or are under development . 
   a substantial improvement in the numerical solution of the oscillation equations , particularly for high - order modes , is the installation of a fourth - order integration scheme , based on the algorithm of @xcite . 
   this is essentially operational but has so far not been carefully tested . 
   comparisons with the results of the variational expression and the use of richardson extrapolation , of the same formal order , will be particularly interesting .    as discussed by @xcite the use of @xmath169 ( or , as here , @xmath170 ) as one of the integration variables has the disadvantage that the quantity @xmath171 enters into the oscillation equations . 
   in models with a density discontinuity , such as results if the model has a growing convective core and diffusion is neglected , @xmath171 has a delta - function singularity at the point of the discontinuity . in the adipls calculations this 
   is dealt with by replacing the discontinuity by a very steep and well - resolved slope . 
   however , it would obviously be an advantage to avoid this problem altogether . this can be achieved by using instead the lagrangian pressure perturbation @xmath172 as one of the variables . 
   implementing this option would be a relatively straightforward modification to the code and is under consideration . 
   the proper classification of dipolar modes of low order in centrally condensed models has been a long - standing problem in the theory of stellar pulsations , as discussed in section  [ sec : results ] . 
   such a scheme must provide a unique order for each mode , such that the order is invariant under continuous changes of the equilibrium model , e.g. , as a result of stellar evolution . as a major breakthrough , takata in a series of papers has elucidated important properties of these modes and defined a new classification scheme satisfying this requirement @xcite . 
   a preliminary version of this scheme has been implemented and tested ; however , the latest and most convenient form of the takata classification still needs to be installed . 
   a version of the code has been established which computes the first - order rotational splitting for a given rotation profile @xmath173 , in addition to setting up the corresponding kernels . 
   this is being extended by k. burke , sheffield , to cover also second - order effects of rotation , based on the formalism of @xcite . 
   an important motivation for this is the integration , discussed by @xcite , of the pulsation calculation with the astec evolution code to allow full calculation of oscillation frequencies for a model of specified parameters ( mass , age , initial rotation rate , etc . ) as the result of a single subroutine call . 
   i am very grateful to w. dziembowski and d.  o.  gough for illuminating discussions of the properties of stellar oscillations , and to a. moya and m. j. p. f. g. monteiro for organizing the comparisons of stellar oscillation and model calculations within the esta collaboration . 
   i thank the referee for useful comments which , i hope , have helped improving the presentation . 
   this project is being supported by the danish natural science research council and by the european helio- and asteroseismology network ( helas ) , a major international collaboration funded by the european commission s sixth framework programme . 
   christensen - dalsgaard , j. , berthomieu , g. : theory of solar oscillations . in : cox , a. n. , livingston , w. c. , matthews , m. ( eds ) , _ solar interior and atmosphere _ , p. 401 
   space science series , university of arizona press ( 1991 )                          moya , a. , christensen - dalsgaard , j. , charpinet , s. , lebreton , y. , miglio , a. , montalbn , j. , monteiro , m. j. p. f. g. , provost , j. , roxburgh , i. , scuflaire , r. , surez , j. c. , suran , m. : inter - comparison of the g- , f- and p - modes calculated using different oscillation codes for a given stellar model . 
   apss , this volume ( 2007 )
summary: development of the aarhus adiabatic pulsation code started around 1978 . 
   although the main features have been stable for more than a decade , development of the code is continuing , concerning numerical properties and output . 
   the code has been provided as a generally available package and has seen substantial use at a number of installations . 
   further development of the package , including bringing the documentation closer to being up to date , is planned as part of the helas coordination action .