Merge pull request #128 from LudwigBoess/RM_fix

Rm fix

Project.toml

@@ -1,7 +1,7 @@
 name = "SPHtoGrid"
 uuid = "8a806948-229a-4117-b79b-c12ac366026f"
 authors = ["LudwigBoess <[email protected]>"]
-version = "0.4.6"
+version = "0.4.7"
 
 [deps]
 Cosmology = "76746363-e552-5dba-9a5a-cef6fa9cc5ab"

docs/src/effects.md

@@ -100,6 +100,16 @@ fo_image = map_path * "B.fits"
 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
+```@docs
+rotation_measure
+```
+For the latter we provide the simple helper function
+```@docs
+get_dz
+```
 
 ## X-Ray emission
 

src/effects/RM.jl

@@ -1,32 +1,65 @@
-# """
-#     rotation_measure(n_cm3::Real, B_par::Real; ν_obs=nothing)
-
-# Computes the rotation measure of the parallel magnetic field along the LOS.
-# """
-# function rotation_measure(n_cm3::Real, B_par::Real; ν_obs=nothing)
-#     RM = 812 * n_cm3 * 1.e3 * B_par * 1.e6
-#     if !isnothing(ν_obs)
-#         RM *= c_light^2 * 1.e-4 / ν_obs
-#     end
-#     return mod(RM, π)
-# end
-
 """
     get_dz(M, rho, hsml)
 
-Helper function to get the particle depth along the LOS.
+Helper function to get the particle depth along the LOS in the given units.
+All additional unit conversion must be performed by hand.
 """
 get_dz(M, rho, hsml) = M / rho / (4 * hsml^2)
 
 """
     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`](@ref)
+- reduce image: `false`
 """
-function rotation_measure(n_cm3::Real, B_los::Real, dz::Real; ν_obs=nothing)
+function rotation_measure(n_cm3::Real, B_los::Real)
+    # RM in rad/m^2
+    return faraday_prefac * n_cm3 * B_los * 100.0^2
+end
+
+
+"""
+    rotation_measure(n_cm3::Real, B_los::Real, dz::Real, ν_obs::Real)
+
+Computes the rotation measure of the parallel magnetic field along the LOS at a given frequency.
+To be used for continuous rotation of polarized emission along the LOS.
+
+## Arguments
+- `n_cm3::Real`: Electron number density in [cm^-3].
+- `B_los::Real`: Magnetic field strength along the LOS in [G].
+- `dz::Real`: Depth along the LOS in [cm]. See [`get_dz`](@ref) for a convenient helper function.
+- `ν_obs::Real`: Observing frequency in [Hz].
+
+## Returns
+- `RM::Real`: Rotation measure in [rad/cm^2].
+
+## Mapping settings
+- weight function: [`part_weight_physical`](@ref)
+- reduce image: `false`
+- stokes: `true`
+- sort_z: `true`
+- `RM`: use output of this function.
+"""
+function rotation_measure(n_cm3::Real, B_los::Real, dz::Real, ν_obs::Real)
+    # RM in rad/cm^2
     RM = faraday_prefac * n_cm3 * B_los
-    if !isnothing(ν_obs)
-        RM *= c_light^2 * dz / ν_obs
-    end
-    return mod(RM, π)
+
+    # for Stokes parameters we need rotation at given frequency
+    RM *= c_light^2 * dz / ν_obs^2
+
+    # store sign of rotation 
+    _sign = sign(RM)
+
+    # reduce to minimum angle
+    return _sign * mod(RM, π)
 end

src/effects/constants.jl

@@ -10,6 +10,7 @@ global const q_e        = 1.602176487e-20 * c_light
 global const yPrefac    = σ_T * k_B / m_e / c_light^2
 global const eV2cgs     = 1.60218e-12
 global const cgs2eV     = 1.0/eV2cgs
+global const kpc        = 3.085678e21
 
 # synchrotron
 const global C_crit = 3q_e / (4π * m_e^3 * c_light^5) # Donnert+16, MNRAS 462, 2014–2032 (2016), Eg. 20 

src/effects/synchrotron_Hoeft.jl

@@ -118,7 +118,7 @@ function analytic_synchrotron_HB07(rho_cgs::Array{<:Real}, m_cgs::Array{<:Real},
 
         if isnothing(dsa_model)
             # interpolate Ψ
-            η_tot = ηB * Ψ[Mach[i], T_keV[i]]
+            η_tot = ηB * Ψ(Mach[i], T_keV[i])
         else
             # evaluate DSA model
             η_tot = ηB * η_Ms_acc(η_model, Mach[i])

test/runtests.jl

@@ -625,6 +625,12 @@ addprocs(2)
             # check if the maximum density error is within 20% of the Gadget density
             @test maximum(L1) < 0.2
         end
+
+        @testset "Faraday Rotation" begin 
+            @test rotation_measure(1.0, 1.0) ≈ SPHtoGrid.faraday_prefac * 1.e4
+
+            @test rotation_measure(1.0, 1.e-6, SPHtoGrid.kpc, 1.4e9) ≈ 2.6705855715515447
+        end
     end
 
     @testset "Image Functions" begin