Merge pull request #86 from kateharborne/dev-thermal

dev-Thermal

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SimSpin
 Type: Package
 Title: SimSpin - A package for the kinematic analysis of galaxy simulations
-Version: 2.4.8
+Version: 2.4.9
 Author: Katherine Harborne
 Co-author: Alice Serene
 Maintainer: <katherine.harborne@icrar.org>

NEWS.md

@@ -1,6 +1,6 @@
 # SimSpin v2.4.8 News
 
-### Last edit: 09/06/23
+### Last edit: 13/06/23
 
 Below is a table containing a summary of all changes made to SimSpin, since the date this file was created on 26/08/2021.
 
@@ -16,7 +16,8 @@ All changes are noted in the changelog table below.
 
 | Date     	| Summary of change                                                                                                                                                                                                                                                                                                                                                                                                        	| Version 	| Commit                                   	| Author            |
 |----------	|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------	|---------	|------------------------------------------	| ----------------- |
-| 09/06/23 | Fixing the SFR maps produced in method = `gas` or `sf gas` mode. Previous values returned were mass-weighted mean SFR along the line-of-sight, now adjusted to give the sum of spawned fraction SFR along the LOS. | 2.4.8 | | Kate Harborne | 
+| 13/06/23 | Adding the thermal dispersion component to the observed gas maps. | 2.4.9 | | Kate Harborne |
+| 09/06/23 | Fixing the SFR maps produced in method = `gas` or `sf gas` mode. Previous values returned were mass-weighted mean SFR along the line-of-sight, now adjusted to give the sum of spawned fraction SFR along the LOS. | 2.4.8 | e1a6dd49d0aa74510c27aba5ccfeac8ab8a5e920 | Kate Harborne | 
 | 26/04/23 | Inverse variance cubes have large values that do not save to FITS successfully. Dividing by 1e40 to remove this issue and incorporating this with the units. Also modifying S/N description to give the MINIMUM S/N in the image, not the maximum. | 2.4.7 | fb904cd302cc6b8edda2032bb68fa8cd6b233f8a | Kate Harborne |
 | 23/04/23 | Fixing the noise implementation in methods spectral and velocity. Adding inverse variance cubes to the output when a S/N value is specified. Incorporating tests to check proper incorporation and behaviour. | 2.4.6 | 8eef4e6dea8f9901b83ced8f2b43cee1438fe034 | Kate Harborne |
 | 18/04/23 | Updating the methodology for method = "velocity" (details below) and fixing dependencies (resolving the Rfits and ProFit versions clash). Added header information to the SimSpin file made such that the whole file can be recreated using the information in the file rather than having to retain the code used to build the file. | 2.4.5 | cdf5522fb6040b7a5fdcf5ea4a54df7e576f57d5 | Kate Harborne |

R/make_simspin_file.R

@@ -236,9 +236,6 @@ make_simspin_file = function(filename, cores=1, disk_age=5, bulge_age=10,
       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 := rep(0, galaxy_data$head$NumPart_Total[1]*sph_spawn_n),]
       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),]

R/utilities.R

@@ -10,6 +10,8 @@
 .cm_to_kpc          = 3.24078e-22
 .cms_to_kms         = 1e-5
 .g_to_msol          = 5.02785e-34
+.gcm1_to_msolkm1    = 5.02785e-29
+.gcm3_to_msolkm3    = 5.02785e-28
 .gcm3_to_msolkpc3   = 1.477e+31
 .g_constant_cgs     = 6.67430e-11
 .g_in_kpcMsolkms2   = 4.3009e-6
@@ -23,7 +25,7 @@
 globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen",
                   "hcl.colors", "ID", "Initial_Mass", "luminosity", "Mass",
                   "Metallicity", "Oxygen", "par", "SFR", "SFT", "SmoothingLength",
-                  "sed_id", "Temperature", "text", "vx", "vy", "vz", "x", "y",
+                  "sed_id", "Temperature", "ThermalDispersion", "text", "vx", "vy", "vz", "x", "y",
                   "z"))
 
 # Functions for computing weighted means
