Implement File writing

Author: Aman-Pandey-afk

src/gti.jl

@@ -24,8 +24,13 @@ function get_gti_from_hdu(gtihdu::TableHDU)
     gtistart = FITSIO.read(gtihdu,startstr)
     gtistop = FITSIO.read(gtihdu,stopstr)
 
-    return mapreduce(permutedims, vcat, 
-    [[a, b] for (a,b) in zip(gtistart, gtistop)])
+    final_gti = [[a, b] for (a,b) in zip(gtistart, gtistop)]
+
+    if isempty(final_gti)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, final_gti)
+    end
 end
 
 function check_gtis(gti::AbstractMatrix)
@@ -96,7 +101,13 @@ function create_gti_mask(times::AbstractVector{<:Real},gtis::AbstractMatrix{<:Re
         end
     end
 
-    return mask, mapreduce(permutedims, vcat, keepat!(new_gtis,new_gti_mask))
+    keepat!(new_gtis,new_gti_mask)
+
+    if isempty(new_gtis)
+        return mask, reshape(Float64[],0,2)
+    else
+        return mask, mapreduce(permutedims, vcat, new_gtis)
+    end
 end
 
 function create_gti_from_condition(time::AbstractVector{<:Real}, condition::AbstractVector{Bool};
@@ -124,7 +135,11 @@ function create_gti_from_condition(time::AbstractVector{<:Real}, condition::Abst
         end
         push!(gtis,[t0, t1])
     end
-    return mapreduce(permutedims, vcat, gtis)
+    if isempty(gtis)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, gtis)
+    end
 end
 
 function operations_on_gtis(gti_list::AbstractVector{<:AbstractMatrix{T}}, 
@@ -189,7 +204,11 @@ function get_btis(gtis::AbstractMatrix{T}, start_time, stop_time) where {T<:Real
         push!(btis, [first(interval), last(interval)])
     end
 
-    return mapreduce(permutedims, vcat, btis)
+    if isempty(btis)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, btis)
+    end
 end
 
 function time_intervals_from_gtis(gtis::AbstractMatrix{<:Real}, segment_size::Real;

src/io.jl

@@ -19,32 +19,49 @@ function load_events_from_fits(filename::String)
         end
         ev.gti = load_gtis(filename)
         header = FITSIO.read_header(eventHDU)
-        if "MJDREFI" in header.keys
+        if haskey(header,"MJDREFI") 
             ev.mjdref += header["MJDREFI"]
         end
-        if "MJDREFF" in header.keys
+        if haskey(header,"MJDREFF") 
             ev.mjdref += header["MJDREFF"]
         end
-        if "INSTRUME" in header.keys
+        if haskey(header,"INSTRUME") 
             ev.instr = header["INSTRUME"]
         end
-        if "TIMEREF" in header.keys
+        if haskey(header,"TIMEREF") 
             ev.timeref = header["TIMEREF"]
         end
-        if "TIMESYS" in header.keys
+        if haskey(header,"TIMESYS") 
             ev.timesys = header["TIMESYS"]
         end
-        if "PLEPHEM" in header.keys
+        if haskey(header,"PLEPHEM") 
             ev.ephem = header["PLEPHEM"]
         end
         ev
     end
     return a
 end
 
-# TODO
 function write_events_to_fits(filename::String, ev::EventList)
-    # FITS(filename,"w") do event
-    #      FITSIO.write(event, ev)
-    # end
+    FITS(filename,"w") do hduList
+        
+        tkeys = String[]
+        values = []
+
+        for field in fieldnames(EventList)
+            fval = getfield(ev, field)
+            if !(fval isa AbstractVecOrMat)
+                push!(tkeys, String(field))
+                push!(values, fval)
+            end
+        end
+
+        header = FITSHeader(tkeys, values,fill("",length(tkeys)))
+
+        eventData = Dict("TIME"=>ev.time, "ENERGY"=>ev.energy, "PI"=>ev.PI);
+        FITSIO.write(hduList, eventData, header = header, name = "EVENTS")
+
+        gtiData = Dict("START"=>ev.gti[:,1], "STOP"=>ev.gti[:,2])
+        FITSIO.write(hduList, gtiData, name = "GTI")
+    end
 end

src/lightcurve.jl

@@ -19,7 +19,7 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0
     sort(toa)
     if isnothing(tstart)
         # if tstart is not set, assume light curve starts with first photon
-        tstart = toa[0]
+        tstart = toa[begin]
     end
 
     # compute the number of bins in the light curve
@@ -28,7 +28,7 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0
     # are not throwing away good events.
 
     if isnothing(tseg)
-        tseg = toa[-1] - tstart
+        tseg = toa[end-1] - tstart
     end
 
     timebin = tseg ÷ dt
@@ -38,15 +38,16 @@ function make_lightcurves(toa::AbstractVector{<:Real}, dt::Float64; tseg::Real=0
     end
 
     tend = tstart + timebin * dt
-    good = (tstart <= toa) & (toa < tend)
+    good = tstart <= toa < tend
     if !use_hist
-        binned_toas = (toa[good] - tstart) ÷ dt
+        binned_toas = (keepat!(toa,good) .- tstart) .÷ dt
         counts = bincount(binned_toas, minlength=timebin)
         time = tstart + range(0.5, 0.5 + len(counts)) * dt
     else
         histbins = range(tstart, tend + dt, dt)
-        counts, histbins = histogram(toa[good], bins=histbins)
-        time = histbins[begin:end] + 0.5 * dt
+        counts, histbins = histogram(keepat!(toa,good), bins=histbins)
+        time = @view(histbins[begin:end-1]) + 0.5 * dt
+    end
 
     return Lightcurve(time, counts; gti=gti, mjdref=mjdref, dt=dt,
                       skip_checks=True, err_dist="poisson")

test/runtests.jl

@@ -2,6 +2,6 @@ using Stingray
 using Test
 using FFTW, Distributions, Statistics, StatsBase, Metadata, HDF5
 
-# include("test_fourier.jl")
-# include("test_gti.jl")
+include("test_fourier.jl")
+include("test_gti.jl")
 include("test_events.jl")

test/test_events.jl

@@ -1,5 +1,5 @@
 @testset "test_events" begin
-    times = [0.5, 1.5, 2.5, 3.5]
+    time = [0.5, 1.5, 2.5, 3.5]
     counts = [3000, 2000, 2200, 3600]
     counts_flat = [3000, 3000, 3000, 3000]
     spectrum = [[1, 2, 3, 4, 5, 6], [1000, 2040, 1000, 3000, 4020, 2070]]
@@ -96,6 +96,14 @@
             ev = Stingray.read(fname, "fits")
             @test ev.mjdref == 55197.00076601852
         end
+
+        @testset "test_io_with_fits" begin
+            ev = EventList(time=time, mjdref=54000)
+            Stingray.write(ev, "ev.fits", "fits")
+            new_ev = Stingray.read("ev.fits", "fits")
+            @test new_ev.time == time
+            rm("ev.fits")
+        end
     end
 end