Skip to content

Commit

Permalink
Merge pull request #41 from kateharborne/dev-LSF_bugfix
Browse files Browse the repository at this point in the history 
dev-LSF_bugfix
  • Loading branch information
@kateharborne
kateharborne committed Jun 13, 2022
2 parents bb371d9 + c4f3b80 commit 62afd85
Show file tree
Hide file tree
Showing 15 changed files with 298 additions and 190 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: SimSpin
Type: Package
Title: SimSpin - A package for the kinematic analysis of galaxy simulations
Version: 2.2.0
Version: 2.3.0
Author: Katherine Harborne
Maintainer: <katherine.harborne@icrar.org>
Description: The purpose of this package is to provide a set of tools to "observe" simulations of galaxies as if you were an observer using an integral field unit (IFU). The galaxy model can be observed at a chosen inclination and an IFU observation can be produced in a selection of different spatial shapes (square, circular, and hexagonal) to account for the variety of CCD arrangements used by current telescopes. This data cube can be used to calculate the value an observer would measure for the spin parameter compared to the value that is measured directly from the simulation.
Expand Down
79 changes: 41 additions & 38 deletions NEWS.md
View file

Large diffs are not rendered by default.

43 changes: 36 additions & 7 deletions R/build_datacube.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,18 +169,47 @@ build_datacube = function(simspin_file, telescope, observing_strategy,

if (observation$method == "spectral"){

wavelength = (observation$z * simspin_data$wave) + simspin_data$wave # applying a shift due to redshift, z, to original wavelength
lsf_fwhm = (observation$z * observation$lsf_fwhm) + observation$lsf_fwhm # adjusting the LSF for the resolution at z
# read original wavelengths of the template spectra
wavelength = simspin_data$wave * (observation$z + 1) # and then applying a shift to those spectra due to redshift, z

# If the requested wavelength resolution of the telescope is a smaller number than the intrinsic wavelength
# resolution of the template spectra, the interpolation onto a finer grid can cause errors that pPXF cannot
# account for in extreme cases. Issue a warning to the user in this case.

if (observation$wave_res < min(diff(wavelength))){
warning(cat("WARNING! - Wavelength resolution of provided template spectra at this redshift is too coarse for the requested telescope wavelength resolution.\n"))
cat("Dlambda_telescope = ", observation$wave_res, " A < Dlambda_templates ", min(diff(wavelength)), " A. \n")
cat("This will cause some interpolation that may make spectral fitting techniques fail. \n")
}

# Similarly, the template spectra for each particle have some intrinsic spectral resolution
# These vary dependent on the template. If the requested spectral resolution of the telescope
# is lower than the spectral resolution of the templates, we can't convolve them.

lsf_fwhm = observation$lsf_fwhm
lsf_fwhm_temp = simspin_data$header$Template_LSF * (observation$z + 1)
# applying a shift to that intrinsic template LSF due to redshift, z

spec_res_sigma_sq = ((lsf_fwhm^2) - (lsf_fwhm_temp^2))

spec_res_sigma_sq = lsf_fwhm^2 - (min(diff(wavelength)))^2
if (spec_res_sigma_sq < 0){ # if the lsf is smaller than the wavelength resolution of the spectra
warning(cat("WARNING! - Wavelength resolution of provided spectra is lower than the requested telescope resolution.\n"))
cat("LSF = ", observation$lsf_fwhm, " A < wavelength resolution ", min(diff(wavelength)), " A. \n")
cat("No LSF will be applied in this case.\n")
warning(cat("WARNING! - Spectral resolution of provided template spectra is greater than the requested telescope spectral resolution.\n"))
cat("LSF_telescope = ", lsf_fwhm, " A < LSF_templates (at redshift z) ", lsf_fwhm_temp, " A. \n")
cat("No LSF convolution will be applied in this case. \n")
cat("Intrinsic LSF of observation = ", lsf_fwhm_temp, " A for comparison with kinematic cubes. \n")
observation$LSF_conv = FALSE
} else {
observation$LSF_conv = TRUE
observation$lsf_sigma = (sqrt(spec_res_sigma_sq) / (2 * sqrt(2*log(2))))
observation$lsf_sigma = (sqrt(spec_res_sigma_sq) / (2 * sqrt(2*log(2)))) / (simspin_data$header$Template_waveres * (1 + observation$z))
# To get to the telescope's LSF, we only need to convolve with a Gaussian the width of the additional
# difference between the redshifted template and the intrinsic telescope LSF.
# This is the scaled for the wavelength pixel size at redshift "z".

for (spectrum in 1:length(simspin_data$spectra)){
convolved_spectrum = .lsf_convolution(observation, simspin_data$spectra[[spectrum]], observation$lsf_sigma)
simspin_data$spectra[[spectrum]] = convolved_spectrum
} # convolving the intrinsic spectra with the convolution kernel sized for the LSF

}

if (verbose){cat("Generating spectra per spaxel... \n")}
Expand Down
24 changes: 23 additions & 1 deletion R/make_simspin_file.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,15 @@ make_simspin_file = function(filename, cores=1, disk_age=5, bulge_age=10,
write_to_file=TRUE, output, overwrite = F,
centre=NA, half_mass=NA, sph_spawn_n=1){

header = list("InputFile" = filename,
"OutputFile" = NULL,
"Type" = character(1),
"Template" = character(1),
"Template_LSF" = numeric(1),
"Template_waveres" = numeric(1),
"Origin" = paste0("SimSpin_v", packageVersion("SimSpin")),
"Date" = Sys.time())

temp_name = stringr::str_to_upper(template)

if (write_to_file){
Expand All @@ -69,14 +78,24 @@ make_simspin_file = function(filename, cores=1, disk_age=5, bulge_age=10,
stop(cat("FileExists Error:: SimSpin file already exists at: ", output, "\n",
"If you wish to overwrite this file, please specify 'overwrite=T'. \n"))
}
header$OutputFile = output
}

if(temp_name == "BC03LR" | temp_name == "BC03"){
temp = SimSpin::BC03lr
header$Template = "BC03lr"
header$Template_LSF = 3 # as according to Bruzual & Charlot (2003) MNRAS 344, pg 1000-1028
header$Template_waveres = min(diff(temp$Wave))
} else if (temp_name == "BC03HR"){
temp = SimSpin::BC03hr
header$Template = "BC03hr"
header$Template_LSF = 3 # as according to Bruzual & Charlot (2003) MNRAS 344, pg 1000-1028
header$Template_waveres = min(diff(temp$Wave))
} else if (temp_name == "EMILES"){
temp = SimSpin::EMILES
header$Template = "EMILES"
header$Template_LSF = 2.51 # as according to Vazdekis et al (2016) MNRAS 463, pg 3409-3436
header$Template_waveres = min(diff(temp$Wave))
} else {
stop(cat("Error: template specified is unavailable.", "\n",
"Please specify template = 'BC03', 'BC03lr', 'BC03hr' or 'EMILES'"))
Expand All @@ -90,6 +109,8 @@ make_simspin_file = function(filename, cores=1, disk_age=5, bulge_age=10,
galaxy_data = tryCatch(expr = {.read_gadget(filename)},
error = function(e){.read_hdf5(filename, cores)})

header$Type = galaxy_data$head$Type

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
Expand Down Expand Up @@ -208,7 +229,8 @@ make_simspin_file = function(filename, cores=1, disk_age=5, bulge_age=10,

}

simspin_file = list("star_part" = galaxy_data$star_part,
simspin_file = list("header" = header,
"star_part" = galaxy_data$star_part,
"gas_part" = galaxy_data$gas_part,
"spectra" = sed,
"wave" = temp$Wave)
Expand Down
4 changes: 2 additions & 2 deletions R/observation.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,8 @@ observation = function(telescope, observing_strategy, method){
(telescope$sbin * sbin_size) / 2, by=sbin_size) # spatial bin break positions
wave_seq = seq(telescope$wave_range[1], telescope$wave_range[2], by = telescope$wave_res) # wavelength bin break positions
pixel_index = seq(1,telescope$sbin*telescope$sbin, by=1)
vbin_size = (telescope$wave_res / telescope$wave_centre) * (3e8/1e3) # km/s per velocity bin
vbin_error = ((telescope$lsf_fwhm / telescope$wave_centre) * (3e8/1e3)) / (2 * sqrt(2 * log(2))) # velocity uncertainty standard deviation
vbin_size = (telescope$wave_res / telescope$wave_centre) * .speed_of_light # km/s per velocity bin
vbin_error = ((telescope$lsf_fwhm / telescope$wave_centre) * .speed_of_light) / (2 * sqrt(2 * log(2))) # velocity uncertainty standard deviation

if (telescope$aperture_shape == "circular"){
aperture_region = .circular_ap(telescope$sbin)
Expand Down
160 changes: 90 additions & 70 deletions R/utilities.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -835,15 +835,6 @@ globalVariables(c(".N", ":=", "Age", "ID", "Initial_Mass", "Mass", "Metallicity"
return(as.vector(ap_region))
}

.interpolate_spectra = function(shifted_wave, spectra, wave_seq){ # function for interpolating the spectra onto a new grid

shifted_spectra = data.table::as.data.table(matrix(data=0.0, nrow = length(wave_seq), ncol=dim(shifted_wave)[1]))
for(j in 1:dim(shifted_wave)[1]){
data.table::set(x = shifted_spectra, i = NULL, j = j, value = stats::approx(x = shifted_wave[j,], y = spectra[[j]], xout = wave_seq, rule=1)[[2]])
}
return(rowSums(shifted_spectra))
}

.sum_velocities = function(galaxy_sample, observation){
vel_diff = function(lum, vy){diff((lum * pnorm(observation$vbin_edges, mean = vy,
sd = observation$vbin_error)))}
Expand All @@ -867,11 +858,13 @@ globalVariables(c(".N", ":=", "Age", "ID", "Initial_Mass", "Mass", "Metallicity"
# Function to apply LSF to spectra
.lsf_convolution = function(observation, luminosity, lsf_sigma){

scaled_sigma = lsf_sigma / observation$wave_res # scaling the size of the gaussian to match the pixel dimensions
kernel = stats::dnorm(seq(-scaled_sigma*5,scaled_sigma*5,length.out = 25), mean = 0, sd = scaled_sigma)
lsf_gauss = kernel/sum(kernel)
lum = stats::convolve(luminosity, lsf_gauss, type="open")
end = (length(luminosity) + length(lsf_gauss) - 1) - 12
kernel_radius = (4 * lsf_sigma + 0.5)
x = seq(-kernel_radius, kernel_radius, length.out = 25)
phi_x = exp((-0.5 / (lsf_sigma^2)) * (x^2)) / (lsf_sigma * sqrt(2*pi))
phi_x = phi_x / sum(phi_x)

lum = stats::convolve(luminosity, phi_x, type="open")
end = (length(luminosity) + length(phi_x) - 1) - 12

return(lum[13:end])
}
Expand Down Expand Up @@ -940,23 +933,29 @@ globalVariables(c(".N", ":=", "Age", "ID", "Initial_Mass", "Mass", "Metallicity"
part_map[part_in_spaxel$pixel_pos[i]] = num_part

galaxy_sample = galaxy_data[ID %in% part_in_spaxel$val[[i]]]
intrinsic_spectra = simspin_data$spectra[ , galaxy_sample$sed_id, with=FALSE] *
(galaxy_sample$Initial_Mass * 1e10) # reading relavent spectra

# pulling wavelengths and using doppler formula to compute the shift in
# wavelengths caused by LOS velocity
wave = matrix(data = rep(wavelength, num_part), nrow = num_part, byrow=T)
wave_shift = ((galaxy_sample$vy / .speed_of_light) * wave) + wave
luminosity = numeric(length(observation$wave_seq)) # initiallise a spectrum array for this pixel

for (p in 1:num_part){
intrinsic_spectra = simspin_data$spectra[[galaxy_sample$sed_id[p]]] *
(galaxy_sample$Initial_Mass[p] * 1e10) # reading relavent spectra

# interpolate each shifted wavelength to telescope grid of wavelengths
# and sum to one spectra
luminosity = .interpolate_spectra(shifted_wave = wave_shift, spectra = intrinsic_spectra,
wave_seq = observation$wave_seq)
# pulling wavelengths and using doppler formula to compute the shift in
# wavelengths caused by LOS velocity
wave_shift = ((galaxy_sample$vy[p] / .speed_of_light) * wavelength) + wavelength

if (observation$LSF_conv){ # should the spectra be degraded for telescope LSF?
luminosity = .lsf_convolution(observation=observation, luminosity=luminosity, lsf_sigma=observation$lsf_sigma)
# interpolate each shifted wavelength to telescope grid of wavelengths
# and sum to one spectra
part_lum = stats::approx(x = wave_shift, y = intrinsic_spectra, xout = observation$wave_seq, rule=1)[[2]]

luminosity = luminosity + part_lum
}

# if (observation$LSF_conv){ # should the spectra be degraded for telescope LSF?
# luminosity = .lsf_convolution(observation=observation, luminosity=luminosity,
# lsf_sigma=observation$lsf_sigma)
# }

if (!is.na(observation$signal_to_noise)){ # should we add noise?
luminosity = .add_noise(luminosity, observation$signal_to_noise)
}
Expand Down Expand Up @@ -995,24 +994,30 @@ globalVariables(c(".N", ":=", "Age", "ID", "Initial_Mass", "Mass", "Metallicity"
part_map = num_part

galaxy_sample = galaxy_data[ID %in% part_in_spaxel$val[[i]]]
intrinsic_spectra = simspin_data$spectra[ , galaxy_sample$sed_id, with=FALSE] *
(galaxy_sample$Initial_Mass * 1e10) # reading relavent spectra

# pulling wavelengths and using doppler formula to compute the shift in
# wavelengths caused by LOS velocity
wave = matrix(data = rep(wavelength, num_part), nrow = num_part, byrow=T)
wave_shift = ((galaxy_sample$vy / .speed_of_light) * wave) + wave

# interpolate each shifted wavelength to telescope grid of wavelengths
# and sum to one spectra
luminosity = .interpolate_spectra(shifted_wave = wave_shift,
spectra = intrinsic_spectra,
wave_seq = observation$wave_seq)

if (observation$LSF_conv){ # should the spectra be degraded for telescope LSF?
luminosity = .lsf_convolution(observation=observation, luminosity=luminosity,
lsf_sigma=observation$lsf_sigma)

luminosity = numeric(length(observation$wave_seq)) # initialize a spectrum array for this pixel

for (p in 1:num_part){
intrinsic_spectra = simspin_data$spectra[[galaxy_sample$sed_id[p]]] *
(galaxy_sample$Initial_Mass[p] * 1e10) # reading relevant spectra

# pulling wavelengths and using doppler formula to compute the shift in
# wavelengths caused by LOS velocity
wave_shift = ((galaxy_sample$vy[p] / .speed_of_light) * wavelength) + wavelength

# interpolate each shifted wavelength to telescope grid of wavelengths
# and sum to one spectra
part_lum = stats::approx(x = wave_shift, y = intrinsic_spectra, xout = observation$wave_seq, rule=1)[[2]]

luminosity = luminosity + part_lum

}

# if (observation$LSF_conv){ # should the spectra be degraded for telescope LSF?
# luminosity = .lsf_convolution(observation=observation, luminosity=luminosity,
# lsf_sigma=observation$lsf_sigma)
# }

if (!is.na(observation$signal_to_noise)){ # should we add noise?
luminosity = .add_noise(luminosity, observation$signal_to_noise)
}
Expand Down Expand Up @@ -1065,20 +1070,27 @@ globalVariables(c(".N", ":=", "Age", "ID", "Initial_Mass", "Mass", "Metallicity"
galaxy_sample[, luminosity := Mass, ]

} else {
intrinsic_spectra = simspin_data$spectra[ , c(galaxy_sample$sed_id), with=FALSE] *
(galaxy_sample$Initial_Mass * 1e10) # reading relavent spectra

# transform luminosity into flux detected at telescope
# flux in units erg/s/cm^2/Ang
for (j in 1:num_part){
set(intrinsic_spectra, i=NULL, j = j,
value = ((intrinsic_spectra[[j]]*.lsol_to_erg) / (4 * pi * (observation$lum_dist*.mpc_to_cm)^2) / (1 + observation$z)))
}

# computing the r-band luminosity per particle from spectra
galaxy_sample[ , luminosity := .bandpass(wave = simspin_data$wave,
flux = intrinsic_spectra,
filter = filter), ]

band_lum = numeric(num_part)

for (p in 1:num_part){

intrinsic_spectra = simspin_data$spectra[[galaxy_sample$sed_id[p]]] *
(galaxy_sample$Initial_Mass[p] * 1e10) # reading particle spectrum

# transform luminosity into flux detected at telescope
# flux in units erg/s/cm^2/Ang
flux_spectra =
(intrinsic_spectra*.lsol_to_erg) / (4 * pi * (observation$lum_dist*.mpc_to_cm)^2) / (1 + observation$z)

# computing the r-band luminosity per particle from spectra
band_lum[p] = .bandpass(wave = simspin_data$wave,
flux = flux_spectra,
filter = filter)
}

galaxy_sample[ , luminosity := band_lum, ]

}


Expand Down Expand Up @@ -1124,20 +1136,28 @@ globalVariables(c(".N", ":=", "Age", "ID", "Initial_Mass", "Mass", "Metallicity"
galaxy_sample[, luminosity := Mass, ]

} else {
intrinsic_spectra = simspin_data$spectra[ , galaxy_sample$sed_id, with=FALSE] *
(galaxy_sample$Initial_Mass * 1e10) # reading relavent spectra
# transform luminosity into flux detected at telescope
# flux in units erg/s/cm^2/Ang
for (j in 1:num_part){
set(intrinsic_spectra, i=NULL, j = j,
value = ((intrinsic_spectra[[j]]*.lsol_to_erg) / (4 * pi * (observation$lum_dist*.mpc_to_cm)^2) / (1 + observation$z)))

band_lum = numeric(num_part)

for (p in 1:num_part){
intrinsic_spectra = simspin_data$spectra[[galaxy_sample$sed_id[p]]] *
(galaxy_sample$Initial_Mass[p] * 1e10) # reading relavent spectra

# transform luminosity into flux detected at telescope
# flux in units erg/s/cm^2/Ang
flux_spectra =
(intrinsic_spectra*.lsol_to_erg) / (4 * pi * (observation$lum_dist*.mpc_to_cm)^2) / (1 + observation$z)

# computing the r-band luminosity per particle from spectra
band_lum[p] = .bandpass(wave = simspin_data$wave,
flux = flux_spectra,
filter = filter)

}
# computing the r-band luminosity per particle from spectra
# computing the r-band luminosity per particle from spectra
galaxy_sample[ , luminosity := .bandpass(wave = simspin_data$wave,
flux = intrinsic_spectra,
filter = filter), ]
}

galaxy_sample[ , luminosity := band_lum, ]

}

# adding the "gaussians" of each particle to the velocity bins
vel_spec = .sum_velocities(galaxy_sample = galaxy_sample, observation = observation)
Expand Down
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@

<p>&nbsp;</p>

SimSpin v2.2.0 - A package for producing mock observations of particle simulations
SimSpin v2.3.0 - A package for producing mock observations of particle simulations

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.

Expand Down
Binary file modified data/BC03hr.rda
View file
Binary file not shown.
Binary file modified data/BC03lr.rda
View file
Binary file not shown.
Binary file modified data/EMILES.rda
View file
Binary file not shown.
Binary file modified inst/extdata/SimSpin_example_Gadget_spectra.Rdata
View file
Binary file not shown.

0 comments on commit 62afd85

Please sign in to comment.