@@ -381,10 +383,10 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
     }
 
     gas = .check_names(gas)
-    eagle_gas_names = c("SmoothingLength", "Temperature")
+    eagle_gas_names = c("SmoothingLength", "Temperature", "InternalEnergy")
     if (!all(eagle_gas_names %in% names(gas))){
       stop("Error. Missing a necessary dataset for EAGLE PartType0. \n
-           Either `SmoothingLength` or `Temperature`. \n
+           Either `SmoothingLength`, `Temperature`, or `InternalEnergy`. \n
            See https://kateharborne.github.io/SimSpin/examples/generating_hdf5.html#parttype0 for more info.")
     }
 
@@ -403,11 +405,14 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
                                       "Density" = gas$Density*.gcm3_to_msolkpc3, # Density in Msol/kpc^3
                                       "Temperature" = gas$Temperature,
                                       "SmoothingLength" = gas$SmoothingLength*.cm_to_kpc, # Smoothing length in kpc
+                                      "ThermalDispersion" = sqrt((gas$InternalEnergy*.cms_to_kms)*(.adiabatic_index - 1)),
                                       "Metallicity" = gas$Metallicity,
                                       "Carbon" = gas$`ElementAbundance/Carbon`,
                                       "Hydrogen" = gas$`ElementAbundance/Hydrogen`,
                                       "Oxygen" = gas$`ElementAbundance/Oxygen`)
 
+    gas_part$ThermalDispersion[gas_part$Temperature <= 1e4] = 11
+
     remove(gas); remove(PT0_attr)
 
   } else {gas_part=NULL}
@@ -474,10 +479,10 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
     }
 
     gas = .check_names(gas)
-    magneticum_gas_names = c("SmoothingLength", "Temperature")
+    magneticum_gas_names = c("SmoothingLength", "Temperature", "InternalEnergy")
     if (!all(magneticum_gas_names %in% names(gas))){
       stop("Error. Missing a necessary dataset for Magneticum PartType0. \n
-           Either `SmoothingLength` or `Temperature`. \n
+           Either `SmoothingLength`, `Temperature` or `InternalEnergy`. \n
            See https://kateharborne.github.io/SimSpin/examples/generating_hdf5.html#parttype0 for more info.")
     }
 
@@ -497,11 +502,14 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
                                       "Density" = gas$Density*.gcm3_to_msolkpc3, # Density in Msol/kpc^3
                                       "Temperature" = gas$Temperature,
                                       "SmoothingLength" = gas$SmoothingLength*.cm_to_kpc, # Smoothing length in kpc
+                                      "ThermalDispersion" = sqrt((gas$InternalEnergy*.cms_to_kms)*(.adiabatic_index - 1)),
                                       "Metallicity" = gas$Metallicity,
                                       "Carbon" = gas$`ElementAbundance/Carbon`,
                                       "Hydrogen" = gas$`ElementAbundance/Hydrogen`,
                                       "Oxygen" =  gas$`ElementAbundance/Oxygen`)
 
+    gas_part$ThermalDispersion[gas_part$Temperature <= 1e4] = 11
+
     remove(gas); remove(PT0_attr)
 
   } else {gas_part=NULL}
@@ -571,10 +579,10 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
     }
 
     gas = .check_names(gas)
-    horizon_gas_names = c("Temperature")
+    horizon_gas_names = c("Temperature", "Pressure")
     if (!all(horizon_gas_names %in% names(gas))){
       stop("Error. Missing a necessary dataset for HorizonAGN PartType0. \n
-           Missing `Temperature`. \n
+           Missing `Temperature` or `Pressure`. \n
            See https://kateharborne.github.io/SimSpin/examples/generating_hdf5.html#parttype0 for more info.")
     }
 
@@ -592,12 +600,15 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
                                       "SFR" = gas$StarFormationRate*(.g_to_msol/.s_to_yr), #SFR in Msol/yr
                                       "Density" = gas$Density*.gcm3_to_msolkpc3, # Density in Msol/kpc^3
                                       "Temperature" = gas$Temperature,
+                                      "ThermalDispersion" = sqrt((gas$Pressure*.gcm1_to_msolkm1)/(gas$Density*.gcm3_to_msolkm3)),
                                       "SmoothingLength" = 2*(((3/(4*pi))*((gas$Mass*.g_to_msol) / (gas$Density*.gcm3_to_msolkpc3)))^(1/3)), # smoothing length based on mass/density in units of kpc
                                       "Metallicity" = gas$Metallicity,
                                       "Carbon" = gas$`ElementAbundance/Carbon`,
                                       "Hydrogen" = gas$`ElementAbundance/Hydrogen`,
                                       "Oxygen" =  gas$`ElementAbundance/Oxygen`)
 
