diff --git a/Project.toml b/Project.toml
index cd910a8..de20406 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,6 +16,7 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 
 [compat]
 julia = "1.7"
@@ -30,6 +31,7 @@ HDF5 = "0.16"
 FITSIO = "0.16"
 Intervals = "1.8"
 NaNMath = "0.3, 1"
+Parameters = "0.12"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/Stingray.jl b/src/Stingray.jl
index 740e849..414a981 100644
--- a/src/Stingray.jl
+++ b/src/Stingray.jl
@@ -1,7 +1,7 @@
 module Stingray
 
 using ResumableFunctions, StatsBase, Statistics, DataFrames
-using FFTW, Metadata, NaNMath, FITSIO, Intervals
+using FFTW, Metadata, NaNMath, FITSIO, Intervals, Parameters
 using ProgressBars: tqdm as show_progress 
 
 include("fourier.jl")
@@ -30,6 +30,18 @@ export get_btis
 export time_intervals_from_gtis
 export bin_intervals_from_gtis
 
+include("lightcurve.jl")
+export LightCurve
+export rebin
+
+include("events.jl")
+export EventList
+export to_lc
+export from_lc
+
+include("io.jl")
+export load_events_from_fits
+
 include("utils.jl")
 
 end 
diff --git a/src/events.jl b/src/events.jl
new file mode 100644
index 0000000..6624a2a
--- /dev/null
+++ b/src/events.jl
@@ -0,0 +1,150 @@
+import Base: read, write, join, sort, sort!
+
+@with_kw mutable struct EventList{T1<:AbstractVector,T2<:AbstractVector,T3<:AbstractMatrix}
+    time::T1
+    energy::T2=Float64[]
+    ncounts::Int=0
+    mjdref::Float64=0
+    dt::Float64=0
+    notes::String=""
+    gti::T3=reshape(Float64[],0,2)
+    PI::Vector{Int}=Int[]
+    high_precision::Bool=false
+    mission::String=""
+    instr::String=""
+    header::String=""
+    detector_id::Vector{String}=String[]
+    ephem::String=""
+    timeref::String=""
+    timesys::String=""
+end
+
+function read(::Type{EventList},filename::String, format::String)
+    if format=="fits"
+        return load_events_from_fits(filename)
+    else
+        throw(ArgumentError("File format $(format) not supported still."))
+    end
+end
+
+function write(ev::EventList, filename::String, format::String)
+    if format=="fits"
+        write_events_to_fits(filename, ev)
+    else
+        throw(ArgumentError("File format $(format) not supported still."))
+    end
+end
+
+function from_lc(lc::LightCurve)
+    times = [lc.time[i] for i in 1:length(lc.time) for _ in 1:lc.counts[i]]
+    return EventList(time=times,gti=lc.gti)
+end
+
+function to_lc(ev::EventList, dt)
+    if isempty(ev.gti)
+        tstart=0
+        tseg=0
+    else
+        tstart = ev.gti[begin][begin]
+        tseg = ev.gti[end][end]-tstart
+    end
+    return to_lc(ev, dt, tstart, tseg)
+end
+
+to_lc(ev::EventList, dt::Real,
+      tstart::Real, tseg::Real) = make_lightcurves(ev.time, dt, ev.gti;
+                                                   tstart=tstart, tseg=tseg,
+                                                   mjdref=ev.mjdref)
+
+function join(ev1::EventList, ev2::EventList)
+    new_ev = EventList(time=vcat(ev1.time, ev2.time))
+    if ev1.dt!=ev2.dt
+        new_ev.dt = max(ev1.dt,ev2.dt);
+    end
+
+    # Tolerance for MJDREF:1 microsecond
+    if !isapprox(ev1.mjdref, ev2.mjdref, atol=1e-6 / 86400)
+        ev2.mjdref = ev1.mjdref
+    end
+
+    order = sortperm(new_ev.time)
+    new_ev.time = new_ev.time[order]
+
+    if !isempty(ev1.PI) || !isempty(ev2.PI)
+        if isempty(ev1.PI)
+            ev1.PI = zero(ev1.time)
+        end
+        if isempty(ev2.PI)
+            ev2.PI = zero(ev2.time)
+        end
+        new_ev.PI = vcat(ev1.PI, ev2.PI)
+        new_ev.PI = new_ev.PI[order]
+    else
+        new_ev.PI = Float64[]
+    end
+
+    if !isempty(ev1.energy) || !isempty(ev2.energy)
+        if isempty(ev1.energy)
+            ev1.energy = zero(ev1.time)
+        end
+        if isempty(ev2.energy)
+            ev2.energy = zero(ev2.time)
+        end
+        new_ev.energy = vcat(ev1.energy, ev2.energy)
+        new_ev.energy = new_ev.energy[order]
+    else
+        new_ev.energy = Float64[]
+    end
+
+    if !isempty(ev1.gti) || !isempty(ev2.gti)
+        if isempty(ev1.gti)
+            ev1.gti = [ev1.time[begin] - ev1.dt/2 ev1.time[end] + ev1.dt/2]
+        end
+        if isempty(ev2.gti)
+            ev2.gti = [ev2.time[begin] - ev2.dt/2 ev2.time[end] + ev2.dt/2]
+        end
+        new_ev.gti = operations_on_gtis([ev1.gti,ev2.gti], intersect)
+        if isempty(new_ev.gti)
+            @warn """
+            GTIs in these two event lists do not overlap at all.
+            Merging instead of returning an overlap.
+            """
+            new_ev.gti = operations_on_gtis([ev1.gti,ev2.gti], union)
+        end
+    else
+        new_ev.gti = reshape(Float64[],0,2)
+    end
+
+    for attr in (:mission, :instr)
+        if getfield(ev1, attr) != getfield(ev2, attr)
+            setfield!(new_ev, attr, string(getfield(ev1, attr),',',getfield(ev2, attr)))
+        else
+            setfield!(new_ev, attr, getfield(ev1, attr))
+        end
+    end
+
+    new_ev.mjdref = ev1.mjdref
+
+    return new_ev
+end
+                                            
+function sort(ev::EventList)
+    order = sortperm(ev.time)
+    new_ev = deepcopy(ev)
+    apply_mask(new_ev, order)
+    return new_ev
+end
+
+function sort!(ev::EventList)
+    order = sortperm(ev.time)
+    apply_mask(ev, order)
+end
+
+function apply_mask(ev::EventList, mask::AbstractVector{T};) where {T<:Union{Bool,Int}}
+    for field in fieldnames(EventList)
+        fval = getfield(ev, field)
+        if fval isa AbstractVector && !isempty(fval) 
+            setfield!(ev, field, deepcopy(fval)[mask])
+        end
+    end
+end
diff --git a/src/gti.jl b/src/gti.jl
index 1d82065..6b380bd 100644
--- a/src/gti.jl
+++ b/src/gti.jl
@@ -21,11 +21,16 @@ function get_gti_from_hdu(gtihdu::TableHDU)
         stopstr = "Stop"
     end
 
