Quadfit

deps/download_soap_example.jl

@@ -3,7 +3,7 @@
 
 include("download.jl")
 
-download_url = "https://psu.box.com/shared/static/nikva44r7xai9h0ulfj2pxewzv8ez4fo.h5"
+download_url = "https://zenodo.org/record/6077993/files/soap_demo.h5?download=1"
 download_filename = joinpath("..","data","spectra","soap_demo.h5")
 download_md5 = "18a6291039d7684ef9187762e4b2dfd5"
 

src/calc_rv/fit_quadratic_to_ccf.jl

@@ -6,8 +6,6 @@ Based on code by Alex Wise (aw@psu.edu)
 Refactors & optimized by Eric Ford
 """
 
-
-
 """  Functor to estimate RV based on fitting quadratic near minimum of CCF.
 TODO: Revist the logic here and see if need to perform transformation first.
 """
@@ -40,17 +38,9 @@ function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2 ) where {T1<:Real, A
     X[:,2] .= view(vels,inds) .-mean_v
     X[:,3] .= (view(vels,inds) .-mean_v).^2
     (c, b, a)  = (X'*X) \ (X'*view(ccf,inds))
-    #=
-    # do the polyfit
-    pfit = Polynomials.fit(vels[inds], ccf[inds], 2)
-    @assert length(Polynomials.coeffs(pfit)) >= 3   # just in case fails to fit a quadratic
 
-    # get center from coeffs
-    c, b, a = Polynomials.coeffs(pfit)
-    @assert a>0
-    =#
     v_at_min_of_quadratic = -b/(2*a) + mean_v
-    return ( rv=v_at_min_of_quadratic, σ_rv=NaN )
+    return (rv=v_at_min_of_quadratic, σ_rv=NaN)
 end
 
 function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2, ccf_var::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1} }
@@ -68,13 +58,5 @@ function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2, ccf_var::A2 ) where
     v_at_min_of_quadratic = -b/(2*a)
     sigma_rv = abs(v_at_min_of_quadratic) * sqrt(covar_beta[2,2]/b^2+covar_beta[3,3]/a^2)
     v_at_min_of_quadratic += mean_v
-    #=
-    # do the polyfit
-    pfit = Polynomials.fit(vels[inds], ccf[inds], 2)
-    @assert length(Polynomials.coeffs(pfit)) >= 3   # just in case fails to fit a quadratic
-
-    # get center from coeffs
-    c, b, a = Polynomials.coeffs(pfit)
-    =#
-    return ( rv=v_at_min_of_quadratic, σ_rv=sigma_rv  )
+    return (rv=v_at_min_of_quadratic, σ_rv=sigma_rv)
 end

test/Project.toml

@@ -2,4 +2,5 @@
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 RvSpectMLBase = "c48404b2-35ea-40e7-ac7f-06a53de703d6"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/calc_rv.jl

@@ -1,6 +1,6 @@
 using EchelleCCFs
 using Test
-
+import Statistics: mean
 
 @testset "Check CCF accuracy" begin  # TODO repeat for multiple mask shapes
     @testset "Tophat mask w/ 1 line" begin
@@ -135,4 +135,52 @@ using Test
         @test vfit.rv ≈ 0  atol = 500   # TODO: tune tolerance better
         =#
     end
+    @testset "RV injection/extraction" begin
+        # velocities
+        amp = 1.0 # meters per second
+        vel = amp .* sin.(range(0, 2π, length=100))
+        mid = 5434.5
+        dep = 0.75
+        wid = 0.1
+        mid_shift = mid .* EchelleCCFs.calc_doppler_factor.(vel)
+
+        # compute gaussian line spectrum
+        buffer = 0.75
+        resolution = 7e5
+        Δlnλ = (1.0 / resolution)
+        wavs = exp.(range(log(mid - buffer), log(mid + buffer), step=Δlnλ))
+
+        flux = zeros(length(wavs), length(vel))
+        for i in eachindex(vel)
+            flux[:,i] .= @. 1.0 - dep * exp(-((wavs-mid_shift[i])/wid)^2 / 2.0)
+        end
+
+        # measure ccf
+        line_list = EchelleCCFs.BasicLineList([mid], [dep])
+        speed_of_light = EchelleCCFs.speed_of_light_mps
+        mask_width = speed_of_light/resolution
+        mask_shape = TopHatCCFMask(mask_width)
+        ccf_plan = BasicCCFPlan(line_list=line_list, mask_shape=mask_shape)
+        v_grid = EchelleCCFs.calc_ccf_v_grid(ccf_plan)
+        ccf = EchelleCCFs.ccf_1D(wavs, flux, ccf_plan)
+
+        # measure velocities from CCF
+        mrv_gauss = MeasureRvFromCCFGaussian(frac_of_width_to_fit=0.5)
+        mrv_quad = MeasureRvFromCCFQuadratic(frac_of_width_to_fit=0.1)
+
+        rvs_gauss = zeros(length(vel))
+        sig_gauss = zeros(length(vel))
+        rvs_quad = zeros(length(vel))
+        sig_quad = zeros(length(vel))
+        for i in eachindex(vel)
+            rvs_gauss[i], sig_gauss[i] = mrv_gauss(v_grid, view(ccf, :, i))
+            rvs_quad[i], sig_quad[i] = mrv_quad(v_grid, view(ccf, :, i))
+        end
+        err_gauss = abs.(rvs_gauss .- mean(rvs_gauss) .- vel)
+        err_quad = abs.(rvs_quad .- mean(rvs_quad) .- vel)
+
+        # test maximum errors are less than a cm/s
+        @assert maximum(err_gauss) < 1e-2
+        @assert maximum(err_quad) < 1e-2
+    end
 end