-
Notifications
You must be signed in to change notification settings - Fork 3
/
make_simspin_file.R
225 lines (191 loc) · 10.9 KB
/
make_simspin_file.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
# Author: Kate Harborne
# Date: 22/10/2020
# Title: Testing the make_simspin_file.R code
#'Reformating isolated galaxy simulations to contain spectra.
#'
#'The purpose of this function is to construct a SimSpin file containing the
#' mock spectra for each particle contained within the galaxy simulation file.
#' If the snapshot provided is from a cosmological simulation, the SEDs
#' generated will be with respect to the Stellar Formation Time/Age, Metallicity
#' and Initial Mass of each stellar particle. If the system is an N-body model,
#' stellar particles are assumed to have an age and metallicity as provided to
#' the function as \code{disk_age}, \code{bulge_age}, \code{disk_Z} and
#' \code{bulge_Z}. Returned is an .Rdata file in a SimSpin readable format.
#'
#'@param filename The path to the snapshot file.
#'@param cores The number of cores across which to multi-thread the problem.
#'@param disk_age The age of the disk particles in Gyr.
#'@param bulge_age The age of the bulge particles in Gyr.
#'@param disk_Z The metallicity of the disk particles in Gyr.
#'@param bulge_Z The metallicity of the bulge particles in Gyr.
#'@param template The stellar templates from which to derive the SEDs. Options
#' include "BC03lr" (GALEXEV low resolution, Bruzual & Charlot 2003), "BC03hr"
#' (GALEXEV high resolution, Bruzual & Charlot 2003) or "EMILES" (Vazdekis et
#' al, 2016).
#'@param write_to_file Boolean to specify whether the list produced should be
#' written to a ".Rdata" file or output to the environment. Default is TRUE, so
#' that files can be re-observed without having the generate spectra each time.
#'@param output The path at which the output file is written. If not provided,
#' file will be written at the location of the input filename with the addition
#' of "_spectra.Rdata".
#'@param overwrite If true, and the file already exists at the output location,
#' a new file will be written over the old one.
#'@param centre If simulation file contains all particles cutout from a box
#' (rather than just particles from a single galaxy), you can specify the point
#' around which the view should be centred. Numeric length = 3. Default is NA,
#' in which case the system is centred around the median position.
#'@param half_mass If simulation file contains all particles cutout from a box
#' (rather than just particles from a single galaxy), you can the half-mass
#' value at which the alignment function is run. Numeric length = 1. Default is
#' NA, in which case half the total mass of the suplied simulation data is used.
#'@param sph_spawn_n Numeric describing the number of gas particles with which
#' to sample the SPH smoothing length. Default is 1, which will not spawn
#' additional gas particles. Increasing this value increases the number of
#' particles used to model the gas distribution. This value may need to be
#' tested for convergence depending on the resolution of the grid used to image
#' the gas properties at the `build_datacube()` stage.
#'@return Returns an .Rdata file that contains a list of particle positions,
#' velocities, and spectra (or a list containing the same information to the
#' environment without writing to file, when `write_to_file = F`).
#'@examples
#'ss_file = make_simspin_file(filename = system.file("extdata",
#' "SimSpin_example_Gadget",
#' package = "SimSpin"),
#' write_to_file = FALSE)
#'
make_simspin_file = function(filename, cores=1, disk_age=5, bulge_age=10,
disk_Z=0.024, bulge_Z=0.001, template="BC03lr",
write_to_file=TRUE, output, overwrite = F,
centre=NA, half_mass=NA, sph_spawn_n=1){
temp_name = stringr::str_to_upper(template)
if (write_to_file){
if (missing(output)){
output = paste(sub('\\..*', '', filename), "_", temp_name, ".Rdata", sep="")
}
if (file.exists(output) & !overwrite){
stop(cat("FileExists Error:: SimSpin file already exists at: ", output, "\n",
"If you wish to overwrite this file, please specify 'overwrite=T'. \n"))
}
}
if(temp_name == "BC03LR" | temp_name == "BC03"){
temp = SimSpin::BC03lr
} else if (temp_name == "BC03HR"){
temp = SimSpin::BC03hr
} else if (temp_name == "EMILES"){
temp = SimSpin::EMILES
} else {
stop(cat("Error: template specified is unavailable.", "\n",
"Please specify template = 'BC03', 'BC03lr', 'BC03hr' or 'EMILES'"))
}
if (sph_spawn_n%%1!=0){
stop(cat("Error: sph_spawn_n must be a whole number.", "\n",
"Please specify an integer value for sph_spawn_n."))
}
galaxy_data = tryCatch(expr = {.read_gadget(filename)},
error = function(e){.read_hdf5(filename, cores)})
Npart_sum = cumsum(galaxy_data$head$Npart) # Particle indices of each type
galaxy_data = .centre_galaxy(galaxy_data, centre) # centering the galaxy based on stellar particles
galaxy_data = .align_galaxy(galaxy_data, half_mass) # align 3D shape of galaxy
if(!"ssp" %in% names(galaxy_data)){ # if the SSP field does not come from the snapshot file, must be working with N-body
n_disk = galaxy_data$head$Npart[3]; n_bulge = galaxy_data$head$Npart[4] # number of disk and bulge particles
n_stars = n_disk + n_bulge # total number of "stars"
galaxy_data$ssp = data.frame("Initial_Mass"=numeric(n_stars), "Age"=numeric(n_stars),
"Metallicity"=numeric(n_stars))
galaxy_data$ssp$Initial_Mass = galaxy_data$star_part$Mass[Npart_sum[2]+1:Npart_sum[4]]/2 # assuming the initial mass is half of the current mass
if (n_disk > 0 & n_bulge > 0){ # assigning ages and metallities to disk and bulge particles (if present in snap)
galaxy_data$ssp$Age[1:n_disk] = disk_age
galaxy_data$ssp$Age[(n_disk+1):n_stars] = bulge_age
galaxy_data$ssp$Metallicity[1:n_disk] = disk_Z
galaxy_data$ssp$Metallicity[(n_disk+1):n_stars] = bulge_Z
AZ_bins = data.frame("ages" = c(disk_age, bulge_age),
"metallicities" = c(disk_Z, bulge_Z),
"id" = c(1,2),
"sed_id" = c(1,2))
az_pos = c(rep(1, n_disk), rep(2, n_bulge))
} else if (n_disk > 0 & n_bulge == 0){
galaxy_data$ssp$Age = disk_age
galaxy_data$ssp$Metallicity = disk_Z
AZ_bins = data.frame("ages" = disk_age,
"metallicities" = disk_Z,
"id" = 1,
"sed_id" = 1)
az_pos = rep(1, n_disk)
} else if (n_disk == 0 & n_bulge > 0){
galaxy_data$ssp$Age = bulge_age
galaxy_data$ssp$Metallicity = bulge_Z
AZ_bins = data.frame("ages" = bulge_age,
"metallicities" = bulge_Z,
"id" = 1,
"sed_id" = 1)
az_pos = rep(1, n_bulge)
}
}
# If a particle has a metallicity of 0, remove from sample?
Z0_int = which(galaxy_data$ssp$Metallicity == 0)
if (length(Z0_int)!=0){ # remove if there are any Z = 0
galaxy_data$star_part = galaxy_data$star_part[-c(Z0_int),]
galaxy_data$ssp = galaxy_data$ssp[-c(Z0_int),]
}
# now binning stellar particles based on their A/Z position
if (length(unique(galaxy_data$ssp$Age)) > 2){ # if not an N-body simulation
age_grid = 10^(seq(log10(min(galaxy_data$ssp$Age))-0.02, log10(max(galaxy_data$ssp$Age))+0.02, by = 0.02))
Z_grid = 10^(seq(log10(min(galaxy_data$ssp$Metallicity))-0.1, log10(max(galaxy_data$ssp$Metallicity))+0.1, by = 0.1))
age_cen = 10^(seq(log10(min(galaxy_data$ssp$Age))-0.01, log10(max(galaxy_data$ssp$Age))+0.01, by = 0.02))
Z_cen = 10^(seq(log10(min(galaxy_data$ssp$Metallicity))-0.05, log10(max(galaxy_data$ssp$Metallicity))+0.05, by = 0.1))
AZ = data.frame("ages" = rep(age_cen, length = length(age_cen)*length(Z_cen)),
"metallicities" = rep(Z_cen, each = length(age_cen)),
"id" = seq(1, length(age_cen)*length(Z_cen)))
az_pos = cut(galaxy_data$ssp$Age, breaks = age_grid, labels = F) +
(length(age_cen) * cut(galaxy_data$ssp$Metallicity, breaks = Z_grid, labels = F)) - length(age_cen)
AZ_bins = AZ[sort(unique(az_pos)),]
AZ_bins$sed_id = seq(1, length(AZ_bins$id))
}
if (!is.null(galaxy_data$ssp)){ # if there are stellar particles in the file at all
# adding info to stellar_particle data
galaxy_data$star_part[, sed_id := AZ_bins$sed_id[match(az_pos,AZ_bins$id)], ]
galaxy_data$star_part[, Metallicity := galaxy_data$ssp$Metallicity, ]
galaxy_data$star_part[, Age := galaxy_data$ssp$Age, ]
galaxy_data$star_part[, Initial_Mass := galaxy_data$ssp$Initial_Mass, ]
sed = .spectra(Metallicity = AZ_bins$metallicities, Age = AZ_bins$ages, Template = temp, cores = cores) # returns a list
} else {sed = NULL}
if (length(galaxy_data$gas_part$SmoothingLength)>0 & sph_spawn_n>1){ # if we need to spawn gas particles beacuse we are working with SPH models
gas_part_names = names(galaxy_data$gas_part)
new_gas_part = stats::setNames(data.frame(matrix(ncol = length(gas_part_names),
nrow = (sph_spawn_n*galaxy_data$head$NumPart_Total[1]))),
gas_part_names) # generate a new DF containing original gas_part names
new_gas_part$vx = rep(galaxy_data$gas_part$vx, each=sph_spawn_n)
new_gas_part$vy = rep(galaxy_data$gas_part$vy, each=sph_spawn_n)
new_gas_part$vz = rep(galaxy_data$gas_part$vz, each=sph_spawn_n)
new_gas_part$SFR = rep(galaxy_data$gas_part$SFR, each=sph_spawn_n)
new_gas_part$Density = rep(galaxy_data$gas_part$Density, each=sph_spawn_n)
new_gas_part$Temperature = rep(galaxy_data$gas_part$Temperature, each=sph_spawn_n)
new_gas_part$SmoothingLength = 0
new_gas_part$Metallicity = rep(galaxy_data$gas_part$Metallicity, each=sph_spawn_n)
new_gas_part$Carbon = rep(galaxy_data$gas_part$Carbon, each=sph_spawn_n)
new_gas_part$Hydrogen = rep(galaxy_data$gas_part$Hydrogen, each=sph_spawn_n)
new_gas_part$Oxygen = rep(galaxy_data$gas_part$Oxygen, each=sph_spawn_n)
kernel = character(1) # choosing the kernel relevent for the
if (galaxy_data$head$Type == "EAGLE"){
kernel = "WC2"
} else if (galaxy_data$head$Type == "Magneticum"){
kernel = "WC6"
} else {
kernel = "WC2"
}
if (cores > 1){
galaxy_data$gas_part = .sph_spawn_mc(galaxy_data$gas_part, new_gas_part, sph_spawn_n, kernel, cores)
} else {
galaxy_data$gas_part = .sph_spawn(galaxy_data$gas_part, new_gas_part, sph_spawn_n, kernel)
}
}
simspin_file = list("star_part" = galaxy_data$star_part,
"gas_part" = galaxy_data$gas_part,
"spectra" = sed,
"wave" = temp$Wave)
if (write_to_file){
saveRDS(simspin_file, file = output)
return(message("SimSpin file written to: ", output, "\n"))
} else {
return(simspin_file)
}
}