We present the R package stellaR, which is designed to access and manipulate publicly available stellar evolutionary tracks and isochrones from the Pisa low-mass database. The procedures for extracting important stages in the evolution of a star from the database, for constructing isochrones from stellar tracks and for interpolating among tracks are discussed and demonstrated.
Due to the advance in the instrumentation, nowadays astronomers can deal
with a huge amount of high-quality observational data. In the last
decade impressive improvements of spectroscopic and photometric
observational capabilities made available data which stimulated the
research in the globular clusters field. The theoretical effort of
recovering the evolutionary history of the clusters benefits from the
computation of extensive databases of stellar tracks and isochrones,
such as Pietrinferni et al. (2006; Bertelli et al. 2008; Dotter et al. 2008). We recently computed a large
data set of stellar tracks and isochrones, “The Pisa low-mass database”
(Dell’Omodarme et al. 2012), with up to date physical and chemical inputs, and made
available all the calculations to the astrophysical community at the
Centre de Données astronomiques de Strasbourg (CDS)
In most databases, the management of the information and the extraction of the relevant evolutionary properties from libraries of tracks and/or isochrones is the responsibility of the end users. Due to its extensive capabilities of data manipulation and analysis, however, R is an ideal choice for these tasks. Nevertheless R is not yet well known in astrophysics; up to December 2012 only seven astronomical or astrophysical-oriented packages have been published on CRAN (see the CRAN Task View Chemometrics and Computational Physics).
The package stellaR (Dell’Omodarme and Valle 2012) is an effort to make available to the astrophysical community a basic tool set with the following capabilities: retrieve the required calculations from CDS; plot the information in a suitable form; construct by interpolation tracks or isochrones of compositions different to the ones available in the database; construct isochrones for age not included in the database; extract relevant evolutionary points from tracks or isochrones.
The Pisa low-mass database contains computations classified according to
four parameters: the metallicity showComposition()
:
> showComposition()
Mixing-length values:
1.7, 1.8, 1.9
alpha-enhancement values:
0, 1 (i.e. [alpha/Fe] = 0.0 [alpha/Fe] = 0.3)
Chemical compositions:
z y.1 y.2 y.3 y.4 y.5 y.6
1e-04 0.249 0.25 0.27 0.33 0.38 0.42
2e-04 0.249 0.25 0.27 0.33 0.38 0.42
3e-04 0.249 0.25 0.27 0.33 0.38 0.42
4e-04 0.249 0.25 0.27 0.33 0.38 0.42
5e-04 0.250 0.25 0.27 0.33 0.38 0.42
6e-04 0.250 0.25 0.27 0.33 0.38 0.42
7e-04 0.250 0.25 0.27 0.33 0.38 0.42
8e-04 0.250 0.25 0.27 0.33 0.38 0.42
9e-04 0.250 0.25 0.27 0.33 0.38 0.42
1e-03 0.250 0.25 0.27 0.33 0.38 0.42
2e-03 0.252 0.25 0.27 0.33 0.38 0.42
3e-03 0.254 0.25 0.27 0.33 0.38 0.42
4e-03 0.256 0.25 0.27 0.33 0.38 0.42
5e-03 0.258 0.25 0.27 0.33 0.38 0.42
6e-03 0.260 0.25 0.27 0.33 0.38 0.42
7e-03 0.262 0.25 0.27 0.33 0.38 0.42
8e-03 0.264 0.25 0.27 0.33 0.38 0.42
9e-03 0.266 0.25 0.27 0.33 0.38 0.42
1e-02 0.268 0.25 0.27 0.33 0.38 0.42
The table of chemical compositions presents all the
Upon specification of the aforementioned parameters, the stellaR
package can import data from CDS (via anonymous ftp) over an active
Internet connection. The CDS data are stored in ASCII format and include
a header with calculation metadata, such as the metallicity, the initial
helium abundance, and the mixing-length. The import is done via a
read.table()
call, skipping the header of the files.
The following data objects can be downloaded from the database site:
Stellar track: a stellar evolutionary track computed starting from
Pre-Main Sequence (PMS) and ending at the onset of helium flash (for
masses getTrk()
and
getTrkSet()
can be used to access such data; they respectively
return objects of classes "trk"
and "trkset"
.
Stellar ZAHB: Zero-Age Horizontal-Branch models. The function
getZahb()
can be used to access such data; it returns an object of
class "zahb"
.
HB models: computed from ZAHB to the onset of thermal pulses. The
functions getHb()
and getHbgrid()
can be used to access such
data; they respectively return objects of classes "hb"
and
"hbgrid"
.
Stellar isochrones: computed in the age range [8.0 - 15.0] Gyr.
The functions getIso()
and getIsoSet()
can be used to access
such data; they respectively return objects of classes "iso"
and
"isoset"
.
Readers interested in details about the computation procedure are
referred to Dell’Omodarme et al. (2012). The data gathered from CDS are organized into
objects of appropriate classes. The package includes print
and plot
S3 methods for the classes "trk"
, "trkset"
, "zahb"
, "hb"
,
"hbgrid"
, "iso"
, and "isoset"
.
As an example, we illustrate the recovering of the stellar track for a
model of mass
> track <- getTrk(m = 0.80, z = 0.001, y = 0.25, ml = 1.90, afe = 0)
> track
Stellar track
Mass = 0.8 Msun
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
> names(track)
[1] "mass" "z" "y" "ml" "alpha.enh" "data"
> class(track)
[1] "trk" "stellar"
The function getTrk()
returns an object of class "trk"
, which is a
list containing the track metadata, i.e. the star mass, the metallicity,
the initial helium abundance, the mixing-length and the
data
.
Track data contains the values of 15 variables:
> names(track$data)
[1] "mod" "time" "logL" "logTe" "mass" "Hc" "logTc" "logRHOc"
[9] "MHEc" "Lpp" "LCNO" "L3a" "Lg" "radius" "logg"
The included variables are: mod
the progressive model number; time
the logarithm of the stellar age (in yr); logL
the logarithm of the
surface luminosity (in units of solar luminosity); logTe
the logarithm
of the effective temperature (in K); mass
the stellar mass (in units
of solar mass); Hc
the central hydrogen abundance (after hydrogen
exhaustion: central helium abundance); logTc
the logarithm of the
central temperature (in K); logRHOc
the logarithm of the central
density (in g/cmMHEc
the mass of the helium core (in units of
solar mass); Lpp
the luminosity of pp chain (in units of surface
luminosity); LCNO
the luminosity of CNO chain (in units of surface
luminosity); L3a
the luminosity of triple-Lg
luminosity of the gravitational energy (in
units of surface luminosity); radius
the stellar radius (in units of
solar radius); logg
the logarithm of surface gravity (in cm/s
Similarly the part of the track starting from ZAHB and ending at the onset of thermal pulses can be downloaded with the call:
> hbtk <- getHb(m = 0.80, z = 0.001, y = 0.25, ml = 1.90, afe = 0)
> hbtk
Stellar track from ZAHB
Mass = 0.8 Msun
Mass RGB = 0.8 Msun
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
> names(hbtk)
[1] "mass" "massRGB" "z" "y" "ml" "alpha.enh"
[7] "data"
> class(hbtk)
[1] "hb" "stellar"
Function getHb()
returns an object of class "hb"
, which differs from
an object of class "trk"
only for the presence of the variable
massRGB
, i.e. the Red-Giant Branch (RGB) progenitor mass.
Usually a set of tracks with different mass and/or metallicity values
are needed for computations. The package stellaR provides the function
getTrkSet()
, which can download a set of tracks with different values
for mass, metallicity, initial helium abundance, mixing-length and
> mass <- seq(0.3, 1.1, by = 0.05)
> trks <- getTrkSet(m = mass, z = 0.001, y = 0.25, ml = 1.90, afe = 0)
> trks
[[1]]
Stellar track
Mass = 0.3 Msun
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
[[2]]
Stellar track
Mass = 0.35 Msun
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
...
The function getTrkSet()
returns an object of class "trkset"
, a list
containing objects of class "trk"
. The track set can be displayed in
the usual (plot()
:
> plot(trks, lty = 1:2)
The output of the function is shown in Figure 1. The plot
is produced by a call to the function plotAstro()
, which allows the
user to customize several aspects of the plot, such as the axes labels,
the number of minor ticks between two major ticks, the limits of the
axes, the color and type of the lines (as in the example), the type of
the plot (lines, points, both, …).
The use of the plot()
function is further demonstrated in Figure
2, where, for
> plot(track)
> plot(hbtk, add = TRUE, col = "green")
Apart from the plots discussed before, it is easy to display other
relations between the computed variables. In the following example we
get the data for two masses, namely 0.50 and 1.00
> trkr <- getTrkSet(m = c(0.5, 1), z = 0.01, y = 0.25, ml = 1.8, afe = 0)
> mydata <- do.call(rbind, lapply(trkr, "[[", "data"))
> D <- subset(mydata, mod <= 100)
> key <- as.numeric(factor(D$mass))
> plotAstro(D$time, D$radius, type = "p", pch = key, ylab = "Radius (Rsun)",
+ xlab = "log age (yr)")
> legend("topright", c("M=0.50", "M=1.00"), pch = 1:2)
Isochrones can be obtained from CDS and plotted in a similar way. As an
example, we get isochrones of 9 and 12 Gyr for
> isc <- getIsoSet(age = c(9, 12), z = 0.001, y = 0.25, ml = 1.90, afe = 0)
> isc
[[1]]
Stellar isochrone
Age = 9 Gyr
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
[[2]]
Stellar isochrone
Age = 12 Gyr
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
attr(,"class")
[1] "isoset" "stellar"
> names(isc[[1]])
[1] "age" "z" "y" "ml" "alpha.enh" "data"
> names(isc[[1]]$data)
[1] "logL" "logTe" "mass" "radius" "logg"
The function returns an object of class "isoset"
, a list containing
objects of class "iso"
. The latter objects are lists containing
metadata (age, metallicity, initial helium abundance, mixing-length,
data
, which contains the
computed theoretical isochrones data. Figure 4 shows the set
of isochrones plotted with the commands:
> plot(isc, lty = 1:2)
> legend("topleft", c("9 Gyr", "12 Gyr"), lty = 1:2)
Even if the database provides a fine grid of models it is possible that
the specific model needed for a comparison with observational data has
not been calculated. To address this issue, the package stellaR
includes tools for constructing new sets of isochrones. The simplest
case is whenever one desires an isochrone for ages not included in the
database, but for a combination of metallicity, initial helium
abundance, mixing-length and makeIso()
is able to compute the required
isochrones by means of interpolation on a set of tracks. The
interpolation can be performed for age in the range [7.0 - 15.0] Gyr.
The user has the choice to explicitly provide to the function a set of
tracks, previously downloaded from CDS, or to specify the required
composition of the tracks to be downloaded for the interpolation. To
show the usage of the function we use the object trks
downloaded
before to obtain an isochrone of age 9.7 Gyr:
> iso.ip <- makeIso(age = 9.7, tr = trks)
> iso.ip
Stellar isochrone
Age = 9.7 Gyr
Z = 0.001 , Y = 0.25
Mixing length = 1.9
[alpha/Fe] = 0
The call produces a result identical to
makeIso(age = 9.7, z = 0.001, y = 0.25, ml = 1.9, afe = 0)
; in the
latter case the data are taken from CDS before the interpolation
procedure.
The interpolation technique is based upon the fact that all the tracks
contain the same number of points by construction, and that a given
point corresponds to the same evolutionary phase on all the tracks. We
define
The package stellaR provides also a tool for performing a 3D
interpolation on the database to construct a set of tracks for values of
metallicity, initial helium abundance and mixing-length not included in
the computations available at CDS. The function interpTrk()
can be
used for this procedure. A call to this function causes the download
from CDS of the sets of tracks needed for the interpolation.
The new set of tracks is computed by means of a linear interpolation.
The metallicity is log-transformed before the interpolation procedure.
Let
As a demonstration, let us compute a set of tracks with mixing-length
value
> ip.trk <- interpTrk(z = 0.002, y = 0.25, ml = 1.74, afe = 0)
Since the values of
> ip.iso <- makeIso(age = 12, tr = ip.trk)
Important stages in the evolution of a star are defined as
"keypoints", e.g. hydrogen core exhaustion, Turn-Off luminosity, RGB
tip luminosity. To simplify their extraction the package stellaR
provides the function keypoints()
, which operates on an object of
class "trk"
or "iso"
.
The function extracts from the data stored in objects of class "trk"
the rows of the data frame relative to the following evolutionary
stages:
ZAMS. : Zero-Age Main-Sequence, defined as the point for which the central H abundance drops below 99% of its initial value.
TO. : Turn-Off, defined as the point for which the effective temperature reaches its maximum value. If multiple lines satisfy the constraint, the values of all the rows are averaged.
BTO. : Brighter Turn-Off, defined as the point for which the effective temperature drops below the temperature of the TO minus 100 K. This point can not exist for low masses. Details on the advantages of this evolutionary point with respect to the TO can be found in Chaboyer et al. (1996) .
exHc. : Central H exhaustion, defined as the point for which the central H
abundance is zero. For low masses the point can coincide with TO.
This is the last point of the tracks with mass lower or equal to
0.50
Heflash. : Helium flash, the last point of the track for masses higher than
0.50
When the function is called on an object of class "iso"
it returns a
data frame containing only TO and BTO phases.
In both cases the function inserts in the returned data frame the
columns relative to mass (or age for isochrones), metallicity, initial
helium abundance value, mixing-length,
As a demonstration we extract the TO and BTO points from the object
isc
generated in a previous example:
> kp <- keypoints(isc)
> kp
logL logTe mass radius logg age z y ml alpha.enh id
TO 0.4044723 3.816652 0.8396903 1.23558 4.178911 9 0.001 0.25 1.9 0 1
BTO 0.5556971 3.809943 0.8561162 1.51671 4.009265 9 0.001 0.25 1.9 0 2
TO1 0.2917250 3.803611 0.7772973 1.15234 4.205965 12 0.001 0.25 1.9 0 1
BTO1 0.4495597 3.796545 0.7909500 1.42768 4.027426 12 0.001 0.25 1.9 0 2
The points can be easily superimposed to the isochrones in the theoretical plane. The top panel of Figure 5 is obtained with the following commands:
> plot(isc)
> points(logL ~ logTe, data = kp, pch = id, col = id)
> legend("topleft", c("TO", "BTO"), pch = 1:2, col = 1:2)
As a last example we extract a set of tracks for masses in the range
[0.40 - 1.10]
> mT <- seq(0.4, 1.1, by = 0.1)
> zT <- c(0.0001, 0.001, 0.01)
> tr <- getTrkSet(m = mT, z = zT, y = 0.25, ml = 1.9, afe = 1)
> kp <- keypoints(tr)
> kpH <- subset(kp, id == 4)
> symbol <- as.numeric(factor(kpH$z))
> plotAstro(kpH$M, kpH$time, type = "p", pch = symbol, xi = 0.1,
+ xlab = expression(italic(M) ~ (italic(M)[sun])), ylab = "log age (yr)")
> lab <- format(sort(unique(kpH$z)), nsmall = 4, scientific = FALSE)
> legend("topright", lab, pch = 1:3)
This paper demonstrated how the package stellaR can be useful to the astrophysical community as a tool to simplify the access to stellar tracks and isochrone calculations available on-line. A set of tools to access, manipulate and plot data are included in the package and their usage was shown. The interpolation functions included in the package can be used to safely produce tracks or isochrones at compositions not included in the database, without the need that users develop software on their own. A planned extension of the package is the modification of the algorithm of isochrone construction to make the calculation of isochrones of young ages feasible. This step can be useful in view of a possible extension of the Pisa database to higher masses, or to manipulation of data stored in other databases. In fact, while the package is currently developed for accessing data from the Pisa low-mass database, other public databases can be in principle accessed in the same way. This step is however complicated by the fact that no standard for the stellar model output exists in the astrophysical community, requiring the individual adaptation of the interface functions for each data set.
We are grateful to our anonymous referees for many stimulating suggestions that helped to clarify and improve the paper and the package. We thank Steve Shore for a careful reading of the manuscript and many useful suggestions.
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Dell'Omodarme & Valle, "stellaR: A Package to Manage Stellar Evolution Tracks and Isochrones", The R Journal, 2013
BibTeX citation
@article{RJ-2013-011, author = {Dell'Omodarme, Matteo and Valle, Giada}, title = {stellaR: A Package to Manage Stellar Evolution Tracks and Isochrones}, journal = {The R Journal}, year = {2013}, note = {https://rjournal.github.io/}, volume = {5}, issue = {1}, issn = {2073-4859}, pages = {108-116} }