using GWecc
using Test
using FiniteDifferences
using LinearAlgebra
using UnPack
using NumericalIntegration
using Statistics


e_from_Ï„_from_e(ecc::Float64)::Float64 = e_from_Ï„(Ï„_from_e(Eccentricity(ecc))).e

@testset verbose = true "GWecc" begin
    @testset "parameters" begin
        @testset "mass" begin
            @test_throws DomainError Mass(-4e-2, 0.2)
            @test_throws DomainError Mass(1.0, -1.0)
            @test_throws DomainError Mass(1.0, 0.0)
            # @test_throws DomainError Mass(6e5, 0.2)
            @test_throws DomainError Mass(1.0, 0.26)
            @test Mass(1.0, 0.1).m == 1.0
            @test Mass(1.0, 0.25).Î· == 0.25
            @test Mass(1.0, 0.1).m > Mass(1.0, 0.1).Mch
        end

        @testset "time" begin
            @test (-Time(-1.0)).t == (-1.0 * Time(-1.0)).t
            @test (Time(-1.0) * -1).t == (-1.0 * Time(-1.0)).t
        end

        @testset "gammabar" begin
            @test_throws DomainError ScaledPeriastronAngle(0.01, 0.01, -43.0)
            @test_throws DomainError ScaledPeriastronAngle(0.01, 2.50, -44.0)
            @test_throws DomainError ScaledPeriastronAngle(0.10, 0.24, -44.0)
        end

        @testset "init phase params" begin
            @test_throws DomainError InitPhaseParams(-7.0, 1.0)
            @test_throws DomainError InitPhaseParams(1.0, 7.0)
        end

        @testset "projection params" begin
            @test_throws DomainError ProjectionParams(-1e-9, 1.0, 0.3, 1.0, 1.0)
            @test_throws DomainError ProjectionParams(1e-9, 4.0, 0.3, 1.0, 1.0)
            @test_throws DomainError ProjectionParams(1e-9, 1.0, 1.1, 1.0, 1.0)
            @test_throws DomainError ProjectionParams(1e-9, 1.0, 0.3, 4.0, 1.0)
            @test_throws DomainError ProjectionParams(1e-9, 1.0, 0.3, 1.0, 4.0)

            @test_throws DomainError ProjectionParams1psr(-1e-9, 0.1, 0.1)
            @test_throws DomainError ProjectionParams1psr(1e-9, -0.1, 0.1)
            @test_throws DomainError ProjectionParams1psr(1e-9, 0.1, 10.0)
        end

        @testset "sky location params" begin
            @test_throws DomainError SkyLocation(-1.0, 1.0)
            @test_throws DomainError SkyLocation(1.0, 4.1)
        end
    end

    @testset "orbital evlution" begin
        @testset "e and Ï„" begin
            @test_throws DomainError Eccentricity(0.0)
            @test Eccentricity(0.1).e == 0.1
            @test_throws DomainError Eccentricity(1.0)

            @test_throws DomainError ScaledTime(-1.0)
            @test ScaledTime(0.0).Ï„ == 0.0
            @test ScaledTime(10.0).Ï„ == 10.0
            @test_throws DomainError ScaledTime(Inf)

            for e in [eccmin.e / 4, eccmin.e, 0.1, 0.5, 0.9, eccmax.e, (eccmax.e + 1) / 2]
                @test e_from_Ï„_from_e(e) â‰ˆ e
            end
        end

        @testset "coeffs" begin
            mass = Mass(1.0, 0.1)
            n_init = MeanMotion(1e-8)
            e_init = Eccentricity(0.1)
            coeffs = EvolvCoeffs(mass, n_init, e_init)
            @test coeffs.Îº > 0 && isfinite(coeffs.Îº)
            @test coeffs.Î± > 0 && isfinite(coeffs.Î±)
            @test coeffs.Î² > 0 && isfinite(coeffs.Î²)
            @test coeffs.Î²2 > 0 && isfinite(coeffs.Î²2)
            @test coeffs.Î²3 > 0 && isfinite(coeffs.Î²3)
        end

        @testset "lbar and Î³bar" begin
            mass = Mass(1.0, 0.1)
            n = MeanMotion(1e-8)
            e1 = Eccentricity(0.1)
            e2 = Eccentricity(0.2)

            lbar1 = lbar_from_e(e1)
            lbar2 = lbar_from_e(e2)
            @test lbar1.lbar < lbar2.lbar

            Î³bar1 = Î³bar_from_e(e1)
            Î³bar2 = Î³bar_from_e(e2)
            @test Î³bar1 < Î³bar2
        end

        @testset "evolve_orbit" begin
            mass = Mass(5000.0, 0.1)
            n_init = MeanMotion(1e-8)
            e_init = Eccentricity(0.1)
            l_init = Angle(0.1)
            Î³_init = Angle(-1.25)
            coeffs = EvolvCoeffs(mass, n_init, e_init)

            delay0 = Time(0.0)
            n0, e0, l0, Î³0 = evolve_orbit(coeffs, l_init, Î³_init, delay0)
            @test n0.n â‰ˆ n_init.n atol = 1e-9
            @test e0.e â‰ˆ e_init.e atol = 1e-9
            @test l0.Î¸ â‰ˆ l_init.Î¸ atol = 1e-8
            @test Î³0.Î¸ â‰ˆ Î³_init.Î¸ atol = 1e-9

            delay1 = Time(-10000.0)
            n1, e1, l1, Î³1 = evolve_orbit(coeffs, l_init, Î³_init, delay1)
            @test n1.n < n_init.n
            @test e1.e > e_init.e
            @test l1.Î¸ < l_init.Î¸
            @test Î³1.Î¸ < Î³_init.Î¸

            delay2 = Time(10000.0)
            n2, e2, l2, Î³2 = evolve_orbit(coeffs, l_init, Î³_init, delay2)
            @test n2.n > n_init.n
            @test e2.e < e_init.e
            @test l2.Î¸ > l_init.Î¸
            @test Î³2.Î¸ > Î³_init.Î¸
        end
    end

    @testset "derivatives" begin
        numdiff = central_fdm(5, 1)

        mass = Mass(5000.0, 0.1)
        n = MeanMotion(1e-8)
        l = Angle(0.0)
        Î³ = Angle(0.0)

        for e in [0.1, 0.5, 0.9]
            ecc = Eccentricity(e)
            coeffs = EvolvCoeffs(mass, n, ecc)

            dÏ„_de_anl = derivative_dÏ„_de(ecc)
            dÏ„_de_num = numdiff(Ï„_from_e, e)
            @test dÏ„_de_anl â‰ˆ dÏ„_de_num atol = 1e-9

            Ï„ = Ï„_from_e(e)
            de_dÏ„_num = numdiff(e_from_Ï„, Ï„)
            de_dÏ„_anl = 1 / dÏ„_de_anl
            @test de_dÏ„_anl â‰ˆ de_dÏ„_num atol = 1e-9

            Îº = evolv_coeff_Îº(mass, n, ecc)
            de_dt_anl1 = derivative_de_dt(mass, n, ecc)
            de_dt_anl2 = -Îº * de_dÏ„_anl
            @test de_dt_anl1 â‰ˆ de_dt_anl2

            dn_dt_anl = derivative_dn_dt(mass, n, ecc)
            dn_de_anl = dn_dt_anl / de_dt_anl1
            n_from_e_func = e1 -> n_from_e(coeffs, Eccentricity(e1)).n
            dn_de_num = numdiff(n_from_e_func, e)
            @test dn_de_anl â‰ˆ dn_de_num

            dlbar_de_anl = derivative_dlbar_de(ecc)
            dlbar_de_num = numdiff(lbar_from_e, e)
            @test dlbar_de_anl â‰ˆ dlbar_de_num atol = 1e-9

            Î± = evolv_coeff_Î±(mass, n, ecc)
            dlbar_dÏ„_anl = dlbar_de_anl * de_dÏ„_anl
            dl_dÏ„_anl = -Î± * dlbar_dÏ„_anl
            dÏ„_dt_anl = -Îº
            dl_dt_anl2 = dl_dÏ„_anl * dÏ„_dt_anl
            dl_dt_anl1 = n.n
            @test dl_dt_anl2 â‰ˆ dl_dt_anl1 atol = 1e-9

            dÎ³bar_de_anl = derivative_dÎ³bar_de(ecc)
            dÎ³bar_de_num = numdiff(Î³bar_from_e, e)
            @test dÎ³bar_de_anl â‰ˆ dÎ³bar_de_num atol = 1e-9

            Î² = evolv_coeff_Î²(mass, n, ecc)
            dÎ³bar_dÏ„_anl = dÎ³bar_de_anl * de_dÏ„_anl
            dÎ³_dÏ„_anl = -Î² * dÎ³bar_dÏ„_anl
            dÎ³_dt_anl2 = dÎ³_dÏ„_anl * dÏ„_dt_anl
            k1 = advance_of_periastron_1PN(mass, n, ecc).k
            dÎ³_dt_anl1 = k1 * n.n
            @test dÎ³_dt_anl2 â‰ˆ dÎ³_dt_anl1 atol = 1e-9

            dÎ³bar2_de_anl = derivative_dÎ³bar2_de(ecc, mass)
            Î³bar2_from_e_func = e1 -> Î³bar2_from_e(Eccentricity(e1), mass)
            dÎ³bar2_de_num = numdiff(Î³bar2_from_e_func, e)
            @test dÎ³bar2_de_anl â‰ˆ dÎ³bar2_de_num atol = 1e-9

            Î²2 = evolv_coeff_Î²2(mass, n, ecc)
            dÎ³bar2_dÏ„_anl = dÎ³bar2_de_anl * de_dÏ„_anl
            dÎ³2_dÏ„_anl = -Î²2 * dÎ³bar2_dÏ„_anl
            dÎ³2_dt_anl2 = dÎ³2_dÏ„_anl * dÏ„_dt_anl
            k2 = advance_of_periastron_2PN(mass, n, ecc).k
            dÎ³2_dt_anl1 = k2 * n.n
            @test dÎ³2_dt_anl2 â‰ˆ dÎ³2_dt_anl1 atol = 1e-9

            dÎ³bar3_de_anl = derivative_dÎ³bar3_de(ecc, mass)
            Î³bar3_from_e_func = e1 -> Î³bar3_from_e(Eccentricity(e1), mass)
            dÎ³bar3_de_num = numdiff(Î³bar3_from_e_func, e)
            # This comparison is so imprecise because Î³bar3_from_e implements a
            # Pade approximant rather than the exact expression. This is OK because
            # this is a 3PN correction.
            @test dÎ³bar3_de_anl â‰ˆ dÎ³bar3_de_num rtol = 1e-2

            Î²3 = evolv_coeff_Î²3(mass, n, ecc)
            dÎ³bar3_dÏ„_anl = dÎ³bar3_de_anl * de_dÏ„_anl
            dÎ³3_dÏ„_anl = -Î²2 * dÎ³bar3_dÏ„_anl
            dÎ³3_dt_anl2 = dÎ³3_dÏ„_anl * dÏ„_dt_anl
            k3 = advance_of_periastron_3PN(mass, n, ecc).k
            dÎ³3_dt_anl1 = k3 * n.n
            @test dÎ³3_dt_anl2 â‰ˆ dÎ³3_dt_anl1 atol = 1e-9
        end
    end

    @testset "orbital phase" begin
        @testset "sincos" begin
            for Î¸ = -7.0:0.5:7.0
                scÎ¸ = SinCos(Angle(Î¸))
                @test scÎ¸.sinx^2 + scÎ¸.cosx^2 â‰ˆ 1.0
            end
        end

        @testset "mikkola" begin
            e = Eccentricity(0.1)
            l = Angle(0.0)
            @test l == mikkola(e, l)

            for e in [0.1, 0.5, 0.9]
                ecc = Eccentricity(e)
                for _l = -7.0:0.5:7.0
                    l = Angle(_l)
                    u = mikkola(ecc, l)
                    @test l.Î¸ â‰ˆ kepler(ecc, u).Î¸
                end
            end
        end

        @testset "true anomaly" begin
            for e in [0.1, 0.5, 0.9]
                ecc = Eccentricity(e)
                for _l in LinRange(-4 * Ï€, 4 * Ï€, 16)
                    l = Angle(_l)
                    u = mikkola(ecc, l)
                    scu = SinCos(u)
                    v_l = true_anomaly_diff(ecc, scu)
                    v = true_anomaly(ecc, scu)
                    @test v_l.Î¸ â‰ˆ v.Î¸ - l.Î¸ atol = 1e-9
                end
            end
        end

        mass = Mass(5000.0, 0.1)
        n = MeanMotion(1e-8)
        Î³ = Angle(3.5)
        Î³_ = Î³.Î¸

        for l_ in [0.0, 1.3, 2.3, 3.3, 4.3]
            l = Angle(l_)
            for e in [1e-9, 0.1, 0.5, 0.9]
                et = Eccentricity(e)
                orbital_phase = OrbitalPhase(mass, n, et, l, Î³)
                Ï† = orbital_phase.sc2Ï†.x.Î¸ / 2
                u = orbital_phase.scu.x

                if l == 0.0
                    @test Ï† â‰ˆ Î³_
                end

                if e == 1e-9
                    @test Ï† â‰ˆ l_ + Î³_
                end
            end
        end
    end

    @testset "antenna pattern" begin
        ra_p = 1.5
        dec_p = -0.8
        ra_gw = 0.5
        dec_gw = 0.75

        psrpos = SkyLocation(ra_p, dec_p)
        gwpos = SkyLocation(ra_gw, dec_gw)

        nhat = sky_direction_uvec(gwpos)
        @test norm(nhat) â‰ˆ 1.0

        phat = sky_direction_uvec(psrpos)
        @test norm(phat) â‰ˆ 1.0

        ep, ec = gw_polarization_tensors(gwpos)
        @test ep * ec + ec * ep â‰ˆ zeros(3, 3) atol = 1e-9
        @test ep * nhat â‰ˆ zeros(3) atol = 1e-9
        @test ec * nhat â‰ˆ zeros(3) atol = 1e-9

        ap1 = AntennaPattern(psrpos, gwpos)
        ap2 = AntennaPattern(gwpos, psrpos)
        @test ap1.cosÂµ â‰ˆ ap2.cosÂµ
        @test ap1.Fp^2 + ap1.Fx^2 â‰ˆ ap2.Fp^2 + ap2.Fx^2

        @test_throws DomainError AntennaPattern(gwpos, gwpos)

        # @test AzimuthParam(0.5).Î± == 0.5
        # @test_throws DomainError AzimuthParam(1.1)
    end

    @testset "mismatch" begin
        a = [1, 2, 3]
        b = [2, -1, 0]
        K = [1 0 0; 0 1 0; 0 0 0]

        @test mismatch(a, a) â‰ˆ 0.0
        @test mismatch(K, a, a) â‰ˆ 0.0
        @test mismatch(a, -a) â‰ˆ 2.0
        @test mismatch(K, a, -a) â‰ˆ 2.0
        @test mismatch(a, b) â‰ˆ 1.0
        @test mismatch(K, a, b) â‰ˆ 1.0
    end

    @testset "waveform and residuals" begin
        mass = Mass(5000.0, 0.1)
        n_init = MeanMotion(1e-8)
        e_init = Eccentricity(0.1)
        l_init = Angle(0.1)
        Î³_init = Angle(1.25)
        coeffs = EvolvCoeffs(mass, n_init, e_init)
        dl = Distance(1e15)

        h0 = gw_amplitude(mass, n_init, e_init, dl)
        @test h0 > 0

        n1 = MeanMotion(2e-8)
        e1 = Eccentricity(0.09)
        h1 = gw_amplitude(mass, n1, e1, dl)

        c01 = gwres_amplitude_ratio(mass, n_init, e_init, n1, e1)
        @test c01 â‰ˆ (h1 / h0) / (n1.n / n_init.n)

        ra_p = 1.5
        dec_p = -0.8
        ra_gw = 0.5
        dec_gw = 0.75
        psrpos = SkyLocation(ra_p, dec_p)
        gwpos = SkyLocation(ra_gw, dec_gw)
        dp = Distance(1e13)
        ap = AntennaPattern(psrpos, gwpos)
        Î± = azimuth_param(ap)
        Î”p = pulsar_term_delay(ap, dp)
        @test Î”p.t < 0

        Ïˆ = 1.1
        cosÎ¹ = 0.52
        Î³0 = Î³_init.Î¸
        Î³p = Î³0 + 0.2
        S0 = 1e-9
        proj = ProjectionParams(S0, Ïˆ, cosÎ¹, Î³0, Î³p)

        l0p = InitPhaseParams(l_init.Î¸, l_init.Î¸)

        dt = Time(10000.0)
        dtp = dt + Î”p

        tref = Time(5 * 365.25 * 24 * 3600)
        tEs = time_range(Time(0.0), tref, 100)

        @testset "single point functions" begin

            sE = residual(mass, coeffs, l0p, proj, ap, [EARTH], Î”p, dt)
            sP = residual(mass, coeffs, l0p, proj, ap, [PULSAR], Î”p, dt)
            s = residual(mass, coeffs, l0p, proj, ap, [EARTH, PULSAR], Î”p, dt)
            @test s â‰ˆ sP + sE

            spE, sxE = residual_px(mass, coeffs, l0p, proj, false, dt)
            spP, sxP = residual_px(mass, coeffs, l0p, proj, true, dtp)
            @test sP â‰ˆ -(ap.Fp * spP + ap.Fx * sxP)
            @test sE â‰ˆ (ap.Fp * spE + ap.Fx * sxE)

            hE = waveform(mass, coeffs, l0p, proj, ap, [EARTH], Î”p, dt)
            hP = waveform(mass, coeffs, l0p, proj, ap, [PULSAR], Î”p, dt)
            h = waveform(mass, coeffs, l0p, proj, ap, [EARTH, PULSAR], Î”p, dt)
            @test h â‰ˆ hP + hE

            hpE, hxE = waveform_px(mass, coeffs, l0p, proj, false, dt)
            hpP, hxP = waveform_px(mass, coeffs, l0p, proj, true, dtp)
            @test hP â‰ˆ -(ap.Fp * hpP + ap.Fx * hxP)
            @test hE â‰ˆ (ap.Fp * hpE + ap.Fx * hxE)
        end

        @testset "terms and polarizations" begin
            rs = residuals(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [EARTH, PULSAR],
                tref,
                tEs,
            )
            rEs = residuals(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [EARTH],
                tref,
                tEs,
            )
            rPs = residuals(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [PULSAR],
                tref,
                tEs,
            )
            @test all(isfinite.(rs))
            @test all(isapprox.(rs, rEs + rPs))

            sps, sxs = residuals_px(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                EARTH,
                tref,
                tEs,
            )
            @test all(isfinite.(sps)) && all(isfinite.(sxs))

            hs = waveform(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [EARTH, PULSAR],
                tref,
                tEs,
            )
            hEs = waveform(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [EARTH],
                tref,
                tEs,
            )
            hPs = waveform(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [PULSAR],
                tref,
                tEs,
            )
            @test all(isfinite.(hs))
            @test all(isapprox.(hs, hEs + hPs))

            rs1, hs1 = residuals_and_waveform(
                mass,
                n_init,
                e_init,
                l0p,
                proj,
                dp,
                psrpos,
                gwpos,
                [EARTH, PULSAR],
                tref,
                tEs,
            )
            @test all(isapprox.(hs1, hs))
            @test all(isapprox.(rs1, rs))
        end

        @testset "component functions" begin
            for e_init in Eccentricity.([0.1, 0.4, 0.8])
                sE = residuals(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [EARTH],
                    tref,
                    tEs,
                )
                sP = residuals(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [PULSAR],
                    tref,
                    tEs,
                )
                sEc = residuals_from_components(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [EARTH],
                    tref,
                    tEs,
                )
                sPc = residuals_from_components(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [PULSAR],
                    tref,
                    tEs,
                )
                @test all(isapprox.(sEc, sE, atol = 1e-8))
                @test all(isapprox.(sPc, sP, atol = 1e-8))

                ð�’œs = residuals_components_ð�’œ(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    dp,
                    psrpos,
                    gwpos,
                    EARTH,
                    tref,
                    tEs,
                )
                @test all([all(isfinite.(ð�’œ)) for ð�’œ in ð�’œs])
            end
        end

        @testset "1psr functions" begin
            for e_init in Eccentricity.([0.1, 0.4, 0.8])
                dÏˆ = polarization_angle_shift_1psr(ap)
                l0p0 = InitPhaseParams(l_init.Î¸)
                proj0 = ProjectionParams(S0, Ïˆ, cosÎ¹, Î³0, Î³0)
                proj1 = ProjectionParams(S0 * Î±, Ïˆ + dÏˆ, cosÎ¹, Î³0)
                proj1_ = ProjectionParams1psr(proj0, ap)
                Î”p = pulsar_term_delay(ap, dp)
                ss = residuals(
                    mass,
                    n_init,
                    e_init,
                    l0p0,
                    proj0,
                    dp,
                    psrpos,
                    gwpos,
                    [EARTH, PULSAR],
                    tref,
                    tEs,
                )

                hs = waveform(
                    mass,
                    n_init,
                    e_init,
                    l0p0,
                    proj0,
                    dp,
                    psrpos,
                    gwpos,
                    [EARTH, PULSAR],
                    tref,
                    tEs,
                )

                hs_new = waveform_1psr(
                    mass,
                    n_init,
                    e_init,
                    l_init,
                    proj1_,
                    Î”p,
                    [EARTH, PULSAR],
                    tref,
                    tEs,
                )
                @test all(isapprox.(hs, hs_new))

                ss_new = residuals_1psr(
                    mass,
                    n_init,
                    e_init,
                    l_init,
                    proj1_,
                    Î”p,
                    [EARTH, PULSAR],
                    tref,
                    tEs,
                )
                @test all(isapprox.(ss, ss_new))
            end
        end

        @testset "h = ds/dt" begin
            year = 365.25 * 24 * 3600
            _ts = LinRange(0, 10 * year, 10000)
            tEs = Time.(_ts)
            tref = Time(maximum(_ts))

            for term in [EARTH, PULSAR]
                hps, hxs = waveform_px(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    term,
                    tref,
                    tEs,
                )
                rps, rxs = residuals_px(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    term,
                    tref,
                    tEs,
                )

                rps_n = cumul_integrate(_ts, hps)
                rxs_n = cumul_integrate(_ts, hxs)

                rps_n = rps_n .+ (mean(rps) - mean(rps_n))
                rxs_n = rxs_n .+ (mean(rxs) - mean(rxs_n))

                @test mismatch(rps, rps_n) < 1e-3
                @test mismatch(rxs, rxs_n) < 1e-3

                hs = waveform(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [term],
                    tref,
                    tEs,
                )
                rs = residuals(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [term],
                    tref,
                    tEs,
                )

                rs_n = cumul_integrate(_ts, hs)
                rs_n = rs_n .+ (mean(rs) - mean(rs_n))

                @test mismatch(rs, rs_n) < 1e-3
            end
        end

        @testset "spline functions" begin
            year = 365.25 * 24 * 3600
            _ts = LinRange(0, 10 * year, 10000)
            tEs = Time.(_ts)
            tref = Time(maximum(_ts))

            proj1 = ProjectionParams1psr(proj, ap)

            for term in [EARTH, PULSAR]
                rs = residuals(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [term],
                    tref,
                    tEs,
                )
                rs_spl = residuals_spline(
                    mass,
                    n_init,
                    e_init,
                    l0p,
                    proj,
                    dp,
                    psrpos,
                    gwpos,
                    [term],
                    tref,
                    tEs,
                )
                @test mismatch(rs, rs_spl) < 1e-3

                Î”p = pulsar_term_delay(ap, dp)
                rs = residuals_1psr(
                    mass,
                    n_init,
                    e_init,
                    l_init,
                    proj1,
                    Î”p,
                    [term],
                    tref,
                    tEs,
                )
                rs_spl = residuals_1psr_spline(
                    mass,
                    n_init,
                    e_init,
                    l_init,
                    proj1,
                    Î”p,
                    [term],
                    tref,
                    tEs,
                )
                @test mismatch(rs, rs_spl) < 1e-3
            end
        end
    end

    @testset "enterprise functions" begin
        toas = LinRange(0, 1000000, 100)
        pdist = 400.0 # kpc
        # alpha = azimuth_param(ap)
        psi = 1.1
        cos_inc = 0.5
        log10_M = 9.0
        eta = 0.25
        log10_F = -8.0
        e0 = 0.3
        gamma0 = 0.0
        gammap = 0.0
        sigma = 0.2
        rho = 0.15
        l0 = 0.0
        lp = 0.0
        tref = maximum(toas)
        log10_dl = -15.0
        deltap = 100.0
        for spline in [true, false]
            for psrTerm in [true, false]
                res = eccentric_pta_signal_1psr(
                    toas,
                    sigma,
                    rho,
                    log10_M,
                    eta,
                    log10_F,
                    e0,
                    l0,
                    tref,
                    log10_dl,
                    deltap,
                    psrTerm,
                )
                @test all(isfinite.(res))
            end
        end

        ra_p = 1.5
        dec_p = -0.8
        ra_gw = 0.5
        dec_gw = 0.75
        for spline in [true, false]
            for psrTerm in [true, false]
                res = eccentric_pta_signal(
                    toas,
                    Ï€ / 2 - dec_p,
                    ra_p,
                    pdist,
                    sin(dec_gw),
                    ra_gw,
                    psi,
                    cos_inc,
                    log10_M,
                    eta,
                    log10_F,
                    e0,
                    gamma0,
                    gammap,
                    l0,
                    lp,
                    tref,
                    log10_dl,
                    psrTerm,
                    spline,
                )
                @test all(isfinite.(res))
            end
        end

        log10_A = -9.0
        gwdist = 100.0
        for spline in [true, false]
            for psrTerm in [true, false]
                res = eccentric_pta_signal_target(
                    toas,
                    Ï€ / 2 - dec_p,
                    ra_p,
                    pdist,
                    sin(dec_gw),
                    ra_gw,
                    psi,
                    cos_inc,
                    eta,
                    log10_F,
                    e0,
                    gamma0,
                    gammap,
                    l0,
                    lp,
                    tref,
                    log10_A,
                    gwdist,
                    psrTerm,
                    spline,
                )
                @test all(isfinite.(res))
            end
        end

        @test validate_params(log10_M, eta, log10_F, e0, tref, tref)[1]
        @test !validate_params(10.0, eta, -7.0, 0.999, tref, tref)[1]

        gwdist = 100.0
        log10_A = -9.0
        @test validate_params_target(log10_A, eta, log10_F, e0, gwdist, tref, tref)[1]
        @test !validate_params_target(-5.0, eta, -7.0, 0.999, gwdist, tref, tref)[1]

        # Roundtrip test for converting mass to/from amplitude
        # for targetted search.
        eta = 0.23
        m1 = mass_from_log10_mass(9.25, eta)
        log10_F = -7.21
        n0 = mean_motion_from_log10_freq(log10_F)
        e0 = Eccentricity(0.3)
        gwdist = 0.931
        dl = dl_from_gwdist(gwdist)
        ampl = gw_amplitude(m1, n0, e0, dl) / n0.n
        log10_A = log10(ampl)
        m2 = mass_from_gwdist(log10_A, log10_F, e0.e, gwdist, eta)
        @test m1.m â‰ˆ m2.m atol = 1e-6
    end

end
