export SchobelZhuParams, evaluateCharFunc, evaluateLogCharFunc, computeCumulants

struct SchobelZhuParams{T}
    v0::T
    κ::T
    θ::T
    ρ::T
    σ::T
end
DefaultCharFunc(params::SchobelZhuParams{Float64}) = DefaultCharFunc{SchobelZhuParams{Float64},Complex}(params)

cinf(params::SchobelZhuParams{T}, τ::T) where {T} = (params.v0^2/params.σ+τ*params.σ)/2 *sqrt(1-params.ρ^2)

#CC the Nemo arbField or Complex type. z: Nemo complex number or Complex number.
@inline function evaluateCharFunc(p::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T)::CR where {T,CR,CT}
    return exp(evaluateLogCharFunc(p, z, τ))
end

evaluateLogCharFunc(p::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T) where {T,CR,CT} = evaluateLogCharFuncZhu(p, z, τ)

function evaluateLogCharFuncLK(cf::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T)::CR where {T,CR,CT}
    #Lord and Kahl formulation
    p = model(cf) #model(cf), abstract cf with concrete types: LK, Zhu.
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    α = -(z * (z + cc1)) / 2
    β = 2 * κ - 2 * σ * ρ * cc1 * z
    γ = 2 * σ^2
    D = sqrt(β^2 - 4 * α * γ)
    G = (β - D) / (β + D)
    Dt = D * τ
    eDT = exp(-Dt)
    oneGeDT = (1 - G * eDT)
    l = log(oneGeDT / (1 - G))
    sqrteDT = exp(-D * τ / 2) #do not use sqrt(eDT) as it may be the other root
    As =
        (β - D) * (κ * θ)^2 / (2 * σ^2 * D^3) * (
            β * (Dt - 4) +
            D * (Dt - 2) +
            4 * sqrteDT * ((D^2 - 2 * β^2) / (β + D) * sqrteDT + 2 * β) / oneGeDT
        )
    A = As + (β - D) * τ / 4 - l / 2
    Bs = 2 * κ * θ * (β - D) * (1 - sqrteDT)^2 / (D * γ * oneGeDT)
    Bv = (1 - eDT) / oneGeDT * (β - D) / (2 * γ)
    return A + Bs * v0 + Bv * v0^2
end

function evaluateLogCharFuncZhu(cf::CharFunc{SchobelZhuParams{T},CR}, z::CT, τ::T)::CR where {T,CR,CT}
    #Zhu formulation: faster but issue with complex log. Here we use Cui log adj.
    p = model(cf)
    v0 = p.v0
    κ = p.κ
    θ = p.θ
    ρ = p.ρ
    σ = p.σ
    cc1 = oneim(cf)
    iu = cc1 * z
    s12 = -iu / 2 * (iu * (1 - ρ^2) - 1 + 2 * ρ * κ / σ)
    s22 = ρ * κ * θ / σ * iu
    s32 = ρ / (2 * σ) * iu
    γ1 = sqrt(2 * σ^2 * s12 + κ^2)
    γ3 = κ^2 * θ - s22 * σ^2
    γ2 = (κ - 2 * σ^2 * s32) / γ1
    ch1 = cosh(γ1 * τ)
    sh1 = sinh(γ1 * τ)
    egt = ch1 + sh1
    γ4 = ch1 + γ2 * sh1
    H3 = (κ - γ1 * (sh1 + γ2 * ch1) / γ4) / σ^2
    H4 = ((κ * θ * γ1 - γ2 * γ3) + γ3 * (sh1 + γ2 * ch1)) / (γ4 * γ1 * σ^2) - κ * θ / σ^2
    #H5 = -(γ1*τ-log(γ1)+ log(γ1* γ4/egt) )/ 2
    H5 = -(γ1 * τ + log(γ4 / egt)) / 2 # here we adjust the log to always use the main branch.
    #original SZ
    #H5 = -log(γ4)/2
    H5 += κ * τ / 2 + (sh1 / γ4 - γ1 * τ) / (2 * γ1^3 * σ^2) * ((κ^2 * θ^2 * γ1^2 - γ3^2))
    H5 += γ3 / (γ1^3 * γ4 * σ^2) * (κ * θ * γ1 - γ2 * γ3) * (ch1 - 1)
    return -s32 * v0^2 - iu * ρ * σ * τ / 2 + H3 / 2 * v0^2 + H4 * v0 + H5
end