+    gas_part$ThermalDispersion[gas_part$Temperature <= 1e4] = 11
+
     remove(gas); remove(PT0_attr)
 
   } else {gas_part=NULL}
@@ -708,13 +719,17 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
                                       "SFR" = gas$StarFormationRate*(.g_to_msol/.s_to_yr),                                        # SFR in Msol/yr
                                       "Density" = gas$Density*.gcm3_to_msolkpc3,                                                  # Density in Msol/kpc^3
                                       "Temperature" = (.adiabatic_index-1)*(gas$InternalEnergy/.Boltzmann_constant)*(4*.mass_of_proton/(1 + gas$`ElementAbundance/Hydrogen`*(3 + 4*gas$ElectronAbundance))),
+                                      "ThermalDispersion" = sqrt((gas$InternalEnergy*.cms_to_kms)*(.adiabatic_index - 1)),
                                       "SmoothingLength" = 2*(((3/(4*pi))*((gas$Mass*.g_to_msol) / (gas$Density*.gcm3_to_msolkpc3)))^(1/3)), # smoothing length based on mass/density in units of kpc
                                       "Metallicity" = gas$Metallicity,
                                       "Carbon" = gas$`ElementAbundance/Carbon`,
                                       "Hydrogen" = gas$`ElementAbundance/Hydrogen`,
                                       "Oxygen" = gas$`ElementAbundance/Oxygen`)
 
+    gas_part$ThermalDispersion[gas_part$Temperature <= 1e4] = 11
+
     remove(gas); remove(PT0_attr)
+
   } else {gas_part=NULL}
 
   if ("PartType4" %in% groups){
@@ -1334,9 +1349,9 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
 }
 
 .sum_gas_velocities = function(galaxy_sample, observation){
-  vel_diff = function(mass, vy){diff((mass * pnorm(observation$vbin_edges, mean = vy, sd = 0)))}
+  vel_diff = function(mass, vy, sigma_t){diff((mass * pnorm(observation$vbin_edges, mean = vy, sd = sigma_t)))}
 
-  bins = mapply(vel_diff, galaxy_sample$Mass, galaxy_sample$vy)
+  bins = mapply(vel_diff, galaxy_sample$Mass, galaxy_sample$vy, (galaxy_sample$ThermalDispersion/sqrt(3)))
 
   return(rowSums(bins))
 
@@ -1872,6 +1887,9 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
       new_gas_part[ind1:ind2, z := part$z+new_xyz[,3],]
       new_gas_part[ind1:ind2, Mass := part$Mass*rand_pos$weight,]
       new_gas_part[ind1:ind2, SFR := part$SFR*rand_pos$weight,]
+      new_gas_part[ind1:ind2, Density := part$Density*rand_pos$weight,]
+      new_gas_part[ind1:ind2, Temperature := part$Temperature*rand_pos$weight,]
+      new_gas_part[ind1:ind2, ThermalDispersion := sqrt((part$ThermalDispersion^2)*rand_pos$weight),]
       # in the new data.frame of particle properties, assign their
       # new positions and masses scaled by the kernel weight.
     }
