diff --git a/Project.toml b/Project.toml index 92e0228..efc9cf0 100644 --- a/Project.toml +++ b/Project.toml @@ -16,8 +16,10 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" ResumableFunctions = "c5292f4c-5179-55e1-98c5-05642aab7184" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -35,8 +37,10 @@ JuliaFormatter = "1.0.62" LinearAlgebra = "1.11.0" Logging = "1.11.0" NaNMath = "0.3, 1" +Plots = "1.40.17" ProgressBars = "1.4" Random = "1.11.0" +RecipesBase = "1.3.4" ResumableFunctions = "0.6" StatsBase = "0.33" julia = "1.11" diff --git a/src/Stingray.jl b/src/Stingray.jl index 533cecd..98a3787 100644 --- a/src/Stingray.jl +++ b/src/Stingray.jl @@ -6,6 +6,7 @@ using ProgressBars: tqdm as show_progress using DocStringExtensions using LinearAlgebra using Random +using RecipesBase include("fourier.jl") export positive_fft_bins @@ -75,6 +76,10 @@ export apply_gtis export fill_bad_time_intervals! export create_filtered_lightcurve export check_gtis -export split_by_gtis +export split_by_gtis, intersect_gtis,get_gti_lengths +include("plotting/plots_recipes_lightcurve.jl") +export create_segments +include("plotting/plots_recipes_gti.jl") +export BTIAnalysisPlot end diff --git a/src/events.jl b/src/events.jl index 143f897..8036722 100644 --- a/src/events.jl +++ b/src/events.jl @@ -43,9 +43,12 @@ mutable struct FITSMetadata{H} end function Base.show(io::IO, ::MIME"text/plain", m::FITSMetadata) - println(io, "FITSMetadata for $(basename(m.filepath))[$(m.hdu)] with $(length(m.extra_columns)) extra column(s)") + println( + io, + "FITSMetadata for $(basename(m.filepath))[$(m.hdu)] with $(length(m.extra_columns)) extra column(s)", + ) if !isnothing(m.gti) - gti_exposure = sum(diff(m.gti; dims=2)) + gti_exposure = sum(diff(m.gti; dims = 2)) println(io, "GTI: $(size(m.gti, 1)) intervals, total exposure: $(gti_exposure) s") end if !isnothing(m.mjd_ref) @@ -124,7 +127,7 @@ function EventList(times::Vector{T}, energies::Union{Nothing,Vector{T}} = nothin nothing, # time_unit nothing, # time_sys nothing, # time_pixr - nothing # time_del + nothing, # time_del ) EventList(times, energies, dummy_meta) end @@ -278,7 +281,7 @@ println("Live time fraction: \$(exposure / time_span)") """ function gti_exposure(ev::EventList) if has_gti(ev) - return sum(diff(ev.meta.gti; dims=2)) + return sum(diff(ev.meta.gti; dims = 2)) else return isempty(ev.times) ? 0.0 : maximum(ev.times) - minimum(ev.times) end @@ -321,9 +324,10 @@ function gti_info(ev::EventList) @warn "No GTI information available" return end - + gti_matrix = ev.meta.gti - @debug "GTI Information" source=ev.meta.gti_source intervals=size(gti_matrix, 1) exposure=gti_exposure(ev) time_range=(minimum(gti_matrix), maximum(gti_matrix)) + @debug "GTI Information" source = ev.meta.gti_source intervals = size(gti_matrix, 1) exposure = + gti_exposure(ev) time_range = (minimum(gti_matrix), maximum(gti_matrix)) end # ============================================================================ @@ -418,11 +422,16 @@ by applying the same filtering mask derived from the source column. function filter_on!(f, src_col::AbstractVector, ev::EventList) @assert size(src_col) == size(ev.times) "Source column size must match times size" + # Handle empty input - return immediately + if isempty(src_col) || isempty(ev.times) + return ev + end + # Modified from Base.filter! implementation for multiple arrays # Use two pointers: i for reading, j for writing j = firstindex(ev.times) - for i in eachindex(ev.times) + for i in eachindex(src_col) # Use src_col indices, not ev.times predicate = f(src_col[i])::Bool if predicate @@ -442,18 +451,18 @@ function filter_on!(f, src_col::AbstractVector, ev::EventList) end end - # Resize all arrays to new length - if j <= lastindex(ev.times) - new_length = j - 1 - resize!(ev.times, new_length) + # Calculate new length correctly + new_length = j - firstindex(ev.times) + + # Resize all arrays to new length (handles 0 length correctly) + resize!(ev.times, new_length) - if !isnothing(ev.energies) - resize!(ev.energies, new_length) - end + if !isnothing(ev.energies) + resize!(ev.energies, new_length) + end - for (_, col) in ev.meta.extra_columns - resize!(col, new_length) - end + for (_, col) in ev.meta.extra_columns + resize!(col, new_length) end ev @@ -673,14 +682,16 @@ end - Errors in individual HDUs are logged but don't stop the search - Auto-detection looks for any HDU with "START" and "STOP" columns """ -function read_gti_from_fits(fits_file::String; - gti_hdu_candidates::Vector{String} = ["GTI", "STDGTI"], - gti_hdu_indices::Union{Vector{Int}, Nothing} = nothing, - combine_gtis::Bool = true)::Tuple{Union{Nothing,Matrix{Float64}}, Union{Nothing,String}} - +function read_gti_from_fits( + fits_file::String; + gti_hdu_candidates::Vector{String} = ["GTI", "STDGTI"], + gti_hdu_indices::Union{Vector{Int},Nothing} = nothing, + combine_gtis::Bool = true, +)::Tuple{Union{Nothing,Matrix{Float64}},Union{Nothing,String}} + all_gtis = Matrix{Float64}[] gti_sources = String[] - + FITS(fits_file, "r") do f for gti_name in gti_hdu_candidates try @@ -700,7 +711,7 @@ function read_gti_from_fits(fits_file::String; else gti_hdu_indices end - + for gti_idx in hdu_indices try if length(f) >= gti_idx @@ -715,21 +726,21 @@ function read_gti_from_fits(fits_file::String; end end end - + if isempty(all_gtis) return nothing, nothing end - + if combine_gtis && length(all_gtis) > 1 # Combine all GTIs into one matrix combined_gti = vcat(all_gtis...) # Sort by start time sort_indices = sortperm(combined_gti[:, 1]) combined_gti = combined_gti[sort_indices, :] - + # Merge overlapping intervals merged_gti = merge_overlapping_gtis(combined_gti) - + return merged_gti, "combined_" * join(gti_sources, "_") else # Return the first GTI found @@ -779,7 +790,7 @@ function extract_timing_keywords(header) return Float64(val) end end - + # Extract MJD reference mjd_ref = nothing if haskey(header, "MJDREF") @@ -1031,6 +1042,9 @@ Times are ready for direct use in timing analysis without additional corrections """ function readevents( path::AbstractString; + mission::Union{String,Nothing} = nothing, + instrument::Union{String,Nothing} = nothing, + epoch::Union{Float64,Nothing} = nothing, hdu::Int = 2, T::Type = Float64, sort::Bool = false, @@ -1038,21 +1052,36 @@ function readevents( energy_alternatives::Vector{String} = ["ENERGY", "PI", "PHA"], load_gti::Bool = true, gti_hdu_candidates::Vector{String} = ["GTI", "STDGTI"], - gti_hdu_indices::Union{Vector{Int}, Nothing} = nothing, + gti_hdu_indices::Union{Vector{Int},Nothing} = nothing, combine_gtis::Bool = true, apply_gti_filter::Bool = false, - convert_to_mjd::Bool = false, # New parameter to control MJD conversion + convert_to_mjd::Bool = false, kwargs..., )::EventList{Vector{T},FITSMetadata{FITSIO.FITSHeader}} + # Get mission support if specified + mission_support = if !isnothing(mission) + ms = get_mission_support(mission, instrument, epoch) + # Use mission-specific energy alternatives if available + energy_alternatives = ms.energy_alternatives + ms + else + nothing + end + # Read GTI first if requested gti_data, gti_source = nothing, nothing if load_gti - gti_data, gti_source = read_gti_from_fits(path; - gti_hdu_candidates, gti_hdu_indices, combine_gtis) - + gti_data, gti_source = read_gti_from_fits( + path; + gti_hdu_candidates, + gti_hdu_indices, + combine_gtis + ) + if !isnothing(gti_data) - @debug "Found GTI data" n_intervals=size(gti_data, 1) time_range=(minimum(gti_data), maximum(gti_data)) + @debug "Found GTI data" n_intervals = size(gti_data, 1) time_range = + (minimum(gti_data), maximum(gti_data)) else @debug "No GTI data found" end @@ -1072,17 +1101,33 @@ function readevents( time_del::Union{Nothing,Float64} = FITS(path, "r") do f selected_hdu = f[hdu] + + # Apply mission-specific FITS interpretation if available + if !isnothing(mission_support) && !isnothing(mission_support.interpretation_func) + interpret_fits_data!(f, mission_support) + end + header = read_header(selected_hdu) all_cols = FITSIO.colnames(selected_hdu) time = convert(Vector{T}, read(selected_hdu, "TIME", case_sensitive = false)) - # Read energy column + # Read energy column using mission-specific alternatives energy_column, energy = read_energy_column( selected_hdu; T = T, energy_alternatives = energy_alternatives, ) + # Apply mission-specific calibration if we have PI data and mission support + if !isnothing(energy) && + !isnothing(mission_support) && + !isnothing(energy_column) && + uppercase(energy_column) == "PI" + energy = convert(Vector{T}, apply_calibration(mission_support, energy)) + # Update the energy column name to reflect that it's now calibrated + energy_column = "ENERGY" + end + # Read extra columns extra_data = Dict{String,Vector}() for col_name in extra_columns @@ -1096,9 +1141,39 @@ function readevents( end # Extract timing keywords - mjd_ref, time_zero, time_unit, time_sys, time_pixr, time_del = extract_timing_keywords(header) + mjd_ref, time_zero, time_unit, time_sys, time_pixr, time_del = + extract_timing_keywords(header) + + ( + time, + energy, + energy_column, + header, + extra_data, + mjd_ref, + time_zero, + time_unit, + time_sys, + time_pixr, + time_del, + ) + end + + # Apply mission-specific header patches if available + if !isnothing(mission_support) + # Convert header to dictionary for patching + header_dict = Dict{String,Any}() + + # Use the proper way to access FITSHeader keys and values + for key in keys(header) + header_dict[key] = header[key] + end - (time, energy, energy_column, header, extra_data, mjd_ref, time_zero, time_unit, time_sys, time_pixr, time_del) + # Apply mission patches + patched_header_dict = patch_mission_info(header_dict, mission) + + # Note: We keep the original header structure but could extend this + # to update the header with patched information if needed end # Apply time corrections and optionally convert to MJD @@ -1110,24 +1185,63 @@ function readevents( if !isnothing(time_del) && time_del > 0.0 effective_timepixr = isnothing(time_pixr) ? 0.0 : time_pixr bin_center_correction = (0.5 - effective_timepixr) * time_del - + corrected_time = time .+ effective_timezero .+ bin_center_correction - @debug "Applied bin-centering correction" time_zero=effective_timezero timedel=time_del timepixr=effective_timepixr bin_correction=bin_center_correction + @debug "Applied bin-centering correction" time_zero = effective_timezero timedel = + time_del timepixr = effective_timepixr bin_correction = bin_center_correction else corrected_time = time .+ effective_timezero - @debug "Applied time zero correction" time_zero=effective_timezero + @debug "Applied time zero correction" time_zero = effective_timezero end - + # Convert to MJD only if requested if convert_to_mjd time = convert(Vector{T}, sec_to_mjd(corrected_time, mjd_ref)) - @debug "Converted to MJD(TT)" mjd_ref=mjd_ref + + # Apply the same time corrections to GTI data + if !isnothing(gti_data) + # Apply time zero correction to GTI + gti_corrected = gti_data .+ effective_timezero + + # Apply bin-centering correction to GTI if applicable + if !isnothing(time_del) && time_del > 0.0 + effective_timepixr = isnothing(time_pixr) ? 0.0 : time_pixr + bin_center_correction = (0.5 - effective_timepixr) * time_del + gti_corrected = gti_corrected .+ bin_center_correction + end + + # Convert GTI to MJD using the same reference + gti_data = sec_to_mjd.(gti_corrected, mjd_ref) + @debug "Applied same time corrections to GTI data" + end + + @debug "Converted to MJD(TT)" mjd_ref = mjd_ref else time = convert(Vector{T}, corrected_time) - @debug "Kept times in seconds (corrected)" + + # Apply time corrections to GTI even when not converting to MJD + if !isnothing(gti_data) + # Apply time zero correction to GTI + gti_corrected = gti_data .+ effective_timezero + + # Apply bin-centering correction to GTI if applicable + if !isnothing(time_del) && time_del > 0.0 + effective_timepixr = isnothing(time_pixr) ? 0.0 : time_pixr + bin_center_correction = (0.5 - effective_timepixr) * time_del + gti_corrected = gti_corrected .+ bin_center_correction + end + + gti_data = gti_corrected + @debug "Applied time corrections to GTI data" + end + + @debug "Kept times in seconds (corrected)" + end + + @debug "Final time range" time_range = (minimum(time), maximum(time)) + if !isnothing(gti_data) + @debug "Final GTI range" gti_range = (minimum(gti_data), maximum(gti_data)) end - - @debug "Final time range" time_range=(minimum(time), maximum(time)) end # Validate data consistency @@ -1137,14 +1251,28 @@ function readevents( # Apply GTI filtering if requested if apply_gti_filter && !isnothing(gti_data) - gti_mask, _ = create_gti_mask(time, gti_data) + @debug "Applying GTI filter" n_events_before = length(time) - time = time[gti_mask] - if !isnothing(energy) - energy = energy[gti_mask] - end - for (col_name, col_data) in extra_data - extra_data[col_name] = col_data[gti_mask] + try + gti_mask, _ = create_gti_mask(time, gti_data) + n_kept = sum(gti_mask) + + @debug "GTI filter results" events_kept = n_kept events_removed = (length(time) - n_kept) + + if n_kept > 0 + time = time[gti_mask] + if !isnothing(energy) + energy = energy[gti_mask] + end + for (col_name, col_data) in extra_data + extra_data[col_name] = col_data[gti_mask] + end + @debug "Successfully applied GTI filter" + else + @warn "GTI filtering removed all events - check time system consistency" event_range = extrema(time) gti_range = extrema(gti_data) + end + catch e + @warn "GTI filtering failed: $e - proceeding without GTI filter" end end @@ -1161,13 +1289,28 @@ function readevents( for (col_name, col_data) in extra_data extra_data[col_name] = col_data[sort_indices] end + @debug "Events sorted by time" else @assert false "Times are not sorted (pass `sort = true` to force sorting)" end end - - meta = FITSMetadata(path, hdu, energy_col, extra_data, header, gti_data, gti_source, - mjd_ref, time_zero, time_unit, time_sys, time_pixr, time_del) + + # Create metadata with all timing information + meta = FITSMetadata( + path, + hdu, + energy_col, + extra_data, + header, + gti_data, + gti_source, + mjd_ref, + time_zero, + time_unit, + time_sys, + time_pixr, + time_del, + ) return EventList(time, energy, meta) end # ============================================================================ @@ -1200,9 +1343,9 @@ println(summary(ev)) function Base.summary(ev::EventList) n_events = length(ev) time_span = isempty(ev.times) ? 0.0 : maximum(ev.times) - minimum(ev.times) - + summary_str = "EventList: $n_events events over $(time_span) time units" - + if has_energies(ev) energy_range = extrema(ev.energies) summary_str *= ", energies: $(energy_range[1]) - $(energy_range[2])" @@ -1227,6 +1370,6 @@ function Base.summary(ev::EventList) effective_extra = max(1, n_extra) # At least 1 when GTI is present summary_str *= ", $effective_extra extra columns" end - + return summary_str end diff --git a/src/gti.jl b/src/gti.jl index 5d74538..99c20aa 100644 --- a/src/gti.jl +++ b/src/gti.jl @@ -1,9 +1,9 @@ -function get_total_gti_length(gti::AbstractMatrix{<:Real}; minlen::Real=0.0) - lengths = diff(gti; dims =2) - return sum(x->x > minlen ? x : zero(x), lengths) +function get_total_gti_length(gti::AbstractMatrix{<:Real}; minlen::Real = 0.0) + lengths = diff(gti; dims = 2) + return sum(x -> x > minlen ? x : zero(x), lengths) end -function load_gtis(fits_file::String, gtistring::String="GTI") +function load_gtis(fits_file::String, gtistring::String = "GTI") gti = FITS(fits_file) do lchdulist gtihdu = lchdulist[gtistring] get_gti_from_hdu(gtihdu) @@ -21,16 +21,19 @@ function get_gti_from_hdu(gtihdu::TableHDU) stopstr = "Stop" end - gtistart = read(gtihdu,startstr) - gtistop = read(gtihdu,stopstr) + gtistart = read(gtihdu, startstr) + gtistop = read(gtihdu, stopstr) if isempty(gtistart) || isempty(gtistop) return reshape(Float64[], 0, 2) end - return mapreduce(permutedims, vcat, - [[a, b] for (a,b) in zip(gtistart, gtistop)], - init=reshape(Float64[], 0, 2)) + return mapreduce( + permutedims, + vcat, + [[a, b] for (a, b) in zip(gtistart, gtistop)], + init = reshape(Float64[], 0, 2), + ) end """ check_gtis(gti::AbstractMatrix) @@ -84,7 +87,7 @@ check_gtis(bad_gtis) # Throws ArgumentError """ function check_gtis(gti::AbstractMatrix) - if ndims(gti) != 2 || size(gti,2) != 2 + if ndims(gti) != 2 || size(gti, 2) != 2 throw(ArgumentError("Please check the formatting of the GTIs. They need to be provided as [[gti00 gti01]; [gti10 gti11]; ...].")) end @@ -98,56 +101,59 @@ function check_gtis(gti::AbstractMatrix) gti_end = @view gti[:, 2] if any(gti_end < gti_start) - throw(ArgumentError( - "The GTI end times must be larger than the GTI start times." - )) + throw(ArgumentError("The GTI end times must be larger than the GTI start times.")) end if any(@view(gti_start[begin+1:end]) < @view(gti_end[begin:end-1])) - throw(ArgumentError( - "This GTI has overlaps" - )) + throw(ArgumentError("This GTI has overlaps")) end end -function create_gti_mask(times::AbstractVector{<:Real},gtis::AbstractMatrix{<:Real}; - safe_interval::AbstractVector{<:Real}=[0,0], min_length::Real=0, - dt::Real = -1, epsilon::Real = 0.001) +function create_gti_mask( + times::AbstractVector{<:Real}, + gtis::AbstractMatrix{<:Real}; + safe_interval::AbstractVector{<:Real} = [0, 0], + min_length::Real = 0, + dt::Real = -1, + epsilon::Real = 0.001, +) if isempty(times) throw(ArgumentError("Passing an empty time array to create_gti_mask")) end check_gtis(gtis) - mask = zeros(Bool,length(times)) + mask = zeros(Bool, length(times)) + + if min_length > 0 + gtis = gtis[min_length.<@view(gtis[:, 2])-@view(gtis[:, 1]), :] - if min_length>0 - gtis = gtis[min_length .< @view(gtis[:,2]) - @view(gtis[:,1]),:] - - if size(gtis,1) < 1 + if size(gtis, 1) < 1 @warn "No GTIs longer than min_length $(min_length)" return mask, reshape(eltype(gtis)[], 0, 2) end - end + end if dt < 0 dt = Statistics.median(diff(times)) end epsilon_times_dt = epsilon * dt - new_gtis = [[0.0, 0.0] for _ in range(1,size(gtis,1))] - new_gti_mask = zeros(Bool, size(gtis,1)) + new_gtis = [[0.0, 0.0] for _ in range(1, size(gtis, 1))] + new_gti_mask = zeros(Bool, size(gtis, 1)) gti_start = @view gtis[:, 1] gti_end = @view gtis[:, 2] - for (ig,(limmin,limmax)) in enumerate(zip(gti_start,gti_end)) + for (ig, (limmin, limmax)) in enumerate(zip(gti_start, gti_end)) limmin += safe_interval[1] limmax -= safe_interval[2] if limmax - limmin >= min_length new_gtis[ig][:] .= limmin, limmax - for (i,t) in enumerate(times) - if (limmin + dt / 2 - epsilon_times_dt) <= t <= (limmax - dt / 2 + epsilon_times_dt) + for (i, t) in enumerate(times) + if (limmin + dt / 2 - epsilon_times_dt) <= + t <= + (limmax - dt / 2 + epsilon_times_dt) mask[i] = true end end @@ -160,16 +166,22 @@ function create_gti_mask(times::AbstractVector{<:Real},gtis::AbstractMatrix{<:Re return mask, reshape(eltype(gtis)[], 0, 2) end - return mask, mapreduce(permutedims, vcat, filtered_gtis, - init=reshape(eltype(gtis)[], 0, 2)) + return mask, + mapreduce(permutedims, vcat, filtered_gtis, init = reshape(eltype(gtis)[], 0, 2)) end -function create_gti_from_condition(time::AbstractVector{<:Real}, condition::AbstractVector{Bool}; - safe_interval::AbstractVector{<:Real}=[0,0], dt::AbstractVector{<:Real}=Float64[]) - +function create_gti_from_condition( + time::AbstractVector{<:Real}, + condition::AbstractVector{Bool}; + safe_interval::AbstractVector{<:Real} = [0, 0], + dt::AbstractVector{<:Real} = Float64[], +) + if length(time) != length(condition) - throw(ArgumentError("The length of the condition and time arrays must be the same.")) + throw( + ArgumentError("The length of the condition and time arrays must be the same."), + ) end idxs = contiguous_regions(condition) @@ -188,22 +200,24 @@ function create_gti_from_condition(time::AbstractVector{<:Real}, condition::Abst if t1 - t0 < 0 continue end - push!(gtis,[t0, t1]) + push!(gtis, [t0, t1]) end - + if isempty(gtis) return reshape(Float64[], 0, 2) end - - return mapreduce(permutedims, vcat, gtis, init=reshape(Float64[], 0, 2)) + + return mapreduce(permutedims, vcat, gtis, init = reshape(Float64[], 0, 2)) end -function operations_on_gtis(gti_list::AbstractVector{<:AbstractMatrix{T}}, - operation::Function) where {T<:Real} - +function operations_on_gtis( + gti_list::AbstractVector{<:AbstractMatrix{T}}, + operation::Function, +) where {T<:Real} + # Convert all GTIs to IntervalSets, handling empty ones interval_sets = IntervalSet[] - + for gti in gti_list if size(gti, 1) == 0 # Empty GTI becomes empty IntervalSet @@ -217,29 +231,34 @@ function operations_on_gtis(gti_list::AbstractVector{<:AbstractMatrix{T}}, push!(interval_sets, IntervalSet(intervals)) end end - + # If no interval sets, return empty if isempty(interval_sets) return reshape(T[], 0, 2) end - + # Apply the operation across all interval sets result_interval = interval_sets[1] - for i in 2:length(interval_sets) + for i = 2:length(interval_sets) result_interval = operation(result_interval, interval_sets[i]) end - + # Convert back to matrix format if isempty(result_interval.items) return reshape(T[], 0, 2) end - + final_gti = Vector{T}[] for interval in result_interval.items push!(final_gti, [first(interval), last(interval)]) end - - return mapreduce(permutedims, vcat, final_gti, init=reshape(T[], 0, 2)) + + return mapreduce(permutedims, vcat, final_gti, init = reshape(T[], 0, 2)) +end +function intersect_gtis(gtis1::AbstractMatrix{<:Real}, gtis2::AbstractMatrix{<:Real}) + check_gtis(gtis1) + check_gtis(gtis2) + return operations_on_gtis([gtis1, gtis2], intersect) end """ @@ -280,7 +299,7 @@ function get_btis(gtis::AbstractMatrix{<:Real}) if isempty(gtis) throw(ArgumentError("Empty GTI and no valid start_time and stop_time")) end - return get_btis(gtis, gtis[1,1], gtis[end,2]) + return get_btis(gtis, gtis[1, 1], gtis[end, 2]) end """ get_btis(gtis::AbstractMatrix{T}, start_time, stop_time) -> Matrix{T} where {T<:Real} @@ -355,18 +374,18 @@ function get_btis(gtis::AbstractMatrix{T}, start_time, stop_time) where {T<:Real if isempty(gtis) return T[start_time stop_time] end - + # Only check GTIs if they're not empty if size(gtis, 1) > 0 check_gtis(gtis) end - total_interval = Interval{T, Closed, Open}[Interval{T, Closed, Open}(start_time, stop_time)] + total_interval = Interval{T,Closed,Open}[Interval{T,Closed,Open}(start_time, stop_time)] total_interval_set = IntervalSet(total_interval) - gti_interval = Interval{T, Closed, Open}[] + gti_interval = Interval{T,Closed,Open}[] for ig in eachrow(gtis) - push!(gti_interval,Interval{T, Closed, Open}(ig[1],ig[2])) + push!(gti_interval, Interval{T,Closed,Open}(ig[1], ig[2])) end gti_interval_set = IntervalSet(gti_interval) @@ -382,29 +401,38 @@ function get_btis(gtis::AbstractMatrix{T}, start_time, stop_time) where {T<:Real return reshape(T[], 0, 2) end - return mapreduce(permutedims, vcat, btis, init=reshape(T[], 0, 2)) + return mapreduce(permutedims, vcat, btis, init = reshape(T[], 0, 2)) end -function time_intervals_from_gtis(gtis::AbstractMatrix{<:Real}, segment_size::Real; - fraction_step::Real=1, epsilon::Real=1e-5) +function time_intervals_from_gtis( + gtis::AbstractMatrix{<:Real}, + segment_size::Real; + fraction_step::Real = 1, + epsilon::Real = 1e-5, +) spectrum_start_times = Float64[] - gti_low = @view gtis[:,1] - gti_up = @view gtis[:,2] + gti_low = @view gtis[:, 1] + gti_up = @view gtis[:, 2] - for (g1,g2) in zip(gti_low,gti_up) + for (g1, g2) in zip(gti_low, gti_up) if g2 - g1 + epsilon < segment_size continue end - newtimes = range(g1, g2 - segment_size + epsilon, step = segment_size* fraction_step) - append!(spectrum_start_times,newtimes) + newtimes = + range(g1, g2 - segment_size + epsilon, step = segment_size * fraction_step) + append!(spectrum_start_times, newtimes) end return spectrum_start_times, spectrum_start_times .+ segment_size end -function calculate_segment_bin_start(startbin::Integer, stopbin::Integer, - nbin::Integer; fraction_step::Real=1) - st = floor.(range(startbin, stopbin, step=Int(nbin * fraction_step))) +function calculate_segment_bin_start( + startbin::Integer, + stopbin::Integer, + nbin::Integer; + fraction_step::Real = 1, +) + st = floor.(range(startbin, stopbin, step = Int(nbin * fraction_step))) if st[end] == stopbin pop!(st) end @@ -414,73 +442,84 @@ function calculate_segment_bin_start(startbin::Integer, stopbin::Integer, return st end -function bin_intervals_from_gtis(gtis::AbstractMatrix{<:Real}, segment_size::Real, - time::AbstractVector{<:Real}; dt=nothing, - fraction_step::Real=1, epsilon::Real=0.001) +function bin_intervals_from_gtis( + gtis::AbstractMatrix{<:Real}, + segment_size::Real, + time::AbstractVector{<:Real}; + dt = nothing, + fraction_step::Real = 1, + epsilon::Real = 0.001, +) if isnothing(dt) dt = Statistics.median(diff(time)) end time_min, time_max = extrema(time) gti_min = minimum(gtis[:, 1]) gti_max = maximum(gtis[:, 2]) - + if gti_max < time_min || gti_min > time_max throw(ArgumentError("GTIs and time array do not overlap")) end - + epsilon_times_dt = epsilon * dt nbin = round(Int, segment_size / dt) spectrum_start_bins = Int[] gti_low = @view(gtis[:, 1]) .+ (dt ./ 2 .- epsilon_times_dt) gti_up = @view(gtis[:, 2]) .- (dt ./ 2 .- epsilon_times_dt) - + for (g0, g1) in zip(gti_low, gti_up) if (g1 - g0 .+ (dt + epsilon_times_dt)) < segment_size continue end startbin, stopbin = searchsortedfirst.(Ref(time), [g0, g1]) startbin -= 1 - + # The issue is here - we need to be more careful with bounds checking if startbin < 0 startbin = 0 end - + if stopbin > length(time) stopbin = length(time) end - + # Only proceed if we have valid indices if startbin >= length(time) continue end - + if startbin + 1 <= length(time) && time[startbin+1] < g0 startbin += 1 end - + # Would be g[1] - dt/2, but stopbin is the end of an interval # so one has to add one bin if stopbin <= length(time) && time[stopbin] > g1 stopbin -= 1 end - + # Make sure we still have valid range after adjustments if startbin >= stopbin || startbin < 0 || stopbin > length(time) continue end - + newbins = calculate_segment_bin_start( - startbin, stopbin, nbin, fraction_step=fraction_step) - + startbin, + stopbin, + nbin, + fraction_step = fraction_step, + ) + append!(spectrum_start_bins, newbins) end return spectrum_start_bins, spectrum_start_bins .+ nbin end -@resumable function generate_indices_of_segment_boundaries_unbinned(times::AbstractVector{<:Real}, - gti::AbstractMatrix{<:Real}, - segment_size::Real) +@resumable function generate_indices_of_segment_boundaries_unbinned( + times::AbstractVector{<:Real}, + gti::AbstractMatrix{<:Real}, + segment_size::Real, +) start, stop = time_intervals_from_gtis(gti, segment_size) startidx = searchsortedfirst.(Ref(times), start) @@ -491,17 +530,22 @@ end end end -@resumable function generate_indices_of_segment_boundaries_binned(times::AbstractVector{<:Real}, - gti::AbstractMatrix{<:Real}, - segment_size::Real; dt=nothing) - startidx, stopidx = bin_intervals_from_gtis(gti, segment_size, times; - dt=dt) +@resumable function generate_indices_of_segment_boundaries_binned( + times::AbstractVector{<:Real}, + gti::AbstractMatrix{<:Real}, + segment_size::Real; + dt = nothing, +) + startidx, stopidx = bin_intervals_from_gtis(gti, segment_size, times; dt = dt) if isnothing(dt) dt = 0 end for (idx0, idx1) in zip(startidx, stopidx) - @yield times[idx0+1] - dt / 2, times[min(idx1, length(times) - 1)] - dt / 2,idx0, idx1 + @yield times[idx0+1] - dt / 2, + times[min(idx1, length(times) - 1)] - dt / 2, + idx0, + idx1 end end """ @@ -548,24 +592,24 @@ end """ function split_by_gtis(el::EventList, gtis::AbstractMatrix{<:Real}) check_gtis(gtis) - + result = EventList[] - - for i in 1:size(gtis, 1) + + for i = 1:size(gtis, 1) gti_start, gti_stop = gtis[i, 1], gtis[i, 2] - + # Create filter function for this specific GTI gti_filter(t) = gti_start ≤ t ≤ gti_stop - + # Apply temporal filtering using existing infrastructure filtered_el = filter_time(gti_filter, el) - + # Only include GTIs that contain events if length(filtered_el.times) > 0 push!(result, filtered_el) end end - + return result end """ @@ -662,26 +706,26 @@ end - [`check_gtis`](@ref): Validate GTI format - [`create_filtered_lightcurve`](@ref): Create filtered light curve segments """ -function apply_gtis(lc::LightCurve{T}, gtis::AbstractMatrix{<:Real}) where T +function apply_gtis(lc::LightCurve{T}, gtis::AbstractMatrix{<:Real}) where {T} check_gtis(gtis) - + result = LightCurve{T}[] - + # Pre-allocate the mask buffer once bin_mask = similar(lc.time, Bool) - - for i in 1:size(gtis, 1) + + for i = 1:size(gtis, 1) gti_start, gti_stop = T(gtis[i, 1]), T(gtis[i, 2]) - + # Use @. to vectorize the mask creation @. bin_mask = (lc.time ≥ gti_start) & (lc.time ≤ gti_stop) - + if any(bin_mask) filtered_lc = create_filtered_lightcurve(lc, bin_mask, gti_start, gti_stop, i) push!(result, filtered_lc) end end - + return result end """ @@ -709,8 +753,13 @@ while preserving all metadata, properties, and statistical characteristics. - Updates metadata to reflect GTI filtering - Recalculates statistical errors for the filtered dataset """ -function create_filtered_lightcurve(lc::LightCurve{T}, mask::AbstractVector{Bool}, - gti_start::T, gti_stop::T, gti_index::Int) where T +function create_filtered_lightcurve( + lc::LightCurve{T}, + mask::AbstractVector{Bool}, + gti_start::T, + gti_stop::T, + gti_index::Int, +) where {T} # Ensure mask is proper boolean vector bool_mask = mask isa BitVector ? mask : BitVector(mask) @@ -718,33 +767,36 @@ function create_filtered_lightcurve(lc::LightCurve{T}, mask::AbstractVector{Bool filtered_time = lc.time[mask] filtered_counts = lc.counts[mask] filtered_exposure = isnothing(lc.exposure) ? nothing : lc.exposure[mask] - + # Filter all computed properties filtered_properties = EventProperty{T}[] for prop in lc.properties filtered_values = prop.values[mask] push!(filtered_properties, EventProperty(prop.name, filtered_values, prop.unit)) end - + # Update metadata with GTI information updated_metadata = LightCurveMetadata( lc.metadata.telescope, - lc.metadata.instrument, + lc.metadata.instrument, lc.metadata.object, lc.metadata.mjdref, (Float64(gti_start), Float64(gti_stop)), # Update time range to GTI bounds lc.metadata.bin_size, lc.metadata.headers, - merge(lc.metadata.extra, Dict{String,Any}( - "gti_applied" => true, - "gti_index" => gti_index, - "gti_bounds" => [Float64(gti_start), Float64(gti_stop)], - "original_time_range" => lc.metadata.time_range, - "filtered_nbins" => length(filtered_time), - "original_nbins" => length(lc.time) - )) + merge( + lc.metadata.extra, + Dict{String,Any}( + "gti_applied" => true, + "gti_index" => gti_index, + "gti_bounds" => [Float64(gti_start), Float64(gti_stop)], + "original_time_range" => lc.metadata.time_range, + "filtered_nbins" => length(filtered_time), + "original_nbins" => length(lc.time), + ), + ), ) - + # Create new LightCurve with filtered data filtered_lc = LightCurve{T}( filtered_time, @@ -754,12 +806,12 @@ function create_filtered_lightcurve(lc::LightCurve{T}, mask::AbstractVector{Bool filtered_exposure, filtered_properties, updated_metadata, - lc.err_method + lc.err_method, ) - + # Recalculate errors for the filtered dataset calculate_errors!(filtered_lc) - + return filtered_lc end """ @@ -852,53 +904,57 @@ end - [`get_btis`](@ref): Function to compute Bad Time Intervals - [`apply_gtis`](@ref): Apply GTIs to segment data """ -function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; - dt::Real=1.0, random_fill_threshold::Real=10.0, - rng::AbstractRNG=Random.GLOBAL_RNG) +function fill_bad_time_intervals!( + el::EventList, + gtis::AbstractMatrix{<:Real}; + dt::Real = 1.0, + random_fill_threshold::Real = 10.0, + rng::AbstractRNG = Random.GLOBAL_RNG, +) check_gtis(gtis) - + # Determine the time range for BTI calculation time_start, time_stop = extrema(el.times) - + # Calculate Bad Time Intervals btis = get_btis(gtis, time_start, time_stop) - + # Track synthetic events for metadata n_synthetic_events = 0 filled_intervals = Float64[] - + # Store all synthetic data before appending to avoid index issues all_synthetic_times = Float64[] all_synthetic_energies = Float64[] - synthetic_extra_columns = Dict{String, Vector}() - + synthetic_extra_columns = Dict{String,Vector}() + # Initialize synthetic extra columns for (col_name, col_data) in el.meta.extra_columns synthetic_extra_columns[col_name] = similar(col_data, 0) end - + # Process each BTI - for i in 1:size(btis, 1) + for i = 1:size(btis, 1) bti_start, bti_stop = btis[i, 1], btis[i, 2] bti_duration = bti_stop - bti_start - + # Skip BTIs that are too long or have zero/negative duration if bti_duration >= random_fill_threshold || bti_duration <= 0 continue end - + # Calculate number of synthetic events based on event rate in nearby GTIs gti_rates = Float64[] - for j in 1:size(gtis, 1) + for j = 1:size(gtis, 1) gti_start, gti_stop = gtis[j, 1], gtis[j, 2] events_in_gti = count(t -> gti_start <= t <= gti_stop, el.times) gti_duration_j = gti_stop - gti_start - + if gti_duration_j > 0 && events_in_gti > 0 push!(gti_rates, events_in_gti / gti_duration_j) end end - + # Determine synthetic event count n_synthetic = 0 if isempty(gti_rates) @@ -915,7 +971,7 @@ function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; median_rate = Statistics.median(gti_rates) n_synthetic = max(1, round(Int, median_rate * bti_duration)) end - + if n_synthetic > 0 # Generate uniformly distributed synthetic times within the BTI # Generate times strictly within (bti_start, bti_stop) exclusive bounds @@ -924,23 +980,25 @@ function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; effective_start = bti_start + epsilon effective_stop = bti_stop - epsilon effective_duration = effective_stop - effective_start - - synthetic_times = sort!(rand(rng, n_synthetic) .* effective_duration .+ effective_start) + + synthetic_times = + sort!(rand(rng, n_synthetic) .* effective_duration .+ effective_start) append!(all_synthetic_times, synthetic_times) - + n_synthetic_events += n_synthetic push!(filled_intervals, bti_duration) - + # Handle energy values if present - sample from nearby GTI events if !isnothing(el.energies) nearby_energies = Float64[] - + # Find energies from GTIs that are adjacent to this BTI - for j in 1:size(gtis, 1) + for j = 1:size(gtis, 1) gti_start, gti_stop = gtis[j, 1], gtis[j, 2] - + # Check if this GTI is adjacent to the BTI - if abs(gti_stop - bti_start) < 1e-10 || abs(gti_start - bti_stop) < 1e-10 + if abs(gti_stop - bti_start) < 1e-10 || + abs(gti_start - bti_stop) < 1e-10 # Find events in this GTI for (idx, t) in enumerate(el.times) if gti_start <= t <= gti_stop @@ -949,7 +1007,7 @@ function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; end end end - + # Sample energies if !isempty(nearby_energies) synthetic_energies = rand(rng, nearby_energies, n_synthetic) @@ -958,20 +1016,21 @@ function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; else synthetic_energies = zeros(eltype(el.energies), n_synthetic) end - + append!(all_synthetic_energies, synthetic_energies) end - + # Handle extra columns with values from nearby GTIs for (col_name, col_data) in el.meta.extra_columns nearby_values = similar(col_data, 0) - + # Find values from GTIs that are adjacent to this BTI - for j in 1:size(gtis, 1) + for j = 1:size(gtis, 1) gti_start, gti_stop = gtis[j, 1], gtis[j, 2] - + # Check if this GTI is adjacent to the BTI - if abs(gti_stop - bti_start) < 1e-10 || abs(gti_start - bti_stop) < 1e-10 + if abs(gti_stop - bti_start) < 1e-10 || + abs(gti_start - bti_stop) < 1e-10 # Find events in this GTI for (idx, t) in enumerate(el.times) if gti_start <= t <= gti_stop @@ -980,7 +1039,7 @@ function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; end end end - + # Sample values if !isempty(nearby_values) synthetic_values = rand(rng, nearby_values, n_synthetic) @@ -989,56 +1048,60 @@ function fill_bad_time_intervals!(el::EventList, gtis::AbstractMatrix{<:Real}; else synthetic_values = zeros(eltype(col_data), n_synthetic) end - + append!(synthetic_extra_columns[col_name], synthetic_values) end end end - + # Now append all synthetic data at once if n_synthetic_events > 0 append!(el.times, all_synthetic_times) - + if !isnothing(el.energies) append!(el.energies, all_synthetic_energies) end - + # Append synthetic extra columns for (col_name, col_data) in el.meta.extra_columns append!(col_data, synthetic_extra_columns[col_name]) end - + # Re-sort all arrays to maintain chronological order sort_indices = sortperm(el.times) el.times[:] = el.times[sort_indices] - + if !isnothing(el.energies) el.energies[:] = el.energies[sort_indices] end - + # Sort extra columns for (col_name, col_data) in el.meta.extra_columns col_data[:] = col_data[sort_indices] end - + # metadata with BTI filling information - merge!(el.meta.headers, Dict{String,Any}( - "BTI_FILLED" => true, - "N_SYNTH_EVENTS" => n_synthetic_events, - "RAND_FILL_THRESH" => random_fill_threshold, - "BTI_FILL_DT" => dt - )) - + merge!( + el.meta.headers, + Dict{String,Any}( + "BTI_FILLED" => true, + "N_SYNTH_EVENTS" => n_synthetic_events, + "RAND_FILL_THRESH" => random_fill_threshold, + "BTI_FILL_DT" => dt, + ), + ) + # Store vector metadata in extra_columns el.meta.extra_columns["filled_bti_durations"] = filled_intervals end - + return el end +get_gti_lengths(gti::AbstractMatrix{<:Real}) = vec(diff(gti; dims=2)) #todo # create a function that fills the bad time intervals in a light curve # in order to maintain optiizational sampling for periodograms # can be intially start like : # fill_bad_time_intervals!(lc::LightCurve{T}, gtis::AbstractMatrix{<:Real}; # dt::Real=1.0, random_fill_threshold::Real=10.0, -# rng::AbstractRNG=Random.GLOBAL_RNG) where T<:Real \ No newline at end of file +# rng::AbstractRNG=Random.GLOBAL_RNG) where T<:Real diff --git a/src/lightcurve.jl b/src/lightcurve.jl index ab76161..0b5c6d7 100644 --- a/src/lightcurve.jl +++ b/src/lightcurve.jl @@ -495,6 +495,77 @@ function apply_filters( return filtered_times, filtered_energies, start_t, stop_t end +""" + apply_filters(times, energies, tstart, tstop, energy_filter) + +Basic event filtering without GTI consideration. + +# Arguments +- `times::AbstractVector{T}`: Event arrival times +- `energies::Union{Nothing,AbstractVector{T}}`: Event energies (optional) +- `tstart::Union{Nothing,Real}`: Start time filter (inclusive) +- `tstop::Union{Nothing,Real}`: Stop time filter (inclusive) +- `energy_filter::Union{Nothing,Tuple{Real,Real}}`: Energy range (emin, emax) + +# Returns +`Tuple{Vector, Union{Nothing,Vector}, Real, Real}`: +- Filtered times +- Filtered energies (if provided) +- Actual start time +- Actual stop time + +# Examples +```julia +# Time filtering only +filtered_times, _, start_t, stop_t = apply_filters(times, nothing, 1000.0, 2000.0, nothing) + +# Energy and time filtering +filtered_times, filtered_energies, start_t, stop_t = apply_filters( + times, energies, 1000.0, 2000.0, (0.5, 10.0) +) +Notes +- Energy filter is applied as: emin ≤ energy < emax +- Time filter is applied as: tstart ≤ time ≤ tstop +""" +function apply_filters( + times::AbstractVector{T}, + energies::Union{Nothing,AbstractVector{T}}, + eventlist::EventList, + tstart::Union{Nothing,Real}, + tstop::Union{Nothing,Real}, + energy_filter::Union{Nothing,Tuple{Real,Real}}, + binsize::Real +) where T + mask = trues(length(times)) + + # Apply energy filter + if !isnothing(energy_filter) && !isnothing(energies) + emin, emax = energy_filter + mask = mask .& (energies .>= emin) .& (energies .< emax) + end + + # Apply time filters + if !isnothing(tstart) + mask = mask .& (times .>= tstart) + end + if !isnothing(tstop) + mask = mask .& (times .<= tstop) + end + # If GTI is present, apply GTI mask + if has_gti(eventlist) + gti_mask, _ = create_gti_mask(times, eventlist.meta.gti, dt=binsize) + mask = mask .& gti_mask + end + + !any(mask) && throw(ArgumentError("No events remain after applying filters")) + + filtered_times = times[mask] + filtered_energies = isnothing(energies) ? nothing : energies[mask] + start_t = isnothing(tstart) ? minimum(filtered_times) : tstart + stop_t = isnothing(tstop) ? maximum(filtered_times) : tstop + + return filtered_times, filtered_energies, start_t, stop_t +end """ calculate_event_properties(times, energies, dt, bin_centers) -> Vector{EventProperty} @@ -605,17 +676,24 @@ function extract_metadata(eventlist::EventList, start_time, stop_time, binsize, n_filtered_events # Assume it's already a count end + # Create extra metadata with GTI information[storing purpose] extra_metadata = Dict{String,Any}( "filtered_nevents" => n_events, "total_nevents" => length(eventlist.times), "energy_filter" => energy_filter, - "binning_method" => "histogram" + "binning_method" => "histogram", + "gti" => eventlist.meta.gti #GTI information[this 'eventlist.meta.gti' need to bw rembered since it is the one where u will all gti information:)] ) if hasfield(typeof(eventlist.meta), :extra) merge!(extra_metadata, eventlist.meta.extra) end + # Add GTI source information if available + if !isnothing(eventlist.meta.gti_source) + extra_metadata["gti_source"] = eventlist.meta.gti_source + end + return LightCurveMetadata( telescope, instrument, object, Float64(mjdref), (Float64(start_time), Float64(stop_time)), Float64(binsize), @@ -691,22 +769,20 @@ function create_lightcurve( !(err_method in [:poisson, :gaussian]) && throw(ArgumentError( "Unsupported error method: $err_method. Use :poisson or :gaussian" )) - binsize_t = convert(T, binsize) - # Apply filters to get filtered times and energies + filtered_times, filtered_energies, start_t, stop_t = apply_filters( eventlist.times, eventlist.energies, + eventlist, tstart, tstop, - energy_filter + energy_filter, + binsize_t ) - # Check if we have any events left after filtering - if isempty(filtered_times) - throw(ArgumentError("No events remain after filtering")) - end + isempty(filtered_times) && throw(ArgumentError("No events remain after filtering")) # Determine time range start_time = minimum(filtered_times) @@ -716,17 +792,17 @@ function create_lightcurve( dt, bin_centers = create_time_bins(start_time, stop_time, binsize_t) counts = bin_events(filtered_times, dt) - @info "Created light curve: $(length(bin_centers)) bins, bin size = $(binsize_t) s" + @debug "Created light curve: $(length(bin_centers)) bins, bin size = $(binsize_t) s" # Calculate exposure and properties exposure = fill(binsize_t, length(bin_centers)) properties = calculate_event_properties(filtered_times, filtered_energies, dt, bin_centers) - # Extract metadata - use filtered time range if time filtering was applied + # Extract metadata with GTI information actual_start = !isnothing(tstart) ? T(tstart) : start_time actual_stop = !isnothing(tstop) ? T(tstop) : stop_time metadata = extract_metadata(eventlist, actual_start, actual_stop, binsize_t, - length(filtered_times), energy_filter) + length(filtered_times), energy_filter) # Create light curve (errors will be calculated when needed) lc = LightCurve{T}( @@ -737,6 +813,11 @@ function create_lightcurve( # Calculate initial errors calculate_errors!(lc) + # Add debug info about GTI + if has_gti(eventlist) + @debug "GTI information preserved" n_intervals=size(eventlist.meta.gti, 1) time_range=extrema(eventlist.meta.gti) + end + return lc end """ diff --git a/src/plotting/plots_recipes_gti.jl b/src/plotting/plots_recipes_gti.jl new file mode 100644 index 0000000..0825d2f --- /dev/null +++ b/src/plotting/plots_recipes_gti.jl @@ -0,0 +1,269 @@ +#wrapper to avoid recpies conflicts +struct BTIAnalysisPlot{T} + eventlist::EventList{T} +end + +""" + btianalysis(events::EventList{T}) -> BTIAnalysisPlot{T} + +Create a BTIAnalysisPlot object that can be plotted to show Bad Time Interval (BTI) length distribution. + +This function creates a plottable object that analyzes the distribution of bad time intervals by extracting +Good Time Intervals (GTIs) from the event metadata and computing the complementary BTIs. When plotted, it +generates a logarithmic histogram showing the frequency distribution of BTI durations. + +# Arguments +- `events::EventList{T}`: The input event list containing timing data and GTI metadata + +# Returns +- `BTIAnalysisPlot{T}`: A plottable object containing the event list data + +# Usage +The returned object can be plotted using the standard `plot()` function with various customization options: + +```julia +# Basic BTI analysis plot +bti_plot = btianalysis(events) +plot(bti_plot, bti_analysis=true) + +# Custom binning and axis limits +plot(bti_plot, bti_analysis=true, nbins=50, xlims_range=(1e-4, 1e4)) + +# Analysis with custom y-axis range +plot(bti_plot, bti_analysis=true, ylims_range=(1, 1000)) +``` + +--- + + plot(bti_plot::BTIAnalysisPlot{T}; bti_analysis=false, nbins=30, min_length=1e-3, max_length=10000.0, xlims_range=nothing, ylims_range=nothing) where T + +Plot a histogram of Bad Time Interval (BTI) lengths from a BTIAnalysisPlot object. + +# Plot Arguments +- `bti_plot::BTIAnalysisPlot{T}`: The BTI analysis object created by `btianalysis()` +- `bti_analysis::Bool=false`: Enable BTI analysis (plot returns `nothing` if `false`) +- `nbins::Int=30`: Number of histogram bins for the distribution plot +- `min_length::Float64=1e-3`: Minimum BTI length threshold (currently unused in filtering) +- `max_length::Float64=10000.0`: Maximum BTI length threshold (currently unused in filtering) +- `xlims_range=nothing`: Custom x-axis limits as a tuple `(min, max)`, or `nothing` for auto-scaling +- `ylims_range=nothing`: Custom y-axis limits as a tuple `(min, max)`, or `nothing` for auto-scaling + +# Returns +- `Vector{Float64}`: Array of BTI lengths in seconds when BTIs are found +- `Vector{Float64}`: Dummy data `[0.5]` when no BTIs exist (creates informational plot) +- `nothing`: When `bti_analysis=false` + +# Behavior +- Returns `nothing` immediately if `bti_analysis=false` +- Creates informational plots with status messages when no GTIs or BTIs are found +- Generates logarithmic histogram of BTI length distribution when BTIs exist +- Prints diagnostic statistics including total exposure and BTI lengths +- Handles various error conditions gracefully with fallback plots + +# Plot Properties +- Uses logarithmic scaling on both axes for better visualization of wide dynamic ranges +- Automatic bin range calculation based on data extent +- Steel blue fill color with transparency +- Grid enabled for easier reading + +# Complete Examples +```julia +using Plots + +# Read event data +events = readevents("your_file.fits") + +# Create BTI analysis object and plot +bti_plot = btianalysis(events) +plot(bti_plot, bti_analysis=true) + +# Customized analysis with specific parameters +plot(bti_plot, + bti_analysis=true, + nbins=50, + xlims_range=(1e-4, 1e4), + ylims_range=(1, 1000)) + +# Save the plot +savefig("bti_analysis.png") +``` + +# Workflow +1. Create EventList with GTI metadata: `events = readevents("file")` +2. Create BTI analysis object: `bti_plot = btianalysis(events)` +3. Generate plot: `plot(bti_plot, bti_analysis=true)` + +# Notes +- Requires GTI information in `events.meta.gti` to function properly +- Short BTIs (< 1.0 second) are tracked separately in diagnostic output +- Function includes extensive error handling for missing or invalid GTI data +- The `bti_analysis=true` parameter must be set to generate the actual analysis plot +""" + +# Constructor function (typically defined elsewhere) +btianalysis(events::EventList{T}) where T = BTIAnalysisPlot(events) + +@recipe function f(events::BTIAnalysisPlot{T}; + bti_analysis=false, + nbins=30, + min_length=1e-3, + max_length=10000.0, + xlims_range=nothing, + ylims_range=nothing) where T + + !bti_analysis && return nothing + + gti_event = events.eventlist + + function create_no_bti_plot() + title --> "Bad Time Interval Analysis: No BTIs Found" + xlabel --> "Observation Quality" + ylabel --> "Coverage Status" + grid --> false + legend --> false + size --> (600, 200) + seriestype --> :scatter + markersize --> 0 + xlims --> (0, 1) + ylims --> (0, 1) + annotations --> [(0.5, 0.5, "No Bad Time Intervals Found\nAll observation time covered by GTIs")] + return [0.5], [0.5] + end + + gtis = gti_event.meta.gti + if isnothing(gtis) || isempty(gtis) + error("No GTI information available in EventList") + end + + total_exposure = gti_exposure(gti_event) + + if isempty(gtis) || size(gtis, 1) == 0 + println("No BTI found - all GTI") + return create_no_bti_plot() + end + + start_time = isempty(gti_event.times) ? 0.0 : minimum(gti_event.times) + stop_time = isempty(gti_event.times) ? 0.0 : maximum(gti_event.times) + + btis = nothing + try + btis = get_btis(gtis, start_time, stop_time) + catch + println("No BTI found - all GTI") + return create_no_bti_plot() + end + + if isnothing(btis) || isempty(btis) || size(btis, 1) == 0 + println("No BTI found - all GTI") + return create_no_bti_plot() + end + + bti_lengths = nothing + try + bti_lengths = get_gti_lengths(btis) + catch + println("No BTI found - all GTI") + return create_no_bti_plot() + end + + if isnothing(bti_lengths) || isempty(bti_lengths) || length(bti_lengths) == 0 + println("No BTI found - all GTI") + return create_no_bti_plot() + end + + total_bti_length = 0.0 + try + total_bti_length = get_total_gti_length(btis) + catch + total_bti_length = 0.0 + end + + # Calculate short BTI length (< 1.0 second) + total_short_bti_length = 0.0 + try + short_bti_mask = bti_lengths .< 1.0 + if any(short_bti_mask) + short_btis = btis[short_bti_mask, :] + total_short_bti_length = get_total_gti_length(short_btis) + end + catch + total_short_bti_length = 0.0 + end + + # Print diagnostic statistics + println("Total exposure: $(total_exposure)") + println("Total BTI length: $(total_bti_length)") + println("Total BTI length (short BTIs): $(total_short_bti_length)") + + data_min, data_max = 0.0, 1.0 + try + data_min = minimum(bti_lengths) + data_max = maximum(bti_lengths) + catch + println("Error calculating data range - no BTI found") + return create_no_bti_plot() + end + + bin_min = max(1e-6, min(1e-3, data_min * 0.9)) # Never go below 1e-6 + bin_max = max(1e-3, max(10000.0, data_max * 1.1)) # Ensure bin_max > bin_min + + # Additional safety check + if bin_min <= 0 || !isfinite(bin_min) + bin_min = 1e-6 + end + if bin_max <= bin_min || !isfinite(bin_max) + bin_max = max(1e-3, bin_min * 1000) + end + + title --> "Distribution of Bad Time Interval Lengths" + xlabel --> "Length of bad time interval (seconds)" + ylabel --> "Number of intervals" + xscale --> :log10 + yscale --> :log10 + + if !isnothing(xlims_range) + xlims --> xlims_range + else + xlims --> (bin_min, bin_max) + end + + if !isnothing(ylims_range) + ylims --> ylims_range + else + ylims --> (0.5, max(100, length(bti_lengths))) + end + + grid --> true + legend --> false + size --> (600, 400) + + # Manual binning for clean bars on log scale + log_min = log10(bin_min) + log_max = log10(bin_max) + bin_edges = 10 .^ range(log_min, log_max, length=nbins+1) + + # Calculate histogram counts + hist_counts = zeros(Int, nbins) + for val in bti_lengths + bin_idx = searchsortedlast(bin_edges, val) + if bin_idx > 0 && bin_idx <= nbins + hist_counts[bin_idx] += 1 + end + end + + # Calculate bin centers (geometric mean for log scale) + bin_centers = sqrt.(bin_edges[1:end-1] .* bin_edges[2:end]) + + # Filter out empty bins + valid_bins = hist_counts .> 0 + plot_x = bin_centers[valid_bins] + plot_y = hist_counts[valid_bins] + + seriestype --> :bar + fillcolor --> :steelblue + fillalpha --> 0.8 + linecolor --> :steelblue + linewidth --> 1 + + return plot_x, plot_y +end \ No newline at end of file diff --git a/src/plotting/plots_recipes_lightcurve.jl b/src/plotting/plots_recipes_lightcurve.jl new file mode 100644 index 0000000..4ea9166 --- /dev/null +++ b/src/plotting/plots_recipes_lightcurve.jl @@ -0,0 +1,1122 @@ +using Plots:text +function create_segments(events::EventList, segment_duration::Real; bin_size::Real = 1.0) + # Get actual time range from the data + t_start = minimum(events.times) + t_stop = maximum(events.times) + total_time = t_stop - t_start + println("Data time range: $(t_start) to $(t_stop) ($(total_time) seconds total)") + + # Calculate number of segments + n_segments = ceil(Int, total_time / segment_duration) + println("Creating $(n_segments) segments of $(segment_duration) seconds each") + + segments = Vector{LightCurve}() + + for i = 1:n_segments + start_time = t_start + (i - 1) * segment_duration + stop_time = min(t_start + i * segment_duration, t_stop) + + # Filter events in this segment using matrix operations + event_times = events.times + mask = (event_times .>= start_time) .& (event_times .<= stop_time) + segment_events = event_times[mask] + + events_in_segment = sum(mask) + println( + "Segment $i: $(events_in_segment) events from $(start_time) to $(stop_time)", + ) + + # Create time bins for this segment + time_edges = collect(start_time:bin_size:stop_time) + if length(time_edges) < 2 + time_edges = [start_time, stop_time] + end + + # Bin centers + time_centers = (time_edges[1:end-1] + time_edges[2:end]) / 2 + n_bins = length(time_centers) + + # Initialize counts + counts = zeros(Int, n_bins) + + # Histogram events into bins using matrix operations + if events_in_segment > 0 + # Use searchsortedlast for efficient binning + for event_time in segment_events + bin_idx = searchsortedlast(time_edges, event_time) + if bin_idx > 0 && bin_idx <= n_bins + counts[bin_idx] += 1 + end + end + end + + # Create light curve data matrix: [time, counts, errors] + # Calculate Poisson errors + count_errors = sqrt.(max.(counts, 1)) # Avoid sqrt(0) + + # Create the data matrix + lc_matrix = hcat(time_centers, counts, count_errors) + + # Create dummy metadata + dummy_metadata = LightCurveMetadata( + "Unknown", + "Unknown", + "Unknown", + 0.0, + (start_time, stop_time), + bin_size, + Vector{Dict{String,Any}}(), + Dict{String,Any}(), + ) + + # Create LightCurve object + lc = LightCurve( + lc_matrix[:, 1], # time + bin_size, # dt + Int.(lc_matrix[:, 2]), # counts + lc_matrix[:, 3], # count_error + nothing, # exposure + Vector{EventProperty{Float64}}(), # properties + dummy_metadata, # metadata + :poisson, # err_method + ) + + push!(segments, lc) + end + + return segments +end +""" + plot(el::EventList{T}, bin_size::Real=1.0; kwargs...) + +Plot a light curve from an EventList with optional Good Time Intervals (GTIs) and Bad Time Intervals (BTIs). + +# Arguments +- `el::EventList{T}`: Event list containing photon arrival times +- `bin_size::Real=1.0`: Time bin size in seconds + +# Keywords +- `tstart=nothing`: Start time for the light curve (defaults to first event) +- `tstop=nothing`: Stop time for the light curve (defaults to last event) +- `energy_filter=nothing`: Energy range filter as (min, max) tuple +- `show_errors=false`: Display error bars using specified error method +- `show_btis=false`: Show Bad Time Intervals as red shaded regions +- `show_bti=false`: Alias for `show_btis` +- `show_gtis=false`: Show Good Time Intervals as green shaded regions +- `show_gti=false`: Alias for `show_gtis` +- `gtis=nothing`: GTI matrix to use (overrides file/metadata GTIs) +- `gti_file=nothing`: FITS file containing GTI extension +- `gti_hdu="GTI"`: HDU name for GTI data +- `bti_alpha=0.3`: Transparency for BTI shading +- `gti_alpha=0.2`: Transparency for GTI shading +- `gap_threshold=10.0`: Minimum gap size to consider as BTI +- `axis_limits=nothing`: Plot limits as `[xmin, xmax]` or `[xmin, xmax, ymin, ymax]` +- `err_method=:poisson`: Error calculation method (`:poisson`, `:gaussian`) + +# Returns +- `Tuple{Vector, Vector}`: Time bins and corresponding count rates + +# Examples +```julia +# Basic light curve +plot(events, 1.0) + +# With error bars and GTIs +plot(events, 0.5, show_errors=true, show_gtis=true) + +# Custom time range with energy filter +plot(events, 2.0, tstart=100.0, tstop=500.0, energy_filter=(2.0, 10.0)) + +# With custom axis limits +plot(events, 1.0, axis_limits=[0, 1000, 0, 100]) +``` + +# Notes +- GTI priority: explicit `gtis` > `gti_file` > `el.meta.gti` +- BTIs are calculated as gaps between GTIs exceeding `gap_threshold` +- Error bars use Poisson statistics by default +""" +@recipe function f( + el::EventList{Vector{T}, M}, + bin_size::Real = 1.0; + tstart = nothing, + tstop = nothing, + energy_filter = nothing, + show_errors = false, + show_btis = false, + show_bti = false, + show_gtis = false, + show_gti = false, + gtis = nothing, + gti_file = nothing, + gti_hdu = "GTI", + bti_alpha = 0.3, + gti_alpha = 0.2, + gap_threshold = 10.0, + axis_limits = nothing, + err_method = :poisson, +) where {T<:Real, M<:FITSMetadata} + + isempty(el.times) && error("EventList is empty") + + show_gtis = show_gtis || show_gti + show_btis = show_btis || show_bti + + # Apply filters directly to EventList data + filtered_el = deepcopy(el) # Work with a copy + + # Apply energy filter if specified + if !isnothing(energy_filter) && has_energies(filtered_el) + e_min, e_max = energy_filter + filter_energy!(e -> e_min ≤ e ≤ e_max, filtered_el) + end + + # Apply time filter if specified + if !isnothing(tstart) || !isnothing(tstop) + t_min = isnothing(tstart) ? minimum(times(filtered_el)) : tstart + t_max = isnothing(tstop) ? maximum(times(filtered_el)) : tstop + filter_time!(t -> t_min ≤ t ≤ t_max, filtered_el) + end + + # Check if we have events after filtering + if isempty(times(filtered_el)) + error("No events remain after applying filters") + end + + # Get filtered times + event_times = times(filtered_el) + time_min = minimum(event_times) + time_max = maximum(event_times) + + # Create histogram bins directly from event times + plot_tstart = isnothing(tstart) ? time_min : max(tstart, time_min) + plot_tstop = isnothing(tstop) ? time_max : min(tstop, time_max) + + # Create time bins + n_bins = max(1, Int(ceil((plot_tstop - plot_tstart) / bin_size))) + bin_edges = range(plot_tstart, plot_tstop, length=n_bins+1) + bin_centers = (bin_edges[1:end-1] + bin_edges[2:end]) / 2 + + # Bin the events + counts = fit(Histogram, event_times, bin_edges).weights + + # Calculate errors if requested + count_errors = if show_errors + if err_method == :poisson + sqrt.(max.(counts, 1)) # Poisson errors, avoid sqrt(0) + else # gaussian + sqrt.(counts) + end + else + nothing + end + + # Basic plot settings + title --> "Event List Light Curve" + xlabel --> "Time" + ylabel --> "Counts per bin" + grid --> true + minorgrid --> true + legend --> :topright + + # Axis limits handling + if !isnothing(axis_limits) + if length(axis_limits) == 4 + xmin, xmax, ymin, ymax = axis_limits + if !isnothing(xmin) || !isnothing(xmax) + xlims --> ( + isnothing(xmin) ? plot_tstart : xmin, + isnothing(xmax) ? plot_tstop : xmax, + ) + end + if !isnothing(ymin) || !isnothing(ymax) + ylims --> ( + isnothing(ymin) ? 0 : ymin, + isnothing(ymax) ? maximum(counts) * 1.1 : ymax, + ) + end + elseif length(axis_limits) == 2 + xmin, xmax = axis_limits + xlims --> ( + isnothing(xmin) ? plot_tstart : xmin, + isnothing(xmax) ? plot_tstop : xmax, + ) + else + @warn "axis_limits should be a vector of length 2 or 4: [xmin, xmax] or [xmin, xmax, ymin, ymax]" + end + end + + # Handle GTI/BTI visualization + if show_btis || show_gtis + effective_gtis = nothing + + # Priority: explicit gtis > gti_file > eventlist.meta.gti + if !isnothing(gtis) + effective_gtis = gtis + elseif !isnothing(gti_file) + try + effective_gtis, _ = read_gti_from_fits(gti_file, gti_hdu_candidates=[gti_hdu]) + catch e + @warn "Could not load GTIs from file: $e" + end + elseif has_gti(el) + effective_gtis = gti(el) + end + + if !isnothing(effective_gtis) + y_min = 0 + y_max = maximum(counts) * 1.1 + + if show_gtis + # Plot GTI intervals as green rectangles + for i in 1:size(effective_gtis, 1) + gti_start, gti_stop = effective_gtis[i, 1], effective_gtis[i, 2] + + # Only show GTI intervals that overlap with plot range + if gti_stop >= plot_tstart && gti_start <= plot_tstop + gti_start_clipped = max(gti_start, plot_tstart) + gti_stop_clipped = min(gti_stop, plot_tstop) + + @series begin + seriestype := :shape + fillcolor := :green + fillalpha := gti_alpha + linecolor := :green + linewidth := 0.5 + label := i == 1 ? "Good Time Intervals" : "" + + # Rectangle coordinates + x_coords = [gti_start_clipped, gti_stop_clipped, gti_stop_clipped, gti_start_clipped, gti_start_clipped] + y_coords = [y_min, y_min, y_max, y_max, y_min] + (x_coords, y_coords) + end + end + end + end + + if show_btis + # Calculate BTI intervals (gaps between GTIs) + btis = [] + + # Add BTI before first GTI if needed + if effective_gtis[1, 1] > plot_tstart + push!(btis, [plot_tstart, effective_gtis[1, 1]]) + end + + # Add BTIs between GTI intervals + for i in 1:(size(effective_gtis, 1) - 1) + gap_start = effective_gtis[i, 2] + gap_end = effective_gtis[i+1, 1] + if gap_end > gap_start && gap_end - gap_start > gap_threshold + push!(btis, [gap_start, gap_end]) + end + end + + # Add BTI after last GTI if needed + if effective_gtis[end, 2] < plot_tstop + push!(btis, [effective_gtis[end, 2], plot_tstop]) + end + + # Plot BTI intervals as red rectangles + for (i, bti) in enumerate(btis) + bti_start, bti_stop = bti[1], bti[2] + + if bti_stop > bti_start + @series begin + seriestype := :shape + fillcolor := :red + fillalpha := bti_alpha + linecolor := :red + linewidth := 0.5 + label := i == 1 ? "Bad Time Intervals" : "" + + # Rectangle coordinates + x_coords = [bti_start, bti_stop, bti_stop, bti_start, bti_start] + y_coords = [y_min, y_min, y_max, y_max, y_min] + (x_coords, y_coords) + end + end + end + end + else + @debug "No GTI data found for visualization" + end + end + + # Main histogram plot + if show_errors && !isnothing(count_errors) + @series begin + seriestype := :scatter + marker := :circle + markersize := 3 + markercolor := :blue + markerstrokewidth := 0.5 + yerror := count_errors + errorbar_color := :black + color := :blue + label := "Event Histogram with $(err_method) errors ($(bin_size)s bins)" + (bin_centers, counts) + end + else + @series begin + seriestype := :steppost + linewidth := 2 + color := :blue + label := "Event Histogram ($(bin_size)s bins, $(length(event_times)) events)" + + # Create step plot coordinates + step_x = Float64[] + step_y = Float64[] + + for i in eachindex(bin_centers) + bin_left = bin_edges[i] + bin_right = bin_edges[i+1] + + push!(step_x, bin_left) + push!(step_y, counts[i]) + push!(step_x, bin_right) + push!(step_y, counts[i]) + end + + (step_x, step_y) + end + end + + # Add summary info as annotation + @series begin + seriestype := :scatter + marker := :none + label := "" + annotations := [(plot_tstart + 0.05 * (plot_tstop - plot_tstart), + maximum(counts) * 0.9, + text("$(length(event_times)) events\nBin size: $(bin_size)s", + :left, 8, :gray))] + ([plot_tstart], [maximum(counts) * 0.9]) + end +end + +""" + plot(lc::LightCurve{T}; kwargs...) + +Plot a pre-computed light curve with optional properties, GTIs, and BTIs. + +# Arguments +- `lc::LightCurve{T}`: Light curve object containing binned time series data + +# Keywords +- `show_errors=false`: Display error bars if available +- `show_properties=false`: Show additional properties on secondary y-axis +- `property_name=:mean_energy`: Which property to display (if `show_properties=true`) +- `show_btis=false`: Show Bad Time Intervals as red shaded regions +- `show_bti=false`: Alias for `show_btis` +- `show_gtis=false`: Show Good Time Intervals as green shaded regions +- `show_gti=false`: Alias for `show_gtis` +- `gtis=nothing`: GTI matrix to use (overrides metadata GTIs) +- `gti_file=nothing`: FITS file containing GTI extension +- `gti_hdu="GTI"`: HDU name for GTI data +- `bti_alpha=0.3`: Transparency for BTI shading +- `gti_alpha=0.2`: Transparency for GTI shading +- `axis_limits=nothing`: Plot limits as `[xmin, xmax]` or `[xmin, xmax, ymin, ymax]` + +# Returns +- `Tuple{Vector, Vector}`: Time bins and corresponding count rates + +# Examples +```julia +# Basic light curve plot +plot(lightcurve) + +# With error bars and mean energy overlay +plot(lightcurve, show_errors=true, show_properties=true, property_name=:mean_energy) + +# Show GTIs with custom transparency +plot(lightcurve, show_gtis=true, gti_alpha=0.4) +``` + +# Notes +- Errors are calculated on-demand if not already present +- Properties must exist in `lc.properties` to be displayed +- GTI priority: explicit `gtis` > `gti_file` > `lc.metadata.extra["gti_bounds"]` +""" +@recipe function f( + lc::LightCurve{T}; + tstart = nothing, + tstop = nothing, + show_errors = false, + show_properties = false, + property_name = :mean_energy, + show_btis = false, + show_bti = false, + show_gtis = false, + show_gti = false, + gtis = nothing, + gti_file = nothing, + gti_hdu = "GTI", + bti_alpha = 0.3, + gti_alpha = 0.2, + axis_limits = nothing, +) where {T} + + show_gtis = show_gtis || show_gti + show_btis = show_btis || show_bti + + if show_errors && isnothing(lc.count_error) + calculate_errors!(lc) + end + + # Apply time filtering if specified + time_mask = trues(length(lc.time)) + + if !isnothing(tstart) + time_mask .&= (lc.time .>= tstart) + end + + if !isnothing(tstop) + time_mask .&= (lc.time .<= tstop) + end + + # Check if any data remains after filtering + if !any(time_mask) + @warn "No data points in specified time range [$tstart, $tstop]" + # Create empty plot + title --> "Light Curve (No Data in Range)" + xlabel --> "Time (s)" + ylabel --> "Counts" + @series begin + seriestype := :scatter + marker := :none + label := "" + ([0], [0]) + end + return + end + + # Filter the data + filtered_time = lc.time[time_mask] + filtered_counts = lc.counts[time_mask] + filtered_errors = isnothing(lc.count_error) ? nothing : lc.count_error[time_mask] + + # Convert to matrix format + lc_matrix = if isnothing(filtered_errors) + hcat(filtered_time, filtered_counts, zeros(length(filtered_time))) + else + hcat(filtered_time, filtered_counts, filtered_errors) + end + + title --> "Light Curve" + xlabel --> "Time (s)" + ylabel --> "Counts" + grid --> true + minorgrid --> true + legend --> :bottomright + + # Axis limits handling + if !isnothing(axis_limits) + if length(axis_limits) == 4 + xmin, xmax, ymin, ymax = axis_limits + + if !isnothing(xmin) || !isnothing(xmax) + xlims --> ( + isnothing(xmin) ? minimum(lc_matrix[:, 1]) : xmin, + isnothing(xmax) ? maximum(lc_matrix[:, 1]) : xmax, + ) + end + + if !isnothing(ymin) || !isnothing(ymax) + ylims --> ( + isnothing(ymin) ? minimum(lc_matrix[:, 2]) : ymin, + isnothing(ymax) ? maximum(lc_matrix[:, 2]) : ymax, + ) + end + elseif length(axis_limits) == 2 + xmin, xmax = axis_limits + xlims --> ( + isnothing(xmin) ? minimum(lc_matrix[:, 1]) : xmin, + isnothing(xmax) ? maximum(lc_matrix[:, 1]) : xmax, + ) + else + @warn "axis_limits should be a vector of length 2 or 4: [xmin, xmax] or [xmin, xmax, ymin, ymax]" + end + else + # Set automatic limits based on filtered data + xlims --> (minimum(filtered_time), maximum(filtered_time)) + end + + # Determine plot time range (use filtered data or tstart/tstop) + plot_tstart = isnothing(tstart) ? minimum(filtered_time) : tstart + plot_tstop = isnothing(tstop) ? maximum(filtered_time) : tstop + + # Handle GTI/BTI visualization + if show_btis || show_gtis + effective_gtis = nothing + + if !isnothing(gtis) + effective_gtis = gtis + elseif !isnothing(gti_file) + try + effective_gtis = load_gtis(gti_file, gti_hdu) + catch e + @warn "Could not load GTIs from file: $e" + end + elseif haskey(lc.metadata.extra, "gti") && !isnothing(lc.metadata.extra["gti"]) + effective_gtis = lc.metadata.extra["gti"] + end + + if !isnothing(effective_gtis) + y_min, y_max = extrema(lc_matrix[:, 2]) + + if show_gtis + n_gtis = size(effective_gtis, 1) + gti_x = Vector{Float64}(undef, n_gtis * 6) + gti_y = Vector{Float64}(undef, n_gtis * 6) + gti_idx = 0 + + @inbounds for i = 1:n_gtis + gti_start, gti_stop = effective_gtis[i, 1], effective_gtis[i, 2] + + if gti_stop >= plot_tstart && gti_start <= plot_tstop + gti_start_clipped = max(gti_start, plot_tstart) + gti_stop_clipped = min(gti_stop, plot_tstop) + + base_idx = gti_idx * 6 + gti_x[base_idx+1] = gti_start_clipped + gti_x[base_idx+2] = gti_stop_clipped + gti_x[base_idx+3] = gti_stop_clipped + gti_x[base_idx+4] = gti_start_clipped + gti_x[base_idx+5] = gti_start_clipped + gti_x[base_idx+6] = NaN + + gti_y[base_idx+1] = y_min + gti_y[base_idx+2] = y_min + gti_y[base_idx+3] = y_max + gti_y[base_idx+4] = y_max + gti_y[base_idx+5] = y_min + gti_y[base_idx+6] = NaN + + gti_idx += 1 + end + end + + if gti_idx > 0 + resize!(gti_x, gti_idx * 6) + resize!(gti_y, gti_idx * 6) + + @series begin + seriestype := :shape + fillcolor := :green + fillalpha := gti_alpha + linecolor := :green + linewidth := 0.5 + label := "Good Time Intervals" + gti_x, gti_y + end + end + end + + if show_btis + btis = get_btis(effective_gtis, plot_tstart, plot_tstop) + + if !isempty(btis) + n_btis = size(btis, 1) + bti_x = Vector{Float64}(undef, n_btis * 6) + bti_y = Vector{Float64}(undef, n_btis * 6) + + @inbounds for i = 1:n_btis + bti_start, bti_stop = btis[i, 1], btis[i, 2] + + base_idx = (i - 1) * 6 + bti_x[base_idx+1] = bti_start + bti_x[base_idx+2] = bti_stop + bti_x[base_idx+3] = bti_stop + bti_x[base_idx+4] = bti_start + bti_x[base_idx+5] = bti_start + bti_x[base_idx+6] = NaN + + bti_y[base_idx+1] = y_min + bti_y[base_idx+2] = y_min + bti_y[base_idx+3] = y_max + bti_y[base_idx+4] = y_max + bti_y[base_idx+5] = y_min + bti_y[base_idx+6] = NaN + end + + @series begin + seriestype := :shape + fillcolor := :red + fillalpha := bti_alpha + linecolor := :red + linewidth := 0.5 + label := "Bad Time Intervals" + bti_x, bti_y + end + end + end + end + end + + # Handle properties display (filter properties too) + if show_properties + prop_idx = findfirst(p -> p.name == property_name, lc.properties) + if !isnothing(prop_idx) + prop = lc.properties[prop_idx] + filtered_prop_values = prop.values[time_mask] + + prop_matrix = hcat(filtered_time, filtered_prop_values) + + @series begin + yaxis := :right + ylabel := "$(prop.name) ($(prop.unit))" + seriestype := :line + color := :red + linewidth := 1.5 + label := String(prop.name) + prop_matrix[:, 1], prop_matrix[:, 2] + end + end + end + + # Main light curve plot + if show_errors && !isnothing(filtered_errors) + seriestype --> :scatter + marker --> :circle + markersize --> 3 + markercolor --> :blue + markerstrokewidth --> 0.5 + yerror --> lc_matrix[:, 3] + errorbar_color --> :black + color --> :blue + label --> "Light Curve with $(lc.err_method) errors" + else + seriestype --> :steppost + linewidth --> 1.5 + color --> :blue + label --> "Light Curve (bin size: $(lc.metadata.bin_size)s)" + end + + return lc_matrix[:, 1], lc_matrix[:, 2] +end + +""" + plot(lc_segments::Vector{<:LightCurve}; kwargs...) + +Plot multiple light curve segments with optional segment boundaries and individual coloring. + +# Arguments +- `lc_segments::Vector{<:LightCurve}`: Vector of light curve segments to plot + +# Keywords +- `show_errors=false`: Display error bars for each segment +- `show_segment_boundaries=true`: Show vertical dashed lines at segment boundaries +- `segment_colors=nothing`: Custom colors for each segment (defaults to standard palette) +- `axis_limits=nothing`: Plot limits as `[xmin, xmax]` or `[xmin, xmax, ymin, ymax]` + +# Returns +- Multiple plot series, one per segment + +# Examples +```julia +# Basic segmented light curve +plot(segments) + +# With custom colors and no boundaries +plot(segments, segment_colors=[:red, :blue, :green], show_segment_boundaries=false) + +# With error bars and custom limits +plot(segments, show_errors=true, axis_limits=[0, 1000, 0, 50]) +``` + +# Notes +- Default color palette cycles through 8 colors for segments +- Segment boundaries are drawn at the start of each segment after the first +- Each segment can have different bin sizes and error methods +""" +@recipe function f( + lc_segments::Vector{<:LightCurve}; + show_errors = false, + show_segment_boundaries = true, + segment_colors = nothing, + axis_limits = nothing, +) + + title --> "Segmented Light Curve" + xlabel --> "Time (s)" + ylabel --> "Counts" + grid --> true + minorgrid --> true + legend --> :bottomright + + # Axis limits handling + if !isnothing(axis_limits) + if length(axis_limits) == 4 + xmin, xmax, ymin, ymax = axis_limits + + # Get overall data bounds + all_times = vcat([lc.time for lc in lc_segments]...) + all_counts = vcat([lc.counts for lc in lc_segments]...) + + if !isnothing(xmin) || !isnothing(xmax) + xlims --> ( + isnothing(xmin) ? minimum(all_times) : xmin, + isnothing(xmax) ? maximum(all_times) : xmax, + ) + end + + if !isnothing(ymin) || !isnothing(ymax) + ylims --> ( + isnothing(ymin) ? minimum(all_counts) : ymin, + isnothing(ymax) ? maximum(all_counts) : ymax, + ) + end + elseif length(axis_limits) == 2 + xmin, xmax = axis_limits + all_times = vcat([lc.time for lc in lc_segments]...) + xlims --> ( + isnothing(xmin) ? minimum(all_times) : xmin, + isnothing(xmax) ? maximum(all_times) : xmax, + ) + else + @warn "axis_limits should be a vector of length 2 or 4: [xmin, xmax] or [xmin, xmax, ymin, ymax]" + end + end + + default_colors = [:blue, :red, :green, :orange, :purple, :brown, :pink, :gray] + colors = isnothing(segment_colors) ? default_colors : segment_colors + n_colors = length(colors) + + boundaries = Vector{Float64}() + + for (i, lc) in enumerate(lc_segments) + color = colors[((i-1)%n_colors)+1] + + if show_errors && isnothing(lc.count_error) + calculate_errors!(lc) + end + + lc_matrix = hcat(lc.time, lc.counts, lc.count_error) + + @series begin + if show_errors + seriestype := :scatter + marker := :circle + markersize := 3 + markerstrokewidth := 0.5 + yerror := lc_matrix[:, 3] + errorbar_color := color + else + seriestype := :steppost + linewidth := 1.5 + end + color := color + label := "Segment $i" + + lc_matrix[:, 1], lc_matrix[:, 2] + end + + if show_segment_boundaries && i > 1 + push!(boundaries, minimum(lc.time)) + end + end + + if show_segment_boundaries && !isempty(boundaries) + @series begin + seriestype := :vline + color := :black + linestyle := :dash + linewidth := 1 + alpha := 0.7 + label := "Segment boundaries" + boundaries + end + end +end +""" + plot(lc::LightCurve{T}, new_binsize::Real; kwargs...) + +Plot a light curve after rebinning it to a new bin size. + +# Arguments +- `lc::LightCurve{T}`: Original light curve to rebin and plot +- `new_binsize::Real`: New bin size in seconds (must be larger than current bin size) + +# Keywords +- `show_errors=false`: Display error bars using rebinned error estimates +- `show_properties=false`: Show additional properties on secondary y-axis +- `property_name=:mean_energy`: Which property to display (if `show_properties=true`) +- `show_btis=false`: Show Bad Time Intervals as red shaded regions +- `show_bti=false`: Alias for `show_btis` +- `show_gtis=false`: Show Good Time Intervals as green shaded regions +- `show_gti=false`: Alias for `show_gtis` +- `gtis=nothing`: GTI matrix to use (overrides metadata GTIs) +- `gti_file=nothing`: FITS file containing GTI extension +- `gti_hdu="GTI"`: HDU name for GTI data +- `bti_alpha=0.3`: Transparency for BTI shading +- `gti_alpha=0.2`: Transparency for GTI shading +- `axis_limits=nothing`: Plot limits as `[xmin, xmax]` or `[xmin, xmax, ymin, ymax]` +- `show_original=false`: Overlay original light curve for comparison +- `original_alpha=0.3`: Transparency for original light curve overlay + +# Returns +- `Tuple{Vector, Vector}`: Rebinned time bins and corresponding count rates + +# Examples +```julia +# Basic rebinned light curve +plot(lc, 100.0) # Rebin to 100s + +# With error bars and GTIs +plot(lc, 50.0, show_errors=true, show_gtis=true) + +# Compare with original +plot(lc, 200.0, show_original=true, original_alpha=0.4) + +# With properties overlay +plot(lc, 30.0, show_properties=true, property_name=:mean_energy) +``` + +# Notes +- The new bin size must be larger than the current bin size +- Original light curve data is preserved; only the plot shows rebinned data +- Error bars are recalculated for the rebinned light curve +- GTI/BTI regions are preserved from the original light curve metadata +""" +@recipe function f( + lc::LightCurve{T}, + new_binsize::Real; + show_errors = false, + show_properties = false, + property_name = :mean_energy, + show_btis = false, + show_bti = false, + show_gtis = false, + show_gti = false, + gtis = nothing, + gti_file = nothing, + gti_hdu = "GTI", + bti_alpha = 0.3, + gti_alpha = 0.2, + axis_limits = nothing, + show_original = false, + original_alpha = 0.3, +) where {T} + + # Validate new bin size + new_binsize <= lc.metadata.bin_size && throw( + ArgumentError( + "New bin size ($new_binsize s) must be larger than current bin size ($(lc.metadata.bin_size) s)", + ), + ) + + # Create rebinned light curve + rebinned_lc = rebin(lc, new_binsize) + + show_gtis = show_gtis || show_gti + show_btis = show_btis || show_bti + + if show_errors && isnothing(rebinned_lc.count_error) + calculate_errors!(rebinned_lc) + end + + # Convert to matrix format + lc_matrix = hcat(rebinned_lc.time, rebinned_lc.counts, rebinned_lc.count_error) + + title --> "Rebinned Light Curve ($(lc.metadata.bin_size)s → $(new_binsize)s)" + xlabel --> "Time (s)" + ylabel --> "Counts" + grid --> true + minorgrid --> true + legend --> :bottomright + + # Axis limits handling + if !isnothing(axis_limits) + if length(axis_limits) == 4 + xmin, xmax, ymin, ymax = axis_limits + + if !isnothing(xmin) || !isnothing(xmax) + xlims --> ( + isnothing(xmin) ? minimum(lc_matrix[:, 1]) : xmin, + isnothing(xmax) ? maximum(lc_matrix[:, 1]) : xmax, + ) + end + + if !isnothing(ymin) || !isnothing(ymax) + ylims --> ( + isnothing(ymin) ? minimum(lc_matrix[:, 2]) : ymin, + isnothing(ymax) ? maximum(lc_matrix[:, 2]) : ymax, + ) + end + elseif length(axis_limits) == 2 + xmin, xmax = axis_limits + xlims --> ( + isnothing(xmin) ? minimum(lc_matrix[:, 1]) : xmin, + isnothing(xmax) ? maximum(lc_matrix[:, 1]) : xmax, + ) + else + @warn "axis_limits should be a vector of length 2 or 4: [xmin, xmax] or [xmin, xmax, ymin, ymax]" + end + end + + plot_tstart, plot_tstop = rebinned_lc.metadata.time_range + + # Show original light curve for comparison + if show_original + original_matrix = hcat(lc.time, lc.counts, lc.count_error) + + @series begin + seriestype := :steppost + linewidth := 1 + color := :gray + alpha := original_alpha + label := "Original ($(lc.metadata.bin_size)s bins)" + original_matrix[:, 1], original_matrix[:, 2] + end + end + + # Handle GTI/BTI visualization (use original metadata) + if show_btis || show_gtis + effective_gtis = nothing + + if !isnothing(gtis) + effective_gtis = gtis + elseif !isnothing(gti_file) + try + effective_gtis = load_gtis(gti_file, gti_hdu) + catch e + @warn "Could not load GTIs from file: $e" + end + elseif haskey(lc.metadata.extra, "gti_applied") && + haskey(lc.metadata.extra, "gti_bounds") + gti_bounds = lc.metadata.extra["gti_bounds"] + effective_gtis = reshape(gti_bounds, 1, 2) + end + + if !isnothing(effective_gtis) + y_min, y_max = extrema(lc_matrix[:, 2]) + + if show_gtis + n_gtis = size(effective_gtis, 1) + gti_x = Vector{Float64}(undef, n_gtis * 6) + gti_y = Vector{Float64}(undef, n_gtis * 6) + gti_idx = 0 + + @inbounds for i = 1:n_gtis + gti_start, gti_stop = effective_gtis[i, 1], effective_gtis[i, 2] + + if gti_stop >= plot_tstart && gti_start <= plot_tstop + gti_start = max(gti_start, plot_tstart) + gti_stop = min(gti_stop, plot_tstop) + + base_idx = gti_idx * 6 + gti_x[base_idx+1] = gti_start + gti_x[base_idx+2] = gti_stop + gti_x[base_idx+3] = gti_stop + gti_x[base_idx+4] = gti_start + gti_x[base_idx+5] = gti_start + gti_x[base_idx+6] = NaN + + gti_y[base_idx+1] = y_min + gti_y[base_idx+2] = y_min + gti_y[base_idx+3] = y_max + gti_y[base_idx+4] = y_max + gti_y[base_idx+5] = y_min + gti_y[base_idx+6] = NaN + + gti_idx += 1 + end + end + + if gti_idx > 0 + resize!(gti_x, gti_idx * 6) + resize!(gti_y, gti_idx * 6) + + @series begin + seriestype := :shape + fillcolor := :green + fillalpha := gti_alpha + linecolor := :green + linewidth := 0.5 + label := "Good Time Intervals" + gti_x, gti_y + end + end + end + + if show_btis + btis = get_btis(effective_gtis, plot_tstart, plot_tstop) + + if !isempty(btis) + n_btis = size(btis, 1) + bti_x = Vector{Float64}(undef, n_btis * 6) + bti_y = Vector{Float64}(undef, n_btis * 6) + + @inbounds for i = 1:n_btis + bti_start, bti_stop = btis[i, 1], btis[i, 2] + + base_idx = (i - 1) * 6 + bti_x[base_idx+1] = bti_start + bti_x[base_idx+2] = bti_stop + bti_x[base_idx+3] = bti_stop + bti_x[base_idx+4] = bti_start + bti_x[base_idx+5] = bti_start + bti_x[base_idx+6] = NaN + + bti_y[base_idx+1] = y_min + bti_y[base_idx+2] = y_min + bti_y[base_idx+3] = y_max + bti_y[base_idx+4] = y_max + bti_y[base_idx+5] = y_min + bti_y[base_idx+6] = NaN + end + + @series begin + seriestype := :shape + fillcolor := :red + fillalpha := bti_alpha + linecolor := :red + linewidth := 0.5 + label := "Bad Time Intervals" + bti_x, bti_y + end + end + end + end + end + + # Handle properties display + if show_properties + prop_idx = findfirst(p -> p.name == property_name, rebinned_lc.properties) + if !isnothing(prop_idx) + prop = rebinned_lc.properties[prop_idx] + + prop_matrix = hcat(rebinned_lc.time, prop.values) + + @series begin + yaxis := :right + ylabel := "$(prop.name) ($(prop.unit))" + seriestype := :line + color := :red + linewidth := 1.5 + label := String(prop.name) + prop_matrix[:, 1], prop_matrix[:, 2] + end + end + end + + if show_errors + seriestype --> :scatter + marker --> :circle + markersize --> 3 + markercolor --> :blue + markerstrokewidth --> 0.5 + yerror --> lc_matrix[:, 3] + errorbar_color --> :black + color --> :blue + label --> "Rebinned LC ($(new_binsize)s) with $(rebinned_lc.err_method) errors" + else + seriestype --> :steppost + linewidth --> 1.5 + color --> :blue + label --> "Rebinned LC ($(new_binsize)s bins)" + end + + return lc_matrix[:, 1], lc_matrix[:, 2] +end diff --git a/test/runtests.jl b/test/runtests.jl index 12ca46a..2c58ac4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ using FFTW, Distributions, Statistics, StatsBase, HDF5, FITSIO using Logging ,LinearAlgebra using CFITSIO using Random +using Plots include("test_fourier.jl") @testset "GTI" begin @@ -15,4 +16,9 @@ end @testset "lightcurve" begin include("test_lightcurve.jl") -end \ No newline at end of file +end + +@testset "recipes" begin + include("test_plotting/test_plots_recipes_lightcurve.jl") + include("test_plotting/test_plots_recipes_gti.jl") +end diff --git a/test/test_plotting/test_plots_recipes_gti.jl b/test/test_plotting/test_plots_recipes_gti.jl new file mode 100644 index 0000000..2909d37 --- /dev/null +++ b/test/test_plotting/test_plots_recipes_gti.jl @@ -0,0 +1,83 @@ +# BTI Histogram Recipe Tests +let + # Suppress prints output during tests(prints to avoid length limit) + original_stdout = stdout + redirect_stdout(devnull) + + try + # Test 1: No BTIs + times = collect(0.0:1.0:100.0) + energies = rand(length(times)) .* 100 + gti_full = reshape([-1.0, 101.0], 1, 2) + metadata = FITSMetadata("test_file.fits", 1, "keV", Dict{String,Vector}(), + Dict{String,Any}(), gti_full, "GTI",nothing, nothing, nothing, nothing, nothing, nothing) + el_no_btis = EventList(times, energies, metadata) + p1 = plot(BTIAnalysisPlot(el_no_btis), bti_analysis=true) + @test p1 isa Plots.Plot + + # Test 2: With BTIs + times2 = collect(0.0:0.5:50.0) + energies2 = rand(length(times2)) .* 100 + gti_gaps = [0.0 10.0; 20.0 30.0; 40.0 50.0] # Creates BTIs at [10-20], [30-40] + metadata2 = FITSMetadata("test_file2.fits", 1, "keV", Dict{String,Vector}(), + Dict{String,Any}(), gti_gaps, "GTI",nothing, nothing, nothing, nothing, nothing, nothing) + el_with_btis = EventList(times2, energies2, metadata2) + p2 = plot(BTIAnalysisPlot(el_with_btis), bti_analysis=true) + @test p2 isa Plots.Plot + + # Test 3: Parameter variations with BTIs + p3 = plot(BTIAnalysisPlot(el_with_btis), bti_analysis=true, nbins=50) + @test p3 isa Plots.Plot + p4 = plot(BTIAnalysisPlot(el_with_btis), bti_analysis=true, nbins=20, xlims_range=(5, 15)) + @test p4 isa Plots.Plot + p5 = plot(BTIAnalysisPlot(el_with_btis), bti_analysis=true, ylims_range=(0.5, 5)) + @test p5 isa Plots.Plot + + # Test 4: Complex GTI pattern with multiple BTIs + times3 = collect(0.0:0.1:30.0) + energies3 = rand(length(times3)) .* 100 + # Creates multiple BTIs: [5-6], [10-12], [15-18], [22-25] + gti_complex = [0.0 5.0; 6.0 10.0; 12.0 15.0; 18.0 22.0; 25.0 30.0] + metadata3 = FITSMetadata("test_file3.fits", 1, "keV", Dict{String,Vector}(), + Dict{String,Any}(), gti_complex, "GTI",nothing, nothing, nothing, nothing, nothing, nothing) + el_complex = EventList(times3, energies3, metadata3) + p6 = plot(BTIAnalysisPlot(el_complex), bti_analysis=true) + @test p6 isa Plots.Plot + + # Test 5: Single very large BTI + times4 = [0.0, 1.0, 100.0, 101.0] + energies4 = [10.0, 20.0, 30.0, 40.0] + gti_large_gap = [0.0 1.0; 100.0 101.0] # Creates 99s BTI gap + metadata4 = FITSMetadata("test_file4.fits", 1, "keV", Dict{String,Vector}(), + Dict{String,Any}(), gti_large_gap, "GTI",nothing, nothing, nothing, nothing, nothing, nothing) + el_large_gap = EventList(times4, energies4, metadata4) + p8 = plot(BTIAnalysisPlot(el_large_gap), bti_analysis=true) + @test p8 isa Plots.Plot + + # Test 6: Edge case - very short BTIs only + times5 = collect(0.0:0.01:10.0) + energies5 = rand(length(times5)) .* 100 + # Creates very short BTIs: [1-1.1], [2-2.05], [3-3.02] + gti_short = [0.0 1.0; 1.1 2.0; 2.05 3.0; 3.02 10.0] + metadata5 = FITSMetadata("test_file5.fits", 1, "keV", Dict{String,Vector}(), + Dict{String,Any}(), gti_short, "GTI",nothing, nothing, nothing, nothing, nothing, nothing) + el_short_btis = EventList(times5, energies5, metadata5) + p9 = plot(BTIAnalysisPlot(el_short_btis), bti_analysis=true) + @test p9 isa Plots.Plot + + # Test 7: Mixed BTI lengths (short and long) + times6 = collect(0.0:0.1:100.0) + energies6 = rand(length(times6)) .* 100 + # Creates mixed BTIs: [5-5.1] (0.1s), [10-20] (10s), [30-30.5] (0.5s), [50-80] (30s) + gti_mixed = [0.0 5.0; 5.1 10.0; 20.0 30.0; 30.5 50.0; 80.0 100.0] + metadata6 = FITSMetadata("test_file6.fits", 1, "keV", Dict{String,Vector}(), + Dict{String,Any}(), gti_mixed, "GTI",nothing, nothing, nothing, nothing, nothing, nothing) + el_mixed_btis = EventList(times6, energies6, metadata6) + p10 = plot(BTIAnalysisPlot(el_mixed_btis), bti_analysis=true) + @test p10 isa Plots.Plot + + finally + # Always restore stdout + redirect_stdout(original_stdout) + end +end \ No newline at end of file diff --git a/test/test_plotting/test_plots_recipes_lightcurve.jl b/test/test_plotting/test_plots_recipes_lightcurve.jl new file mode 100644 index 0000000..dbcdc70 --- /dev/null +++ b/test/test_plotting/test_plots_recipes_lightcurve.jl @@ -0,0 +1,345 @@ +# EventList plotting tests +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + + @test plot(el, 1.0) isa Plots.Plot + @test plot(el, 2.0) isa Plots.Plot +end + +# Empty EventList error handling +let + empty_el = EventList(Float64[], Float64[]) + @test_throws ErrorException plot(empty_el, 1.0) +end + +# Time range filtering +let + times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0] + el = EventList(times, energies) + + @test plot(el, 1.0, tstart=2.0, tstop=5.0) isa Plots.Plot + @test plot(el, 1.0, tstart=2.0) isa Plots.Plot + @test plot(el, 1.0, tstop=4.0) isa Plots.Plot +end + +# Energy filtering +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + + @test plot(el, 1.0, energy_filter=(15.0, 45.0)) isa Plots.Plot + @test plot(el, 1.0, energy_filter=(10.0, 30.0)) isa Plots.Plot +end + +# Error display options +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + + @test plot(el, 1.0, show_errors=true) isa Plots.Plot + @test plot(el, 1.0, show_errors=false) isa Plots.Plot + @test plot(el, 1.0, show_errors=true, err_method=:poisson) isa Plots.Plot +end + +# GTI visualization +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + gti_matrix = [1.0 3.0; 4.0 6.0] + meta = FITSMetadata{Dict{String,Any}}( + "test.fits", 2, "ENERGY", Dict{String,Vector}(), Dict{String,Any}(), + gti_matrix, "GTI",nothing, nothing, nothing, nothing, nothing, nothing + ) + el = EventList(times, energies, meta) + + @test plot(el, 1.0, show_gtis=true) isa Plots.Plot + @test plot(el, 1.0, show_gti=true) isa Plots.Plot + @test plot(el, 1.0, show_gtis=true, gti_alpha=0.5) isa Plots.Plot +end + +# BTI visualization +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + gti_matrix = [1.0 3.0; 4.0 6.0] + meta = FITSMetadata{Dict{String,Any}}( + "test.fits", 2, "ENERGY", Dict{String,Vector}(), Dict{String,Any}(), + gti_matrix, "GTI",nothing, nothing, nothing, nothing, nothing, nothing + ) + el = EventList(times, energies, meta) + + @test plot(el, 1.0, show_btis=true) isa Plots.Plot + @test plot(el, 1.0, show_bti=true) isa Plots.Plot + @test plot(el, 1.0, show_btis=true, bti_alpha=0.4) isa Plots.Plot +end + +# External GTI support +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + external_gtis = [1.0 2.5; 3.5 5.0] + + @test plot(el, 1.0, show_gtis=true, gtis=external_gtis) isa Plots.Plot + @test plot(el, 1.0, show_btis=true, gtis=external_gtis) isa Plots.Plot +end + +# Axis limits configuration +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + + @test plot(el, 1.0, axis_limits=[1.5, 4.5, 0, 10]) isa Plots.Plot + @test plot(el, 1.0, axis_limits=[1.5, 4.5]) isa Plots.Plot + @test plot(el, 1.0, axis_limits=[nothing, 4.5, 0, nothing]) isa Plots.Plot +end + +# Invalid axis limits warning +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + + @test plot(el, 1.0, axis_limits=[1.5, 4.5, 0]) isa Plots.Plot +end + +# LightCurve plotting +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + counts = [10, 15, 20, 12, 8] + errors = [3.0, 4.0, 5.0, 3.5, 2.8] + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 5.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc = LightCurve(times, 1.0, counts, errors, nothing, EventProperty{Float64}[], metadata, :poisson) + + @test plot(lc) isa Plots.Plot + @test plot(lc, show_errors=true) isa Plots.Plot + @test plot(lc, show_errors=false) isa Plots.Plot +end + +# LightCurve without errors +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + counts = [10, 15, 20, 12, 8] + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 5.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc = LightCurve(times, 1.0, counts, nothing, nothing, EventProperty{Float64}[], metadata, :poisson) + + @test plot(lc, show_errors=true) isa Plots.Plot +end + +# LightCurve with event properties +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + counts = [10, 15, 20, 12, 8] + errors = [3.0, 4.0, 5.0, 3.5, 2.8] + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 5.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + properties = [EventProperty(:mean_energy, [25.0, 30.0, 35.0, 28.0, 22.0], "keV")] + lc = LightCurve(times, 1.0, counts, errors, nothing, properties, metadata, :poisson) + + @test plot(lc, show_properties=true) isa Plots.Plot + @test plot(lc, show_properties=true, property_name=:mean_energy) isa Plots.Plot +end + +# Non-existent property handling +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + counts = [10, 15, 20, 12, 8] + errors = [3.0, 4.0, 5.0, 3.5, 2.8] + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 5.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc = LightCurve(times, 1.0, counts, errors, nothing, EventProperty{Float64}[], metadata, :poisson) + + @test plot(lc, show_properties=true, property_name=:nonexistent) isa Plots.Plot +end + +# LightCurve with GTI metadata +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + counts = [10, 15, 20, 12, 8] + errors = [3.0, 4.0, 5.0, 3.5, 2.8] + extra_data = Dict{String,Any}("gti_applied" => true, "gti_bounds" => [1.0, 5.0]) + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 5.0), 1.0, Dict{String,Any}[], extra_data) + lc = LightCurve(times, 1.0, counts, errors, nothing, EventProperty{Float64}[], metadata, :poisson) + + @test plot(lc, show_gtis=true) isa Plots.Plot + @test plot(lc, show_btis=true) isa Plots.Plot +end + +# LightCurve axis limits +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + counts = [10, 15, 20, 12, 8] + errors = [3.0, 4.0, 5.0, 3.5, 2.8] + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 5.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc = LightCurve(times, 1.0, counts, errors, nothing, EventProperty{Float64}[], metadata, :poisson) + + @test plot(lc, axis_limits=[1.5, 4.5, 5, 25]) isa Plots.Plot + @test plot(lc, axis_limits=[1.5, 4.5]) isa Plots.Plot +end + +# Segmented LightCurve plotting +let + times1 = [1.0, 2.0, 3.0] + counts1 = [10, 15, 20] + errors1 = [3.0, 4.0, 5.0] + metadata1 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 3.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc1 = LightCurve(times1, 1.0, counts1, errors1, nothing, EventProperty{Float64}[], metadata1, :poisson) + + times2 = [4.0, 5.0, 6.0] + counts2 = [12, 8, 18] + errors2 = [3.5, 2.8, 4.2] + metadata2 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (4.0, 6.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc2 = LightCurve(times2, 1.0, counts2, errors2, nothing, EventProperty{Float64}[], metadata2, :poisson) + + segments = [lc1, lc2] + + @test plot(segments) isa Plots.Plot + @test plot(segments, show_errors=true) isa Plots.Plot + @test plot(segments, show_segment_boundaries=true) isa Plots.Plot + @test plot(segments, show_segment_boundaries=false) isa Plots.Plot +end + +# Segmented LightCurve with custom colors +let + times1 = [1.0, 2.0, 3.0] + counts1 = [10, 15, 20] + errors1 = [3.0, 4.0, 5.0] + metadata1 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 3.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc1 = LightCurve(times1, 1.0, counts1, errors1, nothing, EventProperty{Float64}[], metadata1, :poisson) + + times2 = [4.0, 5.0, 6.0] + counts2 = [12, 8, 18] + errors2 = [3.5, 2.8, 4.2] + metadata2 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (4.0, 6.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc2 = LightCurve(times2, 1.0, counts2, errors2, nothing, EventProperty{Float64}[], metadata2, :poisson) + + segments = [lc1, lc2] + custom_colors = [:red, :green] + + @test plot(segments, segment_colors=custom_colors) isa Plots.Plot +end + +# Segmented LightCurve axis limits +let + times1 = [1.0, 2.0, 3.0] + counts1 = [10, 15, 20] + errors1 = [3.0, 4.0, 5.0] + metadata1 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 3.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc1 = LightCurve(times1, 1.0, counts1, errors1, nothing, EventProperty{Float64}[], metadata1, :poisson) + + times2 = [4.0, 5.0, 6.0] + counts2 = [12, 8, 18] + errors2 = [3.5, 2.8, 4.2] + metadata2 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (4.0, 6.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc2 = LightCurve(times2, 1.0, counts2, errors2, nothing, EventProperty{Float64}[], metadata2, :poisson) + + segments = [lc1, lc2] + + @test plot(segments, axis_limits=[1.5, 5.5, 5, 25]) isa Plots.Plot + @test plot(segments, axis_limits=[1.5, 5.5]) isa Plots.Plot +end + +# Single segment handling +let + times1 = [1.0, 2.0, 3.0] + counts1 = [10, 15, 20] + errors1 = [3.0, 4.0, 5.0] + metadata1 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 3.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc1 = LightCurve(times1, 1.0, counts1, errors1, nothing, EventProperty{Float64}[], metadata1, :poisson) + + segments = [lc1] + + @test plot(segments) isa Plots.Plot + @test plot(segments, show_segment_boundaries=true) isa Plots.Plot +end + +# Segmented LightCurve without errors +let + times1 = [1.0, 2.0, 3.0] + counts1 = [10, 15, 20] + metadata1 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (1.0, 3.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc1 = LightCurve(times1, 1.0, counts1, nothing, nothing, EventProperty{Float64}[], metadata1, :poisson) + + times2 = [4.0, 5.0, 6.0] + counts2 = [12, 8, 18] + metadata2 = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (4.0, 6.0), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc2 = LightCurve(times2, 1.0, counts2, nothing, nothing, EventProperty{Float64}[], metadata2, :poisson) + + segments = [lc1, lc2] + + @test plot(segments, show_errors=true) isa Plots.Plot +end + +# Color cycling for multiple segments +let + segments = LightCurve[] + for i in 1:10 + times = [Float64(i), Float64(i+1)] + counts = [10, 15] + errors = [3.0, 4.0] + metadata = LightCurveMetadata("TEST", "TEST_INST", "TEST_OBJ", 58000.0, (Float64(i), Float64(i+1)), 1.0, Dict{String,Any}[], Dict{String,Any}()) + lc = LightCurve(times, 1.0, counts, errors, nothing, EventProperty{Float64}[], metadata, :poisson) + push!(segments, lc) + end + + @test plot(segments) isa Plots.Plot +end + +# GTI file loading error handling +let + times = [1.0, 2.0, 3.0, 4.0, 5.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0] + el = EventList(times, energies) + + @test plot(el, 1.0, show_gtis=true, gti_file="nonexistent.fits") isa Plots.Plot +end + +# Complex GTI/BTI visualization +let + times = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] + energies = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0] + gti_matrix = [1.0 2.5; 3.5 5.0; 6.5 8.0] + meta = FITSMetadata{Dict{String,Any}}( + "test.fits", 2, "ENERGY", Dict{String,Vector}(), Dict{String,Any}(), + gti_matrix, "GTI",nothing, nothing, nothing, nothing, nothing, nothing + ) + el = EventList(times, energies, meta) + + @test plot(el, 1.0, show_gtis=true, show_btis=true, tstart=0.5, tstop=8.5) isa Plots.Plot +end + +# GTI boundary edge cases +let + times = [2.0, 3.0, 4.0] + energies = [10.0, 20.0, 30.0] + gti_matrix = [1.0 2.5; 3.5 5.0] + meta = FITSMetadata{Dict{String,Any}}( + "test.fits", 2, "ENERGY", Dict{String,Vector}(), Dict{String,Any}(), + gti_matrix, "GTI",nothing, nothing, nothing, nothing, nothing, nothing + ) + el = EventList(times, energies, meta) + + @test plot(el, 1.0, show_gtis=true, tstart=2.0, tstop=4.0) isa Plots.Plot + @test plot(el, 1.0, show_btis=true, tstart=2.0, tstop=4.0) isa Plots.Plot +end +# Basic rebinning functionality test +let + times = collect(1.0:0.1:10.0) + counts = rand(1:10, length(times)) + + # Create mock light curve + metadata = LightCurveMetadata( + "TEST", "TEST", "TEST", 0.0, (1.0, 10.0), 0.1, + [Dict{String,Any}()], Dict{String,Any}() + ) + lc = LightCurve(times, 0.1, counts, nothing, nothing, + EventProperty{Float64}[], metadata, :poisson) + + @test plot(lc, 1.0) isa Plots.Plot + @test plot(lc, 0.5) isa Plots.Plot + @test plot(lc, 2.0) isa Plots.Plot +end