Merge pull request #54 from kateharborne/dev-ADACS

dev-ADACS

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SimSpin
 Type: Package
 Title: SimSpin - A package for the kinematic analysis of galaxy simulations
-Version: 2.3.7
+Version: 2.3.8
 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.

NEWS.md

@@ -1,14 +1,19 @@
-# SimSpin v2.3.7 News
+# SimSpin v2.3.8 News
 
 ### Author: Kate Harborne
 
-### Last edit: 19/07/22
+### Last edit: 21/07/22
 
 Below is a table containing a summary of all changes made to SimSpin, since the date this file was created on 26/08/2021.
 
 | Date     	| Summary of change                                                                                                                                                                                                                                                                                                                                                                                                    	| Version 	| Commit                                   	|
 |----------	|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------	|---------	|------------------------------------------	|
-| 19/07/22 	| **New feature!** Update to `build_datacube`. Observed images produced with method = 'velocity' now include measured h3 and h4 higher order kinematics. This is done by fitting the LOSVD within the output velocity cube with a Gauss-Hermite polynomial. These images are now included and written to single or split FITS.                                                                                         	| 2.3.7   	|                                          	|
+| 21/07/22 	| *Minor bugfix*. With update to v2.3.5, flux conservation was implemented in mode="velocity". During this upgrade, filters were properly implemented meaning that re-run observations will now produce flux maps = 0 where previously they had flux. This occurs when the filter wavelength range does not overlap with the observed telescope wavelength range.                                                      	| 2.3.8   	|                                          	|
+|          	| (1) Errors now issued by `telescope` if you request a filter that will not allow any flux through the telescope wavelength range.                                                                                                                                                                                                                                                                                    	|         	|                                          	|
+|          	| (2) Default filter has been changed in `telescope` to `filter = "g"` so that there is overlap between the default parameters.                                                                                                                                                                                                                                                                                        	|         	|                                          	|
+|          	| (3) Errors issued if the requested `half_mass` in `make_simspin_file` is greater than the total mass contained in the input simulation file OR if smaller than the mass of a single particle.                                                                                                                                                                                                                        	|         	|                                          	|
+|          	| (4) All images are now written to FITS (not just a selection dependent on the method selected). Default names for files now include a prefix "raw_" or "obs_" to denote which list the image comes from within the build_datacube output.                                                                                                                                                                            	|         	|                                          	|
+| 19/07/22 	| **New feature!** Update to `build_datacube`. Observed images produced with method = 'velocity' now include measured h3 and h4 higher order kinematics. This is done by fitting the LOSVD within the output velocity cube with a Gauss-Hermite polynomial. These images are now included and written to single or split FITS.                                                                                         	| 2.3.7   	| 2dbca740c29637f8802f1c92ce48ca0b215f5458 	|
 | 18/07/22 	| *Minor bugfix*. Update to `make_simspin_file`. Modifying how gas particles are smoothed for HorizonAGN simulations - rather than a uniform smoothing across the cell cube, we now assign a smoothing length dependent on the mass/density and smooth with a cubic spline.                                                                                                                                            	| 2.3.6   	| fcd915985e6ad477538c7f4be3c08428d0fd7fe1 	|
 | 18/07/22 	| *Minor bugfix*. Update to `build_datacube`. Implementing enforced flux conservation across the telescope wavelength range between input template and output observed spectrum. Fix for both spectral and velocity cubes. Tests added to ensure that outputs in parallel and series are identical.                                                                                                                    	| 2.3.5   	| bf4556234127c7352604863b24b4649c9fda5341 	|
 | 12/07/22 	| *Minor bugfix*. Update to `build_datacube`. Returning a meaningful error if the requested `centre` position is outside the range of the provided input file. Tests added to ensure behaviour now works as expected.                                                                                                                                                                                                  	| 2.3.4   	| 30a30eb4758a87edbeb7345f1d62d3b256ad5ff2 	|

R/telescope.R