@@ -1891,7 +1909,7 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
   i = integer()
 
   output = foreach(i = 1:no_gas, .combine='.comb', .multicombine=TRUE,
-                   .init=list(list(), list(), list(), list(), list())) %dopar% {
+                   .init=list(list(), list(), list(), list(), list(), list(), list(), list())) %dopar% {
 
                        part = gas_part[i,]
 
@@ -1904,8 +1922,11 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
                        z = part$z+new_xyz[,3]
                        Mass = part$Mass*rand_pos$weight
                        SFR  = part$SFR*rand_pos$weight
+                       Density = part$Density*rand_pos$weight
+                       Temperature = part$Temperature*rand_pos$weight
+                       ThermalDispersion = sqrt((part$ThermalDispersion^2)*rand_pos$weight)
 
-                       return(list(x, y, z, Mass, SFR))
+                       return(list(x, y, z, Mass, SFR, Density, Temperature, ThermalDispersion))
                        closeAllConnections()
                      }
 
@@ -1915,6 +1936,9 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
   new_gas_part[, z := as.numeric(unlist(output[[3]])),]
   new_gas_part[, Mass := as.numeric(unlist(output[[4]])),]
   new_gas_part[, SFR := as.numeric(unlist(output[[5]])),]
+  new_gas_part[, Density := as.numeric(unlist(output[[6]])),]
+  new_gas_part[, Temperature := as.numeric(unlist(output[[7]])),]
+  new_gas_part[, ThermalDispersion := as.numeric(unlist(output[[8]])),]
 
   return(new_gas_part)
 }

tests/testthat/test_build_datacube.R

@@ -1420,3 +1420,34 @@ test_that("Noise increases as expected in spectral cubes", {
                 sd(gadget_spectra_sn5$spectral_cube[15,15,] - gadget_spectra_nonoise$spectral_cube[15,15,]))
 
 })
+
+# Testing that spawning conserves different properties -----------------------
+test_that("Gas spawning creates even images despite different n-values",{
+
+  ss_eagle_500  = make_simspin_file(ss_pd_eagle, write_to_file = FALSE,
+                                       template = "BC03lr", sph_spawn_n = 500)
+
+  ss_eagle_100 = ss_illustris  = make_simspin_file(ss_pd_eagle,
+                                                       write_to_file = FALSE,
+                                                       template = "BC03lr",
+                                                       sph_spawn_n = 100)
+
+  n500 =  build_datacube(simspin_file = ss_eagle_500,
+                        telescope = telescope(type="IFU", lsf_fwhm = 3.6, signal_to_noise = NA, fov = 25),
+                        observing_strategy = observing_strategy(dist_z = 0.1, inc_deg = 45, blur = F),
+                        method = "gas")
+
+  n100 = build_datacube(simspin_file = ss_eagle_100,
+                        telescope = telescope(type="IFU", lsf_fwhm = 3.6, signal_to_noise = NA, fov = 25),
+                        observing_strategy = observing_strategy(dist_z = 0.1, inc_deg = 45, blur = F),
+                        method = "gas")
+
+
+  mask100 = n100$raw_images$particle_image; mask100[mask100 >= 1] = 1; mask100[mask100!=1]=NA
+  mask500 = n500$raw_images$particle_image; mask500[mask500 >= 1] = 1; mask500[mask500!=1]=NA
+
+  expect_true(sum(n500$raw_images$mass_image*mask500, na.rm=T) == sum(n100$raw_images$mass_image*mask100, na.rm=T))
+  expect_true(sum(n500$raw_images$SFR_image) == sum(n100$raw_images$SFR_image))
+
+})
+

tests/testthat/test_make_simspin_file.R

