text:
# Reference journal article:
# Nikkhoo M. and Walter T.R., 2015. Triangular dislocation: An analytical,
# artefact-free solution.
# Submitted to Geophysical Journal International
#
# Copyright (c) 2014 Mehdi Nikkhoo
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files
# (the "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to permit
# persons to whom the Software is furnished to do so, subject to the
# following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
# NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
# DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
# OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
# USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# I appreciate any comments or bug reports.
#
# Mehdi Nikkhoo
# created: 2013.1.24
# Last modified: 2014.7.30
#
# VolcanoTectonics Research Group
# Section 2.1, Physics of Earthquakes and Volcanoes
# Department 2, Physics of the Earth
# Helmholtz Centre Potsdam
# German Research Centre for Geosciences (GFZ)
#
# email:
# mehdi.nikkhoo@gfz-potsdam.de
# mehdi.nikkhoo@gmail.com
# -------------------------------------------------------
# Translated by Pengcheng Shi (shipengcheng1230@gmail.com) 08/2018

export disp_tri3_fs, disp_tri3_hs, strain_tri3_fs, strain_tri3_hs, stress_tri3_fs, stress_tri3_hs

const _ey = [0, 1, 0]
const _ez = [0, 0, 1]

"""
    disp_tri3_hs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, nu::T) where {T, V}

Compute displacement risen from triangular dislocation in elastic *halfspace*.
    Please see [original version (in supporting information)](https://academic.oup.com/gji/article/201/2/1119/572006#86405752)
    for details, especially the **coordinate system** used here.

## Arguments
- `X`, `Y`, `Z`: observational coordinates
- `P1`, `P2`, `P3`: three triangular vertices coordinates respectively
- `Ss`, `Ds`, `Ts`: triangular dislocation vector, Strike-slip, Dip-slip, Tensile-slip respectively
- `nu`: poisson ratio

## Output
By order: ``u_x``, ``u_y``, ``u_z``
"""
function disp_tri3_hs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, nu::T) where {T, V}
    @assert (Z ≤ zero(T) && P1[3] ≤ zero(T) && P2[3] ≤ zero(T) && P3[3] ≤ zero(T))  "Half-space solution: Z coordinates must be negative!"
    A = transform_matrix(P1, P2, P3)
    ueMS, unMS, uvMS = _disp_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, nu, A')
    ueFSC, unFSC, uvFSC = _td_disp_harmonic_func(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, nu, A)
    P1[3] *= -one(T)
    P2[3] *= -one(T)
    P3[3] *= -one(T)
    allatsurface = P1[3] == zero(T) && P2[3] == zero(T) && P3[3] == zero(T)
    if !allatsurface
        A *= -one(T)
        A[3,1] *= -one(T)
        A[3,2] *= -one(T)
        A[3,3] *= -one(T)
    end
    ueIS, unIS, uvIS = _disp_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, nu, A')
    uvIS = ifelse(allatsurface, -uvIS, uvIS)
    ue = ueMS + ueIS + ueFSC
    un = unMS + unIS + unFSC
    uv = uvMS + uvIS + uvFSC
    allatsurface && begin ue *= -one(T); un *= -one(T); uv *= -one(T) end
    P1[3] *= -one(T)
    P2[3] *= -one(T)
    P3[3] *= -one(T)
    return ue, un, uv
end

"""
    disp_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, nu::T) where {T, V}

Compute displacement risen from triangular dislocation in elastic *fullspace*.
    Please see [original version (in supporting information)](https://academic.oup.com/gji/article/201/2/1119/572006#86405752)
    for details, especially the **coordinate system** used here.

## Arguments
The same as [`disp_tri3_hs`](@ref)
"""
function disp_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, nu::T) where {T, V}
    A = transform_matrix(P1, P2, P3)
    _disp_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, nu, A')
end

@inline function transform_matrix(P1::V, P2::V, P3::V) where V<:AbstractVector{T} where T
    Vnorm = cross(P2 - P1, P3 - P1)
    normalize!(Vnorm)
    Vstrike = cross(_ez, Vnorm)

    if norm(Vstrike) == zero(T)
        Vstrike = _ey * Vnorm[3]
        if P1[3] > zero(T)
            Vstrike *= -one(T)
        end
    end
    normalize!(Vstrike)
    Vdip = Vnorm × Vstrike
    hcat(Vnorm, Vstrike, Vdip)
end