function computeCumulants(p::SchobelZhuParams{T}, τ::T) where {T}
    κ = p.κ
    θ = p.θ
    v0 = p.v0
    σ = p.σ
    ρ = p.ρ

    shk = sinh(κ * τ)
    chk = cosh(κ * τ)
    c1H3 = (chk * (chk - shk) * σ + ρ * κ - σ) / (σ * κ)
    c1H4 = (chk - 1) * θ * (shk - chk) / κ
    c1H5 =
        (
            (2 * κ * θ^2 - σ^2) * chk^2 +
            ((-2 * κ * θ^2 + σ^2) * shk - 4 * κ * θ^2) * chk +
            4 * shk * κ * θ^2 +
            2 * τ * (ρ * σ - θ^2) * κ^2 +
            (-σ^2 * τ + 2 * θ^2) * κ +
            σ^2
        ) / (4 * κ^2)
    c1 = -ρ / (2 * σ) * v0^2 - ρ * σ * τ / 2 + c1H3 / 2 * v0^2 + c1H4 * v0 + c1H5
    c2H3 =
        (
            (-2 * κ * ρ * σ + 2 * κ^2) * chk^2 -
            2 * (κ * ρ * σ - 1 / 2 * σ^2 - κ^2) * shk * chk +
            (2 * ρ * σ * τ - 2) * κ^2 +
            (-σ^2 * τ + 2 * ρ * σ) * κ
        ) / (κ^3 * (2 * chk * shk + 2 * chk^2 - 1))
    c2H4 =
        2 *
        (
            (-2 * κ * ρ * σ + κ^2 + σ^2) * chk^2 +
            (
                (-2 * κ * ρ * σ + 1 / 2 * σ^2 + κ^2) * shk +
                (ρ * σ * τ - 1) * κ^2 +
                (-1 / 2 * σ^2 * τ + 2 * ρ * σ) * κ - σ^2
            ) * chk +
            ((ρ * σ * τ - 1) * κ^2 + (-1 / 2 * σ^2 * τ + 2 * ρ * σ) * κ - σ^2 / 2) * shk -
            κ * σ * (ρ * κ - σ / 2) * τ
        ) *
        θ / (κ^3 * (2 * chk * shk + 2 * chk^2 - 1))
    c2H5 =
        (
            (-8 * κ * σ^2 * θ^2 + 2 * σ^4) * chk^4 -
            2 * σ^2 * ((-4 * κ * θ^2 + σ^2) * shk - 8 * κ * θ^2) * chk^3 -
            8 *
            κ *
            (
                2 * σ^2 * θ^2 * shk +
                (-2 * ρ * σ * τ * θ^2 + θ^2) * κ^2 +
                σ * (ρ * τ * σ^2 + (τ * θ^2 - 1 / 2) * σ - 3 * ρ * θ^2) * κ +
                σ^2 * (-1 / 2 * σ^2 * τ + ρ * σ + 2 * θ^2)
            ) *
            chk^2 +
            (
                (
                    (-16 * ρ * σ * τ * θ^2 + 8 * θ^2) * κ^3 +
                    (8 * ρ * τ * σ^3 + (8 * τ * θ^2 - 4) * σ^2 - 24 * ρ * σ * θ^2) * κ^2 +
                    (-4 * σ^4 * τ + 8 * ρ * σ^3 + 20 * σ^2 * θ^2) * κ - σ^4
                ) * shk -
                16 *
                ((ρ * σ * τ - 1) * κ^2 + (-1 / 2 * σ^2 * τ + 3 * ρ * σ) * κ - σ^2) *
                κ *
                θ^2
            ) * chk +
            16 *
            κ *
            θ^2 *
            ((ρ * σ * τ - 1) * κ^2 + (-1 / 2 * σ^2 * τ + 3 * ρ * σ) * κ - (3 * σ^2) / 2) *
            shk +
            8 * κ^4 * τ * θ^2 +
            (-24 * ρ * σ * τ * θ^2 + 4 * σ^2 * τ - 8 * θ^2) * κ^3 +
            ((12 * τ * θ^2 - 4) * σ^2 + 24 * ρ * σ * θ^2) * κ^2 +
            (-σ^4 * τ + 8 * ρ * σ^3 - 8 * σ^2 * θ^2) * κ - 2 * σ^4
        ) / (8 * κ^4)
    c2 = c2H3 / 2 * v0^2 + c2H4 * v0 + c2H5

    c4H3 =
        32 *
        σ^2 *
        (
            (
                (3 * ρ^2 + 3 / 4) * κ^4 - 3 * σ * ρ * (ρ^2 + 3 / 2) * κ^3 +
                9 * (ρ^2 + 1 / 4) * σ^2 * κ^2 / 2 - (15 * κ * ρ * σ^3) / 8 + (3 * σ^4) / 16
            ) * chk^4 -
            3 *
            (
                (-ρ^2 - 1 / 4) * κ^4 + σ * ρ * (ρ^2 + 3 / 2) * κ^3 -
                3 * (ρ^2 + 1 / 4) * σ^2 * κ^2 / 2 + (5 * κ * ρ * σ^3) / 8 - (3 * σ^4) / 32
            ) *
            shk *
            chk^3 +
            (
                ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^6 +
                (3 * τ * (-ρ^2 * σ^2 * τ^2 + ρ * τ * (ρ^2 + 2) * σ - ρ^2 - 1 / 2) * κ^5) /
                2 +
                (
                    (3 * ρ * σ^3 * τ^3) / 4 - 3 * τ^2 * (ρ^2 + 3 / 8) * σ^2 +
                    (3 * ρ * τ * (ρ^2 + 2) * σ) / 2 - 3 / 4 - (15 * ρ^2) / 4
                ) * κ^4 +
                (
                    15 *
                    σ *
                    (
                        -σ^3 * τ^3 / 30 +
                        (2 * ρ * τ^2 * σ^2) / 5 +
                        (-3 / 5 * ρ^2 * τ - 3 / 20 * τ) * σ +
                        ρ^3 +
                        (6 * ρ) / 5
                    ) *
                    κ^3
                ) / 4 -
                (9 * σ^2 * (1 / 24 * σ^2 * τ^2 - 5 / 24 * ρ * σ * τ + ρ^2 + 1 / 4) * κ^2) /
                2 + 15 * (-σ * τ / 10 + ρ) * σ^3 * κ / 8 - (3 * σ^4) / 16
            ) * chk^2 +
            (
                ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^6 +
                (3 * τ * (-ρ^2 * σ^2 * τ^2 + ρ * τ * (ρ^2 + 2) * σ - ρ^2 - 1 / 2) * κ^5) /
                2 +
                (
                    (3 * ρ * σ^3 * τ^3) / 4 - 3 * τ^2 * (ρ^2 + 3 / 8) * σ^2 +
                    (3 * ρ * τ * (ρ^2 + 2) * σ) / 2 - 3 / 8 - (9 * ρ^2) / 4
                ) * κ^4 +
                (
                    9 *
                    σ *
                    (
                        -σ^3 * τ^3 / 18 +
                        (2 * ρ * τ^2 * σ^2) / 3 +
                        (-ρ^2 * τ - 1 / 4 * τ) * σ +
                        ρ^3 +
                        ρ
                    ) *
                    κ^3
                ) / 4 -
                9 * (1 / 12 * σ^2 * τ^2 - 1 / 6 * ρ * σ * τ + ρ^2 + 1 / 8) * σ^2 * κ^2 / 4 +
                (15 * σ^3 * (σ * τ / 5 + ρ) * κ) / 32 - (3 * σ^4) / 64
            ) *
            shk *
            chk -
            (
                (
                    ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^5 +
                    (
                        3 *
                        τ *
                        (-ρ^2 * σ^2 * τ^2 + ρ * τ * (ρ^2 + 2) * σ - ρ^2 - 1 / 2) *
                        κ^4
                    ) / 2 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 - (3 * τ^2 * (ρ^2 + 3 / 4) * σ^2) / 2 +
                        (3 * ρ * τ * (ρ^2 + 1) * σ) / 2 - (3 * ρ^2) / 2
                    ) * κ^3 +
                    (-1 / 8 * σ^4 * τ^3 + 3 / 2 * ρ^3 * σ + 9 / 16 * σ^2 * τ) * κ^2 -
                    (15 * σ^3 * (-σ * τ / 5 + ρ) * τ * κ) / 16 + (3 * σ^4 * τ) / 32
                ) * κ
            ) / 2
        ) / (κ^7 * (8 * chk^4 + 8 * chk^3 * shk - 8 * chk^2 - 4 * chk * shk + 1))
    c4H4 =
        -32 *
        σ^2 *
        θ *
        (
            -(3 * chk^8 * σ^4) / 4 +
            (3 * σ^4 * (shk + 1) * chk^7) / 4 +
            (
                9 *
                σ^2 *
                (-σ^2 * shk / 6 + κ * ((ρ * σ * τ - 1 / 2) * κ + σ * (-σ * τ / 2 + ρ))) *
                chk^6
            ) / 2 -
            9 *
            (
                ((ρ * σ * τ - 1 / 2) * κ^2 + σ * (-σ * τ / 2 + ρ) * κ - σ^2 / 12) * shk +
                5 * ((ρ * σ * τ - 3 / 5) * κ - σ^2 * τ / 2 + (6 * ρ * σ) / 5) * κ / 6
            ) *
            σ^2 *
            chk^5 / 2 +
            (
                15 *
                (
                    (ρ * σ * τ - 3 / 5) * κ^2 + (-1 / 2 * σ^2 * τ + 6 / 5 * ρ * σ) * κ -
                    σ^2 / 10
                ) *
                σ^2 *
                shk / 4 -
                6 *
                κ *
                (
                    (ρ^2 * σ^2 * τ^2 - ρ * σ * τ + 1 / 8) * κ^3 +
                    2 * σ * (-ρ * τ^2 * σ^2 / 2 + τ * (ρ^2 + 3 / 8) * σ - ρ / 2) * κ^2 +
                    (1 / 4 * τ^2 * σ^4 - 7 / 8 * ρ * τ * σ^3 + 5 / 4 * ρ^2 * σ^2) * κ -
                    σ^4 * τ / 8
                )
            ) * chk^4 +
            (
                (
                    (6 * ρ^2 * σ^2 * τ^2 - 6 * ρ * σ * τ + 3 / 4) * κ^4 +
                    12 * σ * (-ρ * τ^2 * σ^2 / 2 + τ * (ρ^2 + 3 / 8) * σ - ρ / 2) * κ^3 +
                    (15 * σ^2 * (1 / 5 * σ^2 * τ^2 - ρ * σ * τ + ρ^2 + 3 / 20) * κ^2) / 2 -
                    (9 * σ^3 * (-σ * τ / 6 + ρ) * κ) / 4 + (9 * σ^4) / 32
                ) * shk +
                (
                    27 *
                    κ *
                    (
                        (ρ^2 * σ^2 * τ^2 - 4 / 3 * ρ * σ * τ + 2 / 9) * κ^3 +
                        (
                            8 *
                            σ *
                            (
                                -(3 * ρ * τ^2 * σ^2) / 8 + τ * (ρ^2 + 3 / 8) * σ -
                                (2 * ρ) / 3
                            ) *
                            κ^2
                        ) / 3 +
                        (1 / 4 * τ^2 * σ^4 - 10 / 9 * ρ * τ * σ^3 + 20 / 9 * ρ^2 * σ^2) * κ -
                        (7 * σ^4 * τ) / 36
                    )
                ) / 8
            ) * chk^3 +
            (
                (
                    (-27 / 8 * ρ^2 * σ^2 * τ^2 + 9 / 2 * ρ * σ * τ - 3 / 4) * κ^4 -
                    9 *
                    σ *
                    (-(3 * ρ * τ^2 * σ^2) / 8 + τ * (ρ^2 + 3 / 8) * σ - (2 * ρ) / 3) *
                    κ^3 -
                    (
                        15 *
                        σ^2 *
                        (9 / 80 * σ^2 * τ^2 - 3 / 4 * ρ * σ * τ + ρ^2 + 3 / 20) *
                        κ^2
                    ) / 2 + 9 * (-σ * τ / 8 + ρ) * σ^3 * κ / 4 - (9 * σ^4) / 32
                ) * shk +
                κ * (
                    ρ^2 * τ^2 * (-3 / 2 + ρ * σ * τ) * κ^5 +
                    3 *
                    τ *
                    (-ρ^2 * σ^2 * τ^2 / 2 + ρ * τ * (ρ^2 + 1) * σ - ρ^2 - 1 / 4) *
                    κ^4 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 - (9 * σ^2 * τ^2) / 8 +
                        9 * (ρ^2 + 1 / 3) * τ * ρ * σ / 2 - (9 * ρ^2) / 4
                    ) * κ^3 +
                    3 *
                    σ *
                    (
                        -σ^3 * τ^3 / 24 - (7 * ρ * τ^2 * σ^2) / 8 +
                        τ * (ρ^2 / 2 + 9 / 16) * σ +
                        ρ^3
                    ) *
                    κ^2 - 27 * (-(5 * σ * τ) / 18 + ρ) * σ^3 * τ * κ / 8 +
                    (3 * σ^4 * τ) / 8
                )
            ) * chk^2 +
            (
                (
                    (-ρ^3 * σ * τ^3 + 3 / 2 * ρ^2 * τ^2) * κ^6 -
                    3 *
                    τ *
                    (-ρ^2 * σ^2 * τ^2 / 2 + ρ * τ * (ρ^2 + 1) * σ - ρ^2 - 1 / 4) *
                    κ^5 +
                    (
                        -(3 * ρ * σ^3 * τ^3) / 4 + 3 * τ^2 * (ρ^2 + 3 / 8) * σ^2 -
                        (9 * ρ * τ * (ρ^2 + 1) * σ) / 2 +
                        3 / 8 +
                        (9 * ρ^2) / 4
                    ) * κ^4 -
                    3 *
                    (
                        -σ^3 * τ^3 / 24 +
                        ρ * τ^2 * σ^2 / 8 +
                        τ * (-(3 * ρ^2) / 2 - 3 / 16) * σ +
                        ρ^3 +
                        ρ
                    ) *
                    σ *
                    κ^3 +
                    15 *
                    (-1 / 20 * σ^2 * τ^2 - 1 / 4 * ρ * σ * τ + ρ^2 + 9 / 40) *
                    σ^2 *
                    κ^2 / 4 - (27 * σ^3 * (-σ * τ / 18 + ρ) * κ) / 16 + (15 * σ^4) / 64
                ) * shk -
                κ * (
                    (ρ^3 * σ * τ^3 - 3 * ρ^2 * τ^2) * κ^5 +
                    6 *
                    (-ρ^2 * σ^2 * τ^2 / 4 + ρ * τ * (ρ^2 + 1) * σ - 2 * ρ^2 - 1 / 2) *
                    τ *
                    κ^4 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 +
                        τ^2 * (9 * ρ^2 - 9 / 4) * σ^2 +
                        18 * ρ^3 * σ * τ - 18 * ρ^2
                    ) * κ^3 +
                    (
                        -σ^4 * τ^3 / 8 - (57 * ρ * τ^2 * σ^3) / 4 +
                        τ * (18 * ρ^2 + 45 / 4) * σ^2 +
                        24 * ρ^3 * σ
                    ) * κ^2 +
                    (33 / 8 * τ^2 * σ^4 - 45 / 2 * ρ * τ * σ^3) * κ +
                    (21 * σ^4 * τ) / 8
                ) / 8
            ) * chk +
            (
                ρ^2 * τ^2 * (ρ * σ * τ - 3) * κ^6 / 8 +
                3 *
                (-ρ^2 * σ^2 * τ^2 / 4 + ρ * τ * (ρ^2 + 1) * σ - 2 * ρ^2 - 1 / 2) *
                τ *
                κ^5 / 4 +
                (
                    (3 * ρ * σ^3 * τ^3) / 32 - 9 * (ρ^2 + 1 / 2) * τ^2 * σ^2 / 16 +
                    (9 * ρ * τ * (ρ^2 + 1) * σ) / 4 - 3 / 8 - (9 * ρ^2) / 4
                ) * κ^4 +
                3 *
                σ *
                (
                    -σ^3 * τ^3 / 192 - ρ * τ^2 * σ^2 / 32 +
                    τ * (-(3 * ρ^2) / 4 - 3 / 32) * σ +
                    ρ^3 +
                    ρ
                ) *
                κ^3 -
                (
                    15 *
                    σ^2 *
                    (-1 / 40 * σ^2 * τ^2 - 1 / 8 * ρ * σ * τ + ρ^2 + 9 / 40) *
                    κ^2
                ) / 4 + 27 * (-σ * τ / 36 + ρ) * σ^3 * κ / 16 - (15 * σ^4) / 64
            ) * shk -
            κ *
            (
                ρ^2 * τ * (-3 / 2 + ρ * σ * τ) * κ^5 +
                (-(3 * ρ^2 * σ^2 * τ^2) / 2 + 3 * ρ * τ * (ρ^2 + 1) * σ - 3 * ρ^2 - 3 / 4) *
                κ^4 +
                (
                    9 *
                    σ *
                    (ρ * τ^2 * σ^2 / 6 + τ * (-ρ^2 - 1 / 4) * σ + ρ^3 + (4 * ρ) / 3) *
                    κ^3
                ) / 2 -
                15 * (1 / 60 * σ^2 * τ^2 - 1 / 4 * ρ * σ * τ + ρ^2 + 9 / 40) * σ^2 * κ^2 /
                2 + (27 * σ^3 * (-σ * τ / 18 + ρ) * κ) / 8 - (15 * σ^4) / 32
            ) *
            τ / 2
        ) / κ^7
    c4H5 =
        -8 *
        σ^2 *
        (
            -(3 * σ^4 * (-8 * κ * θ^2 + σ^2) * chk^8) / 16 +
            (3 * σ^4 * ((-8 * κ * θ^2 + σ^2) * shk - 16 * κ * θ^2) * chk^7) / 16 +
            (
                3 *
                κ *
                σ^2 *
                (
                    2 * σ^2 * θ^2 * shk +
                    (-6 * ρ * σ * τ * θ^2 + 3 * θ^2) * κ^2 +
                    σ * (ρ * τ * σ^2 + (3 * τ * θ^2 - 1 / 2) * σ - 7 * ρ * θ^2) * κ +
                    σ^2 * (-1 / 2 * σ^2 * τ + ρ * σ + 2 * θ^2)
                ) *
                chk^6
            ) / 2 -
            (
                3 *
                σ^2 *
                (
                    (
                        (-6 * ρ * σ * τ * θ^2 + 3 * θ^2) * κ^3 +
                        σ * (ρ * τ * σ^2 + (3 * τ * θ^2 - 1 / 2) * σ - 7 * ρ * θ^2) * κ^2 +
                        σ^2 * (-1 / 2 * σ^2 * τ + ρ * σ + 5 / 2 * θ^2) * κ - σ^4 / 16
                    ) * shk -
                    10 *
                    (
                        (ρ * σ * τ - 3 / 5) * κ^2 + (-1 / 2 * σ^2 * τ + 7 / 5 * ρ * σ) * κ -
                        σ^2 / 5
                    ) *
                    κ *
                    θ^2
                ) *
                chk^5
            ) / 2 -
            3 *
            κ *
            (
                5 *
                σ^2 *
                (
                    (ρ * σ * τ - 3 / 5) * κ^2 + (-1 / 2 * σ^2 * τ + 7 / 5 * ρ * σ) * κ -
                    (3 * σ^2) / 10
                ) *
                θ^2 *
                shk - 4 * (ρ^2 * σ^2 * τ^2 - ρ * σ * τ + 1 / 8) * θ^2 * κ^4 +
                σ *
                (
                    ρ^2 * σ^3 * τ^2 +
                    (4 * ρ * τ^2 * θ^2 - ρ * τ) * σ^2 +
                    (1 / 8 + (-10 * ρ^2 - 3) * τ * θ^2) * σ +
                    5 * ρ * θ^2
                ) *
                κ^3 +
                2 *
                (
                    -ρ * τ^2 * σ^3 / 2 +
                    τ * (-τ * θ^2 / 2 + ρ^2 + 3 / 8) * σ^2 +
                    (17 / 4 * ρ * τ * θ^2 - 1 / 2 * ρ) * σ - 15 * (ρ^2 + 2 / 5) * θ^2 / 4
                ) *
                σ^2 *
                κ^2 +
                (
                    σ^6 * τ^2 / 4 - ρ * σ^5 * τ +
                    (-(3 * τ * θ^2) / 2 + (5 * ρ^2) / 4) * σ^4 +
                    7 * ρ * σ^3 * θ^2
                ) * κ - θ^2 * σ^4 - σ^6 * τ / 16
            ) *
            chk^4 +
            (
                (
                    -12 * (ρ^2 * σ^2 * τ^2 - ρ * σ * τ + 1 / 8) * θ^2 * κ^5 +
                    3 *
                    σ *
                    (
                        ρ^2 * σ^3 * τ^2 +
                        (4 * ρ * τ^2 * θ^2 - ρ * τ) * σ^2 +
                        (1 / 8 + (-10 * ρ^2 - 3) * τ * θ^2) * σ +
                        5 * ρ * θ^2
                    ) *
                    κ^4 +
                    6 *
                    σ^2 *
                    (
                        -ρ * τ^2 * σ^3 / 2 +
                        τ * (-τ * θ^2 / 2 + ρ^2 + 3 / 8) * σ^2 +
                        (5 * ρ * τ * θ^2 - 1 / 2 * ρ) * σ - 15 * (ρ^2 + 1 / 2) * θ^2 / 4
                    ) *
                    κ^3 +
                    15 *
                    (
                        σ^3 * τ^2 / 5 - ρ * τ * σ^2 +
                        (-(9 * τ * θ^2) / 5 + ρ^2 + 1 / 10) * σ +
                        7 * ρ * θ^2
                    ) *
                    σ^3 *
                    κ^2 / 4 -
                    (3 * σ^4 * (-1 / 4 * σ^2 * τ + ρ * σ + 27 / 4 * θ^2) * κ) / 4 +
                    (9 * σ^6) / 128
                ) * shk -
                (
                    27 *
                    κ *
                    (
                        (ρ^2 * σ^2 * τ^2 - 4 / 3 * ρ * σ * τ + 2 / 9) * κ^4 +
                        (
                            -ρ * τ^2 * σ^3 + τ * ((10 * ρ^2) / 3 + 1) * σ^2 -
                            (20 * ρ * σ) / 9
                        ) * κ^3 +
                        (
                            10 *
                            σ^2 *
                            (3 / 40 * σ^2 * τ^2 - 19 / 30 * ρ * σ * τ + ρ^2 + 1 / 5) *
                            κ^2
                        ) / 3 - (14 * σ^3 * (-(5 * σ * τ) / 56 + ρ) * κ) / 9 +
                        (2 * σ^4) / 9
                    ) *
                    θ^2
                ) / 2
            ) * chk^3 +
            κ *
            (
                (
                    27 *
                    θ^2 *
                    (
                        (ρ^2 * σ^2 * τ^2 - 4 / 3 * ρ * σ * τ + 2 / 9) * κ^4 +
                        (
                            -ρ * τ^2 * σ^3 + τ * ((10 * ρ^2) / 3 + 1) * σ^2 -
                            (20 * ρ * σ) / 9
                        ) * κ^3 +
                        (
                            τ^2 * σ^4 / 4 - (8 * ρ * τ * σ^3) / 3 +
                            ((10 * ρ^2) / 3 + 1) * σ^2
                        ) * κ^2 - (7 * σ^3 * (-(5 * σ * τ) / 28 + ρ) * κ) / 3 +
                        (5 * σ^4) / 12
                    ) *
                    shk
                ) / 2 +
                (-2 * ρ^3 * σ * τ^3 * θ^2 + 3 * ρ^2 * τ^2 * θ^2) * κ^6 +
                (
                    ρ^3 * σ^3 * τ^2 + 3 * τ * ρ^2 * (τ * θ^2 - 1 / 2) * σ^2 -
                    9 * (ρ^2 + 2 / 3) * τ * ρ * θ^2 * σ + 9 * (ρ^2 + 1 / 6) * θ^2
                ) *
                τ *
                κ^5 +
                (
                    -(3 * ρ^2 * σ^4 * τ^3) / 2 +
                    3 * (-τ * θ^2 / 2 + ρ^2 + 1) * τ^2 * ρ * σ^3 +
                    9 * (ρ^2 + 1 / 4) * τ * (τ * θ^2 - 1 / 3) * σ^2 -
                    18 * ρ * τ * θ^2 * (ρ^2 + 1) * σ + (9 * ρ^2 + 3) * θ^2
                ) * κ^4 +
                (
                    9 *
                    σ *
                    (
                        ρ * σ^4 * τ^3 / 6 - ((-τ * θ^2 / 9 + ρ^2 + 1 / 2) * τ^2 * σ^3) / 2 +
                        ρ * τ * (-τ * θ^2 / 3 + ρ^2 + 2 / 3) * σ^2 +
                        (τ * ((20 * ρ^2) / 3 + 5 / 4) * θ^2 - ρ^2 / 2) * σ -
                        (10 * ρ * θ^2 * (ρ^2 + 2)) / 3
                    ) *
                    κ^3
                ) / 2 +
                3 *
                σ^2 *
                (
                    -σ^4 * τ^3 / 24 - ρ * τ^2 * σ^3 / 8 -
                    ((τ * θ^2 / 4 + ρ^2 - 3 / 8) * τ * σ^2) / 2 +
                    (-27 / 8 * ρ * τ * θ^2 + ρ^3) * σ +
                    (15 * ρ^2 + 3) * θ^2
                ) *
                κ^2 +
                (-21 * ρ * σ^3 * θ^2 - 9 / 8 * ρ * σ^5 * τ + 3 / 8 * σ^6 * τ^2) * κ +
                3 * θ^2 * σ^4 +
                (3 * σ^6 * τ) / 32
            ) *
            chk^2 +
            (
                (
                    (2 * ρ^3 * σ * τ^3 * θ^2 - 3 * ρ^2 * τ^2 * θ^2) * κ^7 -
                    (
                        ρ^3 * σ^3 * τ^2 + 3 * τ * ρ^2 * (τ * θ^2 - 1 / 2) * σ^2 -
                        9 * (ρ^2 + 2 / 3) * τ * ρ * θ^2 * σ + 9 * (ρ^2 + 1 / 6) * θ^2
                    ) *
                    τ *
                    κ^6 +
                    (
                        (3 * ρ^2 * σ^4 * τ^3) / 2 -
                        3 * (-τ * θ^2 / 2 + ρ^2 + 1) * τ^2 * ρ * σ^3 +
                        (-15 * τ^2 * (ρ^2 + 3 / 20) * θ^2 + τ * (3 * ρ^2 + 3 / 4)) * σ^2 +
                        18 * (ρ^2 + 4 / 3) * τ * ρ * θ^2 * σ - 9 * (ρ^2 + 5 / 12) * θ^2
                    ) * κ^5 -
                    (
                        9 *
                        σ *
                        (
                            ρ * σ^4 * τ^3 / 6 -
                            (5 * τ^2 * (-τ * θ^2 / 15 + ρ^2 + 3 / 10) * σ^3) / 6 +
                            ρ * τ * (-(5 * τ * θ^2) / 3 + ρ^2 + 1) * σ^2 +
                            (τ * (10 * ρ^2 + 9 / 4) * θ^2 - ρ^2 / 2 - 1 / 24) * σ -
                            10 * (ρ^2 + 5 / 2) * ρ * θ^2 / 3
                        ) *
                        κ^4
                    ) / 2 -
                    3 *
                    σ^2 *
                    (
                        -σ^4 * τ^3 / 24 + (3 * ρ * τ^2 * σ^3) / 8 -
                        3 * (-τ * θ^2 / 4 + ρ^2 + 1 / 8) * τ * σ^2 / 2 +
                        ρ * (-(35 * τ * θ^2) / 4 + ρ^2 + 1 / 2) * σ +
                        75 * (ρ^2 + 27 / 100) * θ^2 / 4
                    ) *
                    κ^3 +
                    (
                        15 *
                        σ^3 *
                        (
                            -ρ * τ * σ^2 / 2 +
                            (-(21 * τ * θ^2) / 10 + ρ^2 + 3 / 20) * σ +
                            (189 * ρ * θ^2) / 10
                        ) *
                        κ^2
                    ) / 8 - (9 * σ^4 * (-1 / 6 * σ^2 * τ + ρ * σ + 65 / 6 * θ^2) * κ) / 16 +
                    (15 * σ^6) / 256
                ) * shk +
                κ *
                (
                    (ρ^3 * σ * τ^3 - 3 * ρ^2 * τ^2) * κ^6 +
                    9 *
                    (-ρ^2 * σ^2 * τ^2 / 6 + ρ * τ * (ρ^2 + 2 / 3) * σ - 2 * ρ^2 - 1 / 3) *
                    τ *
                    κ^5 +
                    (
                        (3 * ρ * σ^3 * τ^3) / 4 +
                        3 * τ^2 * (ρ^2 - 3 / 4) * σ^2 +
                        (36 * ρ^3 + 18 * ρ) * τ * σ - 36 * ρ^2 - 6
                    ) * κ^4 +
                    60 *
                    σ *
                    (
                        -σ^3 * τ^3 / 480 - (7 * ρ * τ^2 * σ^2) / 40 +
                        τ * (-ρ^2 / 4 + 3 / 80) * σ +
                        ρ^3 +
                        ρ
                    ) *
                    κ^3 -
                    90 *
                    σ^2 *
                    (-3 / 80 * σ^2 * τ^2 + 1 / 8 * ρ * σ * τ + ρ^2 + 1 / 5) *
                    κ^2 + (27 / 8 * σ^4 * τ + 42 * ρ * σ^3) * κ - 6 * σ^4
                ) *
                θ^2 / 2
            ) * chk -
            κ *
            (
                (ρ^3 * σ * τ^3 - 3 * ρ^2 * τ^2) * κ^6 +
                9 *
                (-ρ^2 * σ^2 * τ^2 / 6 + ρ * τ * (ρ^2 + 2 / 3) * σ - 2 * ρ^2 - 1 / 3) *
                τ *
                κ^5 +
                (
                    (3 * ρ * σ^3 * τ^3) / 4 - 21 * (ρ^2 + 3 / 14) * τ^2 * σ^2 / 2 +
                    36 * ρ * τ * (ρ^2 + 1) * σ - 36 * ρ^2 - 9
                ) * κ^4 +
                60 *
                (
                    -σ^3 * τ^3 / 480 +
                    ρ * τ^2 * σ^2 / 20 +
                    τ * (-ρ^2 - 3 / 16) * σ +
                    ρ^3 +
                    (3 * ρ) / 2
                ) *
                σ *
                κ^3 - 135 * σ^2 * (-19 / 90 * ρ * σ * τ + ρ^2 + 1 / 4) * κ^2 +
                (-33 / 8 * σ^4 * τ + 315 / 4 * ρ * σ^3) * κ - (105 * σ^4) / 8
            ) *
            θ^2 *
            shk / 2 + (-3 / 2 + ρ * σ * τ) * τ^2 * ρ^2 * θ^2 * κ^7 -
            (
                (
                    ρ^3 * σ^3 * τ^2 + 3 * τ * ρ^2 * (τ * θ^2 - 1 / 2) * σ^2 -
                    9 * (ρ^2 + 2 / 3) * τ * ρ * θ^2 * σ + 18 * (ρ^2 + 1 / 4) * θ^2
                ) *
                τ *
                κ^6
            ) / 2 +
            (
                (3 * ρ^2 * σ^4 * τ^3) / 4 -
                3 * (-τ * θ^2 / 2 + ρ^2 + 1) * τ^2 * ρ * σ^3 / 2 -
                9 * (τ * (ρ^2 + 1 / 8) * θ^2 - ρ^2 / 12 - 1 / 48) * τ * σ^2 +
                15 * (ρ^2 + 17 / 10) * τ * ρ * θ^2 * σ +
                9 * (ρ^2 + 1 / 6) * θ^2
            ) * κ^5 -
            (
                3 *
                σ *
                (
                    ρ * σ^4 * τ^3 / 4 - 3 * (-τ * θ^2 / 18 + ρ^2 + 1 / 4) * τ^2 * σ^3 / 2 +
                    ρ * τ * (-(7 * τ * θ^2) / 2 + ρ^2 + 1) * σ^2 +
                    (τ * ((55 * ρ^2) / 2 + 57 / 8) * θ^2 - (3 * ρ^2) / 2 - 1 / 4) * σ +
                    10 * ρ * θ^2 * (ρ^2 + 1)
                ) *
                κ^4
            ) / 2 -
            3 *
            (
                -σ^4 * τ^3 / 48 + (5 * ρ * τ^2 * σ^3) / 16 -
                5 * (-τ * θ^2 / 2 + ρ^2 + 3 / 20) * τ * σ^2 / 8 +
                ρ * (-(141 * τ * θ^2) / 16 + ρ^2 + 1) * σ - 15 * (ρ^2 + 1 / 5) * θ^2 / 2
            ) *
            σ^2 *
            κ^3 +
            (
                15 *
                σ^3 *
                (
                    σ^3 * τ^2 / 40 - (3 * ρ * τ * σ^2) / 20 +
                    (-(49 * τ * θ^2) / 40 + ρ^2 + 1 / 5) * σ - (14 * ρ * θ^2) / 5
                ) *
                κ^2
            ) / 4 - (3 * σ^4 * (-5 / 128 * σ^2 * τ + ρ * σ - θ^2) * κ) / 2 + (3 * σ^6) / 16
        ) / κ^8
    c4 = c4H3 / 2 * v0^2 + c4H4 * v0 + c4H5
    return c1, c2, c4
end