@@ -34,26 +34,45 @@ test_that("Initial run of each simulation type - HDF5", {
 
 test_that("Initial run of each simulation type - EAGLE", {
   expect_null(make_simspin_file(ss_eagle, output = paste(temp_loc, "/eagle_test", sep=""), centre = c(18318,61583,38667), half_mass = 1483809589))
-  expect_length(readRDS(paste(temp_loc, "/eagle_test", sep="")), ss_file_length)
-  expect_true(length(readRDS(paste(temp_loc, "/eagle_test", sep=""))$gas_part) == 16)
+
+  EAGLE = readRDS(paste(temp_loc, "/eagle_test", sep=""))
+  expect_length(EAGLE, ss_file_length)
+  expect_true(length(EAGLE$gas_part) == 17)
+  expect_false(any(is.na(EAGLE$gas_part$ThermalDispersion)))
+  expect_true(all(EAGLE$gas_part$ThermalDispersion[EAGLE$gas_part$Temperature < 1e4] == 11))
 })
 
 test_that("Initial run of each simulation type - Magneticum", {
   expect_null(make_simspin_file(ss_magneticum, output = paste(temp_loc, "/magneticum_test", sep="")))
-  expect_length(readRDS(paste(temp_loc, "/magneticum_test", sep="")), ss_file_length)
-  expect_true(length(readRDS(paste(temp_loc, "/magneticum_test", sep=""))$gas_part) == 16)
+
+  magneticum = readRDS(paste(temp_loc, "/magneticum_test", sep=""))
+  expect_length(magneticum, ss_file_length)
+  expect_true(length(magneticum$gas_part) == 17)
+  expect_false(any(is.na(magneticum$gas_part$ThermalDispersion)))
+  expect_true(all(magneticum$gas_part$ThermalDispersion[magneticum$gas_part$Temperature < 1e4] == 11))
+
 })
 
 test_that("Initial run of each simulation type - HorizonAGN", {
   expect_null(make_simspin_file(ss_horizon, output = paste(temp_loc, "/horizon_test", sep="")))
-  expect_length(readRDS(paste(temp_loc, "/horizon_test", sep="")), ss_file_length)
-  expect_true(length(readRDS(paste(temp_loc, "/horizon_test", sep=""))$gas_part) == 16)
+
+  HAGN = readRDS(paste(temp_loc, "/horizon_test", sep=""))
+  expect_length(HAGN, ss_file_length)
+  expect_true(length(HAGN$gas_part) == 17)
+  expect_false(any(is.na(HAGN$gas_part$ThermalDispersion)))
+  expect_true(all(HAGN$gas_part$ThermalDispersion[HAGN$gas_part$Temperature < 1e4] == 11))
+
 })
 
 test_that("Initial run of each simulation type - IllustrisTNG", {
   expect_null(make_simspin_file(ss_illustris, output = paste(temp_loc, "/illustris_test", sep="")))
-  expect_length(readRDS(paste(temp_loc, "/illustris_test", sep="")), ss_file_length)
-  expect_true(length(readRDS(paste(temp_loc, "/illustris_test", sep=""))$gas_part) == 16)
+
+  TNG = readRDS(paste(temp_loc, "/illustris_test", sep=""))
+  expect_length(TNG, ss_file_length)
+  expect_true(length(TNG$gas_part) == 17)
+  expect_false(any(is.na(TNG$gas_part$ThermalDispersion)))
+  expect_true(all(TNG$gas_part$ThermalDispersion[TNG$gas_part$Temperature < 1e4] == 11))
+
 })
 
 # Test that the function fails when the file already exists
@@ -305,7 +324,7 @@ test_that("Negative stellar wind particles are removed",{
 # Gas temperature should fall within a reasonable range
 test_that("Temperature does not go outside a reasonable range",{
   illustris = readRDS(paste(temp_loc, "/illustris_test", sep=""))
-  expect_true(min(illustris$gas_part$Temperature)>8000)
+  expect_true(min(illustris$gas_part$Temperature)>0)
   expect_true(max(illustris$gas_part$Temperature)<50000000)
 })
 
@@ -637,3 +656,4 @@ test_that("Expect header information in each of the built simspin files.", {
 
 })
 
+

tests/testthat/test_utilities.R

@@ -122,3 +122,8 @@ test_that("interp_quick fails if given more than a scalar", {
  expect_error(.interp_quick(x = c(1,2,3), params = c(seq(1,10), log = T)))
 })
 
+test_that("kernel weights sum to 1", {
+  expect_true(sum(.generate_uniform_sphere(100, kernel = "WC2")$weight) == 1)
+  expect_true(sum(.generate_uniform_sphere(100, kernel = "WC6")$weight) == 1)
+  expect_true(sum(.generate_uniform_sphere(100, kernel = "M4")$weight) == 1)
+})