@inline function _disp_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, nu::T, At::M) where {T, V, M}
    bx, by, bz = Ts, Ss, Ds
    p1, p2, p3 = [zeros(T, 3) for _ in 1: 3]

    x, y, z = coord_trans(X - P2[1], Y - P2[2], Z - P2[3], At)
    p1[1], p1[2], p1[3] = coord_trans(P1[1] - P2[1], P1[2] - P2[2], P1[3] - P2[3], At)
    p3[1], p3[2], p3[3] = coord_trans(P3[1] - P2[1], P3[2] - P2[2], P3[3] - P2[3], At)

    e12 = normalize(p2 - p1)
    e13 = normalize(p3 - p1)
    e23 = normalize(p3 - p2)

    A = acos(dot(e12, e13))
    B = acos(dot(-e12, e23))
    C = acos(dot(e23, e13))

    Trimode = trimodefinder(y, z, x, p1, p2, p3)

    if Trimode == 1
        u1Tp, v1Tp, w1Tp = TDSetupD(x, y, z, A, bx, by, bz, nu, p1, -e13)
        u2Tp, v2Tp, w2Tp = TDSetupD(x, y, z, B, bx, by, bz, nu, p2, e12)
        u3Tp, v3Tp, w3Tp = TDSetupD(x, y, z, C, bx, by, bz, nu, p3, e23)
        u = u1Tp + u2Tp + u3Tp
        v = v1Tp + v2Tp + v3Tp
        w = w1Tp + w2Tp + w3Tp
    end

    if Trimode == -1
        u1Tn, v1Tn, w1Tn = TDSetupD(x, y, z, A, bx, by, bz, nu, p1, e13)
        u2Tn, v2Tn, w2Tn = TDSetupD(x, y, z, B, bx, by, bz, nu, p2, -e12)
        u3Tn, v3Tn, w3Tn = TDSetupD(x, y, z, C, bx, by, bz, nu, p3, -e23)
        u = u1Tn + u2Tn + u3Tn
        v = v1Tn + v2Tn + v3Tn
        w = w1Tn + w2Tn + w3Tn
    end

    if Trimode == 0
        u = NaN
        v = NaN
        w = NaN
    end

    a1, a2, a3 = -x, p1[2] - y, p1[3] - z
    b1, b2, b3 = -x, -y, -z
    c1, c2, c3 = -x, p3[2] - y, p3[3] - z
    na = hypot(a1, a2, a3)
    nb = hypot(b1, b2, b3)
    nc = hypot(c1, c2, c3)

    Fiy = a1 * (b2 * c3 - b3 * c2)- a2 * (b1 * c3 - b3 * c1) + a3 * (b1 * c2 - b2 * c1)
    Fix = na * nb * nc + (a1 * b1 + a2 * b2 + a3 * b3) * nc + (a1 * c1 + a2 * c2 + a3 * c3) * nb + (b1 * c1 + b2 * c2 + b3 * c3) * na
    Fi = -2 * atan(Fiy, Fix) / 4 / π
    u += bx * Fi
    v += by * Fi
    w += bz * Fi
    ue, un, uv = coord_trans(u, v, w, At')
end

@inline function _td_disp_harmonic_func(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, nu::T, A::M) where {T, V, M}
    bx, by, bz = Ts, Ss, Ds
    bX, bY, bZ = coord_trans(bx, by, bz, A)
    u1, v1, w1 = AngSetupFSC(X, Y, Z, bX, bY, bZ, P1, P2, nu)
    u2, v2, w2 = AngSetupFSC(X, Y, Z, bX, bY, bZ, P2, P3, nu)
    u3, v3, w3 = AngSetupFSC(X, Y, Z, bX, bY, bZ, P3, P1, nu)
    return u1 + u2 + u3, v1 + v2 + v3, w1 + w2 + w3
end

@inline function AngSetupFSC(X::T, Y::T, Z::T, bX::T, bY::T, bZ::T, PA::V, PB::V, nu::T) where {T, V}
    sv1, sv2, sv3 = PB[1] - PA[1], PB[2] - PA[2], PB[3] - PA[3]
    svr = hypot(sv1, sv2, sv3)
    beta = acos(-sv3 / svr)
    if beta ≈ zero(T) || beta ≈ π
        return zero(T), zero(T), zero(T)
    else
        ey1 = [sv1, sv2, zero(T)]
        normalize!(ey1)
        ey3 = -_ez
        ey2 = cross(ey3, ey1)
        A = hcat(ey1, ey2, ey3)
        y1A, y2A, y3A = coord_trans(X - PA[1], Y - PA[2], Z - PA[3], A)
        y1AB, y2AB, y3AB = coord_trans(sv1, sv2, sv3, A)
        y1B = y1A - y1AB
        y2B = y2A - y2AB
        y3B = y3A - y3AB
        b1, b2, b3 = coord_trans(bX, bY, bZ, A)
        angle = ifelse(beta * y1A ≥ zero(T), -π + beta, beta)
        v1A, v2A, v3A = AngDisDispFSC(y1A, y2A, y3A, angle, b1, b2, b3, nu, -PA[3])
        v1B, v2B, v3B = AngDisDispFSC(y1B, y2B, y3B, angle, b1, b2, b3, nu, -PB[3])
        v1 = v1B - v1A
        v2 = v2B - v2A
        v3 = v3B - v3A
        ue, un, uv = coord_trans(v1, v2, v3, A')
        return ue, un, uv
    end
end

@inline function AngDisDispFSC(y1::T, y2::T, y3::T, beta::T, b1::T, b2::T, b3::T, nu::T, a::T) where T
    sinB, cosB = sincos(beta)
    cotB = cot(beta)
    y3b = y3 + 2a
    z1b = y1 * cosB + y3b * sinB
    z3b = -y1 * sinB + y3b * cosB
    rb = hypot(y1, y2, y3b)
    Fib = 2 * atan(-y2 / (-(rb + y3b) * cot(beta/2) + y1))

    v1cb1 = (b1 / 4 / pi / (1 - nu) * (-2 * (1 - nu) * (1 - 2 * nu) * Fib * cotB ^ 2 + (1 - 2 * nu) * y2 /
        (rb + y3b) * ((1 - 2 * nu - a / rb) * cotB - y1 / (rb + y3b) * (nu + a / rb)) + (1 - 2 * nu) *
        y2 * cosB * cotB / (rb + z3b) * (cosB + a / rb) + a * y2 * (y3b - a) * cotB / rb ^ 3 + y2 *
        (y3b - a) / (rb * (rb + y3b)) * (-(1 - 2 * nu) * cotB + y1 / (rb + y3b) * (2 * nu + a / rb) +
        a * y1 / rb ^ 2) +y2 * (y3b - a) / (rb * (rb + z3b)) * (cosB / (rb + z3b) * ((rb *
        cosB + y3b) * ((1 - 2 * nu) * cosB - a / rb) * cotB + 2 * (1 - nu) * (rb * sinB - y1) * cosB) -
        a * y3b * cosB * cotB / rb ^ 2)))

    v2cb1 = ((b1 / 4 / pi / (1 - nu) * ((1 - 2 * nu) * ((2 * (1 - nu) * cotB ^ 2 - nu) * log(rb + y3b) - (2 *
        (1 - nu) * cotB ^ 2 + 1 - 2 * nu)* cosB * log(rb + z3b)) -(1 - 2 * nu) / (rb + y3b) * (y1 *
        cotB * (1 - 2 * nu - a / rb) + nu * y3b - a + y2 ^ 2 / (rb + y3b) * (nu + a / rb)) -(1 - 2 *
        nu) * z1b * cotB / (rb + z3b) * (cosB + a / rb) - a * y1 * (y3b - a) * cotB / rb ^ 3 +
        (y3b - a) / (rb + y3b) * (-2 * nu + 1 / rb * ((1 - 2 * nu) * y1 * cotB - a) + y2 ^ 2 / (rb *
        (rb + y3b)) * (2 * nu + a / rb) + a * y2 ^ 2 / rb ^ 3) +(y3b - a) / (rb + z3b) * (cosB ^ 2 -
        1 / rb * ((1 - 2 * nu) * z1b * cotB + a * cosB) + a * y3b * z1b * cotB / rb ^ 3 - 1 / (rb *
        (rb + z3b)) * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB / rb * (rb * cosB + y3b))))))

    v3cb1 = (b1 / 4 / pi / (1 - nu) * (2 * (1 - nu) * (((1 - 2 * nu) * Fib * cotB) + (y2 / (rb + y3b) * (2 *
        nu + a / rb)) -(y2 * cosB / (rb + z3b) * (cosB + a / rb))) +y2 * (y3b - a) / rb * (2 *
        nu / (rb + y3b) + a / rb ^ 2) +y2 * (y3b - a) * cosB / (rb * (rb + z3b)) * (1 - 2 * nu -
        (rb * cosB + y3b) / (rb + z3b) * (cosB + a / rb) - a * y3b / rb ^ 2)))

    v1cb2 = (b2 / 4 / pi / (1 - nu) * ((1 - 2 * nu) * ((2 * (1 - nu) * cotB ^ 2 + nu) * log(rb + y3b) - (2 *
        (1 - nu) * cotB ^ 2 + 1)* cosB * log(rb + z3b)) +(1 - 2 * nu) / (rb + y3b) * (-(1 - 2 * nu) *
        y1 * cotB + nu * y3b - a + a * y1 * cotB / rb + y1 ^ 2 / (rb + y3b) * (nu + a / rb)) -(1 - 2 *
        nu)* cotB / (rb + z3b) * (z1b * cosB - a * (rb * sinB - y1) / (rb * cosB)) - a * y1 *
        (y3b - a) * cotB / rb ^ 3 + (y3b - a) / (rb + y3b) * (2 * nu + 1 / rb * ((1 - 2 * nu) * y1 *
        cotB + a) -y1 ^ 2 / (rb * (rb + y3b)) * (2 * nu + a / rb) - a * y1 ^ 2 / rb ^ 3) +(y3b - a) *
        cotB / (rb + z3b) * (-cosB * sinB + a * y1 * y3b / (rb ^ 3 * cosB) + (rb * sinB - y1) /
        rb * (2 * (1 - nu) * cosB - (rb * cosB + y3b) / (rb + z3b) * (1 + a / (rb * cosB))))))

    v2cb2 = (b2 / 4 / pi / (1 - nu) * (2 * (1 - nu) * (1 - 2 * nu) * Fib * cotB ^ 2 + (1 - 2 * nu) * y2 /
        (rb + y3b) * (-(1 - 2 * nu - a / rb) * cotB + y1 / (rb + y3b) * (nu + a / rb)) - (1 - 2 * nu) *
        y2 * cotB / (rb + z3b) * (1 + a / (rb * cosB)) - a * y2 * (y3b - a) * cotB / rb ^ 3 + y2 *
        (y3b - a) / (rb * (rb + y3b)) * ((1 - 2 * nu) * cotB - 2 * nu * y1 / (rb + y3b) - a * y1 / rb *
        (1 / rb + 1 / (rb + y3b))) +y2 * (y3b - a) * cotB / (rb * (rb + z3b)) * (-2 * (1 - nu) *
        cosB + (rb * cosB + y3b) / (rb + z3b) * (1 + a / (rb * cosB)) + a * y3b / (rb ^ 2 * cosB))))

    v3cb2 = (b2 / 4 / pi / (1 - nu) * (-2 * (1 - nu) * (1 - 2 * nu) * cotB * (log(rb + y3b) - cosB *
        log(rb + z3b)) -2 * (1 - nu) * y1 / (rb + y3b) * (2 * nu + a / rb) + 2 * (1 - nu) * z1b / (rb +
        z3b) * (cosB + a / rb) + (y3b - a) / rb * ((1 - 2 * nu) * cotB - 2 * nu * y1 / (rb + y3b) - a *
        y1 / rb ^ 2) -(y3b - a) / (rb + z3b) * (cosB * sinB + (rb * cosB + y3b) * cotB / rb *
        (2 * (1 - nu) * cosB - (rb * cosB + y3b) / (rb + z3b)) + a / rb * (sinB - y3b * z1b /
        rb ^ 2 - z1b * (rb * cosB + y3b) / (rb * (rb + z3b))))))

    v1cb3 = (b3 / 4 / pi / (1 - nu) * ((1 - 2 * nu) * (y2 / (rb + y3b) * (1 + a / rb) - y2 * cosB / (rb +
        z3b) * (cosB + a / rb)) -y2 * (y3b - a) / rb * (a / rb ^ 2 + 1 / (rb + y3b)) + y2 *
        (y3b - a) * cosB / (rb * (rb + z3b)) * ((rb * cosB + y3b) / (rb + z3b) * (cosB + a /
        rb) +a * y3b / rb ^ 2)))

    v2cb3 = (b3 / 4 / pi / (1 - nu) * ((1 - 2 * nu) * (-sinB * log(rb + z3b) - y1 / (rb + y3b) * (1 + a /
        rb) +z1b / (rb + z3b) * (cosB + a / rb)) +y1 * (y3b - a) / rb * (a / rb ^ 2 + 1 / (rb +
        y3b)) -(y3b - a) / (rb + z3b) * (sinB * (cosB - a / rb) + z1b / rb * (1 + a * y3b /
        rb ^ 2) -1 / (rb * (rb + z3b)) * (y2 ^ 2 * cosB * sinB - a * z1b / rb * (rb * cosB + y3b)))))

    v3cb3 = (b3 / 4 / pi / (1 - nu) * (2 * (1 - nu) * Fib + 2 * (1 - nu) * (y2 * sinB / (rb + z3b) * (cosB +
        a / rb)) +y2 * (y3b - a) * sinB / (rb * (rb + z3b)) * (1 + (rb * cosB + y3b) / (rb +
        z3b) * (cosB + a / rb) + a * y3b / rb ^ 2)))

    return v1cb1 + v1cb2 + v1cb3, v2cb1 + v2cb2 + v2cb3, v3cb1 + v3cb2 + v3cb3
end

@inline function trimodefinder(x::T, y::T, z::T, p1::V, p2::V, p3::V) where {T, V}
    a = ((p2[3] - p3[3]) * (x - p3[2]) + (p3[2] - p2[2]) * (y - p3[3])) /
        ((p2[3] - p3[3]) * (p1[2] - p3[2]) + (p3[2] - p2[2]) * (p1[3] - p3[3]))
    b = ((p3[3] - p1[3]) * (x - p3[2]) + (p1[2] - p3[2]) * (y - p3[3])) /
        ((p2[3] - p3[3]) * (p1[2] - p3[2]) + (p3[2] - p2[2]) * (p1[3] - p3[3]))
    c = 1 - a - b
    trimode = 1
    if ((a ≤ 0) && (b > c) && (c > a)) || ((b ≤ 0) && (c > a) && (a > b)) || ((c ≤ 0) && (a > b) && (b > c))
        trimode = -1
    end
    if ((a == 0) && (b ≥ 0) && (c ≥ 0)) || ((a ≥ 0) && (b == 0) && (c ≥ 0)) || ((a ≥ 0) && (b ≥ 0) && (c == 0))
        trimode = 0
    end
    if (trimode == 0) && (z ≠ 0)
        trimode = 1
    end
    return trimode
end

@inline function TDSetupD(x::T, y::T, z::T, alpha::T, bx::T, by::T, bz::T, nu::T, TriVertex::V, SideVec::V) where {T, V}
    y1 = SideVec[3] * (y - TriVertex[2]) - SideVec[2] * (z - TriVertex[3])
    z1 = SideVec[2] * (y - TriVertex[2]) + SideVec[3] * (z - TriVertex[3])
    by1 = SideVec[3] * by - SideVec[2] * bz
    bz1 = SideVec[2] * by + SideVec[3] * bz
    u, v0, w0 = AngDisDisp(x, y1, z1, -π + alpha, bx, by1, bz1, nu)
    v = SideVec[3] * v0 + SideVec[2] * w0
    w = -SideVec[2] * v0 + SideVec[3] * w0
    return u, v, w
end

@inline function _TDstrain_HarFunc(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T, At::M) where {T, V, M}
    bx, by, bz = Ts, Ss, Ds
    bX, bY, bZ = coord_trans(bx, by, bz, At)
    exx1, eyy1, ezz1, exy1, exz1, eyz1 = AngSetupFSC_S(X, Y, Z, bX, bY, bZ, P1, P2, λ, μ)
    exx2, eyy2, ezz2, exy2, exz2, eyz2 = AngSetupFSC_S(X, Y, Z, bX, bY, bZ, P2, P3, λ, μ)
    exx3, eyy3, ezz3, exy3, exz3, eyz3 = AngSetupFSC_S(X, Y, Z, bX, bY, bZ, P3, P1, λ, μ)
    return exx1 + exx2 + exx3, eyy1 + eyy2 + eyy3, ezz1 + ezz2 + ezz3, exy1 + exy2 + exy3, exz1 + exz2 + exz3, eyz1 + eyz2 + eyz3
end

@inline function AngSetupFSC_S(X::T, Y::T, Z::T, bX::T, bY::T, bZ::T, PA::V, PB::V, λ::T, μ::T) where {T, V}
    nu = λ / (μ + λ) / 2
    sv1, sv2, sv3 = PB[1] - PA[1], PB[2] - PA[2], PB[3] - PA[3]
    svr = hypot(sv1, sv2, sv3)
    beta = acos(-sv3 / svr)
    if beta ≈ zero(T) || beta ≈ π
        return zero(T), zero(T), zero(T), zero(T), zero(T), zero(T)
    else
        ey1 = [sv1, sv2, zero(T)]
        normalize!(ey1)
        ey3 = -_ez
        ey2 = cross(ey3, ey1)
        A = hcat(ey1, ey2, ey3)
        y1A, y2A, y3A = coord_trans(X - PA[1], Y - PA[2], Z - PA[3], A)
        y1AB, y2AB, y3AB = coord_trans(sv1, sv2, sv3, A)
        y1B = y1A - y1AB
        y2B = y2A - y2AB
        y3B = y3A - y3AB
        b1, b2, b3 = coord_trans(bX, bY, bZ, A)
        if beta * y1A ≥ zero(T)
            v11A, v22A, v33A, v12A, v13A, v23A = AngDisStrainFSC(-y1A, -y2A, y3A, π - beta, -b1, -b2, b3, nu, -PA[3])
            v13A *= -one(T)
            v23A *= -one(T)
            v11B, v22B, v33B, v12B, v13B, v23B = AngDisStrainFSC(-y1B, -y2B, y3B, π - beta, -b1, -b2, b3, nu, -PB[3])
            v13B *= -one(T)
            v23B *= -one(T)
        else
            v11A, v22A, v33A, v12A, v13A, v23A = AngDisStrainFSC(y1A, y2A, y3A, beta, b1, b2, b3, nu, -PA[3])
            v11B, v22B, v33B, v12B, v13B, v23B = AngDisStrainFSC(y1B, y2B, y3B, beta, b1, b2, b3, nu, -PB[3])
        end
        v11 = v11B - v11A
        v22 = v22B - v22A
        v33 = v33B - v33A
        v12 = v12B - v12A
        v13 = v13B - v13A
        v23 = v23B - v23A
        Exx, Eyy, Ezz, Exy, Exz, Eyz = TensTrans(v11, v22, v33, v12, v13, v23, A')
        return Exx, Eyy, Ezz, Exy, Exz, Eyz
    end
end

@inline function AngDisStrainFSC(y1::T, y2::T, y3::T, beta::T, b1::T, b2::T, b3::T, nu::T, a::T) where T
    sinB, cosB = sincos(beta)
    cotB = cot(beta)
    y3b = y3 + 2a
    z1b = y1 * cosB + y3b * sinB
    z3b = -y1 * sinB + y3b * cosB
    rb2 = y1^2 + y2^2 + y3b^2
    rb = sqrt(rb2)
    W1 = rb * cosB + y3b
    W2 = cosB + a / rb
    W3 = cosB + y3b / rb
    W4 = nu + a / rb
    W5 = 2nu + a / rb
    W6 = rb + y3b
    W7 = rb + z3b
    W8 = y3 + a
    W9 = 1 + a / rb / cosB
    N1 = 1-2*nu

    rFib_ry2 = z1b / rb / (rb + z3b) - y1 / rb / (rb + y3b) # y2 = x in ADCS
    rFib_ry1 = y2 / rb / (rb + y3b) - cosB * y2 / rb / (rb + z3b) # y1 =y in ADCS
    rFib_ry3 = -sinB * y2 / rb / (rb + z3b) # y3 = z in ADCS

    v11 = (b1 * (1 / 4 * ((-2 + 2 * nu) * N1 * rFib_ry1 * cotB ^ 2 - N1 * y2 / W6 ^ 2 * ((1 - W5) * cotB -
        y1 / W6 * W4) /rb *y1+N1 *y2 /W6 * (a / rb ^ 3 * y1 * cotB - 1 / W6 * W4 + y1 ^ 2 /
        W6 ^ 2 * W4 / rb + y1 ^ 2 / W6 * a / rb ^ 3) -N1 * y2 * cosB * cotB / W7 ^ 2 * W2 * (y1 /
        rb - sinB) -N1 * y2 * cosB * cotB / W7 * a / rb ^ 3 * y1 - 3 * a * y2 * W8 * cotB / rb ^ 5. *
        y1 - y2 * W8 / rb ^ 3 / W6 * (-N1 * cotB + y1 / W6 * W5 + a * y1 / rb2) * y1 - y2 * W8 /
        rb2 / W6 ^ 2 * (-N1 * cotB + y1 / W6 * W5 + a * y1 / rb2) * y1 + y2 * W8 / rb / W6 *
        (1 / W6 * W5 - y1 ^ 2 / W6 ^ 2 * W5 / rb - y1 ^ 2 / W6 * a / rb ^ 3 + a / rb2 - 2 * a * y1 ^
        2 / rb2 ^ 2) -y2 * W8 / rb ^ 3 / W7 * (cosB / W7 * (W1 * (N1 * cosB - a / rb) * cotB +
        (2 - 2 * nu) * (rb * sinB - y1) * cosB) -a * y3b * cosB * cotB / rb2) * y1 - y2 * W8 / rb /
        W7 ^ 2 * (cosB / W7 * (W1 * (N1 * cosB - a / rb) * cotB + (2 - 2 * nu) * (rb * sinB - y1) *
        cosB) -a * y3b * cosB * cotB / rb2) * (y1 / rb - sinB) + y2 * W8 / rb / W7 * (-cosB /
        W7 ^ 2 * (W1 * (N1 * cosB - a / rb) * cotB + (2 - 2 * nu) * (rb * sinB - y1) * cosB) * (y1 /
        rb - sinB) +cosB / W7 * (1 / rb * cosB * y1 * (N1 * cosB - a / rb) * cotB + W1 * a / rb ^
        3 * y1 * cotB + (2 - 2 * nu) * (1 / rb * sinB * y1 - 1) * cosB) +2 * a * y3b * cosB * cotB /
        rb2 ^ 2 * y1)) /π/(1 - nu)) +
        b2 * (1 / 4 * (N1 * (((2 - 2 * nu) * cotB ^ 2 + nu) / rb * y1 / W6 - ((2 - 2 * nu) * cotB ^ 2 + 1) *
        cosB * (y1 / rb - sinB) / W7) -N1 / W6 ^ 2 * (-N1 * y1 * cotB + nu * y3b - a + a * y1 *
        cotB / rb + y1 ^ 2 / W6 * W4) /rb *y1+N1 /W6 * (-N1 * cotB + a * cotB / rb - a *
        y1 ^ 2 * cotB / rb ^ 3 + 2 * y1 / W6 * W4 - y1 ^ 3 / W6 ^ 2 * W4 / rb - y1 ^ 3 / W6 * a /
        rb ^ 3) +N1 * cotB / W7 ^ 2 * (z1b * cosB - a * (rb * sinB - y1) / rb / cosB) * (y1 /
        rb - sinB) -N1 * cotB / W7 * (cosB ^ 2 - a * (1 / rb * sinB * y1 - 1) / rb / cosB + a *
        (rb * sinB - y1) / rb ^ 3 / cosB * y1) -a * W8 * cotB / rb ^ 3 + 3 * a * y1 ^ 2 * W8 *
        cotB / rb ^ 5 - W8 / W6 ^ 2 * (2 * nu + 1 / rb * (N1 * y1 * cotB + a) - y1 ^ 2 / rb / W6 *
        W5 - a * y1 ^ 2 / rb ^ 3) /rb *y1+W8 /W6 * (-1 / rb ^ 3 * (N1 * y1 * cotB + a) * y1 +
        1 / rb * N1 * cotB - 2 * y1 / rb / W6 * W5 + y1 ^ 3 / rb ^ 3 / W6 * W5 + y1 ^ 3 / rb2 /
        W6 ^ 2 * W5 + y1 ^ 3 / rb2 ^ 2 / W6 * a - 2 * a / rb ^ 3 * y1 + 3 * a * y1 ^ 3 / rb ^ 5) -W8 *
        cotB / W7 ^ 2 * (-cosB * sinB + a * y1 * y3b / rb ^ 3 / cosB + (rb * sinB - y1) / rb *
        ((2 - 2 * nu) * cosB - W1 / W7 * W9)) * (y1 / rb - sinB) + W8 * cotB / W7 * (a * y3b /
        rb ^ 3 / cosB - 3 * a * y1 ^ 2 * y3b / rb ^ 5. / cosB + (1 / rb * sinB * y1 - 1) / rb *
        ((2 - 2 * nu) * cosB - W1 / W7 * W9) - (rb * sinB - y1) / rb ^ 3 * ((2 - 2 * nu) * cosB - W1 /
        W7 * W9) * y1 + (rb * sinB - y1) / rb * (-1 / rb * cosB * y1 / W7 * W9 + W1 / W7 ^ 2 *
        W9 * (y1 / rb - sinB) + W1 / W7 * a / rb ^ 3 / cosB * y1))) /π/(1 - nu)) +
        b3 * (1 / 4 * (N1 * (-y2 / W6 ^ 2 * (1 + a / rb) / rb * y1 - y2 / W6 * a / rb ^ 3 * y1 + y2 *
        cosB / W7 ^ 2 * W2 * (y1 / rb - sinB) + y2 * cosB / W7 * a / rb ^ 3 * y1) +y2 * W8 /
        rb ^ 3 * (a / rb2 + 1 / W6) * y1 - y2 * W8 / rb * (-2 * a / rb2 ^ 2 * y1 - 1 / W6 ^ 2 /
        rb * y1) -y2 * W8 * cosB / rb ^ 3 / W7 * (W1 / W7 * W2 + a * y3b / rb2) * y1 - y2 * W8 *
        cosB / rb / W7 ^ 2 * (W1 / W7 * W2 + a * y3b / rb2) * (y1 / rb - sinB) + y2 * W8 *
        cosB / rb / W7 * (1 / rb * cosB * y1 / W7 * W2 - W1 / W7 ^ 2 * W2 * (y1 / rb - sinB) -
        W1 / W7 * a / rb ^ 3 * y1 - 2 * a * y3b / rb2 ^ 2 * y1)) /π/(1 - nu)))

    v22 = (b1 * (1 / 4 * (N1 * (((2 - 2 * nu) * cotB ^ 2 - nu) / rb * y2 / W6 - ((2 - 2 * nu) * cotB ^ 2 + 1 -
        2 * nu)* cosB / rb * y2 / W7) +N1 / W6 ^ 2 * (y1 * cotB * (1 - W5) + nu * y3b - a + y2 ^
        2 / W6 * W4) /rb *y2-N1 /W6 * (a * y1 * cotB / rb ^ 3 * y2 + 2 * y2 / W6 * W4 - y2 ^
        3 / W6 ^ 2 * W4 / rb - y2 ^ 3 / W6 * a / rb ^ 3) +N1 * z1b * cotB / W7 ^ 2 * W2 / rb *
        y2 + N1 * z1b * cotB / W7 * a / rb ^ 3 * y2 + 3 * a * y2 * W8 * cotB / rb ^ 5. * y1 - W8 /
        W6 ^ 2 * (-2 * nu + 1 / rb * (N1 * y1 * cotB - a) + y2 ^ 2 / rb / W6 * W5 + a * y2 ^ 2 /
        rb ^ 3) /rb *y2+W8 /W6 * (-1 / rb ^ 3 * (N1 * y1 * cotB - a) * y2 + 2 * y2 / rb /
        W6 * W5 - y2 ^ 3 / rb ^ 3 / W6 * W5 - y2 ^ 3 / rb2 / W6 ^ 2 * W5 - y2 ^ 3 / rb2 ^ 2 / W6 *
        a + 2 * a / rb ^ 3 * y2 - 3 * a * y2 ^ 3 / rb ^ 5) -W8 / W7 ^ 2 * (cosB ^ 2 - 1 / rb * (N1 *
        z1b * cotB + a * cosB) +a * y3b * z1b * cotB / rb ^ 3 - 1 / rb / W7 * (y2 ^ 2 * cosB ^ 2 -
        a * z1b * cotB / rb * W1)) /rb *y2+W8 /W7 * (1 / rb ^ 3 * (N1 * z1b * cotB + a *
        cosB) * y2 - 3 * a * y3b * z1b * cotB / rb ^ 5. * y2 + 1 / rb ^ 3 / W7 * (y2 ^ 2 * cosB ^ 2 -
        a * z1b * cotB / rb * W1) * y2 + 1 / rb2 / W7 ^ 2 * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB /
        rb * W1) * y2 - 1 / rb / W7 * (2 * y2 * cosB ^ 2 + a * z1b * cotB / rb ^ 3 * W1 * y2 - a *
        z1b * cotB / rb2 * cosB * y2))) /π/(1 - nu)) +
        b2 * (1 / 4 * ((2 - 2 * nu) * N1 * rFib_ry2 * cotB ^ 2 + N1 / W6 * ((W5 - 1) * cotB + y1 / W6 *
        W4) -N1 * y2 ^ 2 / W6 ^ 2 * ((W5 - 1) * cotB + y1 / W6 * W4) / rb + N1 * y2 / W6 * (-a /
        rb ^ 3 * y2 * cotB - y1 / W6 ^ 2 * W4 / rb * y2 - y2 / W6 * a / rb ^ 3 * y1) -N1 * cotB /
        W7 * W9 + N1 * y2 ^ 2 * cotB / W7 ^ 2 * W9 / rb + N1 * y2 ^ 2 * cotB / W7 * a / rb ^ 3 /
        cosB - a * W8 * cotB / rb ^ 3 + 3 * a * y2 ^ 2 * W8 * cotB / rb ^ 5 + W8 / rb / W6 * (N1 *
        cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 / W6)) -y2 ^ 2 * W8 / rb ^ 3 / W6 *
        (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 / W6)) - y2 ^ 2 * W8 / rb2 / W6 ^
        2 * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 / W6)) + y2 * W8 / rb / W6 *
        (2 * nu * y1 / W6 ^ 2 / rb * y2 + a * y1 / rb ^ 3 * (1 / rb + 1 / W6) * y2 - a * y1 / rb *
        (-1 / rb ^ 3 * y2 - 1 / W6 ^ 2 / rb * y2)) +W8 * cotB / rb / W7 * ((-2 + 2 * nu) * cosB +
        W1 / W7 * W9 + a * y3b / rb2 / cosB) -y2 ^ 2 * W8 * cotB / rb ^ 3 / W7 * ((-2 + 2 * nu) *
        cosB + W1 / W7 * W9 + a * y3b / rb2 / cosB) -y2 ^ 2 * W8 * cotB / rb2 / W7 ^ 2 * ((-2 +
        2 * nu)* cosB + W1 / W7 * W9 + a * y3b / rb2 / cosB) +y2 * W8 * cotB / rb / W7 * (1 /
        rb * cosB * y2 / W7 * W9 - W1 / W7 ^ 2 * W9 / rb * y2 - W1 / W7 * a / rb ^ 3 / cosB * y2 -
        2 * a * y3b / rb2 ^ 2 / cosB * y2)) /π/(1 - nu)) +
        b3 * (1 / 4 * (N1 * (-sinB / rb * y2 / W7 + y2 / W6 ^ 2 * (1 + a / rb) / rb * y1 + y2 / W6 *
        a / rb ^ 3 * y1 - z1b / W7 ^ 2 * W2 / rb * y2 - z1b / W7 * a / rb ^ 3 * y2) -y2 * W8 /
        rb ^ 3 * (a / rb2 + 1 / W6) * y1 + y1 * W8 / rb * (-2 * a / rb2 ^ 2 * y2 - 1 / W6 ^ 2 /
        rb * y2) +W8 / W7 ^ 2 * (sinB * (cosB - a / rb) + z1b / rb * (1 + a * y3b / rb2) - 1 /
        rb / W7 * (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1)) /rb *y2-W8 /W7 * (sinB * a /
        rb ^ 3 * y2 - z1b / rb ^ 3 * (1 + a * y3b / rb2) * y2 - 2 * z1b / rb ^ 5 * a * y3b * y2 +
        1 / rb ^ 3 / W7 * (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1) * y2 + 1 / rb2 / W7 ^ 2 *
        (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1) * y2 - 1 / rb / W7 * (2 * y2 * cosB * sinB + a *
        z1b / rb ^ 3 * W1 * y2 - a * z1b / rb2 * cosB * y2))) /π/(1 - nu)))

    v33 = (b1 * (1 / 4 * ((2 - 2 * nu) * (N1 * rFib_ry3 * cotB - y2 / W6 ^ 2 * W5 * (y3b / rb + 1) -
        1 / 2 * y2 / W6 * a / rb ^ 3 * 2 * y3b + y2 * cosB / W7 ^ 2 * W2 * W3 + 1 / 2 * y2 * cosB / W7 *
        a / rb ^ 3 * 2 * y3b) +y2 / rb * (2 * nu / W6 + a / rb2) - 1 / 2 * y2 * W8 / rb ^ 3 * (2 *
        nu / W6 + a / rb2)* 2 * y3b + y2 * W8 / rb * (-2 * nu / W6 ^ 2 * (y3b / rb + 1) - a /
        rb2 ^ 2 * 2 * y3b) +y2 * cosB / rb / W7 * (1 - 2 * nu - W1 / W7 * W2 - a * y3b / rb2) -
        1 / 2 * y2 * W8 * cosB / rb ^ 3 / W7 * (1 - 2 * nu - W1 / W7 * W2 - a * y3b / rb2) * 2 *
        y3b - y2 * W8 * cosB / rb / W7 ^ 2 * (1 - 2 * nu - W1 / W7 * W2 - a * y3b / rb2) * W3 + y2 *
        W8 * cosB / rb / W7 * (-(cosB * y3b / rb + 1) / W7 * W2 + W1 / W7 ^ 2 * W2 * W3 + 1 / 2 *
        W1 / W7 * a / rb ^ 3 * 2 * y3b - a / rb2 + a * y3b / rb2 ^ 2 * 2 * y3b)) /π/(1 - nu)) +
        b2 * (1 / 4 * ((-2 + 2 * nu) * N1 * cotB * ((y3b / rb + 1) / W6 - cosB * W3 / W7) + (2 - 2 * nu) *
        y1 / W6 ^ 2 * W5 * (y3b / rb + 1) + 1 / 2 * (2 - 2 * nu) * y1 / W6 * a / rb ^ 3 * 2 * y3b + (2 -
        2 * nu)* sinB / W7 * W2 - (2 - 2 * nu) * z1b / W7 ^ 2 * W2 * W3 - 1 / 2 * (2 - 2 * nu) * z1b /
        W7 * a / rb ^ 3 * 2 * y3b + 1 / rb * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb2) - 1 / 2 *
        W8 / rb ^ 3 * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb2) * 2 * y3b + W8 / rb * (2 * nu *
        y1 / W6 ^ 2 * (y3b / rb + 1) + a * y1 / rb2 ^ 2 * 2 * y3b) -1 / W7 * (cosB * sinB + W1 *
        cotB / rb * ((2 - 2 * nu) * cosB - W1 / W7) + a / rb * (sinB - y3b * z1b / rb2 - z1b *
        W1 / rb / W7)) +W8 / W7 ^ 2 * (cosB * sinB + W1 * cotB / rb * ((2 - 2 * nu) * cosB - W1 /
        W7) +a / rb * (sinB - y3b * z1b / rb2 - z1b * W1 / rb / W7)) * W3 - W8 / W7 * ((cosB *
        y3b / rb + 1)* cotB / rb * ((2 - 2 * nu) * cosB - W1 / W7) - 1 / 2 * W1 * cotB / rb ^ 3 *
        ((2 - 2 * nu) * cosB - W1 / W7) * 2 * y3b + W1 * cotB / rb * (-(cosB * y3b / rb + 1) / W7 +
        W1 / W7 ^ 2 * W3) -1 / 2 * a / rb ^ 3 * (sinB - y3b * z1b / rb2 - z1b * W1 / rb / W7) *
        2 * y3b + a / rb * (-z1b / rb2 - y3b * sinB / rb2 + y3b * z1b / rb2 ^ 2 * 2 * y3b -
        sinB * W1 / rb / W7 - z1b * (cosB * y3b / rb + 1) / rb / W7 + 1 / 2 * z1b * W1 / rb ^ 3 /
        W7 * 2 * y3b + z1b * W1 / rb / W7 ^ 2 * W3))) /π/(1 - nu)) +
        b3 * (1 / 4 * ((2 - 2 * nu) * rFib_ry3 - (2 - 2 * nu) * y2 * sinB / W7 ^ 2 * W2 * W3 - 1 / 2 *
        (2 - 2 * nu) * y2 * sinB / W7 * a / rb ^ 3 * 2 * y3b + y2 * sinB / rb / W7 * (1 + W1 / W7 *
        W2 + a * y3b / rb2) -1 / 2 * y2 * W8 * sinB / rb ^ 3 / W7 * (1 + W1 / W7 * W2 + a * y3b /
        rb2)* 2 * y3b - y2 * W8 * sinB / rb / W7 ^ 2 * (1 + W1 / W7 * W2 + a * y3b / rb2) * W3 +
        y2 * W8 * sinB / rb / W7 * ((cosB * y3b / rb + 1) / W7 * W2 - W1 / W7 ^ 2 * W2 * W3 -
        1 / 2 * W1 / W7 * a / rb ^ 3 * 2 * y3b + a / rb2 - a * y3b / rb2 ^ 2 * 2 * y3b)) /π/(1 - nu)))

    v12 = (b1 / 2 * (1 / 4 * ((-2 + 2 * nu) * N1 * rFib_ry2 * cotB ^ 2 + N1 / W6 * ((1 - W5) * cotB - y1 /
        W6 * W4) -N1 * y2 ^ 2 / W6 ^ 2 * ((1 - W5) * cotB - y1 / W6 * W4) / rb + N1 * y2 / W6 *
        (a / rb ^ 3 * y2 * cotB + y1 / W6 ^ 2 * W4 / rb * y2 + y2 / W6 * a / rb ^ 3 * y1) + N1 *
        cosB * cotB / W7 * W2 - N1 * y2 ^ 2 * cosB * cotB / W7 ^ 2 * W2 / rb - N1 * y2 ^ 2 * cosB *
        cotB / W7 * a / rb ^ 3 + a * W8 * cotB / rb ^ 3 - 3 * a * y2 ^ 2 * W8 * cotB / rb ^ 5 + W8 /
        rb / W6 * (-N1 * cotB + y1 / W6 * W5 + a * y1 / rb2) - y2 ^ 2 * W8 / rb ^ 3 / W6 * (-N1 *
        cotB + y1 / W6 * W5 + a * y1 / rb2) -y2 ^ 2 * W8 / rb2 / W6 ^ 2 * (-N1 * cotB + y1 /
        W6 * W5 + a * y1 / rb2) +y2 * W8 / rb / W6 * (-y1 / W6 ^ 2 * W5 / rb * y2 - y2 / W6 *
        a / rb ^ 3 * y1 - 2 * a * y1 / rb2 ^ 2 * y2) +W8 / rb / W7 * (cosB / W7 * (W1 * (N1 *
        cosB - a / rb)* cotB + (2 - 2 * nu) * (rb * sinB - y1) * cosB) -a * y3b * cosB * cotB /
        rb2) -y2 ^ 2 * W8 / rb ^ 3 / W7 * (cosB / W7 * (W1 * (N1 * cosB - a / rb) * cotB + (2 -
        2 * nu) * (rb * sinB - y1) * cosB) -a * y3b * cosB * cotB / rb2) -y2 ^ 2 * W8 / rb2 /
        W7 ^ 2 * (cosB / W7 * (W1 * (N1 * cosB - a / rb) * cotB + (2 - 2 * nu) * (rb * sinB - y1) *
        cosB) -a * y3b * cosB * cotB / rb2) +y2 * W8 / rb / W7 * (-cosB / W7 ^ 2 * (W1 *
        (N1 * cosB - a / rb) * cotB + (2 - 2 * nu) * (rb * sinB - y1) * cosB) /rb *y2+cosB /
        W7 * (1 / rb * cosB * y2 * (N1 * cosB - a / rb) * cotB + W1 * a / rb ^ 3 * y2 * cotB + (2 - 2 *
        nu) /rb*sinB *y2*cosB)+2*a *y3b*cosB*cotB /rb2 ^ 2 * y2)) /π/(1 - nu)) +
        b2 / 2 * (1 / 4 * (N1 * (((2 - 2 * nu) * cotB ^ 2 + nu) / rb * y2 / W6 - ((2 - 2 * nu) * cotB ^ 2 + 1) *
        cosB / rb * y2 / W7) -N1 / W6 ^ 2 * (-N1 * y1 * cotB + nu * y3b - a + a * y1 * cotB / rb +
        y1 ^ 2 / W6 * W4) /rb *y2+N1 /W6 * (-a * y1 * cotB / rb ^ 3 * y2 - y1 ^ 2 / W6 ^
        2 * W4 / rb * y2 - y1 ^ 2 / W6 * a / rb ^ 3 * y2) +N1 * cotB / W7 ^ 2 * (z1b * cosB - a *
        (rb * sinB - y1) / rb / cosB) /rb *y2-N1*cotB /W7 * (-a / rb2 * sinB * y2 /
        cosB + a * (rb * sinB - y1) / rb ^ 3 / cosB * y2) +3 * a * y2 * W8 * cotB / rb ^ 5. * y1 -
        W8 / W6 ^ 2 * (2 * nu + 1 / rb * (N1 * y1 * cotB + a) - y1 ^ 2 / rb / W6 * W5 - a * y1 ^ 2 /
        rb ^ 3) /rb *y2+W8 /W6 * (-1 / rb ^ 3 * (N1 * y1 * cotB + a) * y2 + y1 ^ 2 / rb ^
        3 / W6 * W5 * y2 + y1 ^ 2 / rb2 / W6 ^ 2 * W5 * y2 + y1 ^ 2 / rb2 ^ 2 / W6 * a * y2 + 3 *
        a * y1 ^ 2 / rb ^ 5. * y2) -W8 * cotB / W7 ^ 2 * (-cosB * sinB + a * y1 * y3b / rb ^ 3 /
        cosB + (rb * sinB - y1) / rb * ((2 - 2 * nu) * cosB - W1 / W7 * W9)) /rb *y2+W8*cotB /
        W7 * (-3 * a * y1 * y3b / rb ^ 5. / cosB * y2 + 1 / rb2 * sinB * y2 * ((2 - 2 * nu) * cosB -
        W1 / W7 * W9) -(rb * sinB - y1) / rb ^ 3 * ((2 - 2 * nu) * cosB - W1 / W7 * W9) * y2 + (rb *
        sinB - y1) /rb *(-1 /rb * cosB * y2 / W7 * W9 + W1 / W7 ^ 2 * W9 / rb * y2 + W1 / W7 *
        a / rb ^ 3 / cosB * y2))) /π/(1 - nu)) +
        b3 / 2 * (1 / 4 * (N1 * (1 / W6 * (1 + a / rb) - y2 ^ 2 / W6 ^ 2 * (1 + a / rb) / rb - y2 ^ 2 /
        W6 * a / rb ^ 3 - cosB / W7 * W2 + y2 ^ 2 * cosB / W7 ^ 2 * W2 / rb + y2 ^ 2 * cosB / W7 *
        a / rb ^ 3) -W8 / rb * (a / rb2 + 1 / W6) + y2 ^ 2 * W8 / rb ^ 3 * (a / rb2 + 1 / W6) -
        y2 * W8 / rb * (-2 * a / rb2 ^ 2 * y2 - 1 / W6 ^ 2 / rb * y2) + W8 * cosB / rb / W7 *
        (W1 / W7 * W2 + a * y3b / rb2) - y2 ^ 2 * W8 * cosB / rb ^ 3 / W7 * (W1 / W7 * W2 + a *
        y3b / rb2) -y2 ^ 2 * W8 * cosB / rb2 / W7 ^ 2 * (W1 / W7 * W2 + a * y3b / rb2) + y2 *
        W8 * cosB / rb / W7 * (1 / rb * cosB * y2 / W7 * W2 - W1 / W7 ^ 2 * W2 / rb * y2 - W1 /
        W7 * a / rb ^ 3 * y2 - 2 * a * y3b / rb2 ^ 2 * y2)) /π/(1 - nu)) +
        b1 / 2 * (1 / 4 * (N1 * (((2 - 2 * nu) * cotB ^ 2 - nu) / rb * y1 / W6 - ((2 - 2 * nu) * cotB ^ 2 + 1 -
        2 * nu)* cosB * (y1 / rb - sinB) / W7) +N1 / W6 ^ 2 * (y1 * cotB * (1 - W5) + nu * y3b -
        a + y2 ^ 2 / W6 * W4) /rb *y1-N1 /W6 * ((1 - W5) * cotB + a * y1 ^ 2 * cotB / rb ^ 3 -
        y2 ^ 2 / W6 ^ 2 * W4 / rb * y1 - y2 ^ 2 / W6 * a / rb ^ 3 * y1) -N1 * cosB * cotB / W7 *
        W2 + N1 * z1b * cotB / W7 ^ 2 * W2 * (y1 / rb - sinB) + N1 * z1b * cotB / W7 * a / rb ^
        3 * y1 - a * W8 * cotB / rb ^ 3 + 3 * a * y1 ^ 2 * W8 * cotB / rb ^ 5 - W8 / W6 ^ 2 * (-2 *
        nu + 1 / rb * (N1 * y1 * cotB - a) + y2 ^ 2 / rb / W6 * W5 + a * y2 ^ 2 / rb ^ 3) /rb *
        y1 + W8 / W6 * (-1 / rb ^ 3 * (N1 * y1 * cotB - a) * y1 + 1 / rb * N1 * cotB - y2 ^ 2 /
        rb ^ 3 / W6 * W5 * y1 - y2 ^ 2 / rb2 / W6 ^ 2 * W5 * y1 - y2 ^ 2 / rb2 ^ 2 / W6 * a * y1 -
        3 * a * y2 ^ 2 / rb ^ 5. * y1) -W8 / W7 ^ 2 * (cosB ^ 2 - 1 / rb * (N1 * z1b * cotB + a *
        cosB) +a * y3b * z1b * cotB / rb ^ 3 - 1 / rb / W7 * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB /
        rb * W1)) * (y1 / rb - sinB) + W8 / W7 * (1 / rb ^ 3 * (N1 * z1b * cotB + a * cosB) *
        y1 - 1 / rb * N1 * cosB * cotB + a * y3b * cosB * cotB / rb ^ 3 - 3 * a * y3b * z1b * cotB /
        rb ^ 5. * y1 + 1 / rb ^ 3 / W7 * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB / rb * W1) * y1 + 1 /
        rb / W7 ^ 2 * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB / rb * W1) * (y1 / rb - sinB) - 1 / rb /
        W7 * (-a * cosB * cotB / rb * W1 + a * z1b * cotB / rb ^ 3 * W1 * y1 - a * z1b * cotB /
        rb2 * cosB * y1))) /π/(1 - nu)) +
        b2 / 2 * (1 / 4 * ((2 - 2 * nu) * N1 * rFib_ry1 * cotB ^ 2 - N1 * y2 / W6 ^ 2 * ((W5 - 1) * cotB +
        y1 / W6 * W4) /rb *y1+N1 *y2 /W6 * (-a / rb ^ 3 * y1 * cotB + 1 / W6 * W4 - y1 ^
        2 / W6 ^ 2 * W4 / rb - y1 ^ 2 / W6 * a / rb ^ 3) +N1 * y2 * cotB / W7 ^ 2 * W9 * (y1 /
        rb - sinB) +N1 * y2 * cotB / W7 * a / rb ^ 3 / cosB * y1 + 3 * a * y2 * W8 * cotB / rb ^
        5. * y1 - y2 * W8 / rb ^ 3 / W6 * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 /
        W6)) * y1 - y2 * W8 / rb2 / W6 ^ 2 * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 /
        rb + 1 / W6)) * y1 + y2 * W8 / rb / W6 * (-2 * nu / W6 + 2 * nu * y1 ^ 2 / W6 ^ 2 / rb - a /
        rb * (1 / rb + 1 / W6) + a * y1 ^ 2 / rb ^ 3 * (1 / rb + 1 / W6) - a * y1 / rb * (-1 /
        rb ^ 3 * y1 - 1 / W6 ^ 2 / rb * y1)) -y2 * W8 * cotB / rb ^ 3 / W7 * ((-2 + 2 * nu) *
        cosB + W1 / W7 * W9 + a * y3b / rb2 / cosB) * y1 - y2 * W8 * cotB / rb / W7 ^ 2 * ((-2 +
        2 * nu)* cosB + W1 / W7 * W9 + a * y3b / rb2 / cosB) * (y1 / rb - sinB) + y2 * W8 *
        cotB / rb / W7 * (1 / rb * cosB * y1 / W7 * W9 - W1 / W7 ^ 2 * W9 * (y1 / rb - sinB) -
        W1 / W7 * a / rb ^ 3 / cosB * y1 - 2 * a * y3b / rb2 ^ 2 / cosB * y1)) /π/(1 - nu)) +
        b3 / 2 * (1 / 4 * (N1 * (-sinB * (y1 / rb - sinB) / W7 - 1 / W6 * (1 + a / rb) + y1 ^ 2 / W6 ^
        2 * (1 + a / rb) / rb + y1 ^ 2 / W6 * a / rb ^ 3 + cosB / W7 * W2 - z1b / W7 ^ 2 * W2 *
        (y1 / rb - sinB) - z1b / W7 * a / rb ^ 3 * y1) +W8 / rb * (a / rb2 + 1 / W6) - y1 ^ 2 *
        W8 / rb ^ 3 * (a / rb2 + 1 / W6) + y1 * W8 / rb * (-2 * a / rb2 ^ 2 * y1 - 1 / W6 ^ 2 /
        rb * y1) +W8 / W7 ^ 2 * (sinB * (cosB - a / rb) + z1b / rb * (1 + a * y3b / rb2) - 1 /
        rb / W7 * (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1)) * (y1 / rb - sinB) - W8 / W7 *
        (sinB * a / rb ^ 3 * y1 + cosB / rb * (1 + a * y3b / rb2) - z1b / rb ^ 3 * (1 + a * y3b /
        rb2) * y1 - 2 * z1b / rb ^ 5 * a * y3b * y1 + 1 / rb ^ 3 / W7 * (y2 ^ 2 * cosB * sinB - a *
        z1b / rb * W1) * y1 + 1 / rb / W7 ^ 2 * (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1) *
        (y1 / rb - sinB) - 1 / rb / W7 * (-a * cosB / rb * W1 + a * z1b / rb ^ 3 * W1 * y1 - a *
        z1b / rb2 * cosB * y1))) /π/(1 - nu)))

    v13 = (b1 / 2 * (1 / 4 * ((-2 + 2 * nu) * N1 * rFib_ry3 * cotB ^ 2 - N1 * y2 / W6 ^ 2 * ((1 - W5) *
        cotB - y1 / W6 * W4) * (y3b / rb + 1) + N1 * y2 / W6 * (1 / 2 * a / rb ^ 3 * 2 * y3b * cotB +
        y1 / W6 ^ 2 * W4 * (y3b / rb + 1) + 1 / 2 * y1 / W6 * a / rb ^ 3 * 2 * y3b) -N1 * y2 * cosB *
        cotB / W7 ^ 2 * W2 * W3 - 1 / 2 * N1 * y2 * cosB * cotB / W7 * a / rb ^ 3 * 2 * y3b + a /
        rb ^ 3 * y2 * cotB - 3 / 2 * a * y2 * W8 * cotB / rb ^ 5 * 2 * y3b + y2 / rb / W6 * (-N1 *
        cotB + y1 / W6 * W5 + a * y1 / rb2) -1 / 2 * y2 * W8 / rb ^ 3 / W6 * (-N1 * cotB + y1 /
        W6 * W5 + a * y1 / rb2)* 2 * y3b - y2 * W8 / rb / W6 ^ 2 * (-N1 * cotB + y1 / W6 * W5 +
        a * y1 / rb2) * (y3b / rb + 1) + y2 * W8 / rb / W6 * (-y1 / W6 ^ 2 * W5 * (y3b / rb +
        1) -1 / 2 * y1 / W6 * a / rb ^ 3 * 2 * y3b - a * y1 / rb2 ^ 2 * 2 * y3b) +y2 / rb / W7 *
        (cosB / W7 * (W1 * (N1 * cosB - a / rb) * cotB + (2 - 2 * nu) * (rb * sinB - y1) * cosB) -
        a * y3b * cosB * cotB / rb2) -1 / 2 * y2 * W8 / rb ^ 3 / W7 * (cosB / W7 * (W1 * (N1 *
        cosB - a / rb)* cotB + (2 - 2 * nu) * (rb * sinB - y1) * cosB) -a * y3b * cosB * cotB /
        rb2)* 2 * y3b - y2 * W8 / rb / W7 ^ 2 * (cosB / W7 * (W1 * (N1 * cosB - a / rb) * cotB +
        (2 - 2 * nu) * (rb * sinB - y1) * cosB) -a * y3b * cosB * cotB / rb2) * W3 + y2 * W8 / rb /
        W7 * (-cosB / W7 ^ 2 * (W1 * (N1 * cosB - a / rb) * cotB + (2 - 2 * nu) * (rb * sinB - y1) *
        cosB) * W3 + cosB / W7 * ((cosB * y3b / rb + 1) * (N1 * cosB - a / rb) * cotB + 1 / 2 * W1 *
        a / rb ^ 3 * 2 * y3b * cotB + 1 / 2 * (2 - 2 * nu) / rb * sinB * 2 * y3b * cosB) -a * cosB *
        cotB / rb2 + a * y3b * cosB * cotB / rb2 ^ 2 * 2 * y3b)) /π/(1 - nu)) +
        b2 / 2 * (1 / 4 * (N1 * (((2 - 2 * nu) * cotB ^ 2 + nu) * (y3b / rb + 1) / W6 - ((2 - 2 * nu) * cotB ^
        2 + 1)* cosB * W3 / W7) -N1 / W6 ^ 2 * (-N1 * y1 * cotB + nu * y3b - a + a * y1 * cotB /
        rb + y1 ^ 2 / W6 * W4) * (y3b / rb + 1) + N1 / W6 * (nu - 1 / 2 * a * y1 * cotB / rb ^ 3 * 2 *
        y3b - y1 ^ 2 / W6 ^ 2 * W4 * (y3b / rb + 1) - 1 / 2 * y1 ^ 2 / W6 * a / rb ^ 3 * 2 * y3b) +
        N1 * cotB / W7 ^ 2 * (z1b * cosB - a * (rb * sinB - y1) / rb / cosB) * W3 - N1 * cotB /
        W7 * (cosB * sinB - 1 / 2 * a / rb2 * sinB * 2 * y3b / cosB + 1 / 2 * a * (rb * sinB - y1) /
        rb ^ 3 / cosB * 2 * y3b) -a / rb ^ 3 * y1 * cotB + 3 / 2 * a * y1 * W8 * cotB / rb ^ 5 * 2 *
        y3b + 1 / W6 * (2 * nu + 1 / rb * (N1 * y1 * cotB + a) - y1 ^ 2 / rb / W6 * W5 - a * y1 ^ 2 /
        rb ^ 3) -W8 / W6 ^ 2 * (2 * nu + 1 / rb * (N1 * y1 * cotB + a) - y1 ^ 2 / rb / W6 * W5 - a *
        y1 ^ 2 / rb ^ 3) * (y3b / rb + 1) + W8 / W6 * (-1 / 2 / rb ^ 3 * (N1 * y1 * cotB + a) * 2 *
        y3b + 1 / 2 * y1 ^ 2 / rb ^ 3 / W6 * W5 * 2 * y3b + y1 ^ 2 / rb / W6 ^ 2 * W5 * (y3b / rb +
        1) +1 / 2 * y1 ^ 2 / rb2 ^ 2 / W6 * a * 2 * y3b + 3 / 2 * a * y1 ^ 2 / rb ^ 5 * 2 * y3b) +
        cotB / W7 * (-cosB * sinB + a * y1 * y3b / rb ^ 3 / cosB + (rb * sinB - y1) / rb * ((2 -
        2 * nu)* cosB - W1 / W7 * W9)) -W8 * cotB / W7 ^ 2 * (-cosB * sinB + a * y1 * y3b / rb ^
        3 / cosB + (rb * sinB - y1) / rb * ((2 - 2 * nu) * cosB - W1 / W7 * W9)) * W3 + W8 * cotB /
        W7 * (a / rb ^ 3 / cosB * y1 - 3 / 2 * a * y1 * y3b / rb ^ 5. / cosB * 2 * y3b + 1 / 2 /
        rb2 * sinB * 2 * y3b * ((2 - 2 * nu) * cosB - W1 / W7 * W9) - 1 / 2 * (rb * sinB - y1) / rb ^
        3 * ((2 - 2 * nu) * cosB - W1 / W7 * W9) * 2 * y3b + (rb * sinB - y1) / rb * (-(cosB * y3b /
        rb + 1) /W7 *W9+W1 /W7 ^ 2 * W9 * W3 + 1 / 2 * W1 / W7 * a / rb ^ 3 / cosB * 2 *
        y3b))) /π/(1 - nu)) +
        b3 / 2 * (1 / 4 * (N1 * (-y2 / W6 ^ 2 * (1 + a / rb) * (y3b / rb + 1) - 1 / 2 * y2 / W6 * a /
        rb ^ 3 * 2 * y3b + y2 * cosB / W7 ^ 2 * W2 * W3 + 1 / 2 * y2 * cosB / W7 * a / rb ^ 3 * 2 *
        y3b) -y2 / rb * (a / rb2 + 1 / W6) + 1 / 2 * y2 * W8 / rb ^ 3 * (a / rb2 + 1 / W6) * 2 *
        y3b - y2 * W8 / rb * (-a / rb2 ^ 2 * 2 * y3b - 1 / W6 ^ 2 * (y3b / rb + 1)) + y2 * cosB /
        rb / W7 * (W1 / W7 * W2 + a * y3b / rb2) - 1 / 2 * y2 * W8 * cosB / rb ^ 3 / W7 * (W1 /
        W7 * W2 + a * y3b / rb2)* 2 * y3b - y2 * W8 * cosB / rb / W7 ^ 2 * (W1 / W7 * W2 + a *
        y3b / rb2) * W3 + y2 * W8 * cosB / rb / W7 * ((cosB * y3b / rb + 1) / W7 * W2 - W1 /
        W7 ^ 2 * W2 * W3 - 1 / 2 * W1 / W7 * a / rb ^ 3 * 2 * y3b + a / rb2 - a * y3b / rb2 ^ 2 * 2 *
        y3b)) /π/(1 - nu)) +
        b1 / 2 * (1 / 4 * ((2 - 2 * nu) * (N1 * rFib_ry1 * cotB - y1 / W6 ^ 2 * W5 / rb * y2 - y2 / W6 *
        a / rb ^ 3 * y1 + y2 * cosB / W7 ^ 2 * W2 * (y1 / rb - sinB) + y2 * cosB / W7 * a / rb ^
        3 * y1) -y2 * W8 / rb ^ 3 * (2 * nu / W6 + a / rb2) * y1 + y2 * W8 / rb * (-2 * nu / W6 ^
        2 / rb * y1 - 2 * a / rb2 ^ 2 * y1) -y2 * W8 * cosB / rb ^ 3 / W7 * (1 - 2 * nu - W1 / W7 *
        W2 - a * y3b / rb2) * y1 - y2 * W8 * cosB / rb / W7 ^ 2 * (1 - 2 * nu - W1 / W7 * W2 - a *
        y3b / rb2) * (y1 / rb - sinB) + y2 * W8 * cosB / rb / W7 * (-1 / rb * cosB * y1 / W7 *
        W2 + W1 / W7 ^ 2 * W2 * (y1 / rb - sinB) + W1 / W7 * a / rb ^ 3 * y1 + 2 * a * y3b / rb2 ^
        2 * y1)) /π/(1 - nu)) +
        b2 / 2 * (1 / 4 * ((-2 + 2 * nu) * N1 * cotB * (1 / rb * y1 / W6 - cosB * (y1 / rb - sinB) / W7) -
        (2 - 2 * nu) / W6 * W5 + (2 - 2 * nu) * y1 ^ 2 / W6 ^ 2 * W5 / rb + (2 - 2 * nu) * y1 ^ 2 / W6 *
        a / rb ^ 3 + (2 - 2 * nu) * cosB / W7 * W2 - (2 - 2 * nu) * z1b / W7 ^ 2 * W2 * (y1 / rb -
        sinB) -(2 - 2 * nu) * z1b / W7 * a / rb ^ 3 * y1 - W8 / rb ^ 3 * (N1 * cotB - 2 * nu * y1 /
        W6 - a * y1 / rb2) * y1 + W8 / rb * (-2 * nu / W6 + 2 * nu * y1 ^ 2 / W6 ^ 2 / rb - a / rb2 +
        2 * a * y1 ^ 2 / rb2 ^ 2) +W8 / W7 ^ 2 * (cosB * sinB + W1 * cotB / rb * ((2 - 2 * nu) *
        cosB - W1 / W7) +a / rb * (sinB - y3b * z1b / rb2 - z1b * W1 / rb / W7)) * (y1 / rb -
        sinB) -W8 / W7 * (1 / rb2 * cosB * y1 * cotB * ((2 - 2 * nu) * cosB - W1 / W7) - W1 *
        cotB / rb ^ 3 * ((2 - 2 * nu) * cosB - W1 / W7) * y1 + W1 * cotB / rb * (-1 / rb * cosB *
        y1 / W7 + W1 / W7 ^ 2 * (y1 / rb - sinB)) -a / rb ^ 3 * (sinB - y3b * z1b / rb2 -
        z1b * W1 / rb / W7) * y1 + a / rb * (-y3b * cosB / rb2 + 2 * y3b * z1b / rb2 ^ 2 * y1 -
        cosB * W1 / rb / W7 - z1b / rb2 * cosB * y1 / W7 + z1b * W1 / rb ^ 3 / W7 * y1 + z1b *
        W1 / rb / W7 ^ 2 * (y1 / rb - sinB)))) /π/(1 - nu)) +
        b3 / 2 * (1 / 4 * ((2 - 2 * nu) * rFib_ry1 - (2 - 2 * nu) * y2 * sinB / W7 ^ 2 * W2 * (y1 / rb -
        sinB) -(2 - 2 * nu) * y2 * sinB / W7 * a / rb ^ 3 * y1 - y2 * W8 * sinB / rb ^ 3 / W7 * (1 +
        W1 / W7 * W2 + a * y3b / rb2) * y1 - y2 * W8 * sinB / rb / W7 ^ 2 * (1 + W1 / W7 * W2 +
        a * y3b / rb2) * (y1 / rb - sinB) + y2 * W8 * sinB / rb / W7 * (1 / rb * cosB * y1 /
        W7 * W2 - W1 / W7 ^ 2 * W2 * (y1 / rb - sinB) - W1 / W7 * a / rb ^ 3 * y1 - 2 * a * y3b /
        rb2 ^ 2 * y1)) /π/(1 - nu)))

    v23 = (b1 / 2 * (1 / 4 * (N1 * (((2 - 2 * nu) * cotB ^ 2 - nu) * (y3b / rb + 1) / W6 - ((2 - 2 * nu) *
        cotB ^ 2 + 1 - 2 * nu)* cosB * W3 / W7) +N1 / W6 ^ 2 * (y1 * cotB * (1 - W5) + nu * y3b - a +
        y2 ^ 2 / W6 * W4) * (y3b / rb + 1) - N1 / W6 * (1 / 2 * a * y1 * cotB / rb ^ 3 * 2 * y3b +
        nu - y2 ^ 2 / W6 ^ 2 * W4 * (y3b / rb + 1) - 1 / 2 * y2 ^ 2 / W6 * a / rb ^ 3 * 2 * y3b) -N1 *
        sinB * cotB / W7 * W2 + N1 * z1b * cotB / W7 ^ 2 * W2 * W3 + 1 / 2 * N1 * z1b * cotB / W7 *
        a / rb ^ 3 * 2 * y3b - a / rb ^ 3 * y1 * cotB + 3 / 2 * a * y1 * W8 * cotB / rb ^ 5 * 2 * y3b +
        1 / W6 * (-2 * nu + 1 / rb * (N1 * y1 * cotB - a) + y2 ^ 2 / rb / W6 * W5 + a * y2 ^ 2 /
        rb ^ 3) -W8 / W6 ^ 2 * (-2 * nu + 1 / rb * (N1 * y1 * cotB - a) + y2 ^ 2 / rb / W6 * W5 +
        a * y2 ^ 2 / rb ^ 3) * (y3b / rb + 1) + W8 / W6 * (-1 / 2 / rb ^ 3 * (N1 * y1 * cotB - a) *
        2 * y3b - 1 / 2 * y2 ^ 2 / rb ^ 3 / W6 * W5 * 2 * y3b - y2 ^ 2 / rb / W6 ^ 2 * W5 * (y3b /
        rb + 1) -1 / 2 * y2 ^ 2 / rb2 ^ 2 / W6 * a * 2 * y3b - 3 / 2 * a * y2 ^ 2 / rb ^ 5 * 2 * y3b) +
        1 / W7 * (cosB ^ 2 - 1 / rb * (N1 * z1b * cotB + a * cosB) + a * y3b * z1b * cotB / rb ^
        3 - 1 / rb / W7 * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB / rb * W1)) -W8 / W7 ^ 2 * (cosB ^ 2 -
        1 / rb * (N1 * z1b * cotB + a * cosB) + a * y3b * z1b * cotB / rb ^ 3 - 1 / rb / W7 *
        (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB / rb * W1)) * W3 + W8 / W7 * (1 / 2 / rb ^ 3 * (N1 *
        z1b * cotB + a * cosB)* 2 * y3b - 1 / rb * N1 * sinB * cotB + a * z1b * cotB / rb ^ 3 + a *
        y3b * sinB * cotB / rb ^ 3 - 3 / 2 * a * y3b * z1b * cotB / rb ^ 5 * 2 * y3b + 1 / 2 / rb ^
        3 / W7 * (y2 ^ 2 * cosB ^ 2 - a * z1b * cotB / rb * W1) * 2 * y3b + 1 / rb / W7 ^ 2 * (y2 ^
        2 * cosB ^ 2 - a * z1b * cotB / rb * W1) * W3 - 1 / rb / W7 * (-a * sinB * cotB / rb * W1 +
        1 / 2 * a * z1b * cotB / rb ^ 3 * W1 * 2 * y3b - a * z1b * cotB / rb * (cosB * y3b / rb +
        1)))) /π/(1 - nu)) +
        b2 / 2 * (1 / 4 * ((2 - 2 * nu) * N1 * rFib_ry3 * cotB ^ 2 - N1 * y2 / W6 ^ 2 * ((W5 - 1) * cotB +
        y1 / W6 * W4) * (y3b / rb + 1) + N1 * y2 / W6 * (-1 / 2 * a / rb ^ 3 * 2 * y3b * cotB - y1 /
        W6 ^ 2 * W4 * (y3b / rb + 1) - 1 / 2 * y1 / W6 * a / rb ^ 3 * 2 * y3b) +N1 * y2 * cotB /
        W7 ^ 2 * W9 * W3 + 1 / 2 * N1 * y2 * cotB / W7 * a / rb ^ 3 / cosB * 2 * y3b - a / rb ^ 3 *
        y2 * cotB + 3 / 2 * a * y2 * W8 * cotB / rb ^ 5 * 2 * y3b + y2 / rb / W6 * (N1 * cotB - 2 *
        nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 / W6)) -1 / 2 * y2 * W8 / rb ^ 3 / W6 * (N1 *
        cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 / W6))* 2 * y3b - y2 * W8 / rb / W6 ^
        2 * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb * (1 / rb + 1 / W6)) * (y3b / rb + 1) + y2 *
        W8 / rb / W6 * (2 * nu * y1 / W6 ^ 2 * (y3b / rb + 1) + 1 / 2 * a * y1 / rb ^ 3 * (1 / rb +
        1 / W6)* 2 * y3b - a * y1 / rb * (-1 / 2 / rb ^ 3 * 2 * y3b - 1 / W6 ^ 2 * (y3b / rb +
        1))) +y2 * cotB / rb / W7 * ((-2 + 2 * nu) * cosB + W1 / W7 * W9 + a * y3b / rb2 / cosB) -
        1 / 2 * y2 * W8 * cotB / rb ^ 3 / W7 * ((-2 + 2 * nu) * cosB + W1 / W7 * W9 + a * y3b /
        rb2 / cosB)* 2 * y3b - y2 * W8 * cotB / rb / W7 ^ 2 * ((-2 + 2 * nu) * cosB + W1 / W7 *
        W9 + a * y3b / rb2 / cosB) * W3 + y2 * W8 * cotB / rb / W7 * ((cosB * y3b / rb + 1) /
        W7 * W9 - W1 / W7 ^ 2 * W9 * W3 - 1 / 2 * W1 / W7 * a / rb ^ 3 / cosB * 2 * y3b + a / rb2 /
        cosB - a * y3b / rb2 ^ 2 / cosB * 2 * y3b)) /π/(1 - nu)) +
        b3 / 2 * (1 / 4 * (N1 * (-sinB * W3 / W7 + y1 / W6 ^ 2 * (1 + a / rb) * (y3b / rb + 1) +
        1 / 2 * y1 / W6 * a / rb ^ 3 * 2 * y3b + sinB / W7 * W2 - z1b / W7 ^ 2 * W2 * W3 - 1 / 2 *
        z1b / W7 * a / rb ^ 3 * 2 * y3b) +y1 / rb * (a / rb2 + 1 / W6) - 1 / 2 * y1 * W8 / rb ^
        3 * (a / rb2 + 1 / W6) * 2 * y3b + y1 * W8 / rb * (-a / rb2 ^ 2 * 2 * y3b - 1 / W6 ^ 2 *
        (y3b / rb + 1)) -1 / W7 * (sinB * (cosB - a / rb) + z1b / rb * (1 + a * y3b / rb2) - 1 /
        rb / W7 * (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1)) +W8 / W7 ^ 2 * (sinB * (cosB -
        a / rb) +z1b / rb * (1 + a * y3b / rb2) - 1 / rb / W7 * (y2 ^ 2 * cosB * sinB - a * z1b /
        rb * W1)) * W3 - W8 / W7 * (1 / 2 * sinB * a / rb ^ 3 * 2 * y3b + sinB / rb * (1 + a * y3b /
        rb2) -1 / 2 * z1b / rb ^ 3 * (1 + a * y3b / rb2) * 2 * y3b + z1b / rb * (a / rb2 - a *
        y3b / rb2 ^ 2 * 2 * y3b) +1 / 2 / rb ^ 3 / W7 * (y2 ^ 2 * cosB * sinB - a * z1b / rb *
        W1)* 2 * y3b + 1 / rb / W7 ^ 2 * (y2 ^ 2 * cosB * sinB - a * z1b / rb * W1) * W3 - 1 /
        rb / W7 * (-a * sinB / rb * W1 + 1 / 2 * a * z1b / rb ^ 3 * W1 * 2 * y3b - a * z1b / rb *
        (cosB * y3b / rb + 1)))) /π/(1 - nu)) +
        b1 / 2 * (1 / 4 * ((2 - 2 * nu) * (N1 * rFib_ry2 * cotB + 1 / W6 * W5 - y2 ^ 2 / W6 ^ 2 * W5 /
        rb - y2 ^ 2 / W6 * a / rb ^ 3 - cosB / W7 * W2 + y2 ^ 2 * cosB / W7 ^ 2 * W2 / rb + y2 ^ 2 *
        cosB / W7 * a / rb ^ 3) +W8 / rb * (2 * nu / W6 + a / rb2) - y2 ^ 2 * W8 / rb ^ 3 * (2 *
        nu / W6 + a / rb2) +y2 * W8 / rb * (-2 * nu / W6 ^ 2 / rb * y2 - 2 * a / rb2 ^ 2 * y2) +
        W8 * cosB / rb / W7 * (1 - 2 * nu - W1 / W7 * W2 - a * y3b / rb2) - y2 ^ 2 * W8 * cosB /
        rb ^ 3 / W7 * (1 - 2 * nu - W1 / W7 * W2 - a * y3b / rb2) - y2 ^ 2 * W8 * cosB / rb2 / W7 ^
        2 * (1 - 2 * nu - W1 / W7 * W2 - a * y3b / rb2) + y2 * W8 * cosB / rb / W7 * (-1 / rb *
        cosB * y2 / W7 * W2 + W1 / W7 ^ 2 * W2 / rb * y2 + W1 / W7 * a / rb ^ 3 * y2 + 2 * a *
        y3b / rb2 ^ 2 * y2)) /π/(1 - nu)) +
        b2 / 2 * (1 / 4 * ((-2 + 2 * nu) * N1 * cotB * (1 / rb * y2 / W6 - cosB / rb * y2 / W7) + (2 -
        2 * nu) * y1 / W6 ^ 2 * W5 / rb * y2 + (2 - 2 * nu) * y1 / W6 * a / rb ^ 3 * y2 - (2 - 2 *
        nu) * z1b / W7 ^ 2 * W2 / rb * y2 - (2 - 2 * nu) * z1b / W7 * a / rb ^ 3 * y2 - W8 / rb ^
        3 * (N1 * cotB - 2 * nu * y1 / W6 - a * y1 / rb2) * y2 + W8 / rb * (2 * nu * y1 / W6 ^ 2 /
        rb * y2 + 2 * a * y1 / rb2 ^ 2 * y2) +W8 / W7 ^ 2 * (cosB * sinB + W1 * cotB / rb * ((2 -
        2 * nu)* cosB - W1 / W7) +a / rb * (sinB - y3b * z1b / rb2 - z1b * W1 / rb / W7)) /
        rb * y2 - W8 / W7 * (1 / rb2 * cosB * y2 * cotB * ((2 - 2 * nu) * cosB - W1 / W7) - W1 *
        cotB / rb ^ 3 * ((2 - 2 * nu) * cosB - W1 / W7) * y2 + W1 * cotB / rb * (-cosB / rb *
        y2 / W7 + W1 / W7 ^ 2 / rb * y2) -a / rb ^ 3 * (sinB - y3b * z1b / rb2 - z1b * W1 /
        rb / W7) * y2 + a / rb * (2 * y3b * z1b / rb2 ^ 2 * y2 - z1b / rb2 * cosB * y2 / W7 +
        z1b * W1 / rb ^ 3 / W7 * y2 + z1b * W1 / rb2 / W7 ^ 2 * y2))) /π/(1 - nu)) +
        b3 / 2 * (1 / 4 * ((2 - 2 * nu) * rFib_ry2 + (2 - 2 * nu) * sinB / W7 * W2 - (2 - 2 * nu) * y2 ^ 2 *
        sinB / W7 ^ 2 * W2 / rb - (2 - 2 * nu) * y2 ^ 2 * sinB / W7 * a / rb ^ 3 + W8 * sinB / rb /
        W7 * (1 + W1 / W7 * W2 + a * y3b / rb2) - y2 ^ 2 * W8 * sinB / rb ^ 3 / W7 * (1 + W1 /
        W7 * W2 + a * y3b / rb2) -y2 ^ 2 * W8 * sinB / rb2 / W7 ^ 2 * (1 + W1 / W7 * W2 + a *
        y3b / rb2) +y2 * W8 * sinB / rb / W7 * (1 / rb * cosB * y2 / W7 * W2 - W1 / W7 ^ 2 *
        W2 / rb * y2 - W1 / W7 * a / rb ^ 3 * y2 - 2 * a * y3b / rb2 ^ 2 * y2)) /π/(1 - nu)))

    return v11, v22, v33, v12, v13, v23
