Fix NaN in radialdblpower

EHTJulia · Feb 14, 2025 · 4d82521 · 4d82521
1 parent c9fc9d4
commit 4d82521
Show file tree

Hide file tree

Showing 2 changed files with 2 additions and 13 deletions.
diff --git a/src/models/imagetemplates/rings.jl b/src/models/imagetemplates/rings.jl
@@ -125,8 +125,8 @@ radialextent(d::RadialDblPower{T}) where {T} = 1 + convert(T, 50^(1 / d.αouter)
 
 @inline @fastmath function radial_profile(d::RadialDblPower, p)
     @unpack_params αinner, αouter = d(p)
-    r = p.r
-    return r^αinner / (1 + r^(αouter + αinner + 1))
+    r = p.r + eps(typeof(p.r)) #avoid singularity at r=0
+    return r^αinner*inv(1 + r^(αinner + αouter + 1)) # we do it like this to avoid underflow
 end
 
 """

diff --git a/src/models/polarized.jl b/src/models/polarized.jl
@@ -322,17 +322,6 @@ where `a,b,c,d` are real numbers with no conditions, and `p=√(a² + b² + c²)
         pimgV[i] = tmp * d[i]
     end
 
-    # p = @.. sqrt(b ^ 2 + c ^ 2 + d ^ 2)
-    # ea = @.. exp(a)
-    # tmp = @.. ea * sinh(p) * inv(p)
-    # pimgI = @.. ea * cosh(p)
-    # pimgQ = @.. tmp * b
-    # pimgU = @.. tmp * c
-    # pimgV = @.. tmp * d
-    # stokesI = pimgI
-    # stokesQ = pimgQ
-    # stokesU = pimgU
-    # stokesV = pimgV
     return stokes_intensitymap(pimgI, pimgQ, pimgU, pimgV, grid)
 end