@@ -45,7 +45,7 @@
 #'
 
 telescope = function(type="IFU", fov, aperture_shape="circular", wave_range=c(3700,5700),
-                     wave_centre, wave_res=1.04, spatial_res, filter="r", lsf_fwhm=2.65,
+                     wave_centre, wave_res=1.04, spatial_res, filter="g", lsf_fwhm=2.65,
                      signal_to_noise = 10, method){
 
   if (!missing(method)){
@@ -396,6 +396,10 @@ telescope = function(type="IFU", fov, aperture_shape="circular", wave_range=c(37
     }
   }
 
+  if( min(output$filter$wave) > max(output$wave_range) | max(output$filter$wave) < min(output$wave_range) ){
+    stop("Error: Requested filter will not overlap with the telescope wavelength range. \n
+         Please select a different filter or extend your telescope wavelength range.")
+  }
 
   return(output)
 }

R/utilities.R

@@ -701,6 +701,7 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
   if (is.na(half_mass)){
     half_mass = sum(galaxy_data$Mass) / 2
   }
+
   ellip_radius = sqrt((x*x) + ((y/p)*(y/p)) + ((z/q)*(z/q)))
 
   int_order = order(ellip_radius) # get the indicies of the radii in order (low to high)
@@ -709,6 +710,7 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
   half_mass_ind = which(abs(cum_mass - half_mass) == min(abs(cum_mass - half_mass)))[1] # at what radius does this half-mass now occur?
 
   return(ordered_galaxy_data[1:half_mass_ind,])
+
 }
 
 .ellipsoid_tensor = function(galaxy_data, p, q){
@@ -825,7 +827,22 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
 
 # Function to align full galaxy based on the stellar particles
 .align_galaxy = function(galaxy_data, half_mass=NA){
+
   if (is.null(galaxy_data$star_part)){ # if there are no stellar particles (just gas), use these
+
+    if(!is.na(half_mass[1])){
+      check_bounds = (sum(galaxy_data$gas_part$Mass) > half_mass) & # check that the total model contains more mass than the requested half-mass
+        (min(galaxy_data$gas_part$Mass) <= half_mass) # check that the requested half mass is greater than a single particle
+
+      if (!check_bounds){
+        stop(paste0("Error: Requested half-mass is outside the range possible within the input simulation. \n",
+                    "Minimum mass of particle = ", min(galaxy_data$gas_part$Mass), "\n",
+                    "Maximum mass of simulation = ", sum(galaxy_data$gas_part$Mass), "\n",
+                    "Requested half mass ", half_mass, " is outside this range. \n",
+                    "Please resubmit your request with an appropriate value OR with NA for self-fit half-mass."))
+      }
+    }
+
     dummy_data = list(star_part = galaxy_data$gas_part,
                       gas_part= galaxy_data$star_part,
                       head = galaxy_data$head,
@@ -837,6 +854,20 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
                             head      = dummy$galaxy_data$head,
                             ssp       = dummy$galaxy_data$ssp)
   } else {
+
+    if(!is.na(half_mass[1])){
+      check_bounds = (sum(galaxy_data$star_part$Mass) > half_mass) & # check that the total model contains more mass than the requested half-mass
+                     (min(galaxy_data$star_part$Mass) < half_mass) # check that the requested half mass is greater than a single particle
+
+      if (!check_bounds){
+        stop(paste0("Error: Requested half-mass is outside the range possible within the input simulation. \n",
+                    "Minimum mass of particle = ", min(galaxy_data$star_part$Mass), "\n",
+                    "Maximum mass of simulation = ", sum(galaxy_data$star_part$Mass), "\n",
+                    "Requested half mass ", half_mass, " is outside this range. \n",
+                    "Please resubmit your request with an appropriate value OR with NA for self-fit half-mass."))
+      }
+    }
+
     data = .measure_pqj(galaxy_data, half_mass)
   }
   return(data$galaxy_data)

R/write_simspin_FITS.R

@@ -232,14 +232,14 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
                              "CUNIT2"="Units of coordinate increment and value",
                              "EXTNAME"="Image extension name")
 
-    extnames = c("FLUX", "RAW_VEL", "RAW_DISP", "NPART")
+    extnames = c("RAW_FLUX", "RAW_VEL", "RAW_DISP", "NPART")
     bunits = c("erg/s/cm**2", "km/s", "km/s", "Particle number")
     image_names = c("flux_image", "velocity_image", "dispersion_image", "particle_image")
     extnum = c(3,4,5,6)
-    output_image_file_names = paste0(output_dir, "/", output_file_root, "_", image_names, ".FITS")
+    output_image_file_names = paste0(output_dir, "/", output_file_root, "_raw_", image_names, ".FITS")
 
     if (split_save){ # if writing each image to a seperate file
-      for (i in 1:4){ # for each image in the build_datacube output,
+      for (i in 1:length(extnum)){ # for each image in the build_datacube output,
 
         # 1. Write the header again to the new file to HDU 1
         Rfits::Rfits_write_header(filename = output_image_file_names[i], keyvalues = header_keyvalues,
@@ -256,7 +256,7 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
                                  create_ext = TRUE, create_file = FALSE, overwrite_file = FALSE)
       }
     } else { # if writing all images to the same file
-      for (i in 1:4){ # for each image in the build_datacube output,
+      for (i in 1:length(extnum)){ # for each image in the build_datacube output,
 
         image_keyvalues$BUNIT = bunits[i]
         image_keyvalues$EXTNAME = extnames[i]
@@ -328,28 +328,31 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
                              "EXTNAME"="Image extension name")
 
     extnames    = if ("flux_image" %in% names(simspin_datacube$raw_images)){
-                     c("FLUX", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RAW_AGE", "RAW_Z", "NPART")
+                     c("OBS_FLUX", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RAW_FLUX", "RAW_VEL", "RAW_DISP", "RAW_AGE", "RAW_Z", "NPART")
                   } else {
-                     c("MASS", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RAW_AGE", "RAW_Z", "NPART")
+                     c("OBS_MASS", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_AGE", "RAW_Z", "NPART")
                   }
 
     bunits      = if ("flux_image" %in% names(simspin_datacube$raw_images)){
-                   c("erg/s/cm**2", "km/s", "km/s", "km/s", "km/s", "Gyr", "Z_solar", "Particle number")
+                   c("erg/s/cm**2", "km/s", "km/s", "unitless", "unitless", "erg/s/cm**2", "km/s", "km/s", "Gyr", "Z_solar", "Particle number")
                   } else {
-                     c("Msol", "km/s", "km/s", "km/s", "km/s", "Gyr", "Z_solar", "Particle number")
+                     c("Msol", "km/s", "km/s", "unitless", "unitless", "Msol", "km/s", "km/s", "Gyr", "Z_solar", "Particle number")
                   }
 
     image_names = if ("flux_image" %in% names(simspin_datacube$raw_images)){
-                    c("flux_image", "velocity_image", "dispersion_image", "h3_image", "h4_image", "age_image", "metallicity_image", "particle_image")
+                    c("flux_image", "velocity_image", "dispersion_image", "h3_image", "h4_image",
+                      "flux_image", "velocity_image", "dispersion_image", "age_image", "metallicity_image", "particle_image")
                   } else {
-                    c("mass_image", "velocity_image", "dispersion_image", "h3_image", "h4_image", "age_image", "metallicity_image", "particle_image")
+                    c("mass_image", "velocity_image", "dispersion_image", "h3_image", "h4_image",
+                      "mass_image", "velocity_image", "dispersion_image", "age_image", "metallicity_image", "particle_image")
                   }
 
-    output_image_file_names = paste0(output_dir, "/", output_file_root, "_", image_names, ".FITS")
-    extnum = c(3,4,5,6,7,8,9,10)
+    rawobs = c("obs", "obs", "obs", "obs", "obs", "raw", "raw", "raw", "raw", "raw", "raw")
+    output_image_file_names = paste0(output_dir, "/", output_file_root, "_", rawobs, "_", image_names, ".FITS")
+    extnum = c(3,4,5,6,7,8,9,10,11,12,13)
 
     if (split_save){ # if writing each image to a seperate file
-      for (i in 1:8){ # for each image in the build_datacube output,
+      for (i in 1:length(extnum)){ # for each image in the build_datacube output,
 
         # 1. Write the header again to the new file to HDU 1
         Rfits::Rfits_write_header(filename = output_image_file_names[i], keyvalues = header_keyvalues,
@@ -374,7 +377,7 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
 
       }
     } else { # if writing all images to the same file
-      for (i in 1:8){ # for each image in the build_datacube output,
+      for (i in 1:length(extnum)){ # for each image in the build_datacube output,
 
         image_keyvalues$BUNIT = bunits[i]
         image_keyvalues$EXTNAME = extnames[i]

README.md

@@ -10,17 +10,15 @@
 
 <p>&nbsp;</p>
 
-SimSpin v2.3.7 - A package for producing mock observations of particle simulations
+v2.3.8 - A package for producing mock observations:
 
-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.
-
-A mock spectral data cube can be produced using the functions in this package. This is a simple process comprised of three steps:
+SimSpin allows you to take a simulation of a galaxy and produce a data cube in the style of an Integral Field Spectroscopy (IFS) instrument. A mock observation is produced using three simple steps:
 
 1.  Read in your particle data and produce the relevant spectra using the `make_simspin_file` function.
 2.  Setup the observation by defining your `telescope` and `observing_strategy`.
 3.  Build your data cube using the `build_datacube`.
 
-From this cube, "observables" can be measured using observational pipelines. Specifically, the observed line-of-sight velocities and dispersions. This package, once installed, is fully documented and tested.
+The data cube produced can be saved as a FITS file for ease of comparison with real observations. This package, once installed, is fully documented and tested.
 
 Another implementation of this code (SimSpin v1.1.3) written in Julia is also available at [SimSpin.jl](https://github.com/kateharborne/SimSpin.jl) developed by [Gerry Gralton](https://github.com/gerrygralton).