astronomR: Astronomy and Cosmology Analysis in R 🌌✨

An R package designed to bridge the gap between data science and the cosmos!

While R is widely used in data science and statistical computing, there has been a lack of tools specifically tailored for astronomical data. {astronomR} aims to fill this gap, offering easy-to-use functions and tools for anyone looking to explore space and cosmology through data. Whether you're an astrophysicist or just a data enthusiast, {astronomR} is here to help you navigate the universe.

Installing this Package

Install the stable release from CRAN:

install.packages("astronomR")

Or install the development version from GitHub:

# Install devtools package if you don't have it
# install.packages("devtools")
devtools::install_github("samrit2442/astronomR")

Introduction

This package was developed for astronomy, cosmological computation, and analysis with R. The source code can be found here https://github.com/samrit2442/astronomR.

Function Overview

Category	Functions
Angular Conversions	`deg_to_hms()`, `hms_to_deg()`, `deg_to_dms()`, `dms_to_deg()`, `deg2rad()`, `rad2deg()`
Physical Constants	`constants_df`, `constant_value()`
Gaia Archive	`get_gaia_data()`
Cosmology	`km_to_Mpc()`, `Mpc_to_km()`, `cosmology_model()`, `age_of_universe()`, `comoving_distance()`, `luminosity_distance()`, `angular_diameter_distance()`
Thermal Physics	`photon_energy_density_fn_T()`, `photon_energy_density_fn_z()`, `photon_number_density_fn_T()`, `photon_number_density_fn_z()`, `Saha_Xe()`, `soln_saha()`
Thermal Cosmology	`hubble_radiation()`, `g_star_eff()`, `entropy_density()`, `equilibrium_number_density()`, `equilibrium_yield()`, `boltzmann_pebble_rhs()`, `solve_relic_abundance()`, `freeze_out_xf()`, `peebles_rhs()`
Drake Equation	`drake_equation()`

Note: the worked examples below also call RA_dec2Alt_azi(), FlatLCDM(), t_as_func_a_in_Year(), and radius_of_curvature(). These aren't in the table above yet -- worth adding (with the correct category) so the overview stays in sync with the examples.

Some Usage and Example Codes

Normally in astronomy, we use a different sort of angular system for location. We can easily do that in our package. Let's see how. Suppose, we have an angular value of d = 177.74208° We want to convert it into an hour-minute-second. There is a very simple function.

library(astronomR)
deg_to_hms(177.74208)
#> [1] 11H50M58.0992S
hms_to_deg(11, 50, 58.09925)
#> [1] 177.7421

What else can be done? Let's say we want to find the path of some star in a particular location for a time interval. This can also be done using our package.

To do that, first, we need to define the RA and Dec value of the star. Let's see how to do this. Also, let's define the time of observation and its location.

ra_hour  <- 16.695      # RA in hours
dec_deg  <- 36.466667   # Dec in degrees
lat_obs  <- 52.5        # Observer's latitude
lon_obs  <- -1.9166667  # Observer's longitude
datetime <- as.POSIXct("1998-08-10 23:10:00", tz = "UTC")  # Observation time

Now, create the star location as seen from that location at the time mentioned. This returns Altitude and Azimuth as normally used.

star_location <- RA_dec2Alt_azi(ra_hour, dec_deg, lat_obs, lon_obs, datetime)
print(paste("Altitude:", star_location$altitude, "degrees"))
#> [1] "Altitude: 49.1688687197424 degrees"
print(paste("Azimuth:", star_location$azimuth, "degrees"))
#> [1] "Azimuth: 269.146669462321 degrees"

This tells us from some location, what is the position of the star so that we can use a telescope to watch it! Using this simple function, we can trace out the path any star travels. Let's see how and also maybe plot it. For that, let's first make a time range for which we want to see the location.

library(ggplot2)
# Plotting Rigel's Motion
start_datetime <- as.POSIXct("2024-10-02 00:00:00", tz = "UTC")
end_datetime   <- as.POSIXct("2024-10-03 00:00:00", tz = "UTC")
# Generate timestamps at 20-minute intervals
timestamps <- seq(from = start_datetime, to = end_datetime, by = "20 mins")
altitude   <- numeric(length(timestamps))
azimuth    <- numeric(length(timestamps))
observer_lat <- 43.1566   # Latitude in degrees
observer_lon <- -77.6088  # Longitude in degrees
# RA and Dec for Rigel (converted RA to hours)
rigel_ra  <- 78.634467 / 15  # RA in hours
rigel_dec <- -8.20164         # Dec in degrees
# Calculate altitude and azimuth for each timestamp
for (i in seq_along(timestamps)) {
  datetime    <- timestamps[i]
  result      <- RA_dec2Alt_azi(rigel_ra, rigel_dec, observer_lat, observer_lon, datetime)
  altitude[i] <- result$altitude
  azimuth[i]  <- result$azimuth
}
# Create a data frame with results
rigel_positions <- data.frame(datetime = timestamps, altitude = altitude, azimuth = azimuth)
# 2D Plot for Altitude over Time
ggplot(rigel_positions, aes(x = datetime, y = altitude)) +
  geom_line() +
  labs(x = "Time", y = "Altitude (degrees)", title = "Altitude of Rigel over 24 hours")

# 2D Plot for Azimuth over Time
ggplot(rigel_positions, aes(x = datetime, y = azimuth)) +
  geom_line() +
  labs(x = "Time", y = "Azimuth (degrees)", title = "Azimuth of Rigel over 24 hours")

