Merge branch 'main' into compathelper/new_version/2020-11-18-00-31-45…

…-861-185838806
RvSpectML · Dec 14, 2020 · d79017a · d79017a
2 parents 5629a57 + e9b3526
commit d79017a
Show file tree

Hide file tree

Showing 4 changed files with 140 additions and 2 deletions.
diff --git a/Project.toml b/Project.toml
@@ -31,7 +31,7 @@ QuadGK = "2.4"
 Query = "1"
 RvSpectMLBase = "0.1"
 SpecialFunctions = "0.10, 1.0"
-StaticArrays = "0.12"
+StaticArrays = "0.12, 1.0"
 StatsBase = "0.33"
 ThreadedIterables = "0.2"
 julia = "1.3"
diff --git a/src/EchelleCCFs.jl b/src/EchelleCCFs.jl
@@ -69,4 +69,7 @@ export RVFromCCF
 export measure_rv_from_ccf, measure_rvs_from_ccf
 export MeasureRvFromMinCCF, MeasureRvFromCCFCentroid, MeasureRvFromCCFQuadratic, MeasureRvFromCCFGaussian, MeasureRvFromCCFTemplate
 
+include("io/write_ccf_fits.jl")
+export write_ccf_fits, write_each_ccf_fits
+
 end
diff --git a/src/ccf/ccf_template.jl b/src/ccf/ccf_template.jl
@@ -5,6 +5,29 @@ Author: Eric Ford
 Created: October 2020
 """
 
+""" calc_normalized_ccfs( ccfs )
+Normalizes each spectrum by fitting a line to the region outside the center
+"""
+function fit_ccf_normalizations( v_grid::A1, ccfs::A2; v_mid::Real, dv_min::Real = RvSpectMLBase.max_bc, dv_max::Real = Inf ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2} }
+    @assert zero(dv_min) < dv_min < Inf
+    @assert dv_min < dv_max
+    idx = findall( x->dv_min < abs(x-v_mid) < dv_max , v_grid )
+    @assert length(idx) >= 2
+    map(col-> mean(ccfs[idx,col] ), 1:size(ccfs,2) )
+end
+
+function calc_normalized_ccfs( v_grid::A1, ccfs::A2; v_mid::Real, dv_min::Real = RvSpectMLBase.max_bc, dv_max::Real = Inf ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2} }
+    norms = fit_ccf_normalizations(v_grid,ccfs,v_mid=v_mid,dv_min=dv_min, dv_max=dv_max)
+    ccfs_out = copy(ccfs)
+    for i in 1:size(ccfs,2)
+        ccfs_out[:,i] ./= norms[i]
+    end
+    return ccfs_out
+end
+
+""" calc_normalized_ccfs( ccfs )
+Normalizes each spectrum by its maximum value.
+"""
 function calc_normalized_ccfs( ccfs::A2 ) where {T2<:Real, A2<:AbstractArray{T2,2} }
     return ccfs ./maximum(ccfs,dims=1)
 end
@@ -16,10 +39,14 @@ function calc_normalized_ccfs( ccfs::A2, ccf_vars::A3 ) where { T2<:Real, A2<:Ab
     return (ccfs=ccfs_norm, ccf_vars=ccf_vars_norm)
 end
 
+""" calc_ccf_template( ccfs, [ccf_vars] ; assume_normalized = false )
+Calculates ccf template
+Warning: uses maximum CCF for normalization, unless you normalize manually.
+"""
 function calc_ccf_template( ccfs::A2; assume_normalized::Bool = false ) where {T2<:Real, A2<:AbstractArray{T2,2} }
     ccfs_norm = assume_normalized ? ccfs : calc_normalized_ccfs(ccfs)
     ccf_template = mean(ccfs_norm,dims=2)
-    est_ccf_vars_norm = var(ccfs_norm.-ccf_template,dims=2)
+    est_ccf_vars_norm = var(ccfs_norm.-ccf_template,dims=1)
     ccf_template = vec(sum(ccfs_norm./est_ccf_vars_norm,dims=2) ./sum(1.0 ./est_ccf_vars_norm,dims=2))
     return ccf_template
 end

diff --git a/src/io/write_ccf_fits.jl b/src/io/write_ccf_fits.jl
@@ -0,0 +1,108 @@
+using FITSIO
+using Pkg
+using Dates
+
+function make_ccf_fits_header(metadata::AbstractDict; line_list_filename::Union{Nothing,String} = "unknown",
+            v::Union{Nothing,Float64} = nothing, e_v::Union{Nothing,Float64} = nothing,
+            blue_ord::Union{Nothing,Int64} = nothing, red_ord::Union{Nothing,Int64} = nothing,
+            mask_shape::Union{Nothing,String} = nothing, wave_cal::Union{Nothing,String} = nothing,
+            blazenorm::Union{Nothing,Bool} = nothing, contnorm::Union{Nothing,Bool} = nothing,
+            div_tell::Union{Nothing,Bool} = nothing, rawnorm::Union{Nothing,Bool} = nothing
+            )
+    header_dict = Dict{String,Tuple{Any,String}}("VERSION"=>(string(Pkg.project().version), string(Pkg.project().name) * " code version"))
+    # Add version info for packages used
+    deps = Pkg.dependencies()
+    for (uuid, dep) in deps
+        dep.is_direct_dep || continue
+        dep.version === nothing && continue
+        header_dict["VERSION_" * dep.name] = (string(dep.version), "version of package dependancy")
+    end
+    # header fields similar to EXPRES CCFs
+    header_dict["DATE-CCF"] = (string(Dates.now()), "Date and time of CCF calculation")
+    if hasproperty(metadata,:Filename)
+        header_dict["FILENAME"] = (last(splitpath(metadata[:Filename])), "Original filename")
+    end
+    if hasproperty(metadata,:bjd)
+        header_dict["MJD"] = (metadata[:bjd], "MJD of geometric midpoint")
+    end
+    if !isnothing(line_list_filename)
+        header_dict["MASK"] = (line_list_filename, "CCF mask linelist used")
+    end
+    # Skipped: EXPCOUNT, SNR, V, E_V, CHI2, BLUE_ORD, RED_ORD
+    if !isnothing(mask_shape)
+        header_dict["WINDOW"] = (mask_shape, "Shape of CCF mask used in CCF")
+    end
+    # Skipped BARYCOR
+    if hasproperty(metadata,:airmass)
+        header_dict["AIRMASS"] = (metadata[:airmass], "Airmass of exposure")
+    end
+    # Skipped EXPTIME
+    if hasproperty(metadata,:expres_epoch)
+        header_dict["RVEPOCH"] = (metadata[:expres_epoch], "RV calibration epoch")
+    end
+    if hasproperty(metadata,:wave_cal)
+        header_dict["WAVE_CAL"] = (wave_cal, "Wavelength calibration used when computing CCF")
+    end
+    # Skipped CCOR
+    if !isnothing(blazenorm)
+        header_dict["BLAZE"] = (blazenorm, "Pixels weighted by blaze in CCF")
+    end
+    if !isnothing(contnorm)
+        header_dict["CONTNORM"] = (contnorm, "Continuum normalized spectrum in CCF")
+    end
+    # Skipped SUB_CONT
+    if !isnothing(div_tell)
+        header_dict["DIV_TELL"] = (div_tell, "Tellurics divided before CCF")
+    end
+    # Skipped PLATESCL
+    # Added some that I thought might be useful to have
+    if hasproperty(metadata,:normalization)
+        header_dict["NORMALIZATION"] = (string(metadata[:normalization]), "How spectrum was normalized before computing CCF")
+    end
+    if hasproperty(metadata,:snr_prelim)
+        header_dict["SNR_PRELIM"] = (string(metadata[:snr_prelim]), "Something related to SNR, but different from what EXPRES provided")
+    end
+    if hasproperty(metadata,:pwv)
+        header_dict["PWV"] = (metadata[:pwv], "precip. water vabor")
+    end
+    if hasproperty(metadata,:sundist)
+        header_dict["SUN_DIST"] = (metadata[:sundist], "Distance to sun")
+    end
+    if hasproperty(metadata,:moondist)
+        header_dict["MOON_DIST"] = (metadata[:moondist], "Distance to moon")
+    end
+    fits_header = FITSHeader(collect(keys(header_dict)),first.(values(header_dict)),last.(values(header_dict)))
+end
+
+function write_ccf_fits(filename::String, v_grid::AbstractArray{T1,1}, ccf_total::AbstractArray{T2,1}, ccf_total_var::AbstractArray{T3,1};
+            orders::AbstractArray{T4,1}=zeros(Int64,0), ccf_orders::AbstractArray{T5,2}=zeros(T2,0,0), ccf_orders_var::AbstractArray{T6,2}=zeros(T3,0,0),
+            fits_hdr::FITSHeader #=, order_vels::AbstractArray{T7,1}, order_sigma_vels::AbstractArray{T8,1} =# ) where {
+                T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, T6<:Real }
+    f = FITS(filename, "w")
+    write(f,UInt8[],header=fits_hdr)
+    t2 = Dict("V_grid"=>collect(v_grid),"ccf"=>ccf_total,"e_ccf"=>sqrt.(ccf_total_var) )
+    write(f,t2,name="ccf_total")
+    if length(orders) > 0 && length(ccf_orders) > 0 && length(ccf_orders_var) > 0
+        t3 = Dict("orders"=>convert.(Int64,orders),"ccfs"=>ccf_orders,"errs"=>sqrt.(ccf_orders_var) )
+        #,"v" => order_vels,"e_v" => order_sigma_vels)
+        write(f,t3,name="ccf_order")
+    end
+    close(f)
+end
+
+function write_each_ccf_fits(metadata::AbstractArray{Dict{Symbol,Any},1}, v_grid::AbstractArray{T1,1}, ccf_total::AbstractArray{T2,2}, ccf_total_var::AbstractArray{T3,2};
+            orders::AbstractArray{T4,1}=zeros(Int64,0), ccf_orders::AbstractArray{T5,3}=zeros(T2,0,0,size(ccf_total,2)), ccf_orders_var::AbstractArray{T6,3}=zeros(T3,0,0,size(ccf_total,2)),
+            fits_hdr::FITSHeader #=, order_vels::AbstractArray{T7,1}, order_sigma_vels::AbstractArray{T8,1} =# ) where {
+                T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, T6<:Real }
+    @assert length(metadata) == size(ccf_total,2) == size(ccf_total_var,2)
+    num_files = length(metadata)
+    for i in 1:num_files
+        fits_fn = last(splitpath(metadata[i][:Filename]))
+        ccf_fn = replace(fits_fn, ".fits" => "_ccf.fits")
+        fits_hdr = make_ccf_fits_header(metadata[i])
+        write_ccf_fits(ccf_fn, v_grid, ccf_total[:,i], ccf_total_var[:,i],
+                        orders=orders, ccf_orders=ccf_orders[:,:,i], ccf_orders_var=ccf_orders_var[:,:,i],
+                        fits_hdr=fits_hdr )
+        #break
+    end
+end