The purpose of this vignette is to allow a user to explore the range of possible outputs that can be produced by SimSpin for a tipsy formatted input file from the Gasoline code. Because the Gasoline code is sparsely documented, we provide this branch of the code (v2.8.5.1) as a seperate branch dev-tipsy that will not be merged with main.

In the example below, you will need to enter the file location of your own tipsy Gasoline output galaxy such that this can be run through the code. The rest of the functions should then run without modification. Set the location of your simulation input file below:

The purpose of the SimSpin R-package is to take a particle simulation of a galaxy and produce a spectral data cube in the style of a specified Integral Field Spectroscopy (IFS) instrument.

In this vignette, we will take you through the process of constructing such a data cube from an input tipsy file. This is a simple process comprised of three steps:

From this spectral data cube, observable properties can be measured using standard observational pipelines, such as pPXF. You can see another long-form example here demonstrating how the kinematics can be fit using the spectral data cube.

Step 0: Installing SimSpin (the compatible version)

The dev-tipsy branch of SimSpin is (currently) the only one that supports the reading of tipsy files. To install a specific branch of the package, we can use the following:

# If you don't already have SimSpin installed for tispy files, uncomment the lines below:
# -------------------------------------------------------- #
# library(devtools) 
# install_github("kateharborne/SimSpin", ref="dev-tipsy")  
# # this points towards the specific branch of the code
# -------------------------------------------------------- #

library(SimSpin)
#> Loading required package: foreach
packageVersion("SimSpin")
#> [1] '2.9.2.1'

Ensure that the package version of the loaded code is “v2.8.7.1” as this is the version that includes support for tipsy output files from the Gasoline code.

Step 1: Importing the Gasoline simulation files

To begin a SimSpin observation, we begin by converting the input simulation into a consistent format. This step creates an R-formatted file that is intended to be universal in structure (to aid in comparing between different simulations) and to allow efficient repeat observations of the same system due to the construction of the file. It is a pre-processing step that allows the user to smooth the gas distribution (in the scenario that the softening lengths of the gas is significantly larger than the spaxel size of the intended observation) and pre-computes the assignment of stellar particles to particular stellar templates.

For Gasoline outputs, you should have your outputs formatted in a directory that contains the output simulation in one named binary file in tipsy format, along with several additional extension files (e.g. required files include *.timeform, *.iord, *.massform, *.OxMassFrac). If you place the parameter file used to build the simulation in the same directory, this information will also be used to format the SimSpin object, else expected defaults for required parameters will be substituted and the user will be warned that this is the case.

If you do not have the parameter file in the output directory, you should see a warning issued when the following code chunk is run. Else, it should proceed without warnings raised.

Be aware that the speed at which this chunk will execute will scale with the size of the simulation. For a file containing ~9x10^6 particles (stellar and gas) will take ~10 min to run this initial make_simspin_file function, though later data cube construction steps should then be faster. I’d suggest going and getting a coffee!

tipsy_output = list.files(dirname(f), full.names = T)

if (any(stringr::str_detect(tipsy_output, ".param"))){
  print("Units as set in the parameter file included in this directory will be used for SimSpin.")
} else {
  print("No parameter file included in this directory. Assumed units (distance unit = 1 kpc and mass unit = 1 Msol and G = 1) will be used for SimSpin.")
}
#> [1] "Units as set in the parameter file included in this directory will be used for SimSpin."

### Making a SimSpin file using this input tipsy file ------------------------
tipsy_simspin_file = 
    SimSpin::make_simspin_file(filename = f,
                               template = "BC03lr",
                               write_to_file = F, 
                               cores = 1)

This SimSpin file can then be used to produce an observation. We can inspect the information included in this file as well:

names(tipsy_simspin_file) # this will tell you what groups are contained within the file
#> [1] "header"           "star_part"        "gas_part"         "spectral_weights"

tipsy_simspin_file$header # inspecting the header information of the 
#> $InputFile
#> [1] "~/Desktop/GLX.01000/GLX.01000"
#> 
#> $OutputFile
#> NULL
#> 
#> $Type
#> [1] "Tipsy"
#> 
#> $Template
#> [1] "BC03lr"
#> 
#> $Template_LSF
#> [1] 3
#> 
#> $Template_waveres
#> [1] 1
#> 
#> $Centre
#> [1] NA
#> 
#> $HalfMass
#> [1] NA
#> 
#> $TotalStellarMass
#> [1] 32803476925
#> 
#> $TotalGasMass
#> [1] 2544862528
#> 
#> $Alignment
#> [1] "Default"
#> 
#> $Npart
#> [1] 7060462
#> 
#> $SmoothingN
#> [1] 1
#> 
#> $Origin
#> [1] "SimSpin_v2.9.2.1"
#> 
#> $Date
#> [1] "2024-10-23 14:29:32 AWST"

