text:
# Exceptions.

"""
An error encountered during the HTTP request.

## Fields
- `response::Union{HTTP.Response, Nothing}`: Set for status exceptions.
- `exception::Exception`
- `stacktrace::StackTrace`
"""
struct HTTPError{E<:Exception} <: ForgeError
    response::Union{HTTP.Response, Nothing}
    exception::E
    stacktrace::StackTrace
end

HTTPError(ex::Exception, st::StackTrace) = HTTPError(nothing, nothing, ex, st)

"""
An error encountered during response postprocessing.

## Fields
- `response::HTTP.Response`
- `exception::Exception`
- `stacktrace::StackTrace`
"""
struct PostProcessorError{E<:Exception} <: ForgeError
    response::HTTP.Response
    exception::E
    stacktrace::StackTrace
end

"""
A signal that a rate limit has been exceeded.

## Fields
- `period::Union{Period, Nothing}` Amount of time until rate limit expiry, if known.
"""
struct RateLimitedError <: ForgeError
    period::Union{Period, Nothing}
end

# Response postprocessing.

"""
Determines the behaviour of [`postprocess`](@ref).
"""
abstract type PostProcessor end

"""
Does nothing and always returns `nothing`.
"""
struct DoNothing <: PostProcessor end

"""
    JSON(f::Function=identity) -> JSON

Parses a JSON response into a given type and runs `f` on that object.
"""
struct JSON <: PostProcessor
    f::Function
end
JSON() = JSON(identity)

"""
    DoSomething(::Function) -> DoSomething

Runs a user-defined postprocessor.
"""
struct DoSomething <: PostProcessor
    f::Function
end

"""
    postprocess(::PostProcessor, ::HTTP.Response, ::Type{T})

Computes a value to be returned from an HTTP response.
"""
postprocess(::DoNothing, ::HTTP.Response, ::Type) = nothing
postprocess(p::JSON, r::HTTP.Response, ::Type{T}) where T =
    p.f(JSON2.read(IOBuffer(r.body), T))
postprocess(p::DoSomething, r::HTTP.Response, ::Type) = p.f(r)

# Requests.

"""
    request(
        f::Forge, fun::Function, ep::Endpoint;
        headers::Vector{<:Pair}=HTTP.Header[],
        query::AbstractDict=Dict(),
        request_opts=Dict(),
        kwargs...,
    ) -> T, HTTP.Response

Make an HTTP request and return `T` and the response, where `T` is determined by [`into`](@ref).

## Arguments
- `f::Forge` A [`Forge`](@ref) subtype.
- `fun::Function`: The API function being called.
- `ep::Endpoint`: The endpoint information.

## Keywords
- `query::AbstractDict=Dict()`: Query string parameters to add to the request.
- `headers::Vector{<:Pair}=HTTP.Header[]`: Headers to add to the request.
- `request_opts=Dict()`: Keywords passed into `HTTP.request`.

Trailing keywords are sent as a JSON body for `PATCH`, `POST`, and `PUT` requests.
For other request types, the keywords are sent as query string parameters.

!!! note
    Every API function passes its keyword arguments into this function.
    Therefore, to customize behaviour for a single request, pass the above keywords to the API function.
"""
function request(
    f::Forge, fun::Function, ep::Endpoint;
    headers::Vector{<:Pair}=HTTP.Header[],
    query::AbstractDict=Dict(),
    request_opts=Dict(),
    kwargs...,
)
    global diag

    rate_limit(f, fun)
    url = base_url(f) * ep.url
    headers = vcat(request_headers(f, fun), ep.headers, headers)
    query = merge(request_query(f, fun), ep.query, query)
    opts = merge(request_kwargs(f, fun), Dict(pairs(request_opts)))
    body = if ep.method in (:PATCH, :POST, :PUT)
        JSON2.write(Dict(kwargs))
    else
        merge!(query, Dict(kwargs))
        HTTP.nobody
    end

    diag && @info "GitForge making request, diag = $diag"
    resp = try
        HTTP.request(
            ep.method, url, headers, body;
            query=query, opts...,
            status_exception=false, # We handle status exceptions ourselvse.
            verbose = diag ? 2 : 0,
        )
    catch e
        rethrow(HTTPError(e, stacktrace(catch_backtrace())))
    end

    has_rate_limits(f, fun) && rate_limit_update!(f, fun, resp)

    resp.status >= 300 && !(resp.status == 404 && ep.allow_404) &&
        throw(HTTPError(resp, HTTP.StatusError(resp.status, resp), stacktrace()))

    return try
        postprocess(postprocessor(f, fun), resp, into(f, fun)), resp
    catch e
        rethrow(PostProcessorError(resp, e, stacktrace(catch_backtrace())))
    end
end

function rate_limit(f::Forge, fun::Function)
    if has_rate_limits(f, fun) && rate_limit_check(f, fun)
        orl = on_rate_limit(f, fun)
        if orl === ORL_THROW
            throw(RateLimitedError(rate_limit_period(f, fun)))
        elseif orl === ORL_WAIT
            rate_limit_wait(f, fun)
        else
            @warn "Ignoring unknown rate limit behaviour $orl"
        end
    end
end


text:
#   This file is part of SyntaxTree.jl. It is licensed under the MIT license
#   Copyright (C) 2018 Michael Reed

"""
    expravg(expr)

Recursively obtains the *number of scalars* in an expression, the *average of those scalars* logarithms, the *number of non-trivial exponents* in the expression, and the *average of the exponents*.
"""
function expravg(expr)
    cs = 0
    s = 0.0
    cp = 0
    p = 0.0
    if typeof(expr) == Expr
        if expr.head == :call && expr.args[1] == :^ &&
            expr.args[3] |> typeof <: Number
            cp += 1
            p  += abs(expr.args[3])
            (cst,st,cpt,pt) = expravg(expr.args[2])
            cs += cst
            s  += cst*st
            cp += cpt
            p  += cpt*pt
        else
            for arg ∈ expr.args
                (cst,st,cpt,pt) = expravg(arg)
                cs += cst
                s  += cst*st
                cp += cpt
                p  += cpt*pt
            end
        end
    elseif typeof(expr) <: Number
        cs += 1
        s  += log(abs(expr))
    end
    return (cs,0 ∈ [s,cs] ? 1.0 : s/cs,cp,cp == 0 ? 1.0 : p/cp)
end

"""
    exprdev(expr)

Returns the standard deviation of the logarithm of the scalars in an expression.
"""
function exprdev(expr,val=expravg(expr)[2],cal=callcount(expr))
    s = 0.0
    if typeof(expr) == Expr
        for arg ∈ expr.args
            s += exprdev(arg,val,cal)
        end
    elseif typeof(expr) <: Number
        s  += ((log(abs(expr))-val)^2)/(cal-1)
    end
    return s
end

"""
    exprval(expr)

Returns the expression value and other characteristics as defined in "Optimal polynomial characteristic methods" by Michael Reed in 2018. This value can be used to order polynomial forms, with lower values being more optimal and efficient for computation.

The **first** value is the *expression value*, the **second** is the *number of operations*, the **third** is the *standard deviation* of the logarithms of the scalars, the **fourth** is the *average* of the logarithms of the scalars, and the **fifth** is the *average of the exponents*.
"""
function exprval(expr)
    val = expravg(expr)
    cal = callcount(expr)
    mal = sqrt(exprdev(expr,val[2],cal))
    cal*sqrt(abs(val[2])*mal)*val[4], cal, mal, val[2], val[4]
end

ext: jl


text:
abstract type AbstractIASCG <: AbstractController end

mutable struct IASCG <: AbstractIASCG
    p::Int
    x_ref::Float64
    h_ref::Float64
    h_des::Float64
    Approx_mode::Int
    n::Float64
    x₀::Float64
    V₀::Float64
    Δx_pred::Float64
    x_f::Float64
    h_f_d::Float64
    γ_f_d::Float64
    V_f_d::Float64
    env_model
end