end

function strain_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T) where {T, V}
    A = transform_matrix(P1, P2, P3)
    _strain_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, λ, μ, A')
end

"""
    stress_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T) where {T, V}

Compute stress risen from triangular dislocation in elastic *fullspace*.
    Please see [original version (in supporting information)](https://academic.oup.com/gji/article/201/2/1119/572006#86405752)
    for details, especially the **coordinate system** used here.

## Arguments
The same as [`stress_tri3_hs`](@ref)
"""
function stress_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T) where {T, V}
    exx, eyy, ezz, exy, exz, eyz = strain_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, λ, μ)
    return _strain2stress(exx, eyy, ezz, exy, exz, eyz, λ, μ)
end

@inline function _strain2stress(exx::T, eyy::T, ezz::T, exy::T, exz::T, eyz::T, λ::T, μ::T) where T
    ekk = exx + eyy + ezz
    σxx = λ * ekk + 2μ * exx
    σyy = λ * ekk + 2μ * eyy
    σzz = λ * ekk + 2μ * ezz
    σxy = 2μ * exy
    σxz = 2μ * exz
    σyz = 2μ * eyz
    return σxx, σyy, σzz, σxy, σxz, σyz
end

@inline function _strain_tri3_fs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T, At::M) where {T, V, M}
    nu = λ / (μ + λ) / 2
    bx, by, bz = Ts, Ss, Ds
    p1, p2, p3 = [zeros(T, 3) for _ in 1: 3]
    x, y, z = coord_trans(X - P2[1], Y - P2[2], Z - P2[3], At)
    p1[1], p1[2], p1[3] = coord_trans(P1[1] - P2[1], P1[2] - P2[2], P1[3] - P2[3], At)
    p3[1], p3[2], p3[3] = coord_trans(P3[1] - P2[1], P3[2] - P2[2], P3[3] - P2[3], At)

    e12 = normalize(p2 - p1)
    e13 = normalize(p3 - p1)
    e23 = normalize(p3 - p2)

    A = acos(dot(e12, e13))
    B = acos(dot(-e12, e23))
    C = acos(dot(e23, e13))

    Trimode = trimodefinder(y, z, x, p1, p2, p3)

    if Trimode == 1
        Exx1Tp, Eyy1Tp, Ezz1Tp, Exy1Tp, Exz1Tp, Eyz1Tp = TDSetupS(x, y, z, A, bx, by, bz, nu, p1, -e13)
        Exx2Tp, Eyy2Tp, Ezz2Tp, Exy2Tp, Exz2Tp, Eyz2Tp = TDSetupS(x, y, z, B, bx, by, bz, nu, p2, e12)
        Exx3Tp, Eyy3Tp, Ezz3Tp, Exy3Tp, Exz3Tp, Eyz3Tp = TDSetupS(x, y, z, C, bx, by, bz, nu, p3, e23)
        exx = Exx1Tp + Exx2Tp + Exx3Tp
        eyy = Eyy1Tp + Eyy2Tp + Eyy3Tp
        ezz = Ezz1Tp + Ezz2Tp + Ezz3Tp
        exy = Exy1Tp + Exy2Tp + Exy3Tp
        exz = Exz1Tp + Exz2Tp + Exz3Tp
        eyz = Eyz1Tp + Eyz2Tp + Eyz3Tp
    end

    if Trimode == -1
        Exx1Tn, Eyy1Tn, Ezz1Tn, Exy1Tn, Exz1Tn, Eyz1Tn = TDSetupS(x, y, z, A, bx, by, bz, nu, p1, e13)
        Exx2Tn, Eyy2Tn, Ezz2Tn, Exy2Tn, Exz2Tn, Eyz2Tn = TDSetupS(x, y, z, B, bx, by, bz, nu, p2, -e12)
        Exx3Tn, Eyy3Tn, Ezz3Tn, Exy3Tn, Exz3Tn, Eyz3Tn = TDSetupS(x, y, z, C, bx, by, bz, nu, p3, -e23)
        exx = Exx1Tn + Exx2Tn + Exx3Tn
        eyy = Eyy1Tn + Eyy2Tn + Eyy3Tn
        ezz = Ezz1Tn + Ezz2Tn + Ezz3Tn
        exy = Exy1Tn + Exy2Tn + Exy3Tn
        exz = Exz1Tn + Exz2Tn + Exz3Tn
        eyz = Eyz1Tn + Eyz2Tn + Eyz3Tn
    end

    if Trimode == 0
        exx = NaN
        eyy = NaN
        ezz = NaN
        exy = NaN
        exz = NaN
        eyz = NaN
    end

    Exx, Eyy, Ezz, Exy, Exz, Eyz = TensTrans(exx, eyy, ezz, exy, exz, eyz, At')
    return Exx, Eyy, Ezz, Exy, Exz, Eyz
end

@inline function TDSetupS(x::T, y::T, z::T, alpha::T, bx::T, by::T, bz::T, nu::T, TriVertex::V, SideVec::V) where {T, V}
    y1 = SideVec[3] * (y - TriVertex[2]) - SideVec[2] * (z - TriVertex[3])
    z1 = SideVec[2] * (y - TriVertex[2]) + SideVec[3] * (z - TriVertex[3])
    by1 = SideVec[3] * by - SideVec[2] * bz
    bz1 = SideVec[2] * by + SideVec[3] * bz

    exx, eyy, ezz, exy, exz, eyz = AngDisStrain(x, y1, z1, -π + alpha, bx, by1, bz1, nu)
    B = zeros(T, 3, 3)
    B[1,1] = one(T)
    B[2,2] = SideVec[3]
    B[2,3] = SideVec[2]
    B[3,2] = -SideVec[2]
    B[3,3] = SideVec[3]
    exx, eyy, ezz, exy, exz, eyz = TensTrans(exx, eyy, ezz, exy, exz, eyz, B)
    return exx, eyy, ezz, exy, exz, eyz
end

function strain_tri3_hs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T) where {T, V}
    A = transform_matrix(P1, P2, P3)
    strms11, strms22, strms33, strms12, strms13, strms23 = _strain_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, λ, μ, A')
    strfsc11, strfsc22, strfsc33, strfsc12, strfsc13, strfsc23 = _TDstrain_HarFunc(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, λ, μ, A)
    P1[3] *= -one(T)
    P2[3] *= -one(T)
    P3[3] *= -one(T)
    allatsurface = P1[3] == zero(T) && P2[3] == zero(T) && P3[3] == zero(T)
    if !allatsurface
        A *= -one(T)
        A[3,1] *= -one(T)
        A[3,2] *= -one(T)
        A[3,3] *= -one(T)
    end
    stris11, stris22, stris33, stris12, stris13, stris23 = _strain_tri3_fs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, λ, μ, A')
    if allatsurface
        stris13 *= -one(T)
        stris23 *= -one(T)
    end
    P1[3] *= -one(T)
    P2[3] *= -one(T)
    P3[3] *= -one(T)
    e11 = strms11 + strfsc11 + stris11
    e22 = strms22 + strfsc22 + stris22
    e33 = strms33 + strfsc33 + stris33
    e12 = strms12 + strfsc12 + stris12
    e13 = strms13 + strfsc13 + stris13
    e23 = strms23 + strfsc23 + stris23
    return e11, e22, e33, e12, e13, e23
