Examples

You can find some examples here for common maps. These examples assume that you also have GadgetIO.jl, GadgetUnits.jl and SPHKernels.jl installed.

In all examples snap_base points to the snapshot you want to map and map_path points to the folder where you want to store the fits files.

Read data and convert units

First we need to read the data and convert it to physical units (unless you want to bother with comoving coordinates and h0 in cosmological simulations).

Here halo_pos is the position of the halo you want to map and rvir is its virial radius.

# read the header
h = read_header(snap_base)
# define a unit conversion struct from `GadgetUnits.jl`
GU = GadgetPhysical(h)

# define the blocks you want to read
blocks = ["POS", "VEL", "HSML", "RHO", "U", "MASS", "BFLD", "MACH"]

# read all particles in a cubic volume around the halo
data = read_particles_in_volume(snap_base, blocks, halo_pos, rvir)

# convert to physical code units for mapping
pos  = data["POS"]  .* GU.x_physical
hsml = data["HSML"] .* GU.x_physical
rho  = data["RHO"]  .* GU.rho_physical
mass = data["MASS"] .* GU.m_physical

# we want to map the entire cube we read
xy_size = 2rvir
z_size  = 2rvir

# define mapping parameters and convert to physical units
par = mappingParameters(center = halo_pos .* GU.x_physical,
                        x_size = xy_size * GU.x_physical,
                        y_size = xy_size * GU.x_physical,
                        z_size = z_size * GU.x_physical,
                        Npixels = 1024,
                        boxsize=h.boxsize)

# define the kernel you want to use for mapping
k = WendlandC4(Float64, 2)

2D surface density

To get the projected 2D surface density from your data you can use

# you need to make a copy of the positions if you plan to re-use them
# the mapping shifts them in place
pos_map = copy(pos)

# convert density to physical cgs units
rho_gcm3 = data["RHO"] .* GU.rho_cgs

# to get the integrated values along the LOS you need physical weights and not reduce the image
weights = part_weight_physical(length(hsml), par, GU.x_cgs)

# actual mapping
map = sphMapping(pos_map, hsml, mass, rho,
                rho_gcm3, weights,
                param = par, kernel = k,
                reduce_image = false)

# filename of the output image
fo_image = map_path * "rho.fits"

# store the fits image
write_fits_image(fo_image, quantitiy_map, par, snap = snap, units = "g/cm^2")

The units in this case are of course $g/cm^2$.

Magnetic Field

To map the mean magnetic field along the LOS you need to use the density as weight (this is the default behaviour, so you can leave the arguemnt empty) and set reduce_image=true

# you need to make a copy of the positions if you plan to re-use them
# the mapping shifts them in place
pos_map = copy(pos)

# compute absolute value of magnetic field in muG
B = @. 1.e6 * √(data["BFLD"][1, :]^2 + data["BFLD"][2, :]^2 + data["BFLD"][3, :]^2)


# actual mapping
map = sphMapping(pos_map, hsml, mass, rho, B,
                param = par, kernel = k,
                reduce_image = true)

# filename of the output image
fo_image = map_path * "B.fits"

# store the fits image
write_fits_image(fo_image, quantitiy_map, par, snap = snap, units = "muG")

Faraday Rotation

You can map the Faraday Rotation due to the LOS magnetic field strength for RM maps, or for continuous rotation of polarized synchrotron emission using

SPHtoGrid.rotation_measure — Function

rotation_measure(n_cm3::Real, B_los::Real, dz::Real; ν_obs=nothing)

Computes the rotation measure of the parallel magnetic field along the LOS.

Arguments

n_cm3::Real: Electron number density in [cm^-3].
B_los::Real: Magnetic field strength along the LOS in [G].

Returns

RM::Real: Rotation measure in [rad/m^2].

Mapping settings

weight function: part_weight_physical
reduce image: false