Merge pull request #15 from JuliaAstro/ml/bug-fixes

fix airy amplitude scaling

Project.toml

@@ -1,7 +1,7 @@
 name = "PSFModels"
 uuid = "9ba017d1-7760-46cd-84a3-1e79e9ae9ddc"
 authors = ["Miles Lucas <mdlucas@hawaii.edu> and contributors"]
-version = "0.7.1"
+version = "0.8.0"
 
 [deps]
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"

src/airy.jl

@@ -43,14 +43,12 @@ function _airydisk(px, py, x, y, fwhm, ratio, amp, theta, background)
     # unnormalized airydisk
     r = sqrt(dx^2 + dy^2)
     # short-circuit
-    iszero(r) && return amp / (1 - ratio^2)^2 + background
+    iszero(r) && return amp + background
     q = AIRY_PRE * r / fwhm
     I1 = 2 * besselj1(q) / q
-    if iszero(ratio)
-        return amp * I1^2 + background
-    end
+    iszero(ratio) && return amp * I1^2 + background
     I2 = 2 * ratio * besselj1(q * ratio) / q
-    return amp * ((I1 - I2) / (1 - ratio^2))^2 + background
+    return amp * (I1 - I2)^2 + background
 end
 
 # bivariate
@@ -66,11 +64,9 @@ function _airydisk(px, py, x, y, fwhm::BivariateLike, ratio, amp, theta, backgro
     fwhmx, fwhmy = fwhm
     q = AIRY_PRE * sqrt((dx / fwhmx)^2 + (dy / fwhmy)^2)
     # short-circuit
-    iszero(q) && return amp / (1 - ratio^2)^2 + background
+    iszero(q) && return amp + background
     I1 = 2 * besselj1(q) / q
-    if iszero(ratio)
-        return amp * I1^2 + background
-    end
+    iszero(ratio) && return amp * I1^2 + background
     I2 = 2 * ratio * besselj1(q * ratio) / q
-    return amp * ((I1 - I2) / (1 - ratio^2))^2 + background
+    return amp * (I1 - I2)^2+ background
 end

src/moffat.jl

@@ -23,11 +23,15 @@ function _moffat(px, py, x, y, fwhm, alpha, amp, theta, background)
     # rotate
     !iszero(theta) && @warn "isotropic moffat is not affected by non-zero rotation angle $theta"
     # unnormalized moffat
-    dist = dx^2 + dy^2
-    dist /= (fwhm / 2)^2
+    gamma = _moffat_fwhm_to_gamma(fwhm, alpha)
+    dist = (dx / gamma)^2 + (dy / gamma)^2
     return amp / (1 + dist)^alpha + background
 end
 
+# http://openafox.com/science/peak-function-derivations.html#moffat
+_moffat_gamma_to_fwhm(gamma, alpha) = 2 * gamma * sqrt(2^(1/alpha) - 1)
+_moffat_fwhm_to_gamma(fwhm, alpha) = fwhm / (2 * sqrt(2^(1/alpha) - 1))
+
 # bivariate
 function _moffat(px, py, x, y, fwhm::BivariateLike, alpha, amp, theta, background)
     # find offset from center
@@ -38,8 +42,8 @@ function _moffat(px, py, x, y, fwhm::BivariateLike, alpha, amp, theta, backgroun
         dx, dy = rotate_point(dx, dy, theta)
     end
     # unnormalized moffat
-    fwhmx, fwhmy = fwhm
+    gammax, gammay = _moffat_fwhm_to_gamma.(fwhm, alpha)
     # unnormalized moffat
-    dist = (dx / (fwhmx/2))^2 + (dy / (fwhmy/2))^2
+    dist = (dx / gammax)^2 + (dy / gammay)^2
     return amp / (1 + dist)^alpha + background
 end

test/runtests.jl

@@ -24,6 +24,12 @@ function test_model_interface(K)
     val_dir = @inferred K(12, 13; x=12, y=13, fwhm=10)
     @test m(12, 13) ≈ val_dir ≈ 1
 
+    # test position off pixel grid
+    m = K(x=12.5, y=13.5, fwhm=10)
+    val_dir = @inferred K(12.5, 13.5; x=12.5, y=13.5, fwhm=10)
+    @test m(12.5, 13.5) ≈ val_dir ≈ 1
+
+
     # test diagonal fwhm
     m = K(x=0, y=0, fwhm=(10, 9))
     @test m(3, 5) ≈ K(3, 5; x=0, y=0, fwhm=(10, 9)) ≈ K(3, 5; x=0, y=0, fwhm=SA[10, 9])
@@ -103,16 +109,16 @@ end
     @test m(0, -r2) ≈ 0 atol=1e-5
 
     # test with ratio
-    mratio = airydisk(x=0, y=0, fwhm=10, ratio=sqrt(0.5))
+    mratio = airydisk(x=0, y=0, fwhm=10, amp=1, ratio=sqrt(0.5))
     # test attenuation
-    @test mratio(0, 0) ≈ 4
+    @test mratio(0, 0) ≈ 1
     @test mratio(radius, 0) > m(radius, 0)
     @test mratio(-radius, 0) > m(-radius, 0)
     @test mratio(0, radius) > m(0, radius)
     @test mratio(0, -radius) > m(0, -radius)
 
 
-    mratio = airydisk(x=0, y=0, fwhm=(10, 6), ratio=sqrt(0.5))
+    mratio = airydisk(x=0, y=0, fwhm=(10, 6), ratio=sqrt(0.5), amp=4)
     r1 = fwhm[1] * 1.18677
     r2 = fwhm[2] * 1.18677
     # test attenuation
@@ -121,6 +127,10 @@ end
     @test mratio(-r1, 0) > m(-r1, 0)
     @test mratio(0, r2) > m(0, r2)
     @test mratio(0, -r2) > m(0, -r2)
+
+    # https://github.com/JuliaAstro/PSFModels.jl/issues/14
+    mratio = airydisk(x = 40.5, y=40.5, fwhm=10, ratio=0.2,amp=1)
+    @test mratio(40.5, 40.5) ≈ 1
 end
 
 @testset "moffat" begin
@@ -135,8 +145,15 @@ end
 
     # different alpha
     m = moffat(x=0, y=0, fwhm=10, alpha=2)
-    expected = inv(1 + sum(abs2, SA[1, 2]) / 25)^2
+    expected = inv(1 + sum(abs2, SA[1, 2] ./  PSFModels._moffat_fwhm_to_gamma(10, 2)))^2
     @test m(1, 2) ≈ expected
+
+    fwhms = randn(rng, 100) .+ 10
+    for alpha in [0.5, 1.0, 1.5]
+        gammas = PSFModels._moffat_fwhm_to_gamma.(fwhms, alpha)
+        fwhm_calc = PSFModels._moffat_gamma_to_fwhm.(gammas, alpha)
+        @test fwhms ≈ fwhm_calc
+    end
 end
 
 include("plotting.jl")