end

"""
    stress_tri3_hs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T) where {T, V}

Compute stress risen from triangular dislocation in elastic *halfspace*.
    Please see [original version (in supporting information)](https://academic.oup.com/gji/article/201/2/1119/572006#86405752)
    for details, especially the **coordinate system** used here.

## Arguments
- `X`, `Y`, `Z`: observational coordinates
- `P1`, `P2`, `P3`: three triangular vertices coordinates respectively
- `Ss`, `Ds`, `Ts`: triangular dislocation vector, Strike-slip, Dip-slip, Tensile-slip respectively
- `λ`: Lamé's first parameter
- `μ`: shear modulus

## Output
By order: ``σ_{xx}``, ``σ_{yy}``, ``σ_{zz}``,
    ``σ_{xy}``, ``σ_{xz}``, ``σ_{yz}``.
"""
function stress_tri3_hs(X::T, Y::T, Z::T, P1::V, P2::V, P3::V, Ss::T, Ds::T, Ts::T, λ::T, μ::T) where {T, V}
    exx, eyy, ezz, exy, exz, eyz = strain_tri3_hs(X, Y, Z, P1, P2, P3, Ss, Ds, Ts, λ, μ)
    return _strain2stress(exx, eyy, ezz, exy, exz, eyz, λ, μ)
end


text:
function EventsLogger()
    return EventsLogger(())
end

function clear_logged_events(obj::EventsLogger)
    return jcall(obj, "clearLoggedEvents", void, ())
end

function get_logged_events(obj::EventsLogger)
    return jcall(obj, "getLoggedEvents", List, ())
end

function monitor_detector(obj::EventsLogger, arg0::EventDetector)
    return jcall(obj, "monitorDetector", EventDetector, (EventDetector,), arg0)
end

ext: jl


text:
## Precedence climbing, for re-arranging infix to prefix

head(x) = nil
peek(x) = isempty(x) ? nil : last(x)

function infixtoprefix(x, c)
    y = reverse(x)
    isempty(y) && return nil
    climb(pop!(y), y, nil, c)
end

function climb(lhs, rest, curop, ctx)

    while true

        nextop = peek(rest)

        isop(nextop, ctx) || break
        curop ≡ nil || curop == nextop || cmp(curop, nextop, ctx) < 0 || break

        op  = pop!(rest)
        rhs = pop!(rest)

        while true

            nextop = peek(rest)

            isop(nextop, ctx) || break
            cmp(op, nextop, ctx) < 0 || break

            rhs = climb(rhs, rest, nextop, ctx)

        end

        if lhs isa Ex && head(lhs) == op
            push!(lhs, rhs)
        else
            lhs = node(op, lhs, rhs, ctx)
        end

    end

    return lhs

end


text:
plot(exp, 0, 1)
plot!(sqrt, 0, 1)
legend()


text:
#=
Convert a vector of pairs of {choice => availaility} to a matrix
if availability is nothing, return nothing
=#
function availability_to_matrix(availability::Union{Nothing, <:AbstractVector{<:Pair{<:Any, <:AbstractVector{Bool}}}}, alt_numbers::Dict{<:Any, Int64})::Union{BitMatrix, Nothing}
    if isnothing(availability)
        return nothing
    else
        avail_mat = BitMatrix(undef, length(availability[1].second), length(alt_numbers))
        for (name, avvec) in availability
            index = alt_numbers[name]
            @assert !isnothing(index)
            avail_mat[:, index] = avvec
        end

        return avail_mat
    end
end

# Parallel functions on Tables.namedtupleiterator. From https://github.com/JuliaData/Tables.jl/pull/187
function SplittablesBase.halve(rows::Tables.NamedTupleIterator{schema}) where schema
    left, right = SplittablesBase.halve(rows.x)
    return (
        Tables.NamedTupleIterator{schema,typeof(left)}(left),
        Tables.NamedTupleIterator{schema,typeof(right)}(right),
    )
end

function consistent_rowcount(cols)
    len = length(cols[1])
    if !all(c -> length(c) == len, cols)
        throw(ArgumentError("`halve` on columns return inconsistent number or rows"))
    end
    return len
end

function SplittablesBase.halve(x::Tables.RowIterator)
    if isempty(Tables.columns(x))
        len = cld(length(x), 2)
        return (Tables.RowIterator(columns(x), len), Tables.RowIterator(columns(x), length(x) - len))
    end
    cs = map(SplittablesBase.halve, Tables.columns(x))
    lefts = map(first, cs)
    rights = map(last, cs)
    return (
        Tables.RowIterator(lefts, consistent_rowcount(lefts)),
        Tables.RowIterator(rights, consistent_rowcount(rights)),
    )
end

# END COPIED CODE

#=
Find unique values of a column
function find_unique_values(table::JuliaDB.AbstractIndexedTable, column)
    reducer(x::NamedTuple, y::NamedTuple) = Set{typeof(x[column])}([x[column], y[column]])
    # treat sets as immutable. not sure if this matters but could imagine data races.
    reducer(x::Set, y::NamedTuple) = y[column] in x ? x : Set{typeof(y[column])}([x..., y[column]])
    reducer(x::Set, y::Set) = union(x, y)

    return reduce(reducer, table)
end
=#

#=
Get column type for a column in an IndexedTable
=#
coltype(table, column) = fieldtype(rowtype(table), column)
#rowtype(table::JuliaDB.AbstractIndexedTable) = eltype(table)
rowtype(table) = typeof(first(Tables.rows(table)))
# special code path for looping over data frames to ensure type stability uses numedtupleiterator
rowtype(table::DataFrame) = typeof(first(Tables.namedtupleiterator(table)))

# Work around https://github.com/JuliaData/Tables.jl/issues/264, using rowtables for dtables improves perf
Tables.columnnames(::Vector{NamedTuple{N, T}}) where {N, T} = N

#=
Evaluate perfect prediction problems, returns true if perfect prediction found

function check_perfect_prediction(table::JuliaDB.AbstractIndexedTable, choice, predictors=collect(filter(c -> c != choice, colnames(table))))
    # first, find all unique values of choice
    unique_choice = find_unique_values(table, choice)

    problems_found = false

    # now, reduce to find min/max of each predictor for each choice
    # TODO this could be more efficient, in theory could call reduce for all cols at once, but I get a
    # JuliaDB error about namedtuples with duplicate names - I think due to using the same selector for
    # a bunch of OnlineStats
    for column in predictors
        ct = coltype(table, column)
        for (i, chc) in enumerate(unique_choice)
            if ct == Bool
                if value(reduce(Sum(Int64), table; select=column)) == 0
                    @error "Column $column has no true values"
                else
                    sums = groupreduce(Series(Sum(Int64), Counter(Int64)), table, choice => x -> x == chc; select=column) |>
                    collect |>
                    DataFrame

                    sum_true = value(sums[sums[:,choice], :Series][1])[1]
                    count_true = value(sums[sums[:,choice], :Series][1])[2]
                    sum_false = value(sums[.!sums[:,choice], :Series][1])[1]
                    count_false = value(sums[.!sums[:,choice], :Series][1])[2]

                    if sum_true == 0
                        @error "Column $column == true implies $chc not selected"
                        problems_found = true
                    elseif sum_true == count_true
                        @error "Column $column == true implies $chc always selected"
                        problems_found = true
                    end
                    if sum_false == 0
                        @error "Column $column == false implies $chc not selected"
                        problems_found = true

                    elseif sum_false == count_false
                        @error "Column $column == false implies $chc always selected"
                        problems_found = true
                    end
                end
            else
                minmax = groupreduce(Extrema(ct), table, choice => x -> x == chc; select=column) |>
                    collect |>
                    DataFrame

                minmax_true = minmax[minmax[:,choice], :Extrema][1]
                minmax_false = minmax[.!minmax[:,choice], :Extrema][1]

                if nobs(minmax_true) == 0
                    @error "Choice $chc not selected"
                    problems_found = true
                elseif nobs(minmax_false) == 0
                    @error "Choice $chc always selected"
                    problems_found = true
                elseif minimum(minmax_true) > maximum(minmax_false)
                    @error "Column $column > $(maximum(minmax_false)) perfectly predicts $choice == $chc"
                    problems_found = true
                elseif maximum(minmax_true) < minimum(minmax_false)
                    @error "Column $column < $(minimum(minmax_false)) perfectly predicts $choice == $chc"
                    problems_found = true
                end
            end
        end

        vrnc = var(reduce(Variance(ct == Bool ? Float64 : ct), table, select=column))
        if vrnc < 1e-5
            @warn "Column $column has near-zero variance"
            problems_found = true
        end
    end

    return problems_found
end
=#

#length(table::DTable) = fetch(reduce(+, map(r -> (unity=1,), table))).unity


text:
@testset "> Smith Normal Form" begin
  @testset ">> choose pivot" begin

    @testset ">>> With zero columns" begin
      smith = [0 4 4; 0 6 12; 0 -4 -16]
      left = eye(eltype(smith), 3)
      jt = Crystals.SNF.choose_pivot!(left, smith, 1, 1)
      @test jt == 2
      @test smith == [0 4 4; 0 6 12; 0 -4 -16]
      @test left == eye(eltype(smith), 3)
    end

    @testset ">>> Switch rows" begin
      smith = [0 0 4; 0 6 12; 0 -4 -16]
      original = deepcopy(smith)
      left = eye(eltype(smith), 3)
      jt = Crystals.SNF.choose_pivot!(left, smith, 1, 1)
      @test jt == 2
      @test smith == [0 6 12; 0 0 4; 0 -4 -16]
      @test left == [0 1 0; 1 0 0; 0 0 1]
      @test left * original == smith
    end
  end

  @testset ">> Improve col pivot" begin
    @testset ">>> All multiples" begin
      smith = [1 0 4; 0 2 12; 0 -4 -16]
      left = eye(eltype(smith), 3)
      Crystals.SNF.improve_col_pivot!(left, smith, 1, 1)
      @test left == eye(eltype(smith), 3)
      @test smith == [1 0 4; 0 2 12; 0 -4 -16]

      Crystals.SNF.improve_col_pivot!(left, smith, 2, 2)
      @test left == eye(eltype(smith), 3)
      @test smith == [1 0 4; 0 2 12; 0 -4 -16]
    end

    @testset ">>> Create multiple" begin
      smith = [1 0 4; 0 3 12; 0 -4 -16]
      original = deepcopy(smith)
      left = eye(eltype(smith), 3)
      Crystals.SNF.improve_col_pivot!(left, smith, 2, 2)
      @test all(smith[:, 2] .% smith[2, 2] .== 0)
      @test left * original == smith
      @test abs(det(left)) ≥ 0
    end
  end

  @testset ">> Improve row pivot" begin
    @testset ">>> All multiples" begin
      smith = transpose([1 0 4; 0 2 12; 0 -4 -16])
      right = eye(eltype(smith), 3)
      Crystals.SNF.improve_row_pivot!(smith, right, 1, 1)
      @test right == eye(eltype(smith), 3)
      @test smith == transpose([1 0 4; 0 2 12; 0 -4 -16])

      Crystals.SNF.improve_row_pivot!(smith, right, 2, 2)
      @test right == eye(eltype(smith), 3)
      @test smith == transpose([1 0 4; 0 2 12; 0 -4 -16]   )
    end

    @testset ">>> Create multiple" begin
      smith = transpose([1 0 4; 0 3 12; 0 -4 -16])
      original = deepcopy(smith)
      right = eye(eltype(smith), 3)
      Crystals.SNF.improve_row_pivot!(smith, right, 2, 2)
      @test all(smith[2, :] .% smith[2, 2] .== 0)
      @test original * right == smith
      @test abs(det(right)) ≥ 0
    end
  end

  @testset ">> Diagnonalize one row-column" begin
    smith = [1 2 4; 3 3 12; 5 -4 -16]
    original = deepcopy(smith)
    right = eye(eltype(smith), 3)
    left = eye(eltype(smith), 3)
    Crystals.SNF.diagonalize_at_entry!(left, smith, right, 2, 2)
    @test smith[2, 2] ≠ 0
    @test countnz(smith[:, 2]) == 1
    @test countnz(smith[2, :]) == 1
    @test left * original * right == smith
    @test abs(det(left)) ≥ 0
    @test abs(det(right)) ≥ 0
  end

  @testset ">> Diaginalize all entries" begin
    original = [1 2 4; 3 3 12; 5 -4 -16]
    smith, left, right = Crystals.SNF.diagonalize_all_entries(original)
    @test isdiag(smith)
    @test left * original * right == smith
    @test abs(det(left)) ≥ 0
    @test abs(det(right)) ≥ 0
  end

  @testset ">> Move zero columns left" begin
    matrices = Any[
      [0 1 0; 0 2 2; 0 3 3],
      [0 1 0; 0 0 1; 0 0 0]
    ]
    for smith in matrices
      original = deepcopy(smith)
      right = eye(eltype(smith), 3)
      Crystals.SNF.move_zero_entries!(smith, right)
      @test smith[:, 1] == original[:, 2]
      @test smith[:, 2] == original[:, 3]
      @test smith[:, 3] == original[:, 1]
      @test original * right == smith
      @test abs(det(right)) ≥ 0
    end
  end

  @testset ">> The whole shebang" begin
    matrices = vcat(
      Any[
        BigInt[1 2; -4 5],
        BigInt[0 -4 -5; 0 0 -4; 0 3 -1],
        BigInt[2 4 4; -6 6 12; 10 -4 -16],
        BigInt[-5 -5 2 -2; -4 -4 0 -1; -1 -1 4 -1; 0 -4 4 -4]
      ],
      Any[rand(-BigInt(5):5, tuple(repeat(rand(2:5, (1, )), inner=[2])...)) for u in 1:10]
    )
    @testset ">>> try $matrix" for matrix in matrices
      smith, left, right = Crystals.SNF.smith_normal_form(matrix)
      @test isdiag(smith)
      @test left * matrix * right == smith
      k = findfirst(x -> x == 0, diag(smith))
      if k ≠ 0
      @test countnz(diag(smith)[k:end]) == 0
      else
        k = size(smith, 2) + 1
      end
      for i in 2:(k - 1)
        @test smith[i, i] % smith[i - 1, i - 1] == 0
      end
      @test all(diag(smith) .≥ 0)
    end
  end

  @testset ">> Units" begin
      matrix = BigInt[1 2; -4 5]
      usmith, uleft, uright = Crystals.SNF.smith_normal_form(matrix * u"nm")
      smith, left, right = Crystals.SNF.smith_normal_form(matrix)
      @test  uleft == left
      @test  uright == right
      @test eltype(usmith) == eltype(matrix * u"nm")
      usmith = reshape([ustrip.(u) for u in usmith], size(usmith))
      @test usmith == smith
  end
end


text:
using Pkg

using Revise

include("00_ingestion.jl")
include("../src/tweet_graphing.jl")

using Chain
using DataFrames
using DataFramesMeta
using Dates

function get_user_data()
    dfs = read_files("Data")
    out = []
    for df in dfs
        push!(out, _load_user_data(df))
    end
    return out
end

function _load_user_data(df)
    user_info = groupby(select(df, :user_username, :retweet_count, :like_count, :quote_count, :reply_count), :user_username)
    user_info = combine(user_info, nrow => :ntweets, :retweet_count => sum => :nretweets, :like_count => sum => :nlikes, :quote_count => sum => :nquotes, :reply_count => sum => :nreplies)
    return user_info
end


text:
md"""
An `S2Cap` is a circular patch on a sphere defined by a center and a angular radius.

The radius is expressed as a chord angle to make computations easier.
"""
struct S2Cap <: S2Region
    center::Vector
    radius::S1ChordAngle
    S2Cap(center::Vector, edge::Vector) = verifyPoint(center) && new(center, S1ChordAngle(center, edge))
    S2Cap(center::Vector, radius::S1ChordAngle) = verifyPoint(center) && new(center, radius)
    S2Cap(center::Vector, angle::S1Angle) = verifyPoint(center) && new(center, S1ChordAngle(angle))
end

equals(a::S2Cap, b::S2Cap) = a.center == b.center && a.radius == b.radius


text:
# quadratic constraints

MOI.supports_constraint(o::Optimizer, ::Type{SQF}, ::Type{<:BOUNDS}) = true

function MOI.add_constraint(o::Optimizer, func::SQF, set::S) where {S <: BOUNDS}
    if func.constant != 0.0
        error("SCIP does not support quadratic constraints with a constant offset.")
    end

    allow_modification(o)

    # affine terms
    linrefs = [VarRef(t.variable.value) for t in func.affine_terms]
    lincoefs = [t.coefficient for t in func.affine_terms]

    # quadratic terms
    quadrefs1 = [VarRef(t.variable_1.value) for t in func.quadratic_terms]
    quadrefs2 = [VarRef(t.variable_2.value) for t in func.quadratic_terms]
    # Divide coefficients by 2 iff they come from the diagonal:
    # Take coef * x * y as-is, but turn coef * x^2 into coef/2 * x^2.
    factor = 1.0 .- 0.5 * (quadrefs1 .== quadrefs2)
    quadcoefs = factor .* [t.coefficient for t in func.quadratic_terms]

    # range
    lhs, rhs = bounds(set)
    lhs = lhs === nothing ? -SCIPinfinity(o) : lhs
    rhs = rhs === nothing ?  SCIPinfinity(o) : rhs

    cr = add_quadratic_constraint(o.inner, linrefs, lincoefs,
                                  quadrefs1, quadrefs2, quadcoefs, lhs, rhs)
    ci = CI{SQF, S}(cr.val)
    register!(o, ci)
    register!(o, cons(o, ci), cr)
    return ci
end

function MOI.set(o::SCIP.Optimizer, ::MOI.ConstraintSet, ci::CI{SQF,S}, set::S) where {S <: BOUNDS}
    allow_modification(o)

    lhs, rhs = bounds(set)
    lhs = lhs === nothing ? -SCIPinfinity(o) : lhs
    rhs = rhs === nothing ?  SCIPinfinity(o) : rhs

    @SCIP_CALL SCIPchgLhsQuadratic(o, cons(o, ci), lhs)
    @SCIP_CALL SCIPchgRhsQuadratic(o, cons(o, ci), rhs)

    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{SQF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    c = cons(o, ci)
    expr_ref = SCIPgetExprNonlinear(c)
    # This call is required to get quaddata computed in the expression
    isq = Ref{UInt32}(100)
    @SCIP_CALL LibSCIP.SCIPcheckExprQuadratic(o, expr_ref, isq)
    if isq[] != 1
        error("Constraint index $ci pointing to a non-quadratic expression $expr_ref")
    end

    constant_ref = Ref{Cdouble}(-1.0)
    n_linear_terms_ref = Ref{Cint}(-1)
    linear_exprs = Ref{Ptr{Ptr{LibSCIP.SCIP_EXPR}}}()
    lincoefs = Ref{Ptr{Cdouble}}()
    n_quad_terms_ref = Ref{Cint}(-1)
    n_bilinear_terms_ref = Ref{Cint}(-1)

    LibSCIP.SCIPexprGetQuadraticData(
        expr_ref,
        constant_ref,
        n_linear_terms_ref,
        linear_exprs,
        lincoefs,
        n_quad_terms_ref,
        n_bilinear_terms_ref,
        C_NULL,
        C_NULL,
    )

    lin_expr_vec = unsafe_wrap(Vector{Ptr{Cvoid}}, linear_exprs[], n_linear_terms_ref[])
    lin_coeff_vec = unsafe_wrap(Vector{Cdouble}, lincoefs[], n_linear_terms_ref[])

    func = SCIP.SQF([], [], constant_ref[])
    for idx in 1:n_linear_terms_ref[]
        var_ptr = LibSCIP.SCIPgetVarExprVar(lin_expr_vec[idx])
        func += lin_coeff_vec[idx] * MOI.VariableIndex(o.reference[var_ptr].val)
    end
    for term_idx in 1:n_quad_terms_ref[]
        var_expr = Ref{Ptr{Cvoid}}()
        lin_coef_ref =  Ref{Cdouble}()
        sqr_coef_ref =  Ref{Cdouble}()
        LibSCIP.SCIPexprGetQuadraticQuadTerm(
            expr_ref,
            term_idx - 1, # 0-indexed terms
            var_expr,
            lin_coef_ref,
            sqr_coef_ref,
            C_NULL,
            C_NULL,
            C_NULL,
        )
        @assert lin_coef_ref[] == 0.0
        if sqr_coef_ref[] != 0.0
            var_ptr = LibSCIP.SCIPgetVarExprVar(var_expr[])
            var_idx = MOI.VariableIndex(o.reference[var_ptr].val)
            MOI.Utilities.operate!(+, Float64, func, sqr_coef_ref[] * var_idx * var_idx)
        end
    end

    for term_idx in 1:n_bilinear_terms_ref[]
        var_expr1 = Ref{Ptr{Cvoid}}()
        var_expr2 = Ref{Ptr{Cvoid}}()
        coef_ref =  Ref{Cdouble}()
        LibSCIP.SCIPexprGetQuadraticBilinTerm(
            expr_ref,
            term_idx - 1,
            var_expr1,
            var_expr2,
            coef_ref,
            C_NULL,
            C_NULL,
        )
        if coef_ref[] != 0.0
            var_ptr1 = LibSCIP.SCIPgetVarExprVar(var_expr1[])
            var_idx1 = MOI.VariableIndex(o.reference[var_ptr1].val)
            var_ptr2 = LibSCIP.SCIPgetVarExprVar(var_expr2[])
            var_idx2 = MOI.VariableIndex(o.reference[var_ptr2].val)
            MOI.Utilities.operate!(+, Float64, func, coef_ref[] * var_idx1 * var_idx2)
        end
    end
    return func
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{SQF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    c = cons(o, ci)
    lhs = SCIPgetLhsNonlinear(c)
    rhs = SCIPgetRhsNonlinear(c)
    return from_bounds(S, lhs, rhs)
end

function MOI.get(o::Optimizer, ::MOI.ConstraintPrimal, ci::CI{SQF, S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    c = cons(o, ci)
    expr_ref = SCIPgetExprNonlinear(c)
    isq = Ref{UInt32}(100)
    @SCIP_CALL LibSCIP.SCIPcheckExprQuadratic(o, expr_ref, isq)
    if isq[] != 1
        error("Constraint index $ci pointing to a non-quadratic expression $expr_ref")
    end
    sol = SCIPgetBestSol(o)
    @SCIP_CALL SCIPevalExpr(o, expr_ref, sol, Clonglong(0))
    return SCIPexprGetEvalValue(expr_ref)
end


text:
################################################################################
#
# The following file implements polynomials (in a hard coded, hacky way). The
# plaintext polynomials included here are hard-coded in order for Flux.jl GPU
# training to work properly without having to offload the operations to the
# CPU.
#
################################################################################

module Polynomials

    using Memoize
    using SEAL

    export PolynomialActivation

    """
        PolynomialActivation(degree::Integer, coeffs::Vector)

    Evaluates the polynomial of the provided degree and coefficients. The
    vector of coefficients should start from the 0:th degree up to the provided
    degree.
    """
    function PolynomialActivation(degree::Integer, coeffs::Vector)
        # TODO: do this in a better way.
        if degree == 2
            return Poly2(tuple(coeffs...))
        elseif degree == 3
            return Poly3(tuple(coeffs...))
        elseif degree == 4
            return Poly4(tuple(coeffs...))
        elseif degree == 5
            return Poly5(tuple(coeffs...))
        elseif degree == 6
            return Poly6(tuple(coeffs...))
        elseif degree == 7
            return Poly7(tuple(coeffs...))
        elseif degree == 8
            return Poly8(tuple(coeffs...))
        elseif degree == 9
            return Poly9(tuple(coeffs...))
        elseif degree == 10
            return Poly10(tuple(coeffs...))
        else
            throw(ArgumentError("""
            Sorry, I hard-coded this for GPU training to work.
            Pull requests with a fix are welcome :)
            """))
        end
    end

    struct Poly2
        coeffs::Tuple{Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly2)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2
    end

    function (p::Poly2)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly2 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly3
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly3)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3
    end

    function (p::Poly3)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly3 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly4
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly4)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4
    end

    function (p::Poly4)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly4 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly5
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly5)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4 .+
        p.coeffs[6] .* x .^5
    end

    function (p::Poly5)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly5 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly6
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly6)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4 .+
        p.coeffs[6] .* x .^5 .+
        p.coeffs[7] .* x .^6
    end

    function (p::Poly6)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly6 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly7
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly7)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4 .+
        p.coeffs[6] .* x .^5 .+
        p.coeffs[7] .* x .^6 .+
        p.coeffs[8] .* x .^7
    end

    function (p::Poly7)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly7 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly8
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly8)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4 .+
        p.coeffs[6] .* x .^5 .+
        p.coeffs[7] .* x .^6 .+
        p.coeffs[8] .* x .^7 .+
        p.coeffs[9] .* x .^8
    end

    function (p::Poly8)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly8 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly9
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly9)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4 .+
        p.coeffs[6] .* x .^5 .+
        p.coeffs[7] .* x .^6 .+
        p.coeffs[8] .* x .^7 .+
        p.coeffs[9] .* x .^8 .+
        p.coeffs[10] .* x .^9
    end

    function (p::Poly9)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly9 layer"

        poly_log_eval(x, p.coeffs)
    end

    struct Poly10
        coeffs::Tuple{Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32,
                      Float32}
    end

    function (p::Poly10)(x::AbstractVecOrMat)
        p.coeffs[1] .+
        p.coeffs[2] .* x .+
        p.coeffs[3] .* x .^2 .+
        p.coeffs[4] .* x .^3 .+
        p.coeffs[5] .* x .^4 .+
        p.coeffs[6] .* x .^5 .+
        p.coeffs[7] .* x .^6 .+
        p.coeffs[8] .* x .^7 .+
        p.coeffs[9] .* x .^8 .+
        p.coeffs[10] .* x .^9 .+
        p.coeffs[11] .* x .^10
    end

    function (p::Poly10)(x::AbstractVecOrMat{T}) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function.")
        else
            SEALParams = Main.SEALParams
        end

        @debug "Performing forward-pass of Poly10 layer"

        poly_log_eval(x, p.coeffs)
    end

    """
        poly_log_eval(x::AbstractVecOrMat, coeffs::Vector)

    Evaluates and returns the polynomial of 'x' defined by the passed in
    coefficients 'coeffs'. The evaluation is performed in a logarithmic
    multiplicative depth fashion (i.e. low number of consecutive multiplications).
    """
    function poly_log_eval(x::AbstractVecOrMat{T}, coeffs) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function")
        else
            SEALParams = Main.SEALParams
        end

        # The polynomial should not _only_ be a constant.
        @assert length(coeffs) > 1

        res = similar(x)

        constant, rest... = coeffs

        res_initalized = map(x -> x = false, Array{Bool}(undef, size(x)))

        # Go through the non-constant coefficients in reverse order:
        for (deg, coeff) in Iterators.reverse(enumerate(rest))
            x_i = power_log_eval(x, deg)

            # Calculate c * x_i
            for i in 1:length(x_i)
                p = Plaintext()
                encode!(p, Float64(coeff), scale(x_i[i]), SEALParams.encoder)
                mod_switch_to_inplace!(p, parms_id(x_i[i]), SEALParams.evaluator)

                c = Ciphertext()

                multiply_plain!(c, x_i[i], p, SEALParams.evaluator)
                relinearize_inplace!(c, SEALParams.relinearization_keys, SEALParams.evaluator)
                rescale_to_next_inplace!(c, SEALParams.evaluator)

                if res_initalized[i]
                    # Make sure that x1 and x2 has the same modulus.
                    c_m = chain_index(get_context_data(SEALParams.context, parms_id(c)))
                    res_m = chain_index(get_context_data(SEALParams.context, parms_id(res[i])))

                    if c_m < res_m
                        mod_switch_to_inplace!(res[i], parms_id(c), SEALParams.evaluator)
                    else
                        mod_switch_to_inplace!(c, parms_id(res[i]), SEALParams.evaluator)
                    end

                    # NOTE: useless scaling, but without it there's a bug in
                    # SEAL that returns NaN...
                    scale!(c, scale(c))

                    # NOTE: Hard-coded rescaling to make addition work :/
                    scale!(c, scale(res[i]))

                    add!(res[i], res[i], c, SEALParams.evaluator)
                else
                    res[i] = c
                    res_initalized[i] = true
                end
            end
        end

        # Add the constant coefficient:
        for i in 1:length(x)
            p = Plaintext()
            encode!(p, Float64(constant), scale(res[i]), SEALParams.encoder)
            mod_switch_to_inplace!(p, parms_id(res[i]), SEALParams.evaluator)

            add_plain!(res[i], res[i], p, SEALParams.evaluator)
        end

        return res
    end

    """
        power_log_eval(x::AbstractVecOrMat, degree::Integer)

    Evaluates and returns 'x' raised to the power of 'degree' in a logarithmic
    multiplicative depth fashion (i.e. low number of consecutive multiplications).
    """
    @memoize function power_log_eval(x::AbstractVecOrMat{T}, degree::Integer) where {T<:Ciphertext}
        if !isdefined(Main, :SEALParams)
            throw("SEALParams needs to be defined to run this function")
        else
            SEALParams = Main.SEALParams
        end

        res = similar(x)

        # Base case (poor man's version of a deep clone since it's not available).
        # We use multiply for this (and thereby consume one multiplication) in
        # order to simplify scale handling for other functions.
        if degree == 1
            for i in 1:length(x)
                p = Plaintext()
                encode!(p, Float64(1), scale(x[i]), SEALParams.encoder)

                x_c = Ciphertext()

                multiply_plain!(x_c, x[i], p, SEALParams.evaluator)
                relinearize_inplace!(x_c, SEALParams.relinearization_keys, SEALParams.evaluator)
                rescale_to_next_inplace!(x_c, SEALParams.evaluator)

                res[i] = x_c
            end

            return res
        end

        # Base case:
        if degree == 2
            for i in 1:length(x)
                x_c = Ciphertext()

                square!(x_c, x[i], SEALParams.evaluator)
                relinearize_inplace!(x_c, SEALParams.relinearization_keys, SEALParams.evaluator)
                rescale_to_next_inplace!(x_c, SEALParams.evaluator)

                res[i] = x_c
            end

            return res
        end

        # Recursive case:
        split = Int(2^(ceil(log2(degree) - 1)))

        x1 = power_log_eval(x, split)
        x2 = power_log_eval(x, degree - split)

        for i in 1:length(x)
            x_c = Ciphertext()

            # Make sure that x1 and x2 has the same modulus.
            x1_m = chain_index(get_context_data(SEALParams.context, parms_id(x1[i])))
            x2_m = chain_index(get_context_data(SEALParams.context, parms_id(x2[i])))

            if x1_m < x2_m
                mod_switch_to_inplace!(x2[i], parms_id(x1[i]), SEALParams.evaluator)
            else
                mod_switch_to_inplace!(x1[i], parms_id(x2[i]), SEALParams.evaluator)
            end

            multiply!(x_c, x1[i], x2[i], SEALParams.evaluator)
            relinearize_inplace!(x_c, SEALParams.relinearization_keys, SEALParams.evaluator)
            rescale_to_next_inplace!(x_c, SEALParams.evaluator)

            res[i] = x_c
        end

        return res
    end
end


text:
module RegexTools

include("hex_escape.jl")

end


text:
module RBM

#using PyPlot

export RBM_t
export serialize
export rbm_copy
export rbm_create
export rbm_create_full
export rbm_init_weights_random!
export rbm_init_visible_bias!
export rbm_init_output_bias_random!, rbm_init_hidden_bias_random!
export rbm_rescale_weights!
#= export rbm_calculate_L1, rbm_calculate_L2 =#
export rbm_write, rbm_read, rbm_read_old
export sigm
export rbm_visualise

type RBM_t
  bins::Int64
  n::Int64            # number of output nodes
  m::Int64            # number of hidden nodes
  k::Int64            # number of input nodes
  uditer::Int64       # number of up-down passes when sampling
  alpha::Float64      # learning rate
  momentum::Float64   # contribution of previous gradient
  weightcost::Float64 # weight cost for L2
  numepochs::Int64    # number of training epochs
  batchsize::Int64    # size of training data batches
  W::Matrix{Float64}  # interaction weights between hidden and outputs
  V::Matrix{Float64}  # interaction weights between hidden and inputs
  b::Vector{Float64}  # bias weights for outputs
  c::Vector{Float64}  # bias weights for hiddens
  vW::Matrix{Float64} # for the momentum
  vV::Matrix{Float64} # for the momentum
  vb::Vector{Float64} # for the momentum
  vc::Vector{Float64} # for the momentum
  dropout::Float64
  perturbation::Float64
end

#= serialize(s, rbm::RBM_t) = serialize(f, rbm) =#

function rbm_create_full(n::Int64, m::Int64, k::Int64, uditer::Int64,
  alpha::Float64, momentum::Float64, weightcost::Float64,
  numepochs::Int64, batchsize::Int64, bins::Int64, dropout::Float64,
  perturbation::Float64)

  W = zeros((m, n))
  V = zeros((m, k))
  b = zeros(n)
  c = zeros(m)

  vW = zeros((m, n))
  vV = zeros((m, k))
  vb = zeros(n)
  vc = zeros(m)

  return RBM_t(bins, n, m, k, uditer, alpha,
               momentum, weightcost, numepochs, batchsize,
               W, V, b, c, vW, vV, vb, vc, dropout, perturbation)
end

rbm_create(n::Int64, m::Int64, k::Int64) = rbm_create_full(n, m, k, 15, 0.01, 0.5, 0.001, 10000, 100, 16, 0.0, 0.1)

function rbm_copy(src::RBM_t)
  copy = rbm_create(src.n, src.m, src.k, src.uditer, src.alpha, src.momentum, src.weightcost, src.numepochs, src.batchsize, src.bins, src.dropout)
  copy.W  = deepcopy(src.W)
  copy.V  = deepcopy(src.V)
  copy.b  = deepcopy(src.b)
  copy.c  = deepcopy(src.c)
  copy.vW = deepcopy(src.vW)
  copy.vV = deepcopy(src.vV)
  copy.vb = deepcopy(src.vb)
  copy.vc = deepcopy(src.vc)
  return copy
end

function rbm_init_weights_random!(rbm::RBM_t)
  rbm.W = randn((rbm.m,rbm.n)) .* 0.01
  rbm.V = randn((rbm.m,rbm.k)) .* 0.01
end

function rbm_init_hidden_bias_random!(rbm::RBM_t)
  rbm.b = rand(rbm.n) .* 0.01
end

function rbm_init_output_bias_random!(rbm::RBM_t)
  rbm.c = rand(rbm.m) .* 0.01
end

function rbm_init_visible_bias!(rbm::RBM_t, data::Array{Int64,2})
  # each row of the data must contain be of type {0,1}^k
  n = size(data)[2]
  N = rbm.n
  @assert (size(data)[2] == rbm.n) "rbm_init_visible_bias!: Each row of the data array must be of type {0,1}^n. It has $n and not $N"

  l = size(data)[1]
  for i=1:rbm.n
    s = sum(data[:,i])
    if s == l # only 1 in data
      p =  (l-1) / l
    elseif s == 0 # no 1 in data
      p =  1 / l
    else
      p = s / l
    end
    rbm.b[i] = log2(p/(1.0 - p))
  end
end

function rbm_rescale_weights!(rbm::RBM_t, abs_maximum::Float64)
  @assert (abs_maximum > 0.0) "rbm_rescale_weights!: abs_maximum must be > 0.0"
  mw = maximum(abs(rbm.W))
  mv = maximum(abs(rbm.V))
  if mw > mv
    m = mw
  else
    m = mv
  end
  if m > abs_maximum
    factor = abs_maximum / m
    rbm.W = rbm.W .* factor
    rbm.V = rbm.V .* factor
  end
end

#= function rbm_calculate_L2(rbm::RBM_t, last::Float64) =#
  #= current = sum(rbm.W^2) + sum(rbm.V^2) =#
  #= diff    = rbm.weightcost * (current - last) * rbm.alpha =#
  #= return diff =#
#= end =#

#= function rbm_calculate_L1(rbm::RBM_t, last::Float64) =#
  #= current = sum(abs(rbm.W)) + sum(abs(rbm.V)) =#
  #diff    = rbm.weightcost * (current - last) * rbm.alpha
  #= return diff =#
#= end =#

function rbm_write(filename, rbm)
  fd = open(filename, "w")
  serialize(fd, rbm)
  close(fd)
end

function rbm_read(filename)
  fd = open(filename, "r")
  rbm = deserialize(fd)
  close(fd)
  return rbm
end

function rbm_read_old(filename)
  f        = open(filename, "r")
  s        = split(strip(readline(f)))

  bins     = int(eval(parse(s[1])))
  n        = int(eval(parse(s[2])))
  m        = int(eval(parse(s[3])))
  k        = int(eval(parse(s[4])))
  uditer   = int(eval(parse(s[5])))
  alpha    = float64(eval(parse(s[6])))
  momentum = float64(eval(parse(s[7])))
  epochs   = int(eval(parse(s[8])))
  batch    = int(eval(parse(s[9])))
  memory   = int(eval(parse(s[10])))

  println("n: ", n, " bins: ", bins, " m: ", m, " k: ", k, " uditer: ", uditer)
  rbm = rbm_create(n, m, k, uditer, alpha, momentum, 0.0, epochs, batch, bins, 0.0, 0.0)

  # reading W
  for i=1:m
    l = split(strip(readline(f)))
    for j=1:n
      rbm.W[i,j] = eval(parse(l[j]))
    end
  end
  # reading V
  for i=1:m
    l = split(strip(readline(f)))
    for j=1:k
      rbm.V[i,j]=eval(parse(l[j]))
    end
  end
  # reading b
  bs = split(strip(readline(f)))
  for i=1:length(bs)
    rbm.b[i] = eval(parse(bs[i]))
  end
  cs = split(strip(readline(f)))
  for i=1:length(cs)
    rbm.c[i] = eval(parse(cs[i]))
  end
  close(f)
  return rbm
end

function rbm_visualise(rbm::RBM_t)

end

end # module


text:
"""
	MMC(μ, λ, c)

Tạo mô hình M/M/c
"""
struct MMC{T} <: AbstractMMCK
	μ::Union{T, Real}
	λ::Union{T, Real}
	c::Union{T, Integer}
	r::Union{T, Real}
	ρ::Union{T, Real}
	function MMC(λ, μ, c)
		r = λ / μ
		ρ = r / c
		T = Union{typeof(μ), typeof(λ)}
		new{T}(μ, λ, c, r, ρ)
	end
end

function pn(m::MMC, n::Integer)
	μ = m.μ
	λ = m.λ
	c = m.c
	r = m.r
	ρ = m.ρ
	if n == 0
		d1 = (c * r^c) / factorial(c) / (c - r)
		d2 = sum(r^i / factorial(i) for i = 0:c-1)
		1 / (d1 + d2)
	elseif n < c
		pn(m, 0) * r^n / factorial(n)
	else pn(m, 0) * r^n / factorial(c) / c^(n-c)
	end
end

function Lq(m::MMC)
	r = m.r
	ρ = m.ρ
	c = m.c
	pn(m, 0) * r^c * ρ / factorial(m.c) / (1 - ρ)^2
end

function Wq(m::MMC)
	Lq(m) / m.λ
end

function W(m::MMC)
	Wq(m) + 1 / m.μ
end

function L(m::MMC)
	Lq(m) + m.r
end


text:
# Test routines releated the the spectral radius

using Test
using LinearAlgebra
using QuantumPropagators
using QuantumControlBase.TestUtils

@testset "Ritz values (non-Herm)" begin

    SHOWVALS = false

    N = 1000
    ρ = 10  # spectral radius
    X = random_complex_matrix(N, ρ)
    Ψ = random_state_vector(N)
    # a non_Hermitian X requires significantly more iterations to converge than
    # the default
    ritzvals = QuantumPropagators.ritzvals(
            X, Ψ, 180, 200; prec=1e-5
    )
    evals = eigvals(X)

    E_min = abs(evals[1])
    Δ_min = abs(evals[1] - ritzvals[1])
    SHOWVALS && @show evals[1], ritzvals[1]
    ϵ_min = Δ_min / E_min
    SHOWVALS && @show E_min, Δ_min, ϵ_min
    @test ϵ_min < 0.01

    E_max = abs(evals[end])
    Δ_max = abs(evals[end] - ritzvals[end])
    SHOWVALS && @show evals[end], ritzvals[end]
    ϵ_max = Δ_max / E_max
    SHOWVALS && @show E_max, Δ_max, ϵ_max
    @test ϵ_max < 0.01

    # Note: E_min and E_max should be roughly equal to ρ (Girko-Ginibri
    # circular law)

end

@testset "Ritz values (Hermitian)" begin

    SHOWVALS = false

    N = 1000
    ρ = 10  # spectral radius
    H = random_hermitian_matrix(N, ρ)
    Ψ = random_state_vector(N)
    # for a Hermitian H, lower default (cf. specrange) should be ok
    ritzvals = QuantumPropagators.ritzvals(H, Ψ, 20, 60; prec=1e-3)
    evals = eigvals(H)

    E_min = abs(evals[1])
    Δ_min = abs(evals[1] - ritzvals[1])
    SHOWVALS && @show evals[1], ritzvals[1]
    ϵ_min = Δ_min / E_min
    SHOWVALS && @show E_min, Δ_min, ϵ_min
    @test ϵ_min < 0.01

    E_max = abs(evals[end])
    Δ_max = abs(evals[end] - ritzvals[end])
    SHOWVALS && @show evals[end], ritzvals[end]
    ϵ_max = Δ_max / E_max
    SHOWVALS && @show E_max, Δ_max, ϵ_max
    @test ϵ_max < 0.01

    # Note: E_min and E_max should be roughly equal to ρ (Girko-Ginibri
    # circular law)

end

@testset "Spectral range" begin

    SHOWVALS = false

    N = 1000
    ρ = 10  # spectral radius
    s = 0.1  # sparsity
    H = random_hermitian_sparse_matrix(N, ρ, s)

    evals = eigvals(Array(H))
    SHOWVALS && @show evals[1], evals[end]

    prec = 1e-4
    E_min, E_max = specrange(H, method=:arnoldi, prec=prec)
    SHOWVALS && @show E_min, E_max

    Δ = evals[end] - evals[1]
    SHOWVALS && @show Δ
    δmin = E_min - evals[1]
    δmax = E_max - evals[end]
    SHOWVALS && @show δmin, δmax
    @test evals[1] - 0.05 * Δ <= E_min <= evals[1]
    @test evals[end] <= E_max < evals[end] + 0.05 * Δ

end