function Path_Generator(s_guidance::IASCG, R, x_t, h_t, γ_t)
    @unpack p, x_ref, h_ref, h_des, x_f, h_f_d, γ_f_d = s_guidance

    order = p-1:-1:0
    ϕ_t  = (x_t / x_ref).^order
    ϕ_f  = (x_f / x_ref).^order
    dϕ_t = order .* [ϕ_t[2:end]; 0]
    dϕ_f = order .* [ϕ_f[2:end]; 0]

    A = [ϕ_t'; ϕ_f'; dϕ_t'; dϕ_f']
    b = [h_t; h_f_d; x_ref * tan(γ_t); x_ref * tan(γ_f_d)] / h_ref

    C = zeros(p, p)
    for i in 1:p
        for j in 1:p
            C[i,j] = 1 / (2 * p - (i + j - 1))
        end
    end
    S = (x_f * ϕ_f * ϕ_f'  - x_t * ϕ_t * ϕ_t') / x_ref .* C
    H = [zeros(p-2, 2) diagm(order[1:end-2] .* (order[1:end-2] .- 1));
         zeros(2,p)]
    P = S + (10.0)^R * H * S * H'
    P = Symmetric((P + P') / 2)
    q = -2 * h_des / h_ref * S[:,end]

    c = Convex.Variable(p)
    problem = minimize(quadform(c, P; assume_psd = true) + q' * c, A * c == b)
    Convex.solve!(problem, Mosek.Optimizer; silent_solver = true)
    # c_opt = evaluate(c)

    return problem.status, problem.optval, c.value
end

function Speed_Dynamics(s_guidance::IASCG, c)
    @unpack p, x_ref, h_ref, n, x₀, V₀, x_f, V_f_d, Approx_mode, env_model = s_guidance
    @unpack m, S, C_D_0, K, ρ₀, H, g = env_model
    return function (V, ps, x)
        order = p-1:-1:0
        ϕ  = (x / x_ref).^order
        ϕ′ = order .* [ϕ[2:end]; 0] / x_ref
        ϕ″ = order .* (order .- 1) .* [ϕ[3:end]; zeros(2)] / x_ref^2
        (h, h′, h″) = h_ref * [ϕ'; ϕ′'; ϕ″'] * c
        # γ  = atan(h′)
        # γ′ = h″ / (1 + h′^2)

        ρ  = ρ₀ * exp(-h / H)
        if Approx_mode == 0
            V̂_guess = V
        elseif Approx_mode == 1
            V̂_guess = V_f_d - (V_f_d - V₀) * ((x_f - x) / (x_f - x₀))^n
        elseif Approx_mode == 2
            V̂_guess = Inf
        end
        C_L = 2 * m / (ρ * S * sqrt(1 + h′^2)) * (h″ / (1 + h′^2) + g / V̂_guess^2)
        C_D = C_D_0 + K * C_L^2

        dV = - ρ * S * C_D * sqrt(1 + h′^2) / ( 2 * m ) * V - g / V * h′

        return dV
    end
end

function Speed_Predictor(s_guidance::IASCG, x, V, c)
    @unpack x_f = s_guidance
    # @unpack m, S, C_D_0, K, ρ₀, H, g = env_model

    # ----------------------------------------------------------------------------
    # Numerical Integration - Solving ODE
    # ----------------------------------------------------------------------------
    prob = ODEProblem(Speed_Dynamics(s_guidance, c), V, (x, x_f))
    sol = OrdinaryDiffEq.solve(prob, Tsit5(), reltol = 1e-6, abstol = 1e-6, dense = false, save_on = false, save_everystep = false, save_start = false, save_end = true)
    V̂_f = sol[1]

    # ----------------------------------------------------------------------------
    # Numerical Integration - Calculation of Definite Integrals
    # ----------------------------------------------------------------------------
    # x_s = x:Δx_pred:x_f

    # order = p-1:-1:0
    # path_deriv_s = x_s |> Map(x -> begin
    #                         ϕ  = (x / x_ref).^order
    #                         ϕ′ = order .* [ϕ[2:end]; 0] / x_ref
    #                         ϕ″ = order .* (order .- 1) .* [ϕ[3:end]; zeros(2)] / x_ref^2
    #                         path_deriv = h_ref * [ϕ'; ϕ′'; ϕ″'] * c
    #                     end) |> tcollect
    # path_deriv_s = hcat(path_deriv_s...)'
    # (h_s, h′_s, h″_s) = path_deriv_s[[1,2,3],:]

    # ρ_s = ρ₀ * exp.( - h_s / H)

    # a = -2 * (γ′_s .+ g ./ V̂_guess_s.^2) .* (C_D_0 .+ K .* C_L_s.^2) ./ C_L_s
    # b = -g * tan.(γ_s)
    # Φ = exp(quadgk(a, x, x_f, rtol=1e-3))

    # V̂_f = sqrt( Φ[1] * V^2 + 2 * ∫Φb)

    return V̂_f
end

function zeroprob_IASCG_Planner(s_guidance::IASCG, x, h, γ, V)
    @unpack V_f_d = s_guidance
    return function (R)
        e_v_f = Speed_Predictor(s_guidance, x, V, Path_Generator(s_guidance, R, x, h, γ)[3]) - V_f_d

        return e_v_f
    end
end

function minbndprob_IASCG_Planner(s_guidance::IASCG, x, h, γ, V)
    @unpack V_f_d = s_guidance
    return function (R)
        e_v_f = Speed_Predictor(s_guidance, x, V, Path_Generator(s_guidance, R, x, h, γ)[3]) - V_f_d

        return e_v_f^2
    end
end

function IASCG_Planner(s_guidance::IASCG, s, params)
    @unpack p, V_f_d   = s_guidance
    @unpack x, h, γ, V = s
    @unpack R, c       = params

    # previous parameters -------------------------------------------
    # R_opt = R
    # c_opt = c

    # Root finding --------------------------------------------------
    # if c == zeros(p)
    #     R_opt = find_zero(zeroprob_IASCG_Planner(s_guidance, x, h, γ, V), (-8.0, 5.0), Bisection())
    # else
    #     R_opt = find_zero(zeroprob_IASCG_Planner(s_guidance, x, h, γ, V), R, Order1())
    # end

    # Minimisation over bounded interval ----------------------------
    res     = Optim.optimize(minbndprob_IASCG_Planner(s_guidance, x, h, γ, V), -8.0, 5.0)
    R_opt   = Optim.minimizer(res)

    c_opt = Path_Generator(s_guidance, R_opt, x, h, γ)[3]
    # V̂_f   = Speed_Predictor(s_guidance, x, V, c_opt)

    return ComponentArray(R = R_opt, c = vec(c_opt))
end

function IASCG_cmd(s_guidance::IASCG)
    @unpack p, x_ref, h_ref, env_model = s_guidance
    @unpack g = env_model
    return function (s, params, t)
        @unpack x, γ, V = s
        @unpack c = params

        order = p-1:-1:0
        ϕ  = (x / x_ref).^order
        # ϕ′ = order .* [ϕ[2:end]; 0] / x_ref
        ϕ″ = order .* (order .- 1) .* [ϕ[3:end]; zeros(2)] / x_ref^2
        h″ = h_ref * ϕ″' * c
        A_L_cmd = ( h″ * (cos(γ))^2 * V^2 + g ) * cos(γ)

        return A_L_cmd
    end
end


text:
using Test
using SemanticModels.ModelTools.PetriModels

# SIR  <- IR  -> SEIR
#  |       |      |
#  v       v      v
# SIRS <- IRS -> SEIRS

# +
using Catlab.WiringDiagrams
using Catlab.Doctrines
import Catlab.Doctrines.⊗
import Catlab.Graphics: to_graphviz
import Catlab.Graphics.Graphviz: run_graphviz
⊗(a::WiringDiagram, b::WiringDiagram) = otimes(a,b)
import Base: ∘
∘(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
⊚(a,b) = b ∘ a
S, E, I, R, D= Ob(FreeSymmetricMonoidalCategory, :S, :E, :I, :R, :D)

inf  = to_wiring_diagram(Hom(:infection, S ⊗ I, I⊗I))
expo = to_wiring_diagram(Hom(:exposure, S⊗I, E⊗I))
rec  = to_wiring_diagram(Hom(:recovery, I,   R))
wan  = to_wiring_diagram(Hom(:waning,   R,   S))

sir_wire  = inf ⊚ (rec ⊗ rec)

sir = model(PetriModel, sir_wire)

dump(sir)

# +
@variables S, E, I, R

dump(S)
# -

sir = model(PetriModel, Petri.Model([S, I, R],
                 [(I, R), (S+I, 2I)],
                 ))

ir = model(PetriModel, Petri.Model([S, I, R],
                 [(I, R)]))

seir = model(PetriModel, Petri.Model([S, I, R],
                 [(I, R), (S+I, I+E), (E, I)],
                 ))

# +
irs = model(PetriModel, Petri.Model([S, I, R],
                 [(I, R), (R, S)],
                 ))

dump(sir)
dump(irs)
# -

sirs = model(PetriModel, Petri.Model([S, I, R],
                 [(I, R), (S+I, 2*I), (R, S)],
                 ))

rule = PetriModels.Span(sir, ir, seir)
sirs′ = PetriModels.pushout(irs, sir)
@test sirs′.model.Δ == sirs.model.Δ
seirs = PetriModels.solve(PetriModels.DPOProblem(rule, irs))
@test all(Set(seirs.model.Δ) .== Set([(S+I, I+E),
                             (E, I),
                             (I, R),
                             (R, S)]))
l = sir
c = ir
r = seir
c′ = irs

l′ = PetriModels.pushout(l, c′)
@test l′.model.Δ == sirs.model.Δ
@test PetriModels.dropdown(l,c,l′).model.Δ == c′.model.Δ
@test PetriModels.pushout(r, c′).model.Δ == seirs.model.Δ

@show "SIR"
Δ(l.model, :state)
# @show "IR"
# funckit(c, :state)
@show "SEIR"
# funckit(r, :state)

@show "SIRS"
# funckit(l′, :state)
# @show "IRS"
# funckit(c′, :state)
@show "SEIRS"
# funckit(seirs, :state)

exprs = Δ(sirs.model, :state)
m = Petri.Model([S, I, R], exprs, [
    quote
    λ_2(state) = state.γ * state.I
    end,
    quote
    λ_1(state) = state.β * state.S * state.I / +(state.S, state.I, state.R)
    end,
    quote
    λ_3(state) = state.μ * state.R
    end],
                [
                    quote b_2(state) = state.I > 0 end,
                    quote b_1(state) = state.S > 0 && state.I > 0 end,
                    quote b_3(state) = state.R > 0 end]
                )

# +
mutable struct SIRState{T,F}
    S::T
    I::T
    R::T
    β::F
    γ::F
    μ::F
end

mutable struct ParamSIR{T, P}
    S::T
    I::T
    R::T
    params::P
end
# -

p = Petri.Problem(eval(m), SIRState(100, 1, 0, 0.5, 0.15, 0.05), 150)
#@show Petri.solve(p)

# @show Petri.funckit(Petri.Problem(l, missing, 10), :state)
function test_1()
    no_transitions = Tuple{Operation, Operation}[]
    @variables A, B, C, D
    states = [A, B, C, D]
    l = model(PetriModel, Petri.Model(states, [(A, B)]))
    c = model(PetriModel, Petri.Model(states, no_transitions))
    r = model(PetriModel, Petri.Model(states, [(A, B + C)]))
    rule = PetriModels.Span(l, c, r)
    c′ = model(PetriModel, Petri.Model(states, [(B, A)]))
    r′ = PetriModels.solve(PetriModels.DPOProblem(rule, c′))
    @test r′.model.Δ == [(A, B+C), (B, A)]

    l′ = PetriModels.pushout(l, c′)
    @test l′.model.Δ == [(A, B), (B, A)]
    @test PetriModels.dropdown(l,c,l′).model.Δ == [(B, A)]
end

function test_2()
    no_transitions = Tuple{Operation, Operation}[]
    @variables A, B, C, D
    states = [A, B, C, D]
    l = model(PetriModel, Petri.Model(states, [(A, B), (B,C)]))
    c = model(PetriModel, Petri.Model(states, no_transitions))
    r = model(PetriModel, Petri.Model(states, [(A, C)]))
    rule = PetriModels.Span(l, c, r)
    c′ = model(PetriModel, Petri.Model(states, [(C, D)]))
    r′ = PetriModels.solve(PetriModels.DPOProblem(rule, c′))
    @test r′.model.Δ == [(A, C), (C, D)]

    l′ = PetriModels.pushout(l, c′)
    @test l′.model.Δ == [(A, B), (B, C), (C, D)]
    @test PetriModels.dropdown(l,c,l′).model.Δ == [(C, D)]
    @test PetriModels.pushout(rule.r, c′).model.Δ == [(A, C), (C, D)]
end
test_1()
test_2()

exprs2 = Δ(Petri.Model([S,I,R], [(2S+I, 3I)]))

# @code_native m′.Δ[1](ParamSIR(100, 1, 0, [ 0.15, 0.55/101, 0.15, 0.1 ]))


text:
module GeoPhyInv

# required packages
using ParallelStencil
using Distances
using TimerOutputs
using LinearMaps
using Ipopt
using Optim, LineSearches
using DistributedArrays
using Calculus
using ProgressMeter
using Distributed
using SharedArrays
using Printf
using DataFrames
using SparseArrays
using Interpolations
using OrderedCollections
using CSV
using Statistics
using LinearAlgebra
using Random
using ImageFiltering
using NamedArrays
using Test
using AxisArrays
using Distributions
using StatsBase
using InteractiveUtils
using RecipesBase
using ColorSchemes
using FFTW
using HDF5
using SpecialFunctions
using DSP
using CUDA
using CUDA.CUSPARSE
CUDA.allowscalar(false)

include("params.jl")

# This type is extensively used to create named arrays with Symbols
NamedStack{T} =
    NamedArray{T,1,Array{T,1},Tuple{OrderedCollections.OrderedDict{Symbol,Int64}}}

# a module with some util methods
include("Utils/Utils.jl")
include("Interpolation/Interpolation.jl")

# some structs used for multiple dispatch throughout this package
include("types.jl")
export Srcs, Recs, SSrcs
export FdtdAcoustic, FdtdElastic, FdtdAcousticBorn, AcousticBorn, ElasticBorn, FdtdAcousticVisco

# mutable data type for seismic medium + related methods
include("media/core.jl")
export Medium

# source-receiver geometry
include("ageom/core.jl")
export AGeom, AGeomss

#
include("database/core.jl")

# store records
include("records/core.jl")
export Records

# mutable data type for storing time-domain source wavelets
include("srcwav/wavelets.jl")
export ricker, ormsby
include("srcwav/core.jl")
export SrcWav

"""
Initialize the package with ParallelStencil, giving access to its main functionality.
```julia
@init_parallel_stencil(ndims, use_gpu, datatype, order)
```

# Arguments

* `ndims::Int`: the number of dimensions used for the stencil computations 2D or 3D
* `use_gpu::Bool` : use GPU for stencil computations or not
* `datatype`: the type of numbers used by field arrays (e.g. Float32 or Float64)
* `order::Int ∈ [2,4,6,8]` : order of the finite-difference stencil
"""
macro init_parallel_stencil(ndims::Int, use_gpu::Bool, datatype, order)
    quote
        # using ParallelStencil
        # ParallelStencil.is_initialized() && ParallelStencil.@reset_parallel_stencil()
        # ParallelStencil.@init_parallel_stencil(CUDA, $datatype, $ndims)
        @eval GeoPhyInv begin
            _fd.ndims = $ndims
            @assert _fd.ndims ∈ [2, 3]
            _fd.order = $order
            @assert _fd.order ∈ [2, 4, 6, 8]
            _fd.use_gpu = $use_gpu
            _fd.datatype = $datatype
            @assert _fd.datatype ∈ [Float32, Float64]
            # extend by more points as we increase order, not sure if that is necessary!!!
            _fd.npml = 20 + ($order - 1)
            _fd.npextend = 20 + ($order - 1) # determines exmedium

            ParallelStencil.is_initialized() && ParallelStencil.@reset_parallel_stencil()
            # check whether 2D or 3D, and initialize ParallelStencils accordingly
            @static if $use_gpu
                # using CUDA
                ParallelStencil.@init_parallel_stencil(CUDA, $datatype, $ndims)
            else
                ParallelStencil.@init_parallel_stencil(Threads, $datatype, $ndims)
            end
            include(
                joinpath(
                    dirname(pathof(GeoPhyInv)),
                    "fdtd",
                    string("diff", $ndims, "D.jl"),
                ),
            )
            # define structs for wavefields in 2D/3D
            include(joinpath(dirname(pathof(GeoPhyInv)), "fields.jl"))
            export Fields
            include(joinpath(dirname(pathof(GeoPhyInv)), "fdtd", "core.jl"))
            include(joinpath(dirname(pathof(GeoPhyInv)), "born", "core.jl"))
            nothing
        end
    end
end
export @init_parallel_stencil
export SeisForwExpt

# export update from GeoPhyInv, as it is commonly used
export update, update!

# # include modules (note: due to dependencies, order is important!)
# include("Operators.jl")
# include("Smooth.jl")
# include("IO.jl")

# include("Coupling.jl")

# include("fwi/fwi.jl")

# include("Poisson/Poisson.jl")
include("plots.jl")

#     padarray,
#     padarray!,
#     SeisInvExpt,
#     LS,
#     LS_prior,
#     Migr,
#     Migr_FD

# # export the Expt for Poisson
# const PoissonExpt=GeoPhyInv.Poisson.ParamExpt
# export PoissonExpt
# mod!(a::PoissonExpt, b, c) = GeoPhyInv.Poisson.mod!(a, b, c)
# mod!(a::PoissonExpt, b) = GeoPhyInv.Poisson.mod!(a, b)
# mod!(a::PoissonExpt) = GeoPhyInv.Poisson.mod!(a)
# operator_Born(a::PoissonExpt, b) = GeoPhyInv.Poisson.operator_Born(a, b)

end # module


text:
# regular prior for Pose3

export PriorPose3, PackedPriorPose3

"""
$(TYPEDEF)

Direct observation information of `Pose3` variable type.
"""
struct PriorPose3{T <: IIF.SamplableBelief, P} <: IncrementalInference.AbstractPrior
  Z::T
  p::P
  # empty constructor
  # PriorPose3{T}() where T = new{T}()
  # regular constructor
  # PriorPose3{T}(z::T) where {T <: IIF.SamplableBelief} = new{T}(z)
end
# convenience and default object helper
PriorPose3(z::IIF.SamplableBelief) = PriorPose3(z, getPointIdentity(Pose3))
PriorPose3{T,P}(z::IIF.SamplableBelief) where {T,P} = PriorPose3{T,P}(z, getPointIdentity(Pose3))
PriorPose3() = PriorPose3(MvNormal(zeros(6), LinearAlgebra.diagm([0.01*ones(3);0.0001*ones(3)])))

DFG.getManifold(::PriorPose3) = SpecialEuclidean(3)

# Standardized sampling function
function getSample(cf::CalcFactor{<:PriorPose3}, N::Int=1)
  Z = cf.factor.Z
  p = cf.factor.p
  M = getManifold(cf.factor)

  Xc = rand(Z)

  # X = get_vector.(Ref(M), Ref(p), Xc, Ref(DefaultOrthogonalBasis()))
  X = hat(M, p, Xc)
  points = exp(M, p, X)

  return points
end

#FIXME Serialization
"""
$(TYPEDEF)

Serialization type for PriorPose3.
"""
mutable struct PackedPriorPose3  <: IncrementalInference.PackedInferenceType
    Zi::String
    PackedPriorPose3() = new()
    PackedPriorPose3(x::AbstractString) = new(x)
end
function convert(::Type{PriorPose3}, packed::PackedPriorPose3)
  return PriorPose3( convert(SamplableBelief, packed.Zi) )
end
function convert(::Type{PackedPriorPose3}, obj::PriorPose3)
  return PackedPriorPose3(convert(PackedSamplableBelief, obj.Z))
end


text:
function check_div(chn; per_chain=false)
    divergent = Array(chn[:numerical_error])

    n_for_chains = sum(divergent, dims=1)

    n_divergent = sum(n_for_chains)
    if n_divergent > 0
        N = length(divergent)
        @info "Warning: $(n_divergent) of $(N) iterations ended with a divergence ($(100 * n_divergent / N) %)"
        if per_chain
            chain_len, num_chains = size(divergent)
            for chain_num in 1:num_chains
                if n_for_chains[chain_num] > 0
                    @info "Chain $(chain_num): $(n_for_chains[chain_num]) of $(chain_len) iterations ended with a divergence ($(100 * n_for_chains[chain_num] / chain_len) %)"
                end
            end
        end

        @info "Try running with larger target acceptance rate to remove the divergences."
        return false
    else
        return true
    end
end

function check_n_eff(chn)
    no_warning = true

    n_steps = size(chn, 1)
    ess_df = DataFrame(Turing.ess(chn))
    for row in eachrow(ess_df)
        param_name = string(row[:parameters])
        n_eff = row[:ess]
        ratio = n_eff / n_steps
        if ratio < 0.1 || isnan(ratio) || isinf(ratio)
            no_warning = false
            @info "n_eff / n_steps for parameter $(param_name) is $(ratio)"
        end
    end

    if no_warning
        @info "n_eff / n_iter looks reasonable for all parameters"
    else
        @info "n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated"
    end

    return no_warning
end

function check_rhat(chn)
    no_warning = true

    ess_df = DataFrame(Turing.ess(chn))
    for row in eachrow(ess_df)
        param_name = string(row[:parameters])
        rhat = row[:rhat]
        if isnan(rhat) || isinf(rhat) || rhat > 1.1 || rhat < 0.9
            no_warning = false
            @info "Rhat for parameter $(param_name) is $(rhat)"
        end
    end

    if no_warning
        @info "Rhat looks reasonable for all parameters"
    end

    return no_warning
end

function check_diagnostics(chn)
    div = check_div(chn, per_chain=true)
    n_eff = check_n_eff(chn)
    rhat = check_rhat(chn)
    return div && n_eff && rhat
end


text:
using FITSIO
using DataFrames
using CSV

crossref_tab = FITS("kepler_dr2_1arcsec.fits")[2];

#gaia_full = DataFrames.readtable("/sagan1/kepler/gaia_kepler_dr25_noblanklines.csv",separator=';',skipstart=243,nastrings=["NOT_AVAILABLE"]);
#deleterows!(gaia_full,1:2)

gaia_full = CSV.read("/sagan1/kepler/gaia_kepler_dr25_noblanklines.csv",delim=';',header=244,datarow=247,missingstring="NOT_AVAILABLE", rows_for_type_detect=1000)

(designation, source_id, gaia_ref_epoch, ra, ra_error, dec, dec_error, parallax, parallax_error, parallax_over_error, pmra, pmra_error, pmdec, pmdec_error, ra_dec_corr, ra_parallax_corr, ra_pmra_corr, ra_pmdec_corr, dec_parallax_corr, dec_pmra_corr, dec_pmdec_corr, parallax_pmra_corr, parallax_pmdec_corr, pmra_pmdec_corr, astrometric_chi2_al, astrometric_excess_noise, astrometric_excess_noise_sig, astrometric_primary_flag, duplicated_source, phot_g_mean_flux, phot_g_mean_flux_error, phot_g_mean_mag, phot_bp_mean_flux, phot_bp_mean_flux_error, phot_bp_mean_mag, phot_rp_mean_flux, phot_rp_mean_flux_error, phot_rp_mean_mag, bp_rp, bp_g, g_rp, radial_velocity, radial_velocity_error, phot_variable_flag, l, b, ecl_lon, ecl_lat, teff_val, teff_percentile_lower, teff_percentile_upper, a_g_val, a_g_percentile_lower, a_g_percentile_upper, e_bp_min_rp_val, e_bp_min_rp_percentile_lower, e_bp_min_rp_percentile_upper, radius_val, radius_percentile_lower, radius_percentile_upper, lum_val, lum_percentile_lower, lum_percentile_upper, r_est, r_lo, r_hi, r_length_prior, r_result_flag, r_modality_flag, kepid, tm_designation, ra_kepler, dec_kepler, kepmag, teff, teff_err1, teff_err2, teff_prov, logg, logg_err1, logg_err2, logg_prov, feh, feh_err1, feh_err2, feh_prov, radius, radius_err1, radius_err2, mass, mass_err1, mass_err2, prov_sec, nconfp, nkoi, ntce, jmag, hmag, kmag, planet_exist, kepler_gaia_ang_dist) = map(x->read(crossref_tab,x), ["designation", "source_id", "gaia_ref_epoch", "ra", "ra_error", "dec", "dec_error", "parallax", "parallax_error", "parallax_over_error", "pmra", "pmra_error", "pmdec", "pmdec_error", "ra_dec_corr", "ra_parallax_corr", "ra_pmra_corr", "ra_pmdec_corr", "dec_parallax_corr", "dec_pmra_corr", "dec_pmdec_corr", "parallax_pmra_corr", "parallax_pmdec_corr", "pmra_pmdec_corr", "astrometric_chi2_al", "astrometric_excess_noise", "astrometric_excess_noise_sig", "astrometric_primary_flag", "duplicated_source", "phot_g_mean_flux", "phot_g_mean_flux_error", "phot_g_mean_mag", "phot_bp_mean_flux", "phot_bp_mean_flux_error", "phot_bp_mean_mag", "phot_rp_mean_flux", "phot_rp_mean_flux_error", "phot_rp_mean_mag", "bp_rp", "bp_g", "g_rp", "radial_velocity", "radial_velocity_error", "phot_variable_flag", "l", "b", "ecl_lon", "ecl_lat", "teff_val", "teff_percentile_lower", "teff_percentile_upper", "a_g_val", "a_g_percentile_lower", "a_g_percentile_upper", "e_bp_min_rp_val", "e_bp_min_rp_percentile_lower", "e_bp_min_rp_percentile_upper", "radius_val", "radius_percentile_lower", "radius_percentile_upper", "lum_val", "lum_percentile_lower", "lum_percentile_upper", "r_est", "r_lo", "r_hi", "r_length_prior", "r_result_flag", "r_modality_flag", "kepid", "tm_designation", "ra_kepler", "dec_kepler", "kepmag", "teff", "teff_err1", "teff_err2", "teff_prov", "logg", "logg_err1", "logg_err2", "logg_prov", "feh", "feh_err1", "feh_err2", "feh_prov", "radius", "radius_err1", "radius_err2", "mass", "mass_err1", "mass_err2", "prov_sec", "nconfp", "nkoi", "ntce", "jmag", "hmag", "kmag", "planet?", "kepler_gaia_ang_dist"])

crossref_df = DataFrame(designation=designation, source_id=source_id, gaia_ref_epoch=gaia_ref_epoch, ra=ra, ra_error=ra_error, dec=dec, dec_error=dec_error, parallax=parallax, parallax_error=parallax_error, parallax_over_error=parallax_over_error, pmra=pmra, pmra_error=pmra_error, pmdec=pmdec, pmdec_error=pmdec_error, ra_dec_corr=ra_dec_corr, ra_parallax_corr=ra_parallax_corr, ra_pmra_corr=ra_pmra_corr, ra_pmdec_corr=ra_pmdec_corr, dec_parallax_corr=dec_parallax_corr, dec_pmra_corr=dec_pmra_corr, dec_pmdec_corr=dec_pmdec_corr, parallax_pmra_corr=parallax_pmra_corr, parallax_pmdec_corr=parallax_pmdec_corr, pmra_pmdec_corr=pmra_pmdec_corr)

dr2name = collect(skipmissing(gaia_full[:Source]))
astrometric_gof_al = collect(skipmissing(gaia_full[:gofAL]))
priam_flags = gaia_full[:fPriam]

crossref_ind = map(x->findfirst(dr2name,x),crossref_df[:source_id])

tmp_gof = Vector{Float64}(length(crossref_ind))
tmp_priam = Vector{Int64}(length(crossref_ind))

for i in 1:length(crossref_ind)
    tmp_gof[i] = astrometric_gof_al[crossref_ind[i]]
    if typeof(priam_flags[crossref_ind[i]]) == Missings.Missing
        tmp_priam[i] = -999
    else
        tmp_priam[i] = priam_flags[crossref_ind[i]]
    end
end

crossref_df[:astrometric_gof_al] = tmp_gof
#crossref_df[:astrometric_gof_al] = recode(crossref_df[:astrometric_gof_al], -999.0=>missing)

new_col = ["astrometric_chi2_al", "astrometric_excess_noise", "astrometric_excess_noise_sig", "astrometric_primary_flag", "duplicated_source", "phot_g_mean_flux", "phot_g_mean_flux_error", "phot_g_mean_mag", "phot_bp_mean_flux", "phot_bp_mean_flux_error", "phot_bp_mean_mag", "phot_rp_mean_flux", "phot_rp_mean_flux_error", "phot_rp_mean_mag", "bp_rp", "bp_g", "g_rp", "radial_velocity", "radial_velocity_error", "phot_variable_flag", "l", "b", "ecl_lon", "ecl_lat"]

tmp_col_names = vcat(names(crossref_df),[Symbol(i) for i in new_col])

crossref_df = hcat(crossref_df, astrometric_chi2_al, astrometric_excess_noise, astrometric_excess_noise_sig, astrometric_primary_flag, duplicated_source, phot_g_mean_flux, phot_g_mean_flux_error, phot_g_mean_mag, phot_bp_mean_flux, phot_bp_mean_flux_error, phot_bp_mean_mag, phot_rp_mean_flux, phot_rp_mean_flux_error, phot_rp_mean_mag, bp_rp, bp_g, g_rp, radial_velocity, radial_velocity_error, phot_variable_flag, l, b, ecl_lon, ecl_lat)

names!(crossref_df, tmp_col_names)

crossref_df[:priam_flags] = tmp_priam
crossref_df[:priam_flags] = recode(crossref_df[:priam_flags], -999=>missing)

new_col = ["teff_val", "teff_percentile_lower", "teff_percentile_upper", "a_g_val", "a_g_percentile_lower", "a_g_percentile_upper", "e_bp_min_rp_val", "e_bp_min_rp_percentile_lower", "e_bp_min_rp_percentile_upper", "radius_val", "radius_percentile_lower", "radius_percentile_upper", "lum_val", "lum_percentile_lower", "lum_percentile_upper", "r_est", "r_lo", "r_hi", "r_length_prior", "r_result_flag", "r_modality_flag", "kepid", "tm_designation", "ra_kepler", "dec_kepler", "kepmag", "teff", "teff_err1", "teff_err2", "teff_prov", "logg", "logg_err1", "logg_err2", "logg_prov", "feh", "feh_err1", "feh_err2", "feh_prov", "radius", "radius_err1", "radius_err2", "mass", "mass_err1", "mass_err2", "prov_sec", "nconfp", "nkoi", "ntce", "jmag", "hmag", "kmag", "planet_cand", "kepler_gaia_ang_dist"]

tmp_col_names = vcat(names(crossref_df),[Symbol(i) for i in new_col])

crossref_df = hcat(crossref_df, teff_val, teff_percentile_lower, teff_percentile_upper, a_g_val, a_g_percentile_lower, a_g_percentile_upper, e_bp_min_rp_val, e_bp_min_rp_percentile_lower, e_bp_min_rp_percentile_upper, radius_val, radius_percentile_lower, radius_percentile_upper, lum_val, lum_percentile_lower, lum_percentile_upper, r_est, r_lo, r_hi, r_length_prior, r_result_flag, r_modality_flag, kepid, tm_designation, ra_kepler, dec_kepler, kepmag, teff, teff_err1, teff_err2, teff_prov, logg, logg_err1, logg_err2, logg_prov, feh, feh_err1, feh_err2, feh_prov, radius, radius_err1, radius_err2, mass, mass_err1, mass_err2, prov_sec, nconfp, nkoi, ntce, jmag, hmag, kmag, planet_exist, kepler_gaia_ang_dist)

names!(crossref_df, tmp_col_names)

CSV.write("gaiadr2_keplerdr25_crossref.csv", crossref_df)


text:
struct ClusterTopKBatchQs <: BatchPQs
    model::SVDD.OCClassifier
    sequentialQs::SequentialPQs
    k::Int

    function ClusterTopKBatchQs(model::SVDD.OCClassifier, sequential_strategy::SequentialPQs; k::Int)::ClusterTopKBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))
        return new(model, sequential_strategy, k)
    end
end

"""
Ensures selection of an informative, representative and diverse batch
by clustering the most informative observations and then using the k
cluster centers as batch query.

Number of most informative observations selected = 10*batch_size,
Number of clusters computed = batch_size
"""
function select_batch(qs::ClusterTopKBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    m = min(num_observations, 10 * qs.k)

    candidate_scores = qs_score(qs.sequentialQs, x, labels)[candidate_indices]
    descending_indices = sortperm(candidate_scores; rev=true)
    best_m = x[:, candidate_indices[descending_indices[1:m]]]

    distances = Distances.pairwise(Distances.Euclidean(), best_m, dims=2)
    clustering = kmedoids(distances, qs.k)

    medoid_indices = clustering.medoids

    return return candidate_indices[descending_indices[1:m]][medoid_indices]
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ f252364e-f2e3-11ea-09b5-ff1a685abd56
using Pkg, DrWatson

# ╔═╡ c126b36a-f2e3-11ea-1b0f-a7d38ecffcd8
md"## Clip-03-01t.jl"

# ╔═╡ f2526f60-f2e3-11ea-2d32-276d0b91d917
@quickactivate "StatisticalRethinkingTuring"

# ╔═╡ f252ea76-f2e3-11ea-093a-01c585ee2ea1
md"### snippet 3.1"

# ╔═╡ f25a9eec-f2e3-11ea-1611-9129817fe7a6
begin
	Pr_Positive_Vampire = 0.95
	Pr_Positive_Mortal = 0.01
	Pr_Vampire = 0.001
	Pr_Positive = Pr_Positive_Vampire * Pr_Vampire + Pr_Positive_Mortal * (1 - Pr_Vampire)
	Pr_Vampire_Positive = Pr_Positive_Vampire * Pr_Vampire / Pr_Positive
	Pr_Vampire_Positive
end

# ╔═╡ f25b24a2-f2e3-11ea-302d-f3e7f49f645b
md"## End of clip-03-01t.jl"

# ╔═╡ Cell order:
# ╟─c126b36a-f2e3-11ea-1b0f-a7d38ecffcd8
# ╠═f252364e-f2e3-11ea-09b5-ff1a685abd56
# ╠═f2526f60-f2e3-11ea-2d32-276d0b91d917
# ╟─f252ea76-f2e3-11ea-093a-01c585ee2ea1
# ╠═f25a9eec-f2e3-11ea-1611-9129817fe7a6
# ╟─f25b24a2-f2e3-11ea-302d-f3e7f49f645b


text:
using PyPlot, ElasticWave

nz = 300; nx = 300; # model size
ext= 20; # thickness of PML layers
iflag = 1; #when iflag=1, four side obsorbing boundary, iflag=2, free surface
dx = 5. ; dz = 5. ; dt = 4e-4; # grid size
tmax = 1.5; f0=20.; # recording time and dominant frequency of source wavelet
vp=2500.*ones(nz,nx); vs=2500./sqrt(3)*ones(nz,nx); rho=2.5*ones(nz,nx);
vp[151:end,:]=3500.;  vs[151:end,:]=3500./sqrt(3.); #physical model
SeisPlot(vp, xlabel="X (km)", ylabel="Z (km)", vmin=1500, vmax=4500, cmap="gray", dx=0.005, dy=0.005);
tight_layout(); savefig("/Users/wenyue/Desktop/v.pdf")
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0); #finite difference matrix

# single source
isz = 5; isx = 150; ot = 0.0; #source location and starting time
flags = vec([false false true true false]); # add source to normal stress, simulate explosive sources
src = InitSource(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
src.Txx = src.Txx * 1e5; src.Tzz = src.Tzz * 1e5;

# # =========finite difference modeling==========
# irx = collect(1:nx); irz = 1*ones(Int64, length(irx)); #location of receivers
# shotv = MultiStepForward(irz, irx, src, fidMtx, tmax=tmax) #generate one shot
root = homedir(); path = join([root "/Desktop/wfd.bin"]);
MultiStepForward(path, src, fidMtx, tmax=tmax);

root = homedir(); path = join([root "/Desktop/decomposed.bin"]);
MultiStepForward(path, src, fidMtx, dz, dx; tmax=tmax)

root = homedir(); path = join([root "/Desktop/wfd.bin"]);
MultiStepForward(path, src, fidMtx, tmax=tmax);

(nz, nx, ext, iflag, dt, nt) = InfoPv(path);
pp = zeros(nt); ps = zeros(nt);
rz = irz[1] + ext; rx = irx[1] + ext; idx = (rx-1)*(nz+2*ext) + rz + ext;
for i = 1151 : nt
    (vxx, vxz, vzx, vzz) = ReadPv(path, i)
    pp[i] = vxx[idx] + vzz[idx]
    ps[i] = vxz[idx] - vzx[idx]
end

path1 = join([root "/Desktop/moive"]);
waveAnim(path1, path, cpt="Vz", cmap="gray", wbox=6, hbox=6, dt1=0.003, pclip=1.0, vlim=0.1)

irz = [1]; irx = [250];
shotv = binWfd2shotv(irz, irx, path)

obs = hcat(shotv.Vz, shotv.Vx); obs[1:1150,:] = 0.0;
spt = hcat(pp, ps);
SeisPlot(obs, style="wiggles", dy=0.0004, xcur=0.8, ylabel="Time (s)", xticks=[]);
tight_layout(); savefig("/Users/wenyue/Desktop/obs.pdf");

SeisPlot(spt, style="wiggles", dy=0.0004, xcur=0.8, ylabel="Time (s)", xticks=[]);
tight_layout(); savefig("/Users/wenyue/Desktop/spt.pdf");

# function MultiStepForward_old(irz::Array{Int64,1}, irx::Array{Int64,1}, src::Source, fidMtx::FidMtx; tmax=0.5)
#     nz  =  src.nz;  nx = src.nx    ;
#     ext = src.ext;  iflag=src.iflag; dt = src.dt;
#     stl = src.ot ;  stu = src.ot+(src.nt-1)*dt;
#     nt  = round(Int64, tmax/dt)+1
#     spt1 = InitSnapShot(nz, nx, ext, iflag, dt, 1)
#     spt2 = InitSnapShot(nz, nx, ext, iflag, dt, 2)
#     AddSource!(spt1, src)
#     shotv = InitShotV(src.isz, src.isx, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
#     Spt2ShotV!(shotv, spt1)
#     for it = 2 : nt
#         OneStepForward!(spt2, spt1, fidMtx)
#         if stl <= (it-1)*dt <= stu
#            AddSource!(spt2, src)
#         end
#         CopySnapShot!(spt1, spt2)
#         Spt2ShotV!(shotv, spt1)
#     end
#     return shotv
# end
# @time shotv1 = MultiStepForward_old(irz, irx, src, fidMtx, tmax=tmax)
# @time shotv1 = MultiStepForward_old(irz, irx, src, fidMtx, tmax=tmax)

# SeisPlot(shotv.Vx, pclip=90)
#
# root = homedir(); path = join([root "/Desktop/wfd.bin"]);
# MultiStepForward(path, src, fidMtx, tmax=tmax);
# shotv1 = binWfd2shotv(irz, irx, path);
#
# path1 = join([root "/Desktop/moive"]);
# waveAnim(path1, path, cpt="Vx")
# fid = open(path, "r")
# lwfd = nz*nx*sizeof(Float32)*5; pre = sizeof(Float32)*5;
# nt = floor(Int64, (filesize(fid)-pre)/lwfd)
# shot = zeros(nt, length(irx))
# for it = 1 : nt
#     pos = pre + (it-1)*lwfd + 0*sizeof(Float32)*nz*nx
#     seek(fid, pos)
#     vx = reshape(convert(Array{Float64}, read(fid, Float32, nz*nx)), nz, nx)
#     SeisPlot(vx)
#     for ir = 1 : length(irx)
#         shot[it, ir] = vx[irz[ir], irx[ir]]
#     end
# end
#
# it = 300;
# wfd = ReadWfd(path, it)
# vx  = reshape(wfd.Vx, nz, nx);
# SeisPlot(vx)
# tmp = 0.0
# for ir = 1 : length(irx)
#     tmp += (shotv1.Vx[it, ir]-vx[irz[ir],irx[ir]])^2
# end

# path = "/Users/wenlei/Desktop/test.bin"
# MultiStepForward(path, src, fidMtx, dz, dx)

# it = 200;
# (vxx, vxz, vzx, vzz) = ReadPv(path, it);
#
# wfd = ReadWfd(path1, it)
# (vxx1, vxz1, vzx1, vzz1) = Dxz(wfd);

# wfd = ReadWfd(path, 300);
# d = wfd.Txx
# imshow(d, cmap="PuOr", vmax=maximum(abs(d)), vmin=-maximum(abs(d)));

# irx = collect(1:5:nx); irz = 5*ones(Int64, length(irx));
# shot = MultiStepForward(irz, irx, src, fidMtx);

max_stars_repo_name: WneleiGao/ElasticWave


text:
using CSV
using Downloads: download
using DataFrames: DataFrame
using DimensionalData: DimArray
using Geotherm: Temp, Press, generate_trace
using Geotherm.Geometry: Point2D

@testset "Test `generate_trace`" begin
    file = download(
        "https://raw.githubusercontent.com/MineralsCloud/geothermpy/master/tests/data/example.csv",
    )
    df = CSV.File(file) |> DataFrame
    ts = 1000:5:4000
    ps = [
        20,
        23,
        26,
        29,
        32,
        35,
        38,
        41,
        44,
        47,
        50,
        53,
        56,
        59,
        62,
        65,
        68,
        71,
        74,
        77,
        80,
        83,
        86,
        89,
        92,
        95,
        98,
        101,
        104,
        107,
        110,
        113,
        116,
        119,
        122,
        125,
        128,
    ]
    gradient = DimArray(Matrix(df[:, 2:end]), (Temp(ts), Press(ps)))
    trace = generate_trace(gradient, Point2D(24, 1876.80005), 0.01, 1000)
    @test map(p -> getfield(p, :x), trace) == [
        24,
        24.01,
        24.020000000000003,
        24.030000000000005,
        24.040000000000006,
        24.050000000000008,
        24.06000000000001,
        24.07000000000001,
        24.080000000000013,
        24.090000000000014,
        24.100000000000016,
        24.110000000000017,
        24.12000000000002,
        24.13000000000002,
        24.140000000000022,
        24.150000000000023,
        24.160000000000025,
        24.170000000000027,
        24.180000000000028,
        24.19000000000003,
        24.20000000000003,
        24.210000000000033,
        24.220000000000034,
        24.230000000000036,
        24.240000000000038,
        24.25000000000004,
        24.26000000000004,
        24.270000000000042,
        24.280000000000044,
        24.290000000000045,
        24.300000000000047,
        24.31000000000005,
        24.32000000000005,
        24.33000000000005,
        24.340000000000053,
        24.350000000000055,
        24.360000000000056,
        24.370000000000058,
        24.38000000000006,
        24.39000000000006,
        24.400000000000063,
        24.410000000000064,
        24.420000000000066,
        24.430000000000067,
        24.44000000000007,
        24.45000000000007,
        24.460000000000072,
        24.470000000000073,
        24.480000000000075,
        24.490000000000077,
        24.500000000000078,
        24.51000000000008,
        24.52000000000008,
        24.530000000000083,
        24.540000000000084,
        24.550000000000086,
        24.560000000000088,
        24.57000000000009,
        24.58000000000009,
        24.590000000000092,
        24.600000000000094,
        24.610000000000095,
        24.620000000000097,
        24.6300000000001,
        24.6400000000001,
        24.6500000000001,
        24.660000000000103,
        24.670000000000105,
        24.680000000000106,
        24.690000000000108,
        24.70000000000011,
        24.71000000000011,
        24.720000000000113,
        24.730000000000114,
        24.740000000000116,
        24.750000000000117,
        24.76000000000012,
        24.77000000000012,
        24.780000000000122,
        24.790000000000123,
        24.800000000000125,
        24.810000000000127,
        24.820000000000128,
        24.83000000000013,
        24.84000000000013,
        24.850000000000133,
        24.860000000000134,
        24.870000000000136,
        24.880000000000138,
        24.89000000000014,
        24.90000000000014,
        24.910000000000142,
        24.920000000000144,
        24.930000000000145,
        24.940000000000147,
        24.95000000000015,
        24.96000000000015,
        24.97000000000015,
        24.980000000000153,
        24.990000000000155,
        25.000000000000156,
        25.010000000000158,
        25.02000000000016,
        25.03000000000016,
        25.040000000000163,
        25.050000000000164,
        25.060000000000166,
        25.070000000000167,
        25.08000000000017,
        25.09000000000017,
        25.100000000000172,
        25.110000000000174,
        25.120000000000175,
        25.130000000000177,
        25.140000000000178,
        25.15000000000018,
        25.16000000000018,
        25.170000000000183,
        25.180000000000184,
        25.190000000000186,
        25.200000000000188,
        25.21000000000019,
        25.22000000000019,
        25.230000000000192,
        25.240000000000194,
        25.250000000000195,
        25.260000000000197,
        25.2700000000002,
        25.2800000000002,
        25.2900000000002,
        25.300000000000203,
        25.310000000000205,
        25.320000000000206,
        25.330000000000208,
        25.34000000000021,
        25.35000000000021,
        25.360000000000213,
        25.370000000000214,
        25.380000000000216,
        25.390000000000217,
        25.40000000000022,
        25.41000000000022,
        25.420000000000222,
        25.430000000000224,
        25.440000000000225,
        25.450000000000227,
        25.46000000000023,
        25.47000000000023,
        25.48000000000023,
        25.490000000000233,
        25.500000000000234,
        25.510000000000236,
        25.520000000000238,
        25.53000000000024,
        25.54000000000024,
        25.550000000000242,
        25.560000000000244,
        25.570000000000245,
        25.580000000000247,
        25.59000000000025,
        25.60000000000025,
        25.61000000000025,
        25.620000000000253,
        25.630000000000255,
        25.640000000000256,
        25.650000000000258,
        25.66000000000026,
        25.67000000000026,
        25.680000000000263,
        25.690000000000264,
        25.700000000000266,
        25.710000000000267,
        25.72000000000027,
        25.73000000000027,
        25.740000000000272,
        25.750000000000274,
        25.760000000000275,
        25.770000000000277,
        25.78000000000028,
        25.79000000000028,
        25.80000000000028,
        25.810000000000283,
        25.820000000000285,
        25.830000000000286,
        25.840000000000288,
        25.85000000000029,
        25.86000000000029,
        25.870000000000292,
        25.880000000000294,
        25.890000000000295,
        25.900000000000297,
        25.9100000000003,
        25.9200000000003,
        25.9300000000003,
        25.940000000000303,
        25.950000000000305,
        25.960000000000306,
        25.970000000000308,
        25.98000000000031,
        25.99000000000031,
        26.000000000000313,
        26.010000000000314,
        26.020000000000316,
        26.030000000000317,
        26.04000000000032,
        26.05000000000032,
        26.060000000000322,
        26.070000000000324,
        26.080000000000325,
        26.090000000000327,
        26.10000000000033,
        26.11000000000033,
        26.12000000000033,
        26.130000000000333,
        26.140000000000335,
        26.150000000000336,
        26.160000000000338,
        26.17000000000034,
        26.18000000000034,
        26.190000000000342,
        26.200000000000344,
        26.210000000000345,
        26.220000000000347,
        26.23000000000035,
        26.24000000000035,
        26.25000000000035,
        26.260000000000353,
        26.270000000000355,
        26.280000000000356,
        26.290000000000358,
        26.30000000000036,
        26.31000000000036,
        26.320000000000363,
        26.330000000000364,
        26.340000000000366,
        26.350000000000367,
        26.36000000000037,
        26.37000000000037,
        26.380000000000372,
        26.390000000000374,
        26.400000000000375,
        26.410000000000377,
        26.42000000000038,
        26.43000000000038,
        26.44000000000038,
        26.450000000000383,
        26.460000000000385,
        26.470000000000386,
        26.480000000000388,
        26.49000000000039,
        26.50000000000039,
        26.510000000000392,
        26.520000000000394,
        26.530000000000395,
        26.540000000000397,
        26.5500000000004,
        26.5600000000004,
        26.5700000000004,
        26.580000000000403,
        26.590000000000405,
        26.600000000000406,
        26.610000000000408,
        26.62000000000041,
        26.63000000000041,
        26.640000000000413,
        26.650000000000414,
        26.660000000000416,
        26.670000000000417,
        26.68000000000042,
        26.69000000000042,
        26.700000000000422,
        26.710000000000424,
        26.720000000000425,
        26.730000000000427,
        26.74000000000043,
        26.75000000000043,
        26.76000000000043,
        26.770000000000433,
        26.780000000000435,
        26.790000000000436,
        26.800000000000438,
        26.81000000000044,
        26.82000000000044,
        26.830000000000442,
        26.840000000000444,
        26.850000000000446,
        26.860000000000447,
        26.87000000000045,
        26.88000000000045,
        26.89000000000045,
        26.900000000000453,
        26.910000000000455,
        26.920000000000456,
        26.930000000000458,
        26.94000000000046,
        26.95000000000046,
        26.960000000000463,
        26.970000000000464,
        26.980000000000466,
        26.990000000000467,
        27.00000000000047,
        27.01000000000047,
        27.020000000000472,
        27.030000000000474,
        27.040000000000475,
        27.050000000000477,
        27.06000000000048,
        27.07000000000048,
        27.08000000000048,
        27.090000000000483,
        27.100000000000485,
        27.110000000000486,
        27.120000000000488,
        27.13000000000049,
        27.14000000000049,
        27.150000000000492,
        27.160000000000494,
        27.170000000000496,
        27.180000000000497,
        27.1900000000005,
        27.2000000000005,
        27.210000000000502,
        27.220000000000503,
        27.230000000000505,
        27.240000000000506,
        27.250000000000508,
        27.26000000000051,
        27.27000000000051,
        27.280000000000513,
        27.290000000000514,
        27.300000000000516,
        27.310000000000517,
        27.32000000000052,
        27.33000000000052,
        27.340000000000522,
        27.350000000000524,
        27.360000000000525,
        27.370000000000527,
        27.38000000000053,
        27.39000000000053,
        27.40000000000053,
        27.410000000000533,
        27.420000000000535,
        27.430000000000536,
        27.440000000000538,
        27.45000000000054,
        27.46000000000054,
        27.470000000000542,
        27.480000000000544,
        27.490000000000546,
        27.500000000000547,
        27.51000000000055,
        27.52000000000055,
        27.530000000000552,
        27.540000000000553,
        27.550000000000555,
        27.560000000000556,
        27.570000000000558,
        27.58000000000056,
        27.59000000000056,
        27.600000000000563,
        27.610000000000564,
        27.620000000000566,
        27.630000000000567,
        27.64000000000057,
        27.65000000000057,
        27.660000000000572,
        27.670000000000574,
        27.680000000000575,
        27.690000000000577,
        27.70000000000058,
        27.71000000000058,
        27.72000000000058,
        27.730000000000583,
        27.740000000000585,
        27.750000000000586,
        27.760000000000588,
        27.77000000000059,
        27.78000000000059,
        27.790000000000592,
        27.800000000000594,
        27.810000000000596,
        27.820000000000597,
        27.8300000000006,
        27.8400000000006,
        27.850000000000602,
        27.860000000000603,
        27.870000000000605,
        27.880000000000607,
        27.890000000000608,
        27.90000000000061,
        27.91000000000061,
        27.920000000000613,
        27.930000000000614,
        27.940000000000616,
        27.950000000000617,
        27.96000000000062,
        27.97000000000062,
        27.980000000000622,
        27.990000000000624,
        28.000000000000625,
        28.010000000000627,
        28.02000000000063,
        28.03000000000063,
        28.04000000000063,
        28.050000000000633,
        28.060000000000635,
        28.070000000000636,
        28.080000000000638,
        28.09000000000064,
        28.10000000000064,
        28.110000000000642,
        28.120000000000644,
        28.130000000000646,
        28.140000000000647,
        28.15000000000065,
        28.16000000000065,
        28.170000000000652,
        28.180000000000653,
        28.190000000000655,
        28.200000000000657,
        28.210000000000658,
        28.22000000000066,
        28.23000000000066,
        28.240000000000663,
        28.250000000000664,
        28.260000000000666,
        28.270000000000667,
        28.28000000000067,
        28.29000000000067,
        28.300000000000672,
        28.310000000000674,
        28.320000000000675,
        28.330000000000677,
        28.34000000000068,
        28.35000000000068,
        28.36000000000068,
        28.370000000000683,
        28.380000000000685,
        28.390000000000686,
        28.400000000000688,
        28.41000000000069,
        28.42000000000069,
        28.430000000000692,
        28.440000000000694,
        28.450000000000696,
        28.460000000000697,
        28.4700000000007,
        28.4800000000007,
        28.490000000000702,
        28.500000000000703,
        28.510000000000705,
        28.520000000000707,
        28.530000000000708,
        28.54000000000071,
        28.55000000000071,
        28.560000000000713,
        28.570000000000714,
        28.580000000000716,
        28.590000000000718,
        28.60000000000072,
        28.61000000000072,
        28.620000000000722,
        28.630000000000724,
        28.640000000000725,
        28.650000000000727,
        28.66000000000073,
        28.67000000000073,
        28.68000000000073,
        28.690000000000733,
        28.700000000000735,
        28.710000000000736,
        28.720000000000738,
        28.73000000000074,
        28.74000000000074,
        28.750000000000743,
        28.760000000000744,
        28.770000000000746,
        28.780000000000747,
        28.79000000000075,
        28.80000000000075,
        28.810000000000752,
        28.820000000000753,
        28.830000000000755,
        28.840000000000757,
        28.850000000000758,
        28.86000000000076,
        28.87000000000076,
        28.880000000000763,
        28.890000000000764,
        28.900000000000766,
        28.910000000000768,
        28.92000000000077,
        28.93000000000077,
        28.940000000000772,
        28.950000000000774,
        28.960000000000775,
        28.970000000000777,
        28.98000000000078,
        28.99000000000078,
        29.00000000000078,
        29.010000000000783,
        29.020000000000785,
        29.030000000000786,
        29.040000000000788,
        29.05000000000079,
        29.06000000000079,
        29.070000000000793,
        29.080000000000794,
        29.090000000000796,
        29.100000000000797,
        29.1100000000008,
        29.1200000000008,
        29.130000000000802,
        29.140000000000803,
        29.150000000000805,
        29.160000000000807,
        29.170000000000808,
        29.18000000000081,
        29.19000000000081,
        29.200000000000813,
        29.210000000000814,
        29.220000000000816,
        29.230000000000818,
        29.24000000000082,
        29.25000000000082,
        29.260000000000822,
        29.270000000000824,
        29.280000000000825,
        29.290000000000827,
        29.30000000000083,
        29.31000000000083,
        29.32000000000083,
        29.330000000000833,
        29.340000000000835,
        29.350000000000836,
        29.360000000000838,
        29.37000000000084,
        29.38000000000084,
        29.390000000000843,
        29.400000000000844,
        29.410000000000846,
        29.420000000000847,
        29.43000000000085,
        29.44000000000085,
        29.450000000000852,
        29.460000000000854,
        29.470000000000855,
        29.480000000000857,
        29.490000000000858,
        29.50000000000086,
        29.51000000000086,
        29.520000000000863,
        29.530000000000864,
        29.540000000000866,
        29.550000000000868,
        29.56000000000087,
        29.57000000000087,
        29.580000000000872,
        29.590000000000874,
        29.600000000000875,
        29.610000000000877,
        29.62000000000088,
        29.63000000000088,
        29.64000000000088,
        29.650000000000883,
        29.660000000000885,
        29.670000000000886,
        29.680000000000888,
        29.69000000000089,
        29.70000000000089,
        29.710000000000893,
        29.720000000000894,
        29.730000000000896,
        29.740000000000897,
        29.7500000000009,
        29.7600000000009,
        29.770000000000902,
        29.780000000000904,
        29.790000000000905,
        29.800000000000907,
        29.810000000000908,
        29.82000000000091,
        29.83000000000091,
        29.840000000000913,
        29.850000000000914,
        29.860000000000916,
        29.870000000000918,
        29.88000000000092,
        29.89000000000092,
        29.900000000000922,
        29.910000000000924,
        29.920000000000925,
        29.930000000000927,
        29.94000000000093,
        29.95000000000093,
        29.96000000000093,
        29.970000000000933,
        29.980000000000935,
        29.990000000000936,
        30.000000000000938,
        30.01000000000094,
        30.02000000000094,
        30.030000000000943,
        30.040000000000944,
        30.050000000000946,
        30.060000000000947,
        30.07000000000095,
        30.08000000000095,
        30.090000000000952,
        30.100000000000954,
        30.110000000000955,
        30.120000000000957,
        30.13000000000096,
        30.14000000000096,
        30.15000000000096,
        30.160000000000963,
        30.170000000000964,
        30.180000000000966,
        30.190000000000968,
        30.20000000000097,
        30.21000000000097,
        30.220000000000972,
        30.230000000000974,
        30.240000000000975,
        30.250000000000977,
        30.26000000000098,
        30.27000000000098,
        30.28000000000098,
        30.290000000000983,
        30.300000000000985,
        30.310000000000986,
        30.320000000000988,
        30.33000000000099,
        30.34000000000099,
        30.350000000000993,
        30.360000000000994,
        30.370000000000996,
        30.380000000000997,
        30.390000000001,
        30.400000000001,
        30.410000000001002,
        30.420000000001004,
        30.430000000001005,
        30.440000000001007,
        30.45000000000101,
        30.46000000000101,
        30.47000000000101,
        30.480000000001013,
        30.490000000001015,
        30.500000000001016,
        30.510000000001018,
        30.52000000000102,
        30.53000000000102,
        30.540000000001022,
        30.550000000001024,
        30.560000000001025,
        30.570000000001027,
        30.58000000000103,
        30.59000000000103,
        30.60000000000103,
        30.610000000001033,
        30.620000000001035,
        30.630000000001036,
        30.640000000001038,
        30.65000000000104,
        30.66000000000104,
        30.670000000001043,
        30.680000000001044,
        30.690000000001046,
        30.700000000001047,
        30.71000000000105,
        30.72000000000105,
        30.730000000001052,
        30.740000000001054,
        30.750000000001055,
        30.760000000001057,
        30.77000000000106,
        30.78000000000106,
        30.79000000000106,
        30.800000000001063,
        30.810000000001065,
        30.820000000001066,
        30.830000000001068,
        30.84000000000107,
        30.85000000000107,
        30.860000000001072,
        30.870000000001074,
        30.880000000001075,
        30.890000000001077,
        30.90000000000108,
        30.91000000000108,
        30.92000000000108,
        30.930000000001083,
        30.940000000001085,
        30.950000000001086,
        30.960000000001088,
        30.97000000000109,
        30.98000000000109,
        30.990000000001093,
        31.000000000001094,
        31.010000000001096,
        31.020000000001097,
        31.0300000000011,
        31.0400000000011,
        31.050000000001102,
        31.060000000001104,
        31.070000000001105,
        31.080000000001107,
        31.09000000000111,
        31.10000000000111,
        31.11000000000111,
        31.120000000001113,
        31.130000000001115,
        31.140000000001116,
        31.150000000001118,
        31.16000000000112,
        31.17000000000112,
        31.180000000001122,
        31.190000000001124,
        31.200000000001125,
        31.210000000001127,
        31.22000000000113,
        31.23000000000113,
        31.24000000000113,
        31.250000000001133,
        31.260000000001135,
        31.270000000001136,
        31.280000000001138,
        31.29000000000114,
        31.30000000000114,
        31.310000000001143,
        31.320000000001144,
        31.330000000001146,
        31.340000000001147,
        31.35000000000115,
        31.36000000000115,
        31.370000000001152,
        31.380000000001154,
        31.390000000001155,
        31.400000000001157,
        31.41000000000116,
        31.42000000000116,
        31.43000000000116,
        31.440000000001163,
        31.450000000001165,
        31.460000000001166,
        31.470000000001168,
        31.48000000000117,
        31.49000000000117,
        31.500000000001172,
        31.510000000001174,
        31.520000000001176,
        31.530000000001177,
        31.54000000000118,
        31.55000000000118,
        31.56000000000118,
        31.570000000001183,
        31.580000000001185,
        31.590000000001186,
        31.600000000001188,
        31.61000000000119,
        31.62000000000119,
        31.630000000001193,
        31.640000000001194,
        31.650000000001196,
        31.660000000001197,
        31.6700000000012,
        31.6800000000012,
        31.690000000001202,
        31.700000000001204,
        31.710000000001205,
        31.720000000001207,
        31.73000000000121,
        31.74000000000121,
        31.75000000000121,
        31.760000000001213,
        31.770000000001215,
        31.780000000001216,
        31.790000000001218,
        31.80000000000122,
        31.81000000000122,
        31.820000000001222,
        31.830000000001224,
        31.840000000001226,
        31.850000000001227,
        31.86000000000123,
        31.87000000000123,
        31.880000000001232,
        31.890000000001233,
        31.900000000001235,
        31.910000000001236,
        31.920000000001238,
        31.93000000000124,
        31.94000000000124,
        31.950000000001243,
        31.960000000001244,
        31.970000000001246,
        31.980000000001247,
        31.99000000000125,
        32.00000000000125,
        32.01000000000125,
        32.02000000000125,
        32.030000000001245,
        32.04000000000124,
        32.05000000000124,
        32.06000000000124,
        32.07000000000124,
        32.080000000001235,
        32.09000000000123,
        32.10000000000123,
        32.11000000000123,
        32.12000000000123,
        32.130000000001225,
        32.14000000000122,
        32.15000000000122,
        32.16000000000122,
        32.17000000000122,
        32.180000000001215,
        32.19000000000121,
        32.20000000000121,
        32.21000000000121,
        32.22000000000121,
        32.230000000001205,
        32.2400000000012,
        32.2500000000012,
        32.2600000000012,
        32.2700000000012,
        32.280000000001195,
        32.29000000000119,
        32.30000000000119,
        32.31000000000119,
        32.32000000000119,
        32.330000000001185,
        32.34000000000118,
        32.35000000000118,
        32.36000000000118,
        32.37000000000118,
        32.380000000001175,
        32.39000000000117,
        32.40000000000117,
        32.41000000000117,
        32.42000000000117,
        32.430000000001165,
        32.44000000000116,
        32.45000000000116,
        32.46000000000116,
        32.47000000000116,
        32.480000000001155,
        32.49000000000115,
        32.50000000000115,
        32.51000000000115,
        32.52000000000115,
        32.530000000001145,
        32.54000000000114,
        32.55000000000114,
        32.56000000000114,
        32.57000000000114,
        32.580000000001135,
        32.59000000000113,
        32.60000000000113,
        32.61000000000113,
        32.62000000000113,
        32.630000000001125,
        32.64000000000112,
        32.65000000000112,
        32.66000000000112,
        32.67000000000112,
        32.680000000001115,
        32.69000000000111,
        32.70000000000111,
        32.71000000000111,
        32.72000000000111,
        32.730000000001105,
        32.7400000000011,
        32.7500000000011,
        32.7600000000011,
        32.7700000000011,
        32.780000000001095,
        32.79000000000109,
        32.80000000000109,
        32.81000000000109,
        32.82000000000109,
        32.830000000001085,
        32.84000000000108,
        32.85000000000108,
        32.86000000000108,
        32.87000000000108,
        32.880000000001075,
        32.89000000000107,
        32.90000000000107,
        32.91000000000107,
        32.92000000000107,
        32.930000000001066,
        32.94000000000106,
        32.95000000000106,
        32.96000000000106,
        32.97000000000106,
        32.980000000001056,
        32.990000000001054,
        33.00000000000105,
        33.01000000000105,
        33.02000000000105,
        33.030000000001046,
        33.040000000001044,
        33.05000000000104,
        33.06000000000104,
        33.07000000000104,
        33.080000000001036,
        33.090000000001034,
        33.10000000000103,
        33.11000000000103,
        33.12000000000103,
        33.130000000001026,
        33.140000000001024,
        33.15000000000102,
        33.16000000000102,
        33.17000000000102,
        33.180000000001016,
        33.190000000001014,
        33.20000000000101,
        33.21000000000101,
        33.22000000000101,
        33.230000000001006,
        33.240000000001004,
        33.250000000001,
        33.260000000001,
        33.270000000001,
        33.280000000000996,
        33.290000000000994,
        33.30000000000099,
        33.31000000000099,
        33.32000000000099,
        33.330000000000986,
        33.340000000000984,
        33.35000000000098,
        33.36000000000098,
        33.37000000000098,
        33.380000000000976,
        33.390000000000974,
        33.40000000000097,
        33.41000000000097,
        33.42000000000097,
        33.430000000000966,
        33.440000000000964,
        33.45000000000096,
        33.46000000000096,
        33.47000000000096,
        33.480000000000956,
        33.490000000000954,
        33.50000000000095,
        33.51000000000095,
        33.52000000000095,
        33.530000000000946,
        33.540000000000944,
        33.55000000000094,
        33.56000000000094,
        33.57000000000094,
        33.580000000000936,
        33.590000000000934,
        33.60000000000093,
        33.61000000000093,
        33.62000000000093,
        33.630000000000926,
        33.640000000000924,
        33.65000000000092,
        33.66000000000092,
        33.67000000000092,
        33.680000000000916,
        33.690000000000914,
        33.70000000000091,
        33.71000000000091,
        33.72000000000091,
        33.730000000000906,
        33.740000000000904,
        33.7500000000009,
        33.7600000000009,
        33.7700000000009,
        33.780000000000896,
        33.790000000000894,
        33.80000000000089,
        33.81000000000089,
        33.82000000000089,
        33.83000000000089,
        33.840000000000884,
        33.85000000000088,
        33.86000000000088,
        33.87000000000088,
        33.88000000000088,
        33.890000000000875,
        33.90000000000087,
        33.91000000000087,
        33.92000000000087,
        33.93000000000087,
        33.940000000000865,
        33.95000000000086,
        33.96000000000086,
        33.97000000000086,
        33.98000000000086,
        33.990000000000855,
    ]
    @test map(p -> getfield(p, :y), trace) == [
        1876.80005,
        1876.8952150604855,
        1876.9903693287204,
        1877.0855128017947,
        1877.1806454767993,
        1877.275767350825,
        1877.3708784209625,
        1877.4659786843033,
        1877.5610681379387,
        1877.6561467789602,
        1877.7512146044598,
        1877.8462716115293,
        1877.9413177972613,
        1878.036353158748,
        1878.131377693082,
        1878.2263913973563,
        1878.321394268664,
        1878.4163863040988,
        1878.5113675007537,
        1878.6063378557226,
        1878.7012973660994,
        1878.7962460289784,
        1878.891183841454,
        1878.9861108006205,
        1879.081026903573,
        1879.1759321474065,
        1879.270826529216,
        1879.365710046097,
        1879.4605826951454,
        1879.5554444734569,
        1879.6502953781276,
        1879.745135406254,
        1879.8399645549323,
        1879.9347828212594,
        1880.0295902020118,
        1880.124386689888,
        1880.219172280366,
        1880.3139469703863,
        1880.408710756891,
        1880.5034636368214,
        1880.5982056071198,
        1880.6929366647282,
        1880.787656806589,
        1880.8823660296448,
        1880.9770643308386,
        1881.0717517071137,
        1881.166428155413,
        1881.2610936726805,
        1881.35574825586,
        1881.4503919018953,
        1881.545024607731,
        1881.6396463703115,
        1881.7342571865815,
        1881.8288570534862,
        1881.9234459679708,
        1882.0180239269805,
        1882.1125909274613,
        1882.2071469663592,
        1882.3016920406203,
        1882.396226147191,
        1882.4907492830182,
        1882.5852614450484,
        1882.679762630229,
        1882.7742528355075,
        1882.8687320578313,
        1882.9632002941482,
        1883.0576575414066,
        1883.1521037965547,
        1883.246539056541,
        1883.3409633183144,
        1883.435376578824,
        1883.529778835019,
        1883.624170083849,
        1883.7185503222638,
        1883.8129195472134,
        1883.907277755648,
        1884.0016249445182,
        1884.0959611107746,
        1884.1902862513682,
        1884.2846003632503,
        1884.3789034433723,
        1884.4731954886859,
        1884.567476496143,
        1884.6617464626959,
        1884.756005385297,
        1884.8502532608989,
        1884.9444900864546,
        1885.03871586027,
        1885.1329305946256,
        1885.2271342897268,
        1885.3213269426763,
        1885.4155085505763,
        1885.5096791105302,
        1885.603838619641,
        1885.697987075012,
        1885.7921244737468,
        1885.8862508129491,
        1885.980366089723,
        1886.0744703011728,
        1886.1685634444032,
        1886.2626455165184,
        1886.3567165146237,
        1886.4507764358243,
        1886.5448252772255,
        1886.638863035933,
        1886.7328897090522,
        1886.8269052936896,
        1886.9209097869514,
        1887.014903185944,
        1887.1088854877742,
        1887.202856689549,
        1887.2968167883755,
        1887.3907657813613,
        1887.484703665614,
        1887.5786304382414,
        1887.6725460963517,
        1887.7664506370531,
        1887.8603440574543,
        1887.954226354664,
        1888.048097525791,
        1888.141957567945,
        1888.2358064782352,
        1888.3296442537714,
        1888.4234708916636,
        1888.5172863890216,
        1888.611090742956,
        1888.7048839505778,
        1888.7986660089973,
        1888.8924369153258,
        1888.9861966666747,
        1889.0799452601552,
        1889.1736826928795,
        1889.2674089619593,
        1889.361124064507,
        1889.454827997635,
        1889.5485207584559,
        1889.6422023440828,
        1889.7358727516287,
        1889.829531978207,
        1889.923180020931,
        1890.0168168769676,
        1890.1104425445444,
        1890.20405702151,
        1890.2976603051168,
        1890.3912523926176,
        1890.4848332812655,
        1890.5784029683136,
        1890.6719614510152,
        1890.765508726624,
        1890.8590447923934,
        1890.9525696455778,
        1891.0460832834312,
        1891.1395857032082,
        1891.2330769021635,
        1891.3265568775519,
        1891.4200256266286,
        1891.5134831466487,
        1891.6069294348679,
        1891.700364488542,
        1891.7937883049267,
        1891.8872008812784,
        1891.9806022148537,
        1892.073992302909,
        1892.1673711427009,
        1892.2607387314868,
        1892.3540950665238,
        1892.4474401450693,
        1892.5407739643813,
        1892.6340965217175,
        1892.7274078143362,
        1892.8207078394955,
        1892.9139965944541,
        1893.007274076471,
        1893.100540282805,
        1893.1937952107153,
        1893.2870388574615,
        1893.3802712203033,
        1893.4734922965004,
        1893.566702083313,
        1893.6599005780015,
        1893.7530877778265,
        1893.8462636800489,
        1893.9394282819294,
        1894.0325815807294,
        1894.12572357371,
        1894.2188542581334,
        1894.311973631261,
        1894.4050816903552,
        1894.4981784326783,
        1894.5912638554926,
        1894.684337956061,
        1894.7774007316466,
        1894.8704521795125,
        1894.9634922969221,
        1895.0565210777756,
        1895.149538504419,
        1895.2425445726008,
        1895.3355392793344,
        1895.4285226216339,
        1895.5214945965133,
        1895.614455200987,
        1895.707405231378,
        1895.800345484368,
        1895.893275957511,
        1895.9861966483606,
        1896.0791075544714,
        1896.1720086733974,
        1896.2649000026931,
        1896.357781539913,
        1896.4506532826122,
        1896.543515228345,
        1896.636367374667,
        1896.7292097191332,
        1896.822042259299,
        1896.9148649927197,
        1897.0076779169513,
        1897.1004810295494,
        1897.1932743280704,
        1897.2860578100701,
        1897.378831473105,
        1897.4715953147315,
        1897.5643493325063,
        1897.6570935239858,
        1897.7498278867276,
        1897.8425524182885,
        1897.9352671162258,
        1898.0279719780967,
        1898.1206670014592,
        1898.2133521838707,
        1898.3060275228895,
        1898.3986930160734,
        1898.4913486609805,
        1898.5839944551694,
        1898.6766303961988,
        1898.769256481627,
        1898.8618727090134,
        1898.9544790759164,
        1899.0470755798956,
        1899.1396622185102,
        1899.23223898932,
        1899.3248058898844,
        1899.4173629177633,
        1899.5099100705165,
        1899.6024473457044,
        1899.6949747408871,
        1899.7874922536255,
        1899.8799998814798,
        1899.972497622011,
        1900.064985478993,
        1900.1574634671238,
        1900.2499315850578,
        1900.3423898302274,
        1900.4348382000649,
        1900.5272766920025,
        1900.6197053034728,
        1900.7121240319086,
        1900.804532874743,
        1900.896931829409,
        1900.9893208933395,
        1901.0817000639684,
        1901.1740693387292,
        1901.2664287150556,
        1901.3587781903816,
        1901.4511177621414,
        1901.5434474277693,
        1901.6357671846997,
        1901.7280770303676,
        1901.8203769622075,
        1901.9126669776547,
        1902.0049470741442,
        1902.0972172491113,
        1902.1894774999917,
        1902.2817278242212,
        1902.3739682192356,
        1902.466198682471,
        1902.5584192113636,
        1902.6506298033498,
        1902.7428304558664,
        1902.8350211663499,
        1902.9272019322375,
        1903.019372750966,
        1903.1115336199732,
        1903.2036845366963,
        1903.295825498573,
        1903.3879565030413,
        1903.480077547539,
        1903.5721886295044,
        1903.664289746376,
        1903.7563808955922,
        1903.8484620745917,
        1903.9405332808135,
        1904.0325945116967,
        1904.1246457646807,
        1904.2166870372048,
        1904.3087183267085,
        1904.4007396306317,
        1904.4927509464146,
        1904.584752271497,
        1904.6767436033197,
        1904.7687249393227,
        1904.8606962769472,
        1904.9526576136338,
        1905.044608947423,
        1905.136550281859,
        1905.2284816150834,
        1905.3204029445283,
        1905.4123142676256,
        1905.5042155818076,
        1905.5961068845065,
        1905.6879881731552,
        1905.7798594451863,
        1905.871720698033,
        1905.9635719291282,
        1906.0554131359052,
        1906.1472443157975,
        1906.239065466239,
        1906.3308765846634,
        1906.4226776685048,
        1906.5144687151974,
        1906.6062497221758,
        1906.6980206868743,
        1906.789781606728,
        1906.8815324791715,
        1906.9732733016401,
        1907.0650040715693,
        1907.1567247863945,
        1907.2484354435514,
        1907.3401360404757,
        1907.4318265746037,
        1907.5235070433716,
        1907.6151774442158,
        1907.7068377745727,
        1907.7984880318793,
        1907.8901282135726,
        1907.9817583170898,
        1908.073378339868,
        1908.164988279345,
        1908.2565881329583,
        1908.348177898146,
        1908.4397575723463,
        1908.5313271529972,
        1908.6228866375372,
        1908.714436023405,
        1908.8059753080393,
        1908.8975044888791,
        1908.9890235633638,
        1909.0805325289327,
        1909.1720313830256,
        1909.263520123082,
        1909.3549987465417,
        1909.446467250845,
        1909.5379256334322,
        1909.6293738917439,
        1909.7208120232206,
        1909.812240025303,
        1909.9036578954326,
        1909.9950656310505,
        1910.0864632315474,
        1910.1778506966639,
        1910.2692280238832,
        1910.3605952106507,
        1910.4519522544117,
        1910.543299152612,
        1910.6346359026973,
        1910.7259625021136,
        1910.817278948307,
        1910.9085852387238,
        1910.9998813708107,
        1911.0911673420144,
        1911.1824431497816,
        1911.2737087915598,
        1911.3649642647958,
        1911.4562095669376,
        1911.5474446954324,
        1911.6386696477284,
        1911.7298844212735,
        1911.8210890135158,
        1911.912283421904,
        1912.0034676438866,
        1912.0946416769123,
        1912.18580551843,
        1912.2769591658891,
        1912.368102616739,
        1912.459235868429,
        1912.550358918409,
        1912.641471764129,
        1912.7325744030388,
        1912.8236668325892,
        1912.9147490502303,
        1913.005821053413,
        1913.096882839588,
        1913.1879344062068,
        1913.27897575072,
        1913.3700068705798,
        1913.4610277632373,
        1913.5520384261445,
        1913.6430388567535,
        1913.7340290525165,
        1913.8250090108857,
        1913.9159787293138,
        1914.0069382052538,
        1914.0978874361585,
        1914.188826419481,
        1914.2797551526749,
        1914.3706736331935,
        1914.4615818584907,
        1914.5524798260205,
        1914.643367533237,
        1914.7342449775945,
        1914.8251121565474,
        1914.9159690675506,
        1915.006815707951,
        1915.097652070475,
        1915.1884781488882,
        1915.2792939407466,
        1915.3700994436072,
        1915.4608946550266,
        1915.5516795725619,
        1915.6424541937704,
        1915.7332185162095,
        1915.8239725374367,
        1915.9147162550098,
        1916.0054496664866,
        1916.0961727694253,
        1916.1868855613843,
        1916.277588039922,
        1916.368280202597,
        1916.458962046968,
        1916.5496335705946,
        1916.6402947710355,
        1916.7309456458504,
        1916.8215861925987,
        1916.91221640884,
        1917.0028362921348,
        1917.0934458400427,
        1917.1840450501243,
        1917.27463391994,
        1917.3652124470507,
        1917.455780629017,
        1917.5463384634002,
        1917.6368859477611,
        1917.7274230796618,
        1917.8179498566633,
        1917.9084662763278,
        1917.998972336217,
        1918.0894680338936,
        1918.1799533669193,
        1918.2704283328571,
        1918.3608929292695,
        1918.4513471537196,
        1918.5417910037704,
        1918.6322244769851,
        1918.7226475709274,
        1918.813060283161,
        1918.9034626112493,
        1918.9938545527568,
        1919.0842361052476,
        1919.1746072662859,
        1919.2649680334366,
        1919.3553184042644,
        1919.445658376334,
        1919.535987947211,
        1919.6263071144604,
        1919.716615875648,
        1919.806914228339,
        1919.8972021701002,
        1919.987479698497,
        1920.0777468174492,
        1920.1680035342526,
        1920.2582498467527,
        1920.3484857524118,
        1920.4387112486927,
        1920.5289263330585,
        1920.6191310029724,
        1920.7093252558975,
        1920.7995090892978,
        1920.8896825006366,
        1920.979845487378,
        1921.0699980469865,
        1921.160140176926,
        1921.2502718746614,
        1921.3403931376572,
        1921.4305039633784,
        1921.5206043492904,
        1921.6106942928584,
        1921.700773791548,
        1921.790842842825,
        1921.8809014441554,
        1921.9709495930051,
        1922.0609872868408,
        1922.151014523129,
        1922.2410312993366,
        1922.3310376129302,
        1922.4210334613774,
        1922.5110188421454,
        1922.6009937527017,
        1922.6909581905143,
        1922.780912153051,
        1922.8708556377803,
        1922.9607886421702,
        1923.0507111636896,
        1923.140623906817,
        1923.2305275761441,
        1923.3204221693106,
        1923.4103076839558,
        1923.5001841177193,
        1923.590051468241,
        1923.6799097331605,
        1923.769758910118,
        1923.8595989967537,
        1923.9494299907078,
        1924.039251889621,
        1924.1290646911339,
        1924.2188683928873,
        1924.3086629925222,
        1924.3984484876796,
        1924.488224876001,
        1924.5779921551275,
        1924.6677503227008,
        1924.7574993763628,
        1924.8472393137554,
        1924.9369701325204,
        1925.0266918306263,
        1925.1164044102984,
        1925.2061078708007,
        1925.2958022097353,
        1925.385487424705,
        1925.4751635133125,
        1925.5648304731606,
        1925.6544883018523,
        1925.744136996991,
        1925.8337765561803,
        1925.9234069770232,
        1926.0130282571238,
        1926.1026403940857,
        1926.192243385513,
        1926.2818372290099,
        1926.3714219221806,
        1926.4609974626298,
        1926.5505638479622,
        1926.6401210757822,
        1926.729669143695,
        1926.8192080493059,
        1926.90873779022,
        1926.9982583640426,
        1927.0877697683798,
        1927.177272000837,
        1927.2667650590204,
        1927.3562489405358,
        1927.4457236429896,
        1927.5351891639882,
        1927.624645501138,
        1927.714092652046,
        1927.803530614319,
        1927.8929593855642,
        1927.9823789633888,
        1928.0717893454,
        1928.1611905292054,
        1928.2505825124128,
        1928.3399652926303,
        1928.4293388674655,
        1928.5187032345268,
        1928.6080583914227,
        1928.6974043357613,
        1928.7867410651518,
        1928.8760685772027,
        1928.9653868695232,
        1929.0546959397225,
        1929.1439957854097,
        1929.2332864041946,
        1929.3225677936866,
        1929.4118399514957,
        1929.501102875232,
        1929.5903565625053,
        1929.6796010109263,
        1929.7688362181052,
        1929.858062181653,
        1929.94727889918,
        1930.0364863688847,
        1930.1256845950647,
        1930.2148735768185,
        1930.3040533118815,
        1930.3932237979898,
        1930.4823850328794,
        1930.5715370142866,
        1930.660679739948,
        1930.7498132075998,
        1930.838937414979,
        1930.9280523598225,
        1931.0171580398671,
        1931.10625445285,
        1931.1953415965086,
        1931.2844194685802,
        1931.3734880668026,
        1931.4625473889134,
        1931.5515974326509,
        1931.6406381957527,
        1931.7296696759572,
        1931.818691871003,
        1931.9077047786284,
        1931.996708396572,
        1932.0857027225732,
        1932.1746877543703,
        1932.2636634897028,
        1932.35262992631,
        1932.4415870619314,
        1932.5305348943066,
        1932.6194734211751,
        1932.7084026402772,
        1932.7973225493529,
        1932.8862331461423,
        1932.975134428386,
        1933.0640263938242,
        1933.1529090401978,
        1933.2417823652474,
        1933.3306463667145,
        1933.41950104234,
        1933.5083463898652,
        1933.5971824070316,
        1933.6860090915807,
        1933.7748264412544,
        1933.8636344537947,
        1933.9524331269436,
        1934.0412224584434,
        1934.1300024460363,
        1934.2187730874653,
        1934.3075343804726,
        1934.3962863228014,
        1934.4850289121946,
        1934.5737621463954,
        1934.6624860231473,
        1934.7512005401936,
        1934.839905695278,
        1934.9286014861443,
        1935.0172879104755,
        1935.105964964822,
        1935.1946326462496,
        1935.283290952452,
        1935.3719398811234,
        1935.4605794299575,
        1935.549209596649,
        1935.6378303788922,
        1935.7264417743816,
        1935.8150437808122,
        1935.9036363958787,
        1935.9922196172765,
        1936.0807934427007,
        1936.1693578698469,
        1936.2579128964103,
        1936.346458520087,
        1936.4349947385729,
        1936.523521549564,
        1936.6120389507566,
        1936.7005469398473,
        1936.7890455145323,
        1936.8775346725085,
        1936.966014411473,
        1937.0544847291228,
        1937.142945623155,
        1937.231397091267,
        1937.3198391311564,
        1937.408271740521,
        1937.496694917059,
        1937.5851086584678,
        1937.6735129624462,
        1937.7619078266923,
        1937.8502932489048,
        1937.9386692267824,
        1938.0270357580239,
        1938.1153928403285,
        1938.2037404713953,
        1938.2920786489237,
        1938.3804073706133,
        1938.4687266341639,
        1938.5570364372752,
        1938.6453367776473,
        1938.7336276529807,
        1938.8219090609755,
        1938.9101809993322,
        1938.9984434657517,
        1939.086696457935,
        1939.1749399735827,
        1939.2631740103964,
        1939.3513985660775,
        1939.4396136383273,
        1939.5278192248477,
        1939.6160153233407,
        1939.7042019315081,
        1939.7923790470525,
        1939.8805466676758,
        1939.9687047910809,
        1940.0568534142067,
        1940.1449925320774,
        1940.2331221423378,
        1940.3212422428587,
        1940.4093528315107,
        1940.4974539061643,
        1940.5855454646905,
        1940.6736275049602,
        1940.7617000248447,
        1940.8497630222153,
        1940.9378164949433,
        1941.0258604409005,
        1941.1138948579585,
        1941.2019197439895,
        1941.2899350968653,
        1941.3779409144583,
        1941.465937194641,
        1941.5539239352856,
        1941.641901134265,
        1941.729868789452,
        1941.8178268987197,
        1941.905775459941,
        1941.9937144709893,
        1942.0816439297382,
        1942.1695638340611,
        1942.2574741818319,
        1942.3453749709242,
        1942.4332661992123,
        1942.5211478645704,
        1942.609019964873,
        1942.6968824979942,
        1942.784735461809,
        1942.8725788541922,
        1942.9604126730185,
        1943.0482369161632,
        1943.1360515815018,
        1943.2238566669093,
        1943.3116521702616,
        1943.3994380894342,
        1943.4872144223032,
        1943.5749811667445,
        1943.6627383206344,
        1943.750485881849,
        1943.8382238482652,
        1943.9259522177595,
        1944.0136709882086,
        1944.1013801574895,
        1944.1890797234794,
        1944.2767696840554,
        1944.3644500370951,
        1944.452120780476,
        1944.5397819120758,
        1944.6274334297725,
        1944.7150753314443,
        1944.802707614969,
        1944.8903302782253,
        1944.9779433190915,
        1945.065546737899,
        1945.1531405380736,
        1945.2407247176686,
        1945.32829927438,
        1945.415864205904,
        1945.503419509937,
        1945.5909651841753,
        1945.6785012263158,
        1945.7660276340553,
        1945.8535444050908,
        1945.9410515371198,
        1946.0285490278395,
        1946.1160368749474,
        1946.2035150761415,
        1946.2909836291194,
        1946.3784425315796,
        1946.46589178122,
        1946.5533313757396,
        1946.6407613128365,
        1946.72818159021,
        1946.8155922055587,
        1946.902993156582,
        1946.990384440979,
        1947.0777660564493,
        1947.1651380006929,
        1947.2525002714094,
        1947.339852866299,
        1947.427195783062,
        1947.5145290193987,
        1947.6018525730096,
        1947.6891664415957,
        1947.7764706228577,
        1947.8637651144968,
        1947.9510499142145,
        1948.0383250197121,
        1948.1255904286913,
        1948.212846138854,
        1948.3000921479022,
        1948.3873284535382,
        1948.4745550534642,
        1948.5617719453828,
        1948.648979126997,
        1948.7361765960093,
        1948.823364350123,
        1948.9105423870415,
        1948.9977107044679,
        1949.0848693001062,
        1949.17201817166,
        1949.2591573168336,
        1949.3462867333308,
        1949.433406418856,
        1949.5205163711141,
        1949.6076165878096,
        1949.69470769888,
        1949.7817903344842,
        1949.868864492659,
        1949.9559301714412,
        1950.0429873703038,
        1950.1300361017554,
        1950.2170763653835,
        1950.3041081592312,
        1950.3911314813417,
        1950.4781463297584,
        1950.5651527025248,
        1950.6521505976846,
        1950.7391400132817,
        1950.8261209473599,
        1950.913093397963,
        1951.0000573631353,
        1951.087012840921,
        1951.1739598293645,
        1951.2608983265102,
        1951.3478283304028,
        1951.4347498390869,
        1951.5216628506075,
        1951.6085673630096,
        1951.695463374338,
        1951.7823508826382,
        1951.8692298859553,
        1951.9561003823349,
        1952.0429623698224,
        1952.1298158464635,
        1952.2166608103041,
        1952.30349725939,
        1952.3903251917675,
        1952.4771446054826,
        1952.5639554985817,
        1952.650757869111,
        1952.7375517151172,
        1952.8243370346468,
        1952.9111138257467,
        1952.9978820864637,
        1953.0846418148449,
        1953.1713930089375,
        1953.2581356667886,
        1953.3448697864455,
        1953.4315953659561,
        1953.5183124033676,
        1953.605020896728,
        1953.6917208440852,
        1953.778412243487,
        1953.8650950929816,
        1953.9517693906173,
        1954.0384351344424,
        1954.1250923225052,
        1954.2117409528546,
        1954.2983810235392,
        1954.3850125326076,
        1954.471635478109,
        1954.5582498580925,
        1954.6448556706073,
        1954.7314529137027,
        1954.8180415854283,
        1954.9046216838335,
        1954.991193206968,
        1955.077756149471,
        1955.1643105046703,
        1955.2508562703842,
        1955.33739344457,
        1955.4239220251843,
        1955.5104420101848,
        1955.5969533975288,
        1955.6834561851738,
        1955.7699503710778,
        1955.8564359531983,
        1955.9429129294936,
        1956.0293812979214,
        1956.1158410564403,
        1956.2022922030085,
        1956.2887347355845,
        1956.375168652127,
        1956.4615939505948,
        1956.5480106289467,
        1956.6344186851418,
        1956.7208181171393,
        1956.8072089228986,
        1956.893591100379,
        1956.9799646475399,
        1957.0663295623413,
        1957.1526858427428,
        1957.2390334867048,
        1957.325372492187,
        1957.4117028571497,
        1957.4980245795537,
        1957.5843376573591,
        1957.6706420885266,
        1957.756937871017,
        1957.8432250027915,
        1957.929503481811,
        1958.0157733060366,
        1958.1020344734297,
        1958.1882869819517,
        1958.2745308295644,
        1958.3607660142293,
        1958.4469925339083,
        1958.5332103865635,
        1958.619419570157,
        1958.7056200826514,
        1958.7918119220087,
        1958.8779950861915,
        1958.9641695731627,
        1959.0503353808851,
        1959.1364925073215,
        1959.222640950435,
        1959.3087807081893,
        1959.394911778547,
        1959.4810341594723,
        1959.5671478489285,
        1959.6532528448795,
        1959.739349145289,
        1959.8254367481213,
        1959.9115156513406,
        1959.9975858529112,
        1960.0836473530912,
        1960.169700152379,
        1960.2557442488194,
        1960.341779640435,
        1960.4278063252489,
        1960.5138243012836,
        1960.5998335665627,
        1960.685834119109,
        1960.7718259569456,
        1960.8578090780964,
        1960.9437834805847,
        1961.0297491624342,
        1961.1157061216688,
        1961.2016543563127,
        1961.2875938643897,
        1961.3735246439242,
        1961.4594466929404,
        1961.545360009463,
        1961.6312645915164,
        1961.7171604371256,
        1961.8030475443154,
        1961.8889259111108,
        1961.974795535537,
        1962.0606564156192,
        1962.146508549383,
        1962.232351934854,
        1962.3181865700576,
        1962.40401245302,
        1962.4898295817668,
        1962.5756379543243,
        1962.6614375687188,
        1962.7472284229766,
        1962.8330105151242,
        1962.918783843188,
        1963.0045484051952,
        1963.0903041991724,
        1963.1760512231467,
        1963.2617894751454,
        1963.3475189531955,
        1963.4332396553248,
        1963.5189515795605,
        1963.6046547239307,
        1963.690349086463,
        1963.7760346651853,
        1963.861711458126,
        1963.947379463313,
        1964.0330386787748,
        1964.1186891025402,
        1964.2043307326376,
        1964.2899635670958,
        1964.3755876039438,
        1964.4612028412105,
        1964.5468092769252,
        1964.6324069091172,
        1964.717995735816,
        1964.803575755051,
        1964.8891469648522,
        1964.9747093632495,
        1965.0602629528835,
        1965.1458077445336,
        1965.2313437372795,
        1965.3168709293086,
        1965.402389318808,
        1965.487898903965,
        1965.5733996829672,
        1965.6588916540022,
        1965.7443748152577,
        1965.8298491649218,
        1965.9153147011823,
        1966.0007714222274,
        1966.0862193262453,
        1966.1716584114242,
        1966.257088675953,
        1966.3425101180196,
        1966.4279227358134,
        1966.5133265275228,
        1966.5987214913368,
        1966.6841076254445,
        1966.7694849280354,
    ]
end