diff --git a/src/Stingray.jl b/src/Stingray.jl index 533cecd..5765474 100644 --- a/src/Stingray.jl +++ b/src/Stingray.jl @@ -77,4 +77,13 @@ export create_filtered_lightcurve export check_gtis export split_by_gtis +include("missionSupport.jl") +export MissionSupport, + get_mission_support, + apply_calibration, + patch_mission_info, + SIMPLE_CALIBRATION_FUNCS, + interpret_fits_data!, + AbstractMissionSupport + end diff --git a/src/events.jl b/src/events.jl index 143f897..190702d 100644 --- a/src/events.jl +++ b/src/events.jl @@ -1031,6 +1031,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, @@ -1045,6 +1048,16 @@ function readevents( 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 @@ -1072,6 +1085,12 @@ 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)) @@ -1083,6 +1102,14 @@ function readevents( 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 @@ -1101,6 +1128,23 @@ function readevents( (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 + + # 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 if !isnothing(mjd_ref) effective_timezero = isnothing(time_zero) ? 0.0 : time_zero @@ -1206,14 +1250,8 @@ function Base.summary(ev::EventList) if has_energies(ev) energy_range = extrema(ev.energies) summary_str *= ", energies: $(energy_range[1]) - $(energy_range[2])" - # Map energy units to display units if !isnothing(ev.meta.energy_units) - display_units = if ev.meta.energy_units in ["ENERGY", "PI", "PHA"] - "keV" - else - ev.meta.energy_units - end - summary_str *= " ($display_units)" + summary_str *= " ($(ev.meta.energy_units))" end end if has_gti(ev) diff --git a/src/gti.jl b/src/gti.jl index 5d74538..d426ead 100644 --- a/src/gti.jl +++ b/src/gti.jl @@ -97,13 +97,13 @@ function check_gtis(gti::AbstractMatrix) gti_start = @view gti[:, 1] gti_end = @view gti[:, 2] - if any(gti_end < gti_start) + if any(gti_end .< gti_start) 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])) + if any(@view(gti_start[begin+1:end]) .< @view(gti_end[begin:end-1])) throw(ArgumentError( "This GTI has overlaps" )) diff --git a/src/missionSupport.jl b/src/missionSupport.jl new file mode 100644 index 0000000..9ac0089 --- /dev/null +++ b/src/missionSupport.jl @@ -0,0 +1,226 @@ +""" +Dictionary of simple conversion functions for different missions. + +This dictionary provides PI (Pulse Invariant) to energy conversion functions +for various X-ray astronomy missions. Each function takes a PI channel value +and returns the corresponding energy in keV. + +Supported missions: +- NuSTAR: Nuclear Spectroscopic Telescope Array +- XMM: X-ray Multi-Mirror Mission +- NICER: Neutron star Interior Composition Explorer +- IXPE: Imaging X-ray Polarimetry Explorer +- AXAF/Chandra: Advanced X-ray Astrophysics Facility +- XTE/RXTE: Rossi X-ray Timing Explorer +""" +const SIMPLE_CALIBRATION_FUNCS = Dict{String, Function}( + "nustar" => (pi) -> pi * 0.04 + 1.62, + "xmm" => (pi) -> pi * 0.001, + "nicer" => (pi) -> pi * 0.01, + "ixpe" => (pi) -> pi / 375 * 15, + "axaf" => (pi) -> (pi - 1) * 14.6e-3, # Chandra/AXAF + "chandra" => (pi) -> (pi - 1) * 14.6e-3, # Explicit chandra entry + "xte" => (pi) -> pi * 0.025 # RXTE/XTE +) + +""" +Abstract type for mission-specific calibration and interpretation. + +This serves as the base type for all mission support implementations, +allowing for extensibility and type safety in mission-specific operations. +""" +abstract type AbstractMissionSupport end + +""" + MissionSupport{T} <: AbstractMissionSupport + +Structure containing mission-specific calibration and interpretation information. + +This structure encapsulates all the necessary information for handling +data from a specific X-ray astronomy mission, including calibration +functions, energy column alternatives, and GTI extension preferences. + +# Fields +- `name::String`: Mission name (normalized to lowercase) +- `instrument::Union{String, Nothing}`: Instrument identifier +- `epoch::Union{T, Nothing}`: Observation epoch in MJD (for time-dependent calibrations) +- `calibration_func::Function`: PI to energy conversion function +- `interpretation_func::Union{Function, Nothing}`: Mission-specific FITS interpretation function +- `energy_alternatives::Vector{String}`: Preferred energy column names in order of preference +- `gti_extensions::Vector{String}`: GTI extension names in order of preference + +# Type Parameters +- `T`: Type of the epoch parameter (typically Float64) +""" +struct MissionSupport{T} <: AbstractMissionSupport + name::String + instrument::Union{String, Nothing} + epoch::Union{T, Nothing} + calibration_func::Function + interpretation_func::Union{Function, Nothing} + energy_alternatives::Vector{String} + gti_extensions::Vector{String} +end + +""" + get_mission_support(mission::String, instrument=nothing, epoch=nothing) -> MissionSupport + +Create mission support object with mission-specific parameters. + +This function creates a MissionSupport object containing all the necessary +information for processing data from a specified X-ray astronomy mission. +It handles mission aliases (e.g., Chandra/AXAF) and provides appropriate +defaults for each mission. + +# Arguments +- `mission::String`: Mission name (case-insensitive) +- `instrument::Union{String, Nothing}=nothing`: Instrument identifier +- `epoch::Union{Float64, Nothing}=nothing`: Observation epoch in MJD + +# Returns +- `MissionSupport{Float64}`: Mission support object + +# Throws +- `ArgumentError`: If mission name is empty + +# Examples +```julia +# Basic usage +ms = get_mission_support("nustar") + +# With instrument specification +ms = get_mission_support("nustar", "FPM_A") + +# With epoch for time-dependent calibrations +ms = get_mission_support("xte", "PCA", 50000.0) +``` +""" +function get_mission_support(mission::String, + instrument::Union{String, Nothing}=nothing, + epoch::Union{Float64, Nothing}=nothing) + + # Check for empty mission string + if isempty(mission) + throw(ArgumentError("Mission name cannot be empty")) + end + + mission_lower = lowercase(mission) + + if mission_lower in ["chandra", "axaf"] + mission_lower = "chandra" + end + + calib_func = if haskey(SIMPLE_CALIBRATION_FUNCS, mission_lower) + SIMPLE_CALIBRATION_FUNCS[mission_lower] + else + @warn "Mission $mission not recognized, using identity function" + identity + end + + energy_alts = if mission_lower in ["chandra", "axaf"] + ["ENERGY", "PI", "PHA"] + elseif mission_lower == "xte" + ["PHA", "PI", "ENERGY"] + elseif mission_lower == "nustar" + ["PI", "ENERGY", "PHA"] + else + ["ENERGY", "PI", "PHA"] + end + + # Mission-specific GTI extensions + gti_exts = if mission_lower == "xmm" + ["GTI", "GTI0", "STDGTI"] + elseif mission_lower in ["chandra", "axaf"] + ["GTI", "GTI0", "GTI1", "GTI2", "GTI3"] + else + ["GTI", "STDGTI"] + end + + MissionSupport{Float64}(mission_lower, instrument, epoch, calib_func, nothing, energy_alts, gti_exts) +end + +""" + apply_calibration(mission_support::MissionSupport, pi_channels::AbstractArray) -> Vector{Float64} + +Apply calibration function to PI channels. + +Converts PI (Pulse Invariant) channel values to energies in keV using +the mission-specific calibration function stored in the MissionSupport object. + +# Arguments +- `mission_support::MissionSupport`: Mission support object containing calibration function +- `pi_channels::AbstractArray{T}`: Array of PI channel values + +# Returns +- `Vector{Float64}`: Array of energy values in keV + +# Examples +```julia +ms = get_mission_support("nustar") +pi_values = [100, 500, 1000] +energies = apply_calibration(ms, pi_values) +``` +""" +function apply_calibration(mission_support::MissionSupport, pi_channels::AbstractArray{T}) where T + if isempty(pi_channels) + return similar(pi_channels, Float64) + end + return mission_support.calibration_func.(pi_channels) +end + +""" + patch_mission_info(info::Dict{String,Any}, mission=nothing) -> Dict{String,Any} + +Apply mission-specific patches to header information. + +This function applies mission-specific modifications to FITS header information +to handle mission-specific quirks and conventions. It's based on the Python +implementation in Stingray's mission interpretation module. + +# Arguments +- `info::Dict{String,Any}`: Dictionary containing header information +- `mission::Union{String,Nothing}=nothing`: Mission name + +# Returns +- `Dict{String,Any}`: Patched header information dictionary + +# Examples +```julia +info = Dict("gti" => "STDGTI", "ecol" => "PHA") +patched = patch_mission_info(info, "xmm") # Adds GTI0 to gti field +``` +""" +function patch_mission_info(info::Dict{String,Any}, mission::Union{String,Nothing}=nothing) + if isnothing(mission) + return info + end + + mission_lower = lowercase(mission) + patched_info = copy(info) + + # Normalize chandra/axaf + if mission_lower in ["chandra", "axaf"] + mission_lower = "chandra" + end + + if mission_lower == "xmm" && haskey(patched_info, "gti") + patched_info["gti"] = string(patched_info["gti"], ",GTI0") + elseif mission_lower == "xte" && haskey(patched_info, "ecol") + patched_info["ecol"] = "PHA" + patched_info["ccol"] = "PCUID" + elseif mission_lower == "chandra" + # Chandra-specific patches + if haskey(patched_info, "DETNAM") + patched_info["detector"] = patched_info["DETNAM"] + end + if haskey(patched_info, "TIMESYS") + patched_info["time_system"] = patched_info["TIMESYS"] + end + end + + return patched_info +end +function interpret_fits_data!(f::FITS, mission_support::MissionSupport) + # Placeholder for mission-specific interpretation + return nothing +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 12ca46a..2882333 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,4 +15,8 @@ end @testset "lightcurve" begin include("test_lightcurve.jl") +end + +@testset "Synthetic Events Tests" begin + include("test_missionSupport.jl") end \ No newline at end of file diff --git a/test/test_events.jl b/test/test_events.jl index cbf10ce..f8fd4a8 100644 --- a/test/test_events.jl +++ b/test/test_events.jl @@ -669,7 +669,7 @@ let @test occursin("4 events", summary_str) @test occursin("over 3.0 time units", summary_str) @test occursin("energies: 100.0 - 400.0", summary_str) - @test occursin("(keV)", summary_str) + @test occursin("(ENERGY)", summary_str) @test occursin("GTI: 2 intervals (3.0 s exposure)", summary_str) @test occursin("1 extra columns", summary_str) diff --git a/test/test_missionSupport.jl b/test/test_missionSupport.jl new file mode 100644 index 0000000..4506324 --- /dev/null +++ b/test/test_missionSupport.jl @@ -0,0 +1,580 @@ +""" + create_synthetic_events(n_events::Int=1000; mission::String="nustar", + time_range::Tuple{Float64,Float64}=(0.0, 1000.0), + pi_range::Tuple{Int,Int}=(1, 4096), + T::Type=Float64) -> EventList{Vector{T}, FITSMetadata} + +Create synthetic X-ray event data for testing purposes. + +This function generates synthetic X-ray event data with proper time ordering, +energy calibration, and mission-specific metadata. The generated events include: +- Random event times within the specified range +- Random PI (Pulse Invariant) channels +- Energy values calculated using mission-specific calibration +- Detector ID assignments +- Realistic FITS-style metadata headers + +# Arguments +- `n_events::Int=1000`: Number of events to generate +- `mission::String="nustar"`: Mission name for calibration and metadata +- `time_range::Tuple{Float64,Float64}=(0.0, 1000.0)`: Time range for events (start, stop) +- `pi_range::Tuple{Int,Int}=(1, 4096)`: PI channel range +- `T::Type=Float64`: Numeric type for the data + +# Returns +- `EventList{Vector{T}, FITSMetadata}`: Synthetic event list with times, energies, and metadata + +# Examples +```julia +# Default NuSTAR events +events = create_synthetic_events() + +# Custom XMM events +events = create_synthetic_events(5000; mission="xmm", time_range=(100.0, 200.0)) + +# NICER events with custom PI range +events = create_synthetic_events(2000; mission="nicer", pi_range=(50, 1000)) +``` +""" +function create_synthetic_events(n_events::Int=1000; + mission::String="nustar", + time_range::Tuple{Float64,Float64}=(0.0, 1000.0), + pi_range::Tuple{Int,Int}=(1, 4096), + T::Type=Float64) + + # Input validation + if n_events <= 0 + throw(ArgumentError("Number of events must be positive")) + end + + if time_range[1] >= time_range[2] + throw(ArgumentError("Time range start must be less than end")) + end + + if pi_range[1] > pi_range[2] + throw(ArgumentError("PI range start must be less than or equal to end")) + end + + # Get mission support for calibration + mission_support = get_mission_support(mission) + + # Generate random event times and sort them + times = sort(rand(n_events) * (time_range[2] - time_range[1]) .+ time_range[1]) + times = convert(Vector{T}, times) + + # Generate random PI channels + pi_channels = rand(pi_range[1]:pi_range[2], n_events) + + # Apply mission-specific calibration to get energies + energies = convert(Vector{T}, apply_calibration(mission_support, pi_channels)) + + # Generate detector IDs (assume 4-detector system like NuSTAR) + detector_ids = rand(0:3, n_events) + + # Create extra columns dictionary + extra_columns = Dict{String, Vector}( + "PI" => pi_channels, + "DETID" => detector_ids + ) + + # Create synthetic filename + filepath = "synthetic_$(mission).fits" + + # Create FITS headers using FITSIO.FITSHeader structure + # Note: This is a simplified header for testing - in real usage, + # headers would come from actual FITS files + header_dict = Dict{String, Any}( + "TELESCOP" => uppercase(mission), + "INSTRUME" => "TEST_INSTRUMENT", + "OBJECT" => "TEST_SOURCE", + "RA_NOM" => 180.0, + "DEC_NOM" => 0.0, + "EQUINOX" => 2000.0, + "RADECSYS" => "FK5", + "TSTART" => time_range[1], + "TSTOP" => time_range[2], + "EXPOSURE" => time_range[2] - time_range[1], + "ONTIME" => time_range[2] - time_range[1], + "LIVETIME" => time_range[2] - time_range[1], + "NAXIS2" => n_events, + "TIMESYS" => "TT", + "TIMEREF" => "LOCAL", + "TIMEUNIT" => "s", + "MJDREFI" => 55197, # Standard MJD reference for many missions + "MJDREFF" => 0.00076601852, + "CLOCKCORR" => "T", + "DATE-OBS" => "2020-01-01T00:00:00", + "DATE-END" => "2020-01-01T01:00:00" + ) + + # Add mission-specific header information + if lowercase(mission) == "nustar" + header_dict["DETNAM"] = "TEST_DET" + header_dict["PIFLTCOR"] = "T" + elseif lowercase(mission) == "xmm" + header_dict["FILTER"] = "Medium" + header_dict["SUBMODE"] = "PrimeFullWindow" + elseif lowercase(mission) == "nicer" + header_dict["DETNAM"] = "TEST_NICER" + header_dict["FILTFILE"] = "NONE" + elseif lowercase(mission) in ["chandra", "axaf"] + header_dict["DETNAM"] = "ACIS-S" + header_dict["GRATING"] = "NONE" + elseif lowercase(mission) == "xte" + header_dict["DETNAM"] = "PCA" + header_dict["LAYERS"] = "ALL" + elseif lowercase(mission) == "ixpe" + header_dict["DETNAM"] = "GPD" + header_dict["POLMODE"] = "ON" + end + + # Create a simple header structure (for testing purposes) + # In real usage, this would be a proper FITSIO.FITSHeader + headers = header_dict + + # Create FITSMetadata + metadata = FITSMetadata( + filepath, # filepath + 2, # hdu (typical for event data) + "ENERGY", # energy_units (after calibration) + extra_columns, # extra_columns + headers, # headers + nothing, # gti + nothing, # gti_source + nothing, # mjd_ref + nothing, # time_zero + nothing, # time_unit + nothing, # time_sys + nothing, # time_pixr + nothing # time_del + ) + + # Create and return EventList + return EventList(times, energies, metadata) +end +# Test: Basic Synthetic Event Creation - Default Parameters +let + events = create_synthetic_events() + + # Test basic structure + @test isa(events, EventList) + @test length(events) == 1000 # Default n_events + @test length(times(events)) == 1000 + @test length(energies(events)) == 1000 + @test events.meta.filepath == "synthetic_nustar.fits" + + # Test times are sorted and in range + @test issorted(times(events)) + @test all(0.0 .<= times(events) .<= 1000.0) + + # Test energies are positive (after calibration) + @test all(energies(events) .> 0) + + # Test extra columns + @test haskey(events.meta.extra_columns, "PI") + @test haskey(events.meta.extra_columns, "DETID") + @test length(events.meta.extra_columns["PI"]) == 1000 + @test length(events.meta.extra_columns["DETID"]) == 1000 + + # Test PI channels are in expected range + @test all(1 .<= events.meta.extra_columns["PI"] .<= 4096) + + # Test detector IDs are in expected range + @test all(0 .<= events.meta.extra_columns["DETID"] .<= 3) + + # Test metadata structure + @test isa(events.meta, FITSMetadata) + @test events.meta.hdu == 2 + @test events.meta.energy_units == "ENERGY" + + # Test metadata content + @test events.meta.headers["TELESCOP"] == "NUSTAR" + @test events.meta.headers["INSTRUME"] == "TEST_INSTRUMENT" + @test events.meta.headers["NAXIS2"] == 1000 + @test events.meta.headers["TSTART"] == 0.0 + @test events.meta.headers["TSTOP"] == 1000.0 +end + +# Test: Basic Synthetic Event Creation - Custom Parameters +let + n_events = 500 + mission = "xmm" + time_range = (100.0, 200.0) + pi_range = (50, 1000) + + events = create_synthetic_events(n_events; + mission=mission, + time_range=time_range, + pi_range=pi_range) + + # Test custom parameters are respected + @test length(events) == n_events + @test length(times(events)) == n_events + @test length(energies(events)) == n_events + @test events.meta.filepath == "synthetic_xmm.fits" + + # Test time range + @test all(time_range[1] .<= times(events) .<= time_range[2]) + @test events.meta.headers["TSTART"] == time_range[1] + @test events.meta.headers["TSTOP"] == time_range[2] + + # Test PI range + @test all(pi_range[1] .<= events.meta.extra_columns["PI"] .<= pi_range[2]) + + # Test mission-specific metadata + @test events.meta.headers["TELESCOP"] == "XMM" +end + +# Test: Different Missions +let + missions = ["nustar", "xmm", "nicer", "ixpe", "axaf", "chandra", "xte"] + + for mission in missions + events = create_synthetic_events(100; mission=mission) + + @test events.meta.filepath == "synthetic_$(mission).fits" + @test events.meta.headers["TELESCOP"] == uppercase(mission) + @test length(events) == 100 + @test length(times(events)) == 100 + @test length(energies(events)) == 100 + + # Test that calibration was applied correctly + ms = get_mission_support(mission) + expected_energies = apply_calibration(ms, events.meta.extra_columns["PI"]) + @test energies(events) ≈ expected_energies + end +end + +# Test: Time Ordering +let + for _ in 1:10 # Test multiple times due to randomness + events = create_synthetic_events(100) + @test issorted(times(events)) + + # Test no duplicate times (very unlikely but possible) + @test length(unique(times(events))) >= 95 # Allow some duplicates due to floating point + end +end + +# Test: Energy Calibration Consistency +let + missions = ["nustar", "xmm", "nicer", "ixpe"] + + for mission in missions + events = create_synthetic_events(200; mission=mission) + + # Manually verify calibration + ms = get_mission_support(mission) + expected_energies = apply_calibration(ms, events.meta.extra_columns["PI"]) + @test energies(events) ≈ expected_energies + + # Test energy ranges are reasonable for each mission + if mission == "nustar" + # NuSTAR: pi * 0.04 + 1.62, PI range 1-4096 + @test minimum(energies(events)) >= 1.66 # 1*0.04 + 1.62 + @test maximum(energies(events)) <= 165.46 # 4096*0.04 + 1.62 + elseif mission == "xmm" + # XMM: pi * 0.001, PI range 1-4096 + @test minimum(energies(events)) >= 0.001 + @test maximum(energies(events)) <= 4.096 + elseif mission == "nicer" + # NICER: pi * 0.01, PI range 1-4096 + @test minimum(energies(events)) >= 0.01 + @test maximum(energies(events)) <= 40.96 + elseif mission == "ixpe" + # IXPE: pi / 375 * 15, PI range 1-4096 + @test minimum(energies(events)) >= 15.0/375 # ≈ 0.04 + @test maximum(energies(events)) <= 4096*15.0/375 # ≈ 163.84 + end + end +end + +# Test: Statistical Properties +let + events = create_synthetic_events(10000) # Large sample for statistics + + # Test time distribution (should be roughly uniform) + time_hist = fit(Histogram, times(events), 0:100:1000) + counts = time_hist.weights + # Expect roughly equal counts in each bin (within statistical fluctuation) + expected_count = 10000 / 10 # 10 bins + @test all(abs.(counts .- expected_count) .< 3 * sqrt(expected_count)) # 3-sigma test + + # Test PI distribution (should be roughly uniform over discrete range) + pi_hist = fit(Histogram, events.meta.extra_columns["PI"], 1:100:4096) + pi_counts = pi_hist.weights + # More lenient test due to discrete uniform distribution + @test std(pi_counts) / mean(pi_counts) < 0.2 # Coefficient of variation < 20% + + # Test detector ID distribution + detid_counts = [count(==(i), events.meta.extra_columns["DETID"]) for i in 0:3] + @test all(abs.(detid_counts .- 2500) .< 3 * sqrt(2500)) # Each detector ~2500 events +end + +# Test: Small Event Counts +let + for n in [1, 2, 5, 10] + events = create_synthetic_events(n) + @test length(events) == n + @test length(times(events)) == n + @test length(energies(events)) == n + @test length(events.meta.extra_columns["PI"]) == n + @test length(events.meta.extra_columns["DETID"]) == n + @test issorted(times(events)) + end +end + +# Test: Large Event Counts +let + events = create_synthetic_events(100000) + @test length(events) == 100000 + @test length(times(events)) == 100000 + @test length(energies(events)) == 100000 + @test issorted(times(events)) + @test events.meta.headers["NAXIS2"] == 100000 +end + +# Test: Extreme Time Ranges +let + # Very short time range + events = create_synthetic_events(100; time_range=(0.0, 0.1)) + @test all(0.0 .<= times(events) .<= 0.1) + @test events.meta.headers["TSTART"] == 0.0 + @test events.meta.headers["TSTOP"] == 0.1 + + # Very long time range + events = create_synthetic_events(100; time_range=(0.0, 1e6)) + @test all(0.0 .<= times(events) .<= 1e6) + @test events.meta.headers["TSTART"] == 0.0 + @test events.meta.headers["TSTOP"] == 1e6 + + # Negative time range + events = create_synthetic_events(100; time_range=(-1000.0, -500.0)) + @test all(-1000.0 .<= times(events) .<= -500.0) + @test issorted(times(events)) +end + +# Test: Extreme PI Ranges +let + # Small PI range + events = create_synthetic_events(100; pi_range=(100, 110)) + @test all(100 .<= events.meta.extra_columns["PI"] .<= 110) + + # Single PI value + events = create_synthetic_events(100; pi_range=(500, 500)) + @test all(events.meta.extra_columns["PI"] .== 500) + @test all(energies(events) .== energies(events)[1]) # All same energy + + # Large PI range + events = create_synthetic_events(100; pi_range=(1, 10000)) + @test all(1 .<= events.meta.extra_columns["PI"] .<= 10000) +end + +# Test: Unknown Mission +let + # Should still work but with warning + @test_logs (:warn, r"Mission unknown_mission not recognized") begin + events = create_synthetic_events(100; mission="unknown_mission") + @test events.meta.filepath == "synthetic_unknown_mission.fits" + @test events.meta.headers["TELESCOP"] == "UNKNOWN_MISSION" + @test length(events) == 100 + # Should use identity calibration + @test energies(events) == Float64.(events.meta.extra_columns["PI"]) + end +end + +# Test: No Missing Data +let + events = create_synthetic_events(1000) + + # Check no NaN or missing values + @test all(isfinite.(times(events))) + @test all(isfinite.(energies(events))) + @test all(isfinite.(events.meta.extra_columns["PI"])) + @test all(isfinite.(events.meta.extra_columns["DETID"])) + + # Check no negative energies (after calibration) + @test all(energies(events) .>= 0) +end + +# Test: Correct Data Types +let + events = create_synthetic_events(100) + + @test eltype(times(events)) == Float64 + @test eltype(energies(events)) == Float64 + @test eltype(events.meta.extra_columns["PI"]) <: Integer + @test eltype(events.meta.extra_columns["DETID"]) <: Integer + @test isa(events.meta, FITSMetadata) + @test isa(events.meta.filepath, String) +end + +# Test: Array Length Consistency +let + for n in [10, 100, 1000, 5000] + events = create_synthetic_events(n) + + @test length(events) == n + @test length(times(events)) == n + @test length(energies(events)) == n + @test length(events.meta.extra_columns["PI"]) == n + @test length(events.meta.extra_columns["DETID"]) == n + @test events.meta.headers["NAXIS2"] == n + end +end + +# Test: Mission Name Handling +let + # Test case sensitivity + missions = ["NUSTAR", "nustar", "NuSTAR", "NuStar"] + for mission in missions + events = create_synthetic_events(100; mission=mission) + @test events.meta.headers["TELESCOP"] == "NUSTAR" + @test events.meta.filepath == "synthetic_$(mission).fits" + end +end + +# Test: Calibration Differences +let + # Same PI values should give different energies for different missions + test_pi_range = (1000, 1000) # Fixed PI value + + results = Dict{String, Float64}() + for mission in ["nustar", "xmm", "nicer", "ixpe"] + events = create_synthetic_events(100; mission=mission, pi_range=test_pi_range) + results[mission] = energies(events)[1] # All should be same since PI is fixed + end + + # Different missions should give different energies + missions = collect(keys(results)) + for i in 1:length(missions) + for j in (i+1):length(missions) + @test results[missions[i]] ≠ results[missions[j]] + end + end +end + +# Test: Metadata Consistency +let + missions = ["nustar", "xmm", "nicer", "ixpe", "axaf", "chandra"] + + for mission in missions + events = create_synthetic_events(100; mission=mission) + + @test events.meta.headers["TELESCOP"] == uppercase(mission) + @test events.meta.headers["INSTRUME"] == "TEST_INSTRUMENT" + @test haskey(events.meta.headers, "NAXIS2") + @test haskey(events.meta.headers, "TSTART") + @test haskey(events.meta.headers, "TSTOP") + end +end + +# Test: EventList Interface Methods +let + events = create_synthetic_events(1000) + + # Test length methods + @test length(events) == 1000 + @test size(events) == (1000,) + + # Test accessor methods + @test times(events) === events.times + @test energies(events) === events.energies + @test has_energies(events) == true + + # Test show methods (should not error) + io = IOBuffer() + show(io, MIME"text/plain"(), events) + @test length(String(take!(io))) > 0 + + show(io, MIME"text/plain"(), events.meta) + @test length(String(take!(io))) > 0 +end + +# Test: Large Dataset Creation +let + # Test that large datasets can be created without issues + @time events = create_synthetic_events(50000) + @test length(events) == 50000 + @test sizeof(times(events)) + sizeof(energies(events)) < 1e6 # Less than 1MB for 50k events +end + +# Test: Memory Efficiency +let + # Test that no unnecessary copies are made + events1 = create_synthetic_events(1000) + events2 = create_synthetic_events(1000) + + # Each should have independent data + @test times(events1) !== times(events2) + @test energies(events1) !== energies(events2) + @test events1.meta.extra_columns["PI"] !== events2.meta.extra_columns["PI"] +end + +# Test: Random Seed Behavior +let + # Test that different calls produce different results (random behavior) + events1 = create_synthetic_events(1000) + events2 = create_synthetic_events(1000) + + # Should be different due to randomness + @test times(events1) != times(events2) + @test events1.meta.extra_columns["PI"] != events2.meta.extra_columns["PI"] + @test events1.meta.extra_columns["DETID"] != events2.meta.extra_columns["DETID"] + + # But same structure + @test length(events1) == length(events2) + @test typeof(events1) == typeof(events2) +end + +# Test: Input Validation - Error Cases +let + # Test negative event count + @test_throws ArgumentError create_synthetic_events(-1) + @test_throws ArgumentError create_synthetic_events(0) + + # Test invalid time range + @test_throws ArgumentError create_synthetic_events(100; time_range=(100.0, 50.0)) + @test_throws ArgumentError create_synthetic_events(100; time_range=(100.0, 100.0)) + + # Test invalid PI range + @test_throws ArgumentError create_synthetic_events(100; pi_range=(1000, 500)) +end + +# Test: FITSMetadata Structure +let + events = create_synthetic_events(100) + meta = events.meta + + # Test all required fields are present + @test isa(meta.filepath, String) + @test isa(meta.hdu, Int) + @test isa(meta.energy_units, Union{String, Nothing}) + @test isa(meta.extra_columns, Dict{String, Vector}) + @test !isnothing(meta.headers) + + # Test metadata consistency + @test meta.hdu == 2 + @test meta.energy_units == "ENERGY" + @test length(meta.extra_columns) == 2 # PI and DETID +end + +# Test: Summary Function +let + events = create_synthetic_events(1000) + summary_str = summary(events) + + @test isa(summary_str, String) + @test occursin("1000 events", summary_str) + # The time span should be close to but less than the full range (1000.0) + actual_time_span = maximum(times(events)) - minimum(times(events)) + expected_pattern = "$(actual_time_span) time units" + @test occursin(expected_pattern, summary_str) + + # Alternative: Use regex to match any floating point number + @test occursin(r"\d+\.\d+ time units", summary_str) + + @test occursin("energies:", summary_str) + @test occursin("(ENERGY)", summary_str) + @test occursin("2 extra columns", summary_str) +end \ No newline at end of file