-    gtistart = read(gtihdu,startstr)
-    gtistop = read(gtihdu,stopstr)
+    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}}, 
@@ -152,7 +167,11 @@ function operations_on_gtis(gti_list::AbstractVector{<:AbstractMatrix{T}},
         push!(final_gti, [first(interval), last(interval)])
     end
 
-    return mapreduce(permutedims, vcat, final_gti)
+    if isempty(final_gti)
+        return reshape(Float64[],0,2)
+    else
+        return mapreduce(permutedims, vcat, final_gti)
+    end
 end
 
 function get_btis(gtis::AbstractMatrix{<:Real})
@@ -185,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;
diff --git a/src/io.jl b/src/io.jl
new file mode 100644
index 0000000..57e2e57
--- /dev/null
+++ b/src/io.jl
@@ -0,0 +1,66 @@
+function load_events_from_fits(filename::String)
+    FITS(filename) do hduList
+        eventHDU = hduList["EVENTS"]
+        cols = FITSIO.colnames(eventHDU)
+        df = DataFrame(eventHDU)
+        if "TIME" in cols
+            time = df[!,"TIME"]
+        else
+            throw(KeyError("Time stamp data not provided or column names not appropriate in the file"))
+        end
+        ev = EventList(time=time)
+        if "PI" in cols
+            PI = df[!,"PI"]
+            ev.PI = PI
+        end
+        if "ENERGY" in cols
+            energy = df[!,"ENERGY"]
+            ev.energy = energy
+        end
+        ev.gti = load_gtis(filename)
+        header = FITSIO.read_header(eventHDU)
+        if haskey(header,"MJDREFI") 
+            ev.mjdref += header["MJDREFI"]
+        end
+        if haskey(header,"MJDREFF") 
+            ev.mjdref += header["MJDREFF"]
+        end
+        if haskey(header,"INSTRUME") 
+            ev.instr = header["INSTRUME"]
+        end
+        if haskey(header,"TIMEREF") 
+            ev.timeref = header["TIMEREF"]
+        end
+        if haskey(header,"TIMESYS") 
+            ev.timesys = header["TIMESYS"]
+        end
+        if haskey(header,"PLEPHEM") 
+            ev.ephem = header["PLEPHEM"]
+        end
+        ev
+    end
+end
+
+function write_events_to_fits(filename::String, ev::EventList)
+    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
diff --git a/src/lightcurve.jl b/src/lightcurve.jl
new file mode 100644
index 0000000..cc57cb1
--- /dev/null
+++ b/src/lightcurve.jl
@@ -0,0 +1,100 @@
+@with_kw mutable struct LightCurve{T1<:AbstractVector,T2<:AbstractMatrix}
+    time::T1
+    counts::Vector{Int}=Int[]
+    err::Vector{Float64}=Float64[]
+    input_counts::Bool=true
+    gti::T2=reshape(Float64[],0,2)
+    err_dist::String=""
+    mjdref::Float64=0
+    dt::Float64=0
+    skip_checks::Bool=false
+    low_memory::Bool=false
+    mission::String=""
+    instr::String=""
+    header::String=""
+end
+
+function make_lightcurves(toa::AbstractVector{<:Real}, dt::Real, 
+                          gti::AbstractMatrix{<:Real}; tseg::Real=-1, 
+                          tstart::Real=-1, mjdref::Real=0)
+    sort!(toa)
+    if tstart == -1
+        # if tstart is not set, assume light curve starts with first photon
+        tstart = toa[begin]
+    end
+
+    # compute the number of bins in the light curve
+    # for cases where tseg/dt is not integer.
+    # TODO: check that this is always consistent and that we
+    # are not throwing away good events.
+
+    if tseg == -1
+        tseg = toa[end-1] - tstart
+    end
+
+    timebin = tseg ÷ dt
+    # If we are missing the next bin by just 1%, let's round up:
+    if tseg / dt - timebin >= 0.99
+        timebin += 1
+    end
+    
+    tend = tstart + timebin * dt
+    good = tstart .<= toa .< tend
+    
+    histbins = range(tstart, tend, step = dt)
+    counts= fit(Histogram, keepat!(toa,good), histbins).weights
+    time = @view(histbins[begin:end-1]) .+ 0.5 .* dt
+
+    return LightCurve(time=time, counts=counts, gti=gti, mjdref=mjdref, dt=dt,
+                      skip_checks=true, err_dist="poisson")
+end
+
+function rebin(lc::LightCurve, dt_new::Real; method::String = "sum")
+    if dt_new<lc.dt
+        throw(ArgumentError("New time resolution must be larger than old time resolution!"))
+    end
+
+    bin_time = Float64[]
+    bin_counts = Int[]
+    bin_err = Float64[]
+    gti_new = Vector{Float64}[]
+
+    for g in eachrow(lc.gti)
+        if g[2] - g[1] < dt_new
+            continue
+        else
+            # find start and end of GTI segment in data
+            start_ind = searchsortedfirst.(Ref(lc.time), g[1])
+            end_ind = searchsortedlast.(Ref(lc.time), g[2])
+
+            t_temp = @view(lc.time[start_ind:end_ind])
+            c_temp = @view(lc.counts[start_ind:end_ind])
+
+            if !isempty(lc.err)
+                e_temp = @view(lc.err[start_ind:end_ind])
+            else
+                e_temp = Float64[]
+            end
+
+            bin_t, bin_c, bin_e, _ = rebin_data(t_temp, c_temp, dt_new,
+                                                      y_err=e_temp, method=method)
+
+            append!(bin_time, bin_t)
+            append!(bin_counts, bin_c)
+            append!(bin_err, bin_e)
+            push!(gti_new, [g[1], g[2]])
+        end
+    end
+
+    if isempty(gti_new) & !isempty(lc.gti)
+        throw(ArgumentError("No valid GTIs after rebin."))
+    elseif isempty(gti_new)
+        gti_new = lc.gti
+    else
+        gti_new = mapreduce(permutedims, vcat, gti_new)
+    end
+
+    return LightCurve(time=bin_time, counts=bin_counts, err=bin_err,
+                      mjdref=lc.mjdref, dt=dt_new, gti=gti_new,
+                      skip_checks=true)
+end
diff --git a/src/utils.jl b/src/utils.jl
index c92520c..8c6e2d7 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -27,3 +27,90 @@ function contiguous_regions(condition::AbstractVector{Bool})
     # Reshape the result into two columns
     return reshape(idx,2,length(idx) ÷ 2)'
 end
+
+function rebin_data(x::AbstractVector{<:Real}, y::AbstractVector{<:Integer}, dx_new::Real;
+                    y_err::AbstractVector{<:Real}=Float64[], method::String = "sum", dx::Real=0)
+    if isempty(y_err)
+        y_err = zero(y)
+    end
+
+    if dx isa AbstractVector
+        dx_old = dx
+    elseif dx!=0
+        dx_old = [dx]
+    else
+        dx_old = diff(x)
+    end
+
+    if any(dx_new .< dx_old)
+        throw(ArgumentError("New frequency resolution must be larger than old frequency resolution."))
+    end
+
+    # left and right bin edges
+    # assumes that the points given in `x` correspond to
+    # the left bin edges
+    xedges = vcat(x, x[end]+dx_old[end])
+
+    # new regularly binned resolution
+    xbin = range(xedges[begin], xedges[end], step = dx_new)
+
+    output = zeros(length(xbin) - 1)
+    outputerr = zeros(length(xbin) - 1)
+    step_size = zeros(length(xbin) - 1)
+
+    all_x = searchsortedfirst.(Ref(xedges), xbin)
+    min_inds = @view(all_x[begin:end-1])
+    max_inds = @view(all_x[begin+1:end])
+    xmins = @view(xbin[begin:end-1])
+    xmaxs = @view(xbin[begin+1:end])
+
+    for (i, (xmin, xmax, min_ind, max_ind)) in enumerate(zip(xmins, xmaxs, min_inds, max_inds))
+        filtered_y = @view(y[min_ind:max_ind-1])
+        filtered_yerr = @view(y_err[min_ind:max_ind-1])
+        output[i] = sum(filtered_y)
+        outputerr[i] = sum(filtered_yerr)
+        step_size[i] = max_ind - 1 - min_ind
+
+        prev_dx = xedges[min_ind+1] - xedges[min_ind]
+        prev_frac = (xedges[min_ind] - xmin)/prev_dx
+        output[i] += y[min_ind]*prev_frac
+        outputerr[i] += y_err[min_ind]*prev_frac
+        step_size[i] += prev_frac
+
+        if !(max_ind == length(xedges))
+            dx_post = xedges[max_ind] - xedges[max_ind-1]
+            post_frac = (xmax-xedges[max_ind-1])/dx_post
+            output[i] += y[max_ind-1]*post_frac
+            outputerr[i] += y_err[max_ind-1]*post_frac
+            step_size[i] += post_frac
+        end
+    end
+
+    if method in ["mean", "avg", "average", "arithmetic mean"]
+        ybin = output / step_size
+        ybinerr = sqrt.(outputerr) / step_size
+    elseif method == "sum"
+        ybin = output
+        ybinerr = sqrt.(outputerr)
+    else
+        throw(ArgumentError("Method for summing or averaging not recognized. Please enter either 'sum' or 'mean'."))
+    end
+    tseg = x[end] - x[begin] + dx_old[end]
+
+    if (tseg / dx_new % 1) > 0
+        pop!(ybin)
+        pop!(ybinerr)
+        pop!(step_size)
+    end
+
+    dx_var = Statistics.var(dx_old) / Statistics.mean(dx_old)
+
+    if size(dx_old) == 1 || dx_var < 1e-6
+        new_step = step_size[1]
+    end
+
+    new_x0 = (x[1] - (0.5 * dx_old[1])) + (0.5 * dx_new)
+    xbin = 1:length(ybin) * dx_new .+ new_x0
+
+    return xbin, ybin, ybinerr, new_step
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a4460f4..a864cae 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,3 +4,5 @@ using FFTW, Distributions, Statistics, StatsBase, Metadata, HDF5
 
 include("test_fourier.jl")
 include("test_gti.jl")
+include("test_events.jl")
+include("test_lightcurve.jl")
diff --git a/test/test_events.jl b/test/test_events.jl
new file mode 100644
index 0000000..a86cd0f
--- /dev/null
+++ b/test/test_events.jl
@@ -0,0 +1,124 @@
+@testset "test_events" begin
+    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]]
+    gti = [0 4]
+
+    @testset "test_initialize_eventList" begin
+        times = sort(rand(Uniform(1e8,1e8+100),101))
+        ev = EventList(time=times, mjdref=54600)
+        @test ev.time ≈ times atol = 1e-15
+        @test ev.mjdref == 54600 
+    end
+
+    @testset "test_to_and_from_lightcurve" begin
+        @testset "test_from_lc" begin
+            lc = LightCurve(time=[0.5, 1.5, 2.5], counts=[2, 1, 2])
+            ev = from_lc(lc)
+            @test ev.time == [0.5, 0.5, 1.5, 2.5, 2.5]
+        end
+
+        @testset "test_to_lc" begin
+            ev = EventList(time=time, gti=gti)
+            lc = to_lc(ev, 1)
+            @test lc.time == [0.5, 1.5, 2.5, 3.5]
+            @test lc.gti == gti
+        end
+    end
+
+    @testset "test_join_events" begin
+        @testset "test_join_different_dt" begin
+            ev = EventList(time=[10, 20, 30], dt=1)
+            ev_other = EventList(time=[40, 50, 60], dt=3)
+            ev_new = Stingray.join(ev,ev_other)
+            @test ev_new.dt == 3
+        end
+        @testset "test_join_different_instr" begin
+            ev = EventList(time=[10, 20, 30], instr="fpma")
+            ev_other = EventList(time=[40, 50, 60], instr="fpmb")
+            ev_new = Stingray.join(ev,ev_other)
+            @test ev_new.instr == "fpma,fpmb"
+        end
+        @testset "test_join_without_energy" begin
+            ev = EventList(time=[1, 2, 3], energy=[3, 3, 3])
+            ev_other = EventList(time=[4, 5])
+            ev_new = Stingray.join(ev, ev_other)
+            @test ev_new.energy == [3, 3, 3, 0, 0]
+        end
+        @testset "test_join_without_pi" begin
+            ev = EventList(time=[1, 2, 3], PI=[3, 3, 3])
+            ev_other = EventList(time=[4, 5])
+            ev_new = Stingray.join(ev, ev_other)
+            @test ev_new.PI == [3, 3, 3, 0, 0]
+        end
+        @testset "test_join_with_gti_none" begin
+            ev = EventList(time=[1, 2, 3])
+            ev_other = EventList(time=[4, 5], gti=[3.5 5.5])
+            ev_new = Stingray.join(ev, ev_other)
+            @test ev_new.gti == [[1 3]; [3.5 5.5];]
+
+            ev = EventList(time=[1, 2, 3], gti=[0.5 3.5])
+            ev_other = EventList(time=[4, 5])
+            ev_new = Stingray.join(ev,ev_other)
+            @test ev_new.gti == [[0.5 3.5]; [4 5];]
+
+            ev = EventList(time=[1, 2, 3])
+            ev_other = EventList(time=[4, 5])
+            ev_new = Stingray.join(ev, ev_other)
+            @test isempty(ev_new.gti)
+        end
+        @testset "test_non_overlapping_join" begin
+            ev = EventList(time=[1, 1, 2, 3, 4],
+                        energy=[3, 4, 7, 4, 3], gti=[[1 2]; [3 4];])
+            ev_other = EventList(time=[5, 6, 6, 7, 10],
+                                energy=[4, 3, 8, 1, 2], gti=[6 7])
+            ev_new = Stingray.join(ev, ev_other)
+
+            @test ev_new.time == [1, 1, 2, 3, 4, 5, 6, 6, 7, 10]
+            @test ev_new.energy == [3, 4, 7, 4, 3, 4, 3, 8, 1, 2]
+            @test ev_new.gti == [[1 2]; [3 4]; [6 7];]
+        end
+        @testset "test_overlapping_join" begin
+            ev = EventList(time=[1, 1, 10, 6, 5],
+                        energy=[10, 6, 3, 11, 2], gti=[[1 3]; [5 6];])
+            ev_other = EventList(time=[5, 7, 6, 6, 10],
+                                energy=[2, 3, 8, 1, 2], gti=[[5 7]; [8 10];])
+            ev_new = Stingray.join(ev, ev_other)
+
+            @test ev_new.time == [1, 1, 5, 5, 6, 6, 6, 7, 10, 10]
+            @test ev_new.energy == [10, 6, 2, 2, 11, 8, 1, 3, 3, 2]
+            @test ev_new.gti == [5 6]
+        end
+    end
+    @testset "test_sort" begin
+        ev = EventList(time=[2, 4, 1, 4, 3], energy=[3, 4, 7, 4, 3])
+        @testset "test_sort_not_inplace" begin
+            new_ev = Stingray.sort(ev);
+            @test new_ev.time == [1, 2, 3, 4, 4]
+            @test new_ev.energy == [7, 3, 3, 4, 4]
+            @test ev.time == [2, 4, 1, 4, 3]
+        end
+        @testset "test_sort_inplace" begin
+            Stingray.sort!(ev);
+            @test ev.time == [1, 2, 3, 4, 4]
+            @test ev.energy == [7, 3, 3, 4, 4] 
+        end
+    end
+    @testset "Input/Output" begin
+        @testset "test_fits_with_standard_file" begin
+            fname = joinpath(@__DIR__ ,"data","monol_testA.evt")
+            ev = read(EventList, fname, "fits")
+            @test ev.mjdref == 55197.00076601852
+        end
+
+        @testset "test_io_with_fits" begin
+            ev = EventList(time=time, mjdref=54000)
+            write(ev, "ev.fits", "fits")
+            new_ev = read(EventList, "ev.fits", "fits")
+            @test new_ev.time == time
+            rm("ev.fits")
+        end
+    end
+end
+
diff --git a/test/test_lightcurve.jl b/test/test_lightcurve.jl
new file mode 100644
index 0000000..f9a0018
--- /dev/null
+++ b/test/test_lightcurve.jl
@@ -0,0 +1,49 @@
+@testset "test_lightcurve" begin
+    @testset "rebin_lightcurve" begin
+        dt = 0.0001220703125
+        n = 1384132
+        mean_counts = 2.0
+        times = range(dt/2, dt/2 + n*dt, step = dt)
+        counts = zero(times) .+ mean_counts
+        lc = LightCurve(time=times, counts=counts)
+
+        @testset "test_rebin_even" begin
+            dt_new = 2.0
+            lc_binned = rebin(lc, dt_new)
+            @test lc_binned.dt == dt_new
+            counts_test = zero(lc_binned.time) .+ lc.counts[1]*dt_new/lc.dt
+            @test lc_binned.counts == counts_test
+        end
+
+        @testset "test_rebin_odd" begin
+            dt_new = 1.5
+            lc_binned = rebin(lc, dt_new)
+            @test lc_binned.dt == dt_new
+            counts_test = zero(lc_binned.time) .+ lc.counts[1]*dt_new/lc.dt
+            @test lc_binned.counts == counts_test
+        end
+
+        @testset "test_rebin_several" for dt in [2,3,pi,5]
+            lc_binned = rebin(lc, dt)
+            @test length(lc_binned.time) == length(lc_binned.counts)
+        end
+
+        @testset "test_rebin_with_gtis" begin
+            times = collect(range(0, 99.9, step=0.1))
+
+            counts = round.(rand(Normal(100, 0.1), length(times)))
+            gti = [[0 40]; [60 100];]
+
+            good, gti = create_gti_mask(times, gti)
+
+            keepat!(times, good)
+            keepat!(counts, good)
+
+            lc = LightCurve(time=times,counts=counts, gti=gti, skip_checks=true, dt=0.1)
+
+            lc_rebin = rebin(lc, 1.0)
+
+            @test (lc_rebin.time[40] - lc_rebin.time[39]) >= 1.0
+        end
+    end
+end