To smooth the gas distribution in an SPH model, we would set the input parameter sph_spawn_n to a value greater than 1. This will spawn n particles for each gas particle scaled in their mass and sampling SPH kernel in random uniform position within the smoothing length of that gas particle. Again, depending on the comparative resolution of the simulation to the spatial resolution of the intended telescope, this smoothing may not be necessary for visualisation. Obviously, this process of spawning particles will increase the size of the SimSpin file produced.

Step 2: Making observations from SimSpin files

First, initialise a telescope and observing strategy:

MUSE = SimSpin::telescope(type="MUSE", fov = 15)
strategy = SimSpin::observing_strategy(dist_z = 0.05, inc_deg=70, blur=F)

Further examples of supported instruments, or how to define your own, can be found online: https://kateharborne.github.io/SimSpin/basic_usage/

An observation can then be made using build_datacube. This will be a stellar kinematic observation, with stellar properties summarised. By changing the method parameter, this can be modified to create an observation of the gas (or star-forming gas where SFR > 0).

observed_tipsy_cube = build_datacube(simspin_file = tipsy_simspin_file,
                                     telescope = MUSE,
                                     observing_strategy = strategy,
                                     method = "velocity", 
                                     write_fits = F, mass_flag = T)

This file will contain a number of images. To construct a spectral data cube (in which you will then need to run pPXF to derive observed kinematics and stellar population parameters), simply switch method = "spectral".

The same kind of observation can also be made for the gas component of the galaxy:

observed_tipsy_cube_gas = build_datacube(simspin_file = tipsy_simspin_file,
                                         telescope = MUSE,
                                         observing_strategy = strategy,
                                         method = "gas", 
                                         write_fits = F)

This file will contain another set of images for the gas component. We summarise some of the other assumptions made for the output Gasoline parameters in the section below and point the user to the lines in the SimSpin package where the specific values have been computed (for bug-fixing in the future)!

Step 3: Considering the results

For a stellar observation

In the read_files.R file, you will find the .read_tipsy file at line 134. This function contains the methods we use to convert the information contained within the Gasoline output files into the output summary images.


mask = observed_tipsy_cube$raw_images$particle_image
mask[mask<200] = NA
mask[!is.na(mask)] = 1

plot_flux(log10(observed_tipsy_cube$observed_images$flux_image*mask), fig = c(0,0.33, 0, 1), labN = 2, titleshift = -7, units = expression("log"[10]*"(Flux, CGS)"))
plot_velocity(observed_tipsy_cube$observed_images$velocity_image*mask, fig = c(0.33, 0.66, 0, 1), new=T, labN = 2, titleshift = -7)
plot_dispersion(observed_tipsy_cube$observed_images$dispersion_image*mask, fig = c(0.66, 0.99, 0, 1), new=T, labN = 2, titleshift = -7)

For computing fluxes, we have used: * stellar ages * stellar metallicities * stellar initial masses

These are used to map to the grid of stellar templates in the chosen template library (in this case, the Bruzal and Charlot 2003 templates).

We’ve assumed that the output stellar ages given in the .timeform file output can be converted to stellar ages in the following way:

# This chunk should not run in the example, but we copy the code used to compute the stellar ages here. This can be found on line 245 in the read_files.R file:

  if (file.exists(paste0(f,".timeform"))){

    age_data = file(paste0(f,".timeform"), "rb")
    stars_formed = readBin(age_data, "numeric", n = nstar, size = 4, endian = endian) # time since the start of the simulation, given in Myr
    close(age_data)

    stars_formed = stars_formed * t_Unit # formation time of stars in Gyrs 
    # t_Unit is determined from the parameter file:
    #    d_Unit = as.numeric(param[which(stringr::str_detect(param$achInFile, "dKpcUnit")), 2])
    #    G = 1 # NB: This is assumed! Can't find where it's specified
    #    v_Unit = sqrt(((.g_constant_cgs * G * 1e-15 / 1e-3) / .kpc_to_km) * .msol_to_kg)
    #    t_Unit = d_Unit/v_Unit * 0.97781311

    stars_age = t0 - stars_formed
    stars_age[which(stars_age < 0)] = 14 # setting any initial condition stars with large ages to the maximum
                                         # option for future users to change to a predetermined age distribution if desired

    ssp$Age = stars_age
  }

Similarly, that initial masses of stars are given in the .massform file output, which we map to initial masses as follows:

# This chunk should not run in the example, but we copy the code used to compute the stellar initial masses here. This can be found on line 255 in the read_files.R file:

  if (file.exists(paste0(f,".massform"))){
    mass_data = file(paste0(f,".massform"), "rb")
    stars_mass_formed = readBin(mass_data, "numeric", n = nstar, size = 4, endian = endian)
    close(mass_data)

    stars_mass_formed = stars_mass_formed*m_Unit*1e10 # initial mass in Msol
    # where m_Unit is taken from the parameter file:
    #   m_Unit = as.numeric(param[which(stringr::str_detect(param$achInFile, "dMsolUnit")), 2])

    ssp$Initial_Mass = stars_mass_formed
  }

For a gas observation

mask = observed_tipsy_cube$raw_images$mass_image
mask[mask<1e4] = NA
mask[!is.na(mask)] = 1

plot_mass(log10(observed_tipsy_cube_gas$observed_images$mass_image*mask), fig = c(0,0.33, 0, 1), labN = 2, titleshift = -7, units = expression("Mass, log"["10"]*"(M"["sol"]*")"))
plot_SFR(log10(observed_tipsy_cube_gas$observed_images$SFR_image)*mask, zlim = c(0, 10), fig = c(0.33,0.66, 0, 1), new=T, labN = 2, titleshift = -7, units = expression("SFR, log"["10"]*"(M"["sol"]*")"))
plot_velocity(observed_tipsy_cube_gas$observed_images$velocity_image*mask, fig = c(0.66, 0.99, 0, 1), new=T, labN = 2, titleshift = -7)

Behaviour of gas in SimSpin is treated in lines 202 - 232, as copied below. We have used treatments from pynbody to compute the mass fractions of different elements, the mean molecular mass given the gas temperature and hence the level of dispersion due to thermal effects in the gas. These are also copied below for consideration:


# This chunk should not run in the example, but we copy the code used to compute gas properties here. This can be found on lines 202-232 in the read_files.R file:

    # Helium mass fraction including correction based on metallicity, from
    # https://pynbody.github.io/pynbody/_modules/pynbody/snapshot/tipsy.html
    hetot = 0.236 + (2.1 * gas_part$Metallicity)
    # Hydrogen mass fraction including correction based on metallicity, from
    # https://pynbody.github.io/pynbody/_modules/pynbody/snapshot/tipsy.html
    gas_part$Hydrogen = 1.0 - gas_part$Metallicity - hetot
    gas_part$ID = 1:ngas

    #mean molecular mass, i.e. the mean atomic mass per particle
    mu = numeric(length = ngas)
    mu[gas_part$Temperature <= 1e4] = 0.59
    mu[gas_part$Temperature > 1e4] = 1.3

    #Gas internal energy derived from temperature
    gas_part$InternalEnergy = gas_part$Temperature / (mu *((5/3) - 1))
    gas_part$ThermalDispersion = sqrt((gas_part$InternalEnergy)*(.adiabatic_index - 1))

    #gridding gas particle positions for SFR computation
    cell_size = 0.1 # 100 pc - can be modified in future
    cell_seq = seq(min(min(gas_part$x), min(gas_part$y), min(gas_part$z))-cell_size,
                   max(max(gas_part$x), max(gas_part$y), max(gas_part$z))+cell_size, by = cell_size)
    cell_x_n = as.numeric(length(unique(cut(gas_part$x, breaks=cell_seq, labels=F))))
    cell_y_n = as.numeric(length(unique(cut(gas_part$y, breaks=cell_seq, labels=F))))
    cell_z_n = as.numeric(length(unique(cut(gas_part$z, breaks=cell_seq, labels=F))))

    gas_part$gas_cell = cut(gas_part$x, breaks=cell_seq, labels=F) +
      (cell_x_n * cut(gas_part$y, breaks=cell_seq, labels=F)) +
      (cell_x_n * cell_y_n * cut(gas_part$z, breaks=cell_seq, labels=F)) -
      (cell_x_n * cell_y_n) - cell_x_n

    gas_part$SFR = 0 # place holder for after reading in stars
    
    # ... while reading in stars, lines 302 - 328:
    if (ngas > 0) { # compute the SFR from stellar formation from coincident cells
      star_part$gas_cell = cut(star_part$x, breaks=cell_seq, labels=F) +
        (cell_x_n * cut(star_part$y, breaks=cell_seq, labels=F)) +
        (cell_x_n * cell_y_n * cut(star_part$z, breaks=cell_seq, labels=F)) -
        (cell_x_n * cell_y_n) - cell_x_n

      sfr_stars = which(stars_formed <= 0.01)

      summed_young = star_part[sfr_stars,
                              list(.N,
                              mass = sum(Mass)),
                              by = "gas_cell"]

      summed_young$SFR = summed_young$mass / 1e6
      summed_young$SFR[is.na(summed_young$gas_cell)] = 0

      part_in_cell = gas_part[, list(val=list(ID), .N), by = "gas_cell"]

      glue = merge(summed_young, part_in_cell, by = "gas_cell")

      vals = unlist(glue$val)

      gas_part$SFR[vals] = rep(glue$SFR, glue$N.y)

      remove(part_in_cell, glue, summed_young, vals, sfr_stars)

      }