R/readgadget.R

In obreschkow/simstar: Reading, processing and visualising astrophysical N-body simulation data

#' Read Gadget simulation data
#' 
#' @importFrom rhdf5 h5ls h5read h5readAttributes
#'
#' @description Reads astrophysical N-body data output by the Gadget code (versions 1/2/3/4, see \url{https://wwwmpa.mpa-garching.mpg.de/gadget4}). Binary and HDF5 formats can be read.
#'
#' @param file filename of file to load
#' @param type character specifying the data format. Must be either of: \code{bin} for binary format, \code{hdf} for HDF5 format or \code{auto} to automatically determine the format from file extension.
#' 
#' @return Returns a list containing the particle data, structured very similarly to the HDF5 format, even for binary files. When reading a binary file, some of the entries in the header are duplicated for compatibility between old (binary) and new (HDF5) parameter names: e.g. \code{NumPart_ThisFile=Npart}, \code{NumPart_Total=Nall}, \code{MassTable=Massarr} (for more, see Table 4 of \url{https://wwwmpa.mpa-garching.mpg.de/gadget/users-guide.pdf}).
#' 
#' @seealso \code{\link{writegadget}}
#' 
#' @author Danail Obreschkow
#'
#' @export

readgadget = function(file, type='auto') {
  
  # check file
  if (!file.exists(file)) stop(paste0('File not found: ',file))
  if (file.access(file,4)!=0) stop(paste0('No permission to read: ',file))
  
  # handle type
  if (type=='binary') type = 'bin'
  if (type=='hdf5' | type=='h5') type = 'hdf'
  if (type=='auto') {
    
    # get file extension
    ext = strsplit(basename(file), split="\\.")[[1]]
    ext = ext[-1]
    
    # check file extension
    if (length(ext)>0) {
      if (ext%in%c('hdf','hdf5','h5')) {
        type = 'hdf'
      } else {
        type = 'bin'
      }
    } else {
      type = 'bin'
    }
    
  }
  
  # load file
  if (type=='bin') {
    
    # initialize
    block.size = NA
    block.start = function(target=NA) {
      if (is.na(block.size)) {
        block.size <<- readBin(data, 'integer', size=4)
        if (target>2^23-1) target=NA
        if (!is.na(target)) {
          if (target!=block.size) {
            close(data)
            print(block.size)
            print(target)
            stop('block opening error')
          }
        }
      } else {
        close(data)
        stop('block.start must be called after block.end')
      }
    }
    block.end = function() {
      if (is.na(block.size)) {
        close(data)
        stop('block.end must be called after block.start')
      } else {
        check = readBin(data, 'integer', size=4)
        if (check!=block.size) {
          close(data)
          stop('block termination error')
        }
        block.size <<- NA
      }
    }
    
    # open file
    data = file(file, "rb")
    dat = list()
    
    # read header
    block.start(256)
    dat$Header = list()
    dat$Header$Npart = readBin(data, 'integer', 6, size = 4)
    dat$Header$Massarr = readBin(data, 'numeric', 6, size = 8)
    dat$Header$Time = readBin(data, 'numeric', 1, size = 8)
    dat$Header$Redshift = readBin(data, 'numeric', 1, size = 8)
    dat$Header$FlagSfr = readBin(data, 'integer', 1, size = 4)
    dat$Header$FlagFeedback = readBin(data, 'integer', 1, size = 4)
    dat$Header$Nall = readBin(data, 'integer', 6, size = 4)
    dat$Header$FlagCooling = readBin(data, 'integer', 1, size = 4)
    dat$Header$NumFiles = readBin(data, 'integer', 1, size = 4)
    dat$Header$BoxSize = readBin(data, 'numeric', 1, size = 8)
    dat$Header$OmegaM = readBin(data, 'numeric', 1, size = 8)
    dat$Header$OmegaL = readBin(data, 'numeric', 1, size = 8)
    dat$Header$h = readBin(data, 'numeric', 1, size = 8)
    dat$Header$FlagAge = readBin(data, 'integer', 1, size = 4)
    dat$Header$FlagMetals = readBin(data, 'integer', 1, size = 4)
    dat$Header$NallHW = readBin(data, 'integer', 6, size = 4)
    dat$Header$flag_entr_ics = readBin(data, 'integer', 1, size = 4)
    dat$Header$unused = readBin(data, 'integer', 15, size = 4)
    block.end()
    
    # duplicate header entries for compatibility with HDF5 names
    dat$Header$NumPart_ThisFile = dat$Header$Npart
    dat$Header$MassTable = dat$Header$Massarr
    dat$Header$Flag_Sfr = dat$Header$FlagSfr
    dat$Header$Flag_Feedback = dat$Header$FlagFeedback
    dat$Header$NumPart_Total = dat$Header$Nall
    dat$Header$Flag_Cooling = dat$Header$FlagCooling
    dat$Header$NumFilesPerSnapshot = dat$Header$NumFiles
    dat$Header$Flag_StellarAge = dat$Header$FlagAge
    dat$Header$Flag_Metals = dat$Header$FlagMetals
    dat$Header$NumPart_Total_HW = dat$Header$NallHW
    dat$Header$Flag_Entropy_ICs = dat$Header$flag_entr_ics
    
    # read positions
    n = sum(dat$Header$Npart)
    block.start(12*n)
    for (i in seq(0,5)) {
      k = dat$Header$Npart[i+1]
      if (k>0) {
        field = sprintf('PartType%d',i)
        dat[[field]] = list()
        dat[[field]]$Coordinates = t(array(readBin(data, 'numeric', 3*k, size = 4),c(3,k)))
      }
    }
    block.end()
    
    # read velocities
    block.start(12*n)
    for (i in seq(0,5)) {
      k = dat$Header$Npart[i+1]
      if (k>0) {
        field = sprintf('PartType%d',i)
        dat[[field]]$Velocities = t(array(readBin(data, 'numeric', 3*k, size = 4),c(3,k)))
      }
    }
    block.end()
    
    # read IDs
    n = sum(dat$Header$Npart)
    block.start(4*n)
    for (i in seq(0,5)) {
      k = dat$Header$Npart[i+1]
      if (k>0) {
        field = sprintf('PartType%d',i)
        dat[[field]]$ParticleIDs = readBin(data, 'integer', k, size = 4)
      }
    }
    block.end()
    
    # read masses
    nmass = dat$Header$Npart[dat$Header$Massarr==0]
    if (sum(nmass)>0) {
      block.start(4*sum(nmass))
      for (i in seq(0,5)) {
        if (nmass[i+1]>0) {
          field = sprintf('PartType%d',i)
          dat[[field]]$Masses = readBin(data, 'numeric', nmass[i+1], size = 4)
        }
      }
      block.end()
    }
    
    # read internal energies
    ngas = dat$Header$Npart[1]
    if (ngas >0) {
      block.start(4*ngas)
      dat$PartType0$InternalEnergy = readBin(data, 'numeric', ngas, size = 4)
      block.end()
    }
    
    # close file
    close(data)
    
    return(dat)
    
  } else if (type=='hdf') {
    
    groups = rhdf5::h5ls(file,recursive=FALSE)$name
    dat = list()
    if (length(groups)<1) stop('something wrong with HDF5 file')
    if ('Config'%in%groups) dat$Config = h5readAttributes(file,'/Config')
    if ('Header'%in%groups) dat$Header = h5readAttributes(file,'/Header')
    if ('Parameters'%in%groups) dat$Parameters = h5readAttributes(file,'/Parameters')
    for (i in seq(0,5)) {
      field = sprintf('PartType%d',i)
      if (field%in%groups) {
        raw = rhdf5::h5read(file,field)
        if (!is.null(raw$Coordinates)) raw$Coordinates = t(raw$Coordinates)
        if (!is.null(raw$Velocities)) raw$Velocities = t(raw$Velocities)
        dat[[field]] = raw
      }
    }
    return(dat)
    
  } else {
    
    stop('Unknown "type".')
    
  }
  
}