# 2D Scatter Plot for Azimuth vs Altitude
ggplot(rigel_positions, aes(x = azimuth, y = altitude)) +
  geom_point(color = "firebrick") +
  labs(x = "Azimuth (degrees)", y = "Altitude (degrees)", title = "Azimuth vs Altitude for Rigel over 24 hours")

The three plots above (altitude vs. time, azimuth vs. time, azimuth vs. altitude) are embedded via GitHub's pasted-image links, which are signed and can expire. If they ever render as broken, re-upload the PNGs (e.g., commit them under man/figures/) and swap in the permanent URLs.

Gaia Data Archive

The Gaia Data Archive is a comprehensive database that houses the data collected by the European Space Agency's Gaia mission. Launched in December 2013, Gaia is designed to create the most accurate three-dimensional map of the Milky Way galaxy by observing and cataloguing the positions, distances, and motions of over a billion stars. In Python, Astropy and Astroquery are used to import data directly. We can do the same thing in our package. Let's see how:

df <- get_gaia_data(vars = "ra, dec, parallax", condition = "parallax > 50")
head(df)
#>         ra       dec  parallax
#> 1 316.7537  38.75607 286.00534
#> 2 316.7485  38.76386 285.99493
#> 3 298.4819  44.41291 214.57451
#> 4 249.3875 -53.69952  53.22840
#> 5 243.4453 -57.57679  73.54824
#> 6 312.2788  37.47123  56.86465

Nice, isn't it? Let's use this for some analysis. Why not create an H-R diagram?

For this, first let's convert parallax to absolute magnitude.

library(ggplot2)
library(dplyr)
# Helper function to convert parallax (mas) and apparent magnitude to absolute magnitude
calculate_absolute_magnitude <- function(parallax, g_mag) {
  distance_pc <- 1 / (parallax / 1000)
  abs_mag     <- g_mag - 5 * (log10(distance_pc) - 1)
  return(abs_mag)
}
# Let's import the data
df <- get_gaia_data(vars = "ra, dec, parallax, phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag",
                     condition = "parallax > 50")
# Data processing
df <- df %>%
  mutate(color_index = phot_bp_mean_mag - phot_rp_mean_mag,
         abs_mag = calculate_absolute_magnitude(parallax, phot_g_mean_mag))
# Filter out any rows with NA values
df <- df %>% filter(!is.na(color_index), !is.na(abs_mag))
# Plot an H-R diagram
ggplot(df, aes(x = color_index, y = abs_mag)) +
  geom_point(alpha = 0.5, color = "blue") +
  labs(x = "Color Index (G_BP - G_RP)", y = "Absolute Magnitude (M)",
       title = "Hertzsprung-Russell Diagram") +
  theme_minimal() +
  scale_y_reverse()  # Reverse the y-axis for magnitude

Cosmological Calculation

We can do many things related to cosmology using this package. Let's see a few of them. Before that, it's worth running this code. As we know, there are many cosmological models, and you can define these in our package.

cosmo <- FlatLCDM(0.6774, 0.6911, 0.3089, 8.4e-5)
cosmo
#> $hubble_constant_fact
#> [1] 0.6774
#>
#> $dark_matter_crit
#> [1] 0.6911
#>
#> $matter_crit
#> [1] 0.3089
#>
#> $radiation_crit
#> [1] 8.4e-05
#>
#> $type
#> [1] "FlatLCDM"
#>
#> $h_per_s
#> [1] 2.194948e-20

Now we can directly use this model as input to other functions. As we know, the scale factor of our universe a = 1 today, we can directly use our function to calculate time at any given a value.

t_as_func_a_in_Year(cosmo, 1) # in year
#> [1] 13808979943

We can find the age of our universe using this, or use the function age_of_universe().

age_of_universe(cosmo, unit = "GY") # our universe age !!!!!
#> [1] 13.80898

Our package can also find the radius of curvature of our universe. Let's see:

radius_of_curvature(cosmo)
#> [1] 0

Thermal Cosmology of the Early Universe

The thermal cosmology module covers the complete early-universe thermal history, from the WIMP freeze-out epoch all the way through hydrogen recombination. All functions in this section use natural units (ħ = c = k_B = 1) with temperature and mass in GeV, except peebles_rhs(), which uses SI units.

Hubble Rate & Degrees of Freedom

# Hubble rate H(T) in the radiation-dominated era [GeV]
hubble_radiation(100)                    # T = 100 GeV, full SM  g* = 106.75
#> [1] 1.405077e-14
hubble_radiation(1e-3, g_star = 10.75)   # T = 1 MeV, neutrino era
#> [1] 4.458986e-25
# Effective relativistic degrees of freedom g*(T)
g_star_eff(500)    # > 300 GeV: full Standard Model
#> [1] 106.75
g_star_eff(0.1)    # neutrino era
#> [1] 10.75
g_star_eff(0.01)   # after e+e- annihilation
#> [1] 3.91

Entropy Density & Equilibrium Distributions

# Entropy density s(T) = (2 pi^2 / 45) g*S T^3  [GeV^3]
entropy_density(0.1)                     # T = 100 MeV
#> [1] 0.004715535
# Maxwell-Boltzmann equilibrium number density n_eq(T)
equilibrium_number_density(T_GeV = 5, m_GeV = 100, g_dof = 2)
#> [1] 2.926338e-06
# Equilibrium yield Y_eq(x)  [x = m/T]
equilibrium_yield(x = 20, m_GeV = 100)   # near freeze-out
#> [1] 7.720313e-10

WIMP Freeze-Out: The Boltzmann Pebble Equation

The pebble equation governs the evolution of the comoving yield Y = n/s of a thermally produced dark-matter relic: