text:
#=
#creates the beta_v and beta_T from the vector specified
function gerg_betamatrix_from_vector(v, symmetric_op = (a) -> a)
    N1 = length(v)
    N = Int64((1 + sqrt(1 + 8 * N1)) / 2) #the real size of the matrix, this will fail if x is not a vector with the values of
    #the diagonal
    A = Array{eltype(v),2}(undef, N, N)
    count = 1
    @boundscheck checkbounds(A, N, N)
    @inbounds begin
        for i = 1:N
            A[i, i] = 1.0
            for j = i+1:N
                A[i, j] = v[count]
                A[j, i] = symmetric_op(v[count])
                count += 1
            end
        end
    end
    return A
end

struct GERG2008 <: EmpiricHelmholtzModel
    components::Vector{String} #SingleParam
    acentricfactor::Vector{Float64} #SingleParam
    Mw::Vector{Float64} #SingleParam
    rhoc::Vector{Float64} #SingleParam
    vc::Vector{Float64} #SingleParam
    Tc::Vector{Float64} #SingleParam
    Pc::Vector{Float64} #SingleParam
    lb_v::Vector{Float64} #SingleParam
    ideal_iters::Vector{Vector{Int64}} #SingleParam{Int}
    nr::Matrix{Float64} #SingleParam
    zeta::Matrix{Float64} #SingleParam
    n0ik::Vector{Vector{Float64}} #VectorSingleParam
    t0ik::Vector{Vector{Float64}} #VectorSingleParam
    d0ik::Vector{Vector{Int64}} #VectorSingleParam
    c0ik::Vector{Vector{Int64}} #VectorSingleParam
    k_pol_ik::Vector{Int64}
    k_exp_ik::Vector{Int64}
    gamma_v::Matrix{Float64} #PairParam
    gamma_T::Matrix{Float64} #PairParam
    beta_v::Matrix{Float64} #PairParam
    beta_T::Matrix{Float64} #PairParam
    Aij_indices::Vector{Tuple{Int64,Int64,Int64}} #AssocParam ?
    fij::Vector{Float64} #SparsePairParam
    dijk::Vector{Vector{Int64}} #SparsePairParam
    tijk::Vector{Vector{Float64}} #SparsePairParam
    nijk::Vector{Vector{Float64}} #SparsePairParam
    etaijk::Vector{Vector{Float64}} #SparsePairParam
    epsijk::Vector{Vector{Float64}} #SparsePairParam
    betaijk::Vector{Vector{Float64}}
    gammaijk::Vector{Vector{Float64}}
    k_pol_ijk::Vector{Int64}
    k_exp_ijk::Vector{Int64}
end
function GERG2008_splitter(names::Vector{String})
    symvalues = GERG2008_consts.names
    xsel = Vector{Int}()
    for i in names
        i in symvalues || throw(ArgumentError("invalid name: $i"))
        push!(xsel, findall(isequal(i), symvalues)[1])

    end
    if iszero(length(xsel))
        return ArgumentError("no match found for input components. valid components are: $(GERG2008constsProps.names)")
    else
        return xsel
    end
end
function GERG2008(components::Vector{String})
    xsel = GERG2008_splitter(components)
    comps = components
    pc = GERG2008_consts.pc
    Tc = GERG2008_consts.Tc
    Mw = GERG2008_consts.Mw
    lb_v = GERG2008_consts.lb_volumes
    rhoc = GERG2008_consts.rhoc
    acentricfactor = GERG2008_consts.acentricfactor

    #critical pressure, Pa
    #iterations for the ideal part, for now i can´t design a general function
    #with alternating indices  in this form.
    iter_ideal = [4, 3, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 4, 2, 2, 3, 2, 0, 0]
    #original order
    iter4 = (findall(==(4), iter_ideal[xsel]))
    iter3 = (findall(==(3), iter_ideal[xsel]))
    iter2 = (findall(==(2), iter_ideal[xsel]))
    iter0 = (findall(==(0), iter_ideal[xsel]))
    ideal_iters = [iter4, iter3, iter2, iter0]

    nr = hcat(
        [19.597508817, -83.959667892, 3.00088, 0.76315, 0.00460, 8.74432, -4.46921], #methane
        [11.083407489, -22.202102428, 2.50031, 0.13732, -0.14660, 0.90066, 0.0],     #nitrogen
        [11.925152758, -16.118762264, 2.50002, 2.04452, -1.06044, 2.03366, 0.01393],  #carbon dioxide
        [24.675437527, -77.425313760, 3.00263, 4.33939, 1.23722, 13.19740, -6.01989], #ethane
        [31.602908195, -84.463284382, 3.02939, 6.60569, 3.19700, 19.19210, -8.37267], #propane
        [20.884143364, -91.638478026, 3.33944, 9.44893, 6.89406, 24.46180, 14.78240], #n-butane
        [20.413726078, -94.467620036, 3.06714, 8.97575, 5.25156, 25.14230, 16.13880], #iso-butane
        [14.536611217, -89.919548319, 3.00000, 8.95043, 21.83600, 33.40320, 0.0],     #n-pentane
        [15.449907693, -101.298172792, 3.00000, 11.76180, 20.11010, 33.16880, 0.0],   #isopentane
        [14.345969349, -96.165722367, 3.00000, 11.69770, 26.81420, 38.61640, 0.0],    #n-hexane
        [15.063786601, -97.345252349, 3.00000, 13.72660, 30.47070, 43.55610, 0.0],    #n-heptane
        [15.864687161, -97.370667555, 3.00000, 15.68650, 33.80290, 48.17310, 0.0],   #n-octane
        [16.313913248, -102.160247463, 3.00000, 18.02410, 38.12350, 53.34150, 0.0],  #n-nonane
        [15.870791919, -108.858547525, 3.00000, 21.00690, 43.49310, 58.36570, 0.0],       #n-decane
        [13.796443393, -175.864487294, 1.47906, 0.95806, 0.45444, 1.56039, -1.37560], #hydrogen
        [10.001843586, -14.996095135, 2.50146, 1.07558, 1.01334, 0.0, 0.0],           #oxygen
        [10.813340744, -19.834733959, 2.50055, 1.02865, 0.00493, 0.0, 0.0],           #carbon monoxide
        [8.203520690, -11.996306443, 3.00392, 0.01059, 0.98763, 3.06904, 0.0],            #water
        [9.336197742, -16.266508995, 3.00000, 3.11942, 1.00243, 0.0, 0.0],            #hydrogen sulfide
        [13.628409737, -143.470759602, 1.50000, 0.0, 0.0, 0.0, 0.0],                  #helium
        [8.316631500, -4.946502600, 1.50000, 0.0, 0.0, 0.0, 0.0],                      #argon
    )

    zeta = hcat(
        [4.306474465, 0.936220902, 5.577233895, 5.722644361],                  #methane
        [5.251822620, -5.393067706, 13.788988208, 0],                        #nitrogen
        [3.022758166, -2.844425476, 1.589964364, 1.121596090],                 #carbon dioxide
        [1.831882406, 0.731306621, 3.378007481, 3.508721939],                  #ethane
        [1.297521801, 0.543210978, 2.583146083, 2.777773271],                  #propane
        [1.101487798, 0.431957660, 4.502440459, 2.124516319],                  #n-butane
        [1.074673199, 0.485556021, 4.671261865, 2.191583480],                  #isobutane
        [0.380391739, 1.789520971, 3.777411113, 0],                          #n-pentane
        [0.635392636, 1.977271641, 4.169371131, 0],                          #isopentane
        [0.359036667, 1.691951873, 3.596924107, 0],                          #n-hexane
        [0.314348398, 1.548136560, 3.259326458, 0],                          #n-heptane
        [0.279143540, 1.431644769, 2.973845992, 0],                          #n-octane
        [0.263819696, 1.370586158, 2.848860483, 0],                          #n-nonane
        [0.267034159, 1.353835195, 2.833479035, 0],                          #n-decane
        [6.891654113, 9.847634830, 49.765290750, 50.367279301],                #hydrogen
        [14.461722565, 7.223325463, 0, 0],                                 #oxygen
        [11.669802800, 5.302762306, 0, 0],                                 #carbon monoxide
        [0.415386589, 1.763895929, 3.874803739, 0],                          #water
        [4.914580541, 2.270653980, 0, 0],                                  #hydrogen sulfide
        [0.0, 0, 0.0, 0.0],                                                  #helium
        [0.0, 0, 0.0, 0.0],                                                   #argon
    )

    n0ik = [
        [
            0.57335704239162,
            -0.16760687523730e1,
            0.23405291834916,
            -0.21947376343441,
            0.16369201404128e-1,
            0.15004406389280e-1,
            0.98990489492918e-1,
            0.58382770929055,
            -0.74786867560390,
            0.30033302857974,
            0.20985543806568,
            -0.18590151133061e-1,
            -0.15782558339049,
            0.12716735220791,
            -0.32019743894346e-1,
            -0.68049729364536e-1,
            0.24291412853736e-1,
            0.51440451639444e-2,
            -0.19084949733532e-1,
            0.55229677241291e-2,
            -0.44197392976085e-2,
            0.40061416708429e-1,
            -0.33752085907575e-1,
            -0.25127658213357e-2,
        ],#methane
        [
            0.59889711801201,
            -0.16941557480731e1,
            0.24579736191718,
            -0.23722456755175,
            0.17954918715141e-1,
            0.14592875720215e-1,
            0.10008065936206,
            0.73157115385532,
            -0.88372272336366,
            0.31887660246708,
            0.20766491728799,
            -0.19379315454158e-1,
            -0.16936641554983,
            0.13546846041701,
            -0.33066712095307e-1,
            -0.60690817018557e-1,
            0.12797548292871e-1,
            0.58743664107299e-2,
            -0.18451951971969e-1,
            0.47226622042472e-2,
            -0.52024079680599e-2,
            0.43563505956635e-1,
            -0.36251690750939e-1,
            -0.28974026866543e-2,
        ],#nitrogen
        [
            0.52646564804653,
            -0.14995725042592e1,
            0.27329786733782,
            0.12949500022786,
            0.15404088341841,
            -0.58186950946814,
            -0.18022494838296,
            -0.95389904072812e-1,
            -0.80486819317679e-2,
            -0.35547751273090e-1,
            -0.28079014882405,
            -0.82435890081677e-1,
            0.10832427979006e-1,
            -0.67073993161097e-2,
            -0.46827907600524e-2,
            -0.28359911832177e-1,
            0.19500174744098e-1,
            -0.21609137507166,
            0.43772794926972,
            -0.22130790113593,
            0.15190189957331e-1,
            -0.15380948953300e-1,
        ],#carbon_dioxide
        [
            0.63596780450714,
            -0.17377981785459e1,
            0.28914060926272,
            -0.33714276845694,
            0.22405964699561e-1,
            0.15715424886913e-1,
            0.11450634253745,
            0.10612049379745e1,
            -0.12855224439423e1,
            0.39414630777652,
            0.31390924682041,
            -0.21592277117247e-1,
            -0.21723666564905,
            -0.28999574439489,
            0.42321173025732,
            0.46434100259260e-1,
            -0.13138398329741,
            0.11492850364368e-1,
            -0.33387688429909e-1,
            0.15183171583644e-1,
            -0.47610805647657e-2,
            0.46917166277885e-1,
            -0.39401755804649e-1,
            -0.32569956247611e-2,
        ],#ethane
        [
            0.10403973107358e1,
            -0.28318404081403e1,
            0.84393809606294,
            -0.76559591850023e-1,
            0.94697373057280e-1,
            0.24796475497006e-3,
            0.27743760422870,
            -0.43846000648377e-1,
            -0.26991064784350,
            -0.69313413089860e-1,
            -0.29632145981653e-1,
            0.14040126751380e-1,
        ],#propane
        [
            0.10626277411455e1,
            -0.28620951828350e1,
            0.88738233403777,
            -0.12570581155345,
            0.10286308708106,
            0.25358040602654e-3,
            0.32325200233982,
            -0.37950761057432e-1,
            -0.32534802014452,
            -0.79050969051011e-1,
            -0.20636720547775e-1,
            0.57053809334750e-2,
        ],#n-butane
        [
            0.10429331589100e1,
            -0.28184272548892e1,
            0.86176232397850,
            -0.10613619452487,
            0.98615749302134e-1,
            0.23948208682322e-3,
            0.30330004856950,
            -0.41598156135099e-1,
            -0.29991937470058,
            -0.80369342764109e-1,
            -0.29761373251151e-1,
            0.13059630303140e-1,
        ],#isobutane
        [
            0.10968643098001e1,
            -0.29988888298061e1,
            0.99516886799212,
            -0.16170708558539,
            0.11334460072775,
            0.26760595150748e-3,
            0.40979881986931,
            -0.40876423083075e-1,
            -0.38169482469447,
            -0.10931956843993,
            -0.32073223327990e-1,
            0.16877016216975e-1,
        ],#n-pentane
        [
            0.10963e1,
            -0.30402e1,
            0.10317e1,
            -0.15410,
            0.11535,
            0.29809e-3,
            0.39571,
            -0.45881e-1,
            -0.35804,
            -0.10107,
            -0.35484e-1,
            0.18156e-1,
        ],#isopentane
        [
            0.10553238013661e1,
            -0.26120615890629e1,
            0.76613882967260,
            -0.29770320622459,
            0.11879907733358,
            0.27922861062617e-3,
            0.46347589844105,
            0.11433196980297e-1,
            -0.48256968738131,
            -0.93750558924659e-1,
            -0.67273247155994e-2,
            -0.51141583585428e-2,
        ],#n-hexane
        [
            0.10543747645262e1,
            -0.26500681506144e1,
            0.81730047827543,
            -0.30451391253428,
            0.12253868710800,
            0.27266472743928e-3,
            0.49865825681670,
            -0.71432815084176e-3,
            -0.54236895525450,
            -0.13801821610756,
            -0.61595287380011e-2,
            0.48602510393022e-3,
        ],#n-heptane
        [
            0.10722544875633e1,
            -0.24632951172003e1,
            0.65386674054928,
            -0.36324974085628,
            0.12713269626764,
            0.30713572777930e-3,
            0.52656856987540,
            0.19362862857653e-1,
            -0.58939426849155,
            -0.14069963991934,
            -0.78966330500036e-2,
            0.33036597968109e-2,
        ],#n-octane
        [
            0.11151e1,
            -0.27020e1,
            0.83416,
            -0.38828,
            0.13760,
            0.28185e-3,
            0.62037,
            0.15847e-1,
            -0.61726,
            -0.15043,
            -0.12982e-1,
            0.44325e-2,
        ],#n-nonane
        [
            0.10461e1,
            -0.24807e1,
            0.74372,
            -0.52579,
            0.15315,
            0.32865e-3,
            0.84178,
            0.55424e-1,
            -0.73555,
            -0.18507,
            -0.20775e-1,
            0.12335e-1,
        ],#n-decane
        [
            0.53579928451252e1,
            -0.62050252530595e1,
            0.13830241327086,
            -0.71397954896129e-1,
            0.15474053959733e-1,
            -0.14976806405771,
            -0.26368723988451e-1,
            0.56681303156066e-1,
            -0.60063958030436e-1,
            -0.45043942027132,
            0.42478840244500,
            -0.21997640827139e-1,
            -0.10499521374530e-1,
            -0.28955902866816e-2,
        ],#hydrogen
        [
            0.88878286369701,
            -0.24879433312148e1,
            0.59750190775886,
            0.96501817061881e-2,
            0.71970428712770e-1,
            0.22337443000195e-3,
            0.18558686391474,
            -0.38129368035760e-1,
            -0.15352245383006,
            -0.26726814910919e-1,
            -0.25675298677127e-1,
            0.95714302123668e-2,
        ],#oxygen
        [
            0.90554,
            -0.24515e1,
            0.53149,
            0.24173e-1,
            0.72156e-1,
            0.18818e-3,
            0.19405,
            -0.43268e-1,
            -0.12778,
            -0.27896e-1,
            -0.34154e-1,
            0.16329e-1,
        ],#carbon_monoxide
        [
            0.82728408749586,
            -0.18602220416584e1,
            -0.11199009613744e1,
            0.15635753976056,
            0.87375844859025,
            -0.36674403715731,
            0.53987893432436e-1,
            0.10957690214499e1,
            0.53213037828563e-1,
            0.13050533930825e-1,
            -0.41079520434476,
            0.14637443344120,
            -0.55726838623719e-1,
            -0.11201774143800e-1,
            -0.66062758068099e-2,
            0.46918522004538e-2,
        ],#water
        [
            0.87641,
            -0.20367e1,
            0.21634,
            -0.50199e-1,
            0.66994e-1,
            0.19076e-3,
            0.20227,
            -0.45348e-2,
            -0.22230,
            -0.34714e-1,
            -0.14885e-1,
            0.74154e-2,
        ],#hydrogen_sulfide
        [
            -0.45579024006737,
            0.12516390754925e1,
            -0.15438231650621e1,
            0.20467489707221e-1,
            -0.34476212380781,
            -0.20858459512787e-1,
            0.16227414711778e-1,
            -0.57471818200892e-1,
            0.19462416430715e-1,
            -0.33295680123020e-1,
            -0.10863577372367e-1,
            -0.22173365245954e-1,
        ],#helium
        [
            0.85095714803969,
            -0.24003222943480e1,
            0.54127841476466,
            0.16919770692538e-1,
            0.68825965019035e-1,
            0.21428032815338e-3,
            0.17429895321992,
            -0.33654495604194e-1,
            -0.13526799857691,
            -0.16387350791552e-1,
            -0.24987666851475e-1,
            0.88769204815709e-2,
        ],#argon
    ]

    #most_compounds
    c0ik1 = [1, 1, 2, 2, 3, 3]
    d0ik1 = [1, 1, 1, 2, 3, 7, 2, 5, 1, 4, 3, 4]
    t0ik1 = [
        0.250,
        1.125,
        1.500,
        1.375,
        0.250,
        0.875,
        0.625,
        1.750,
        3.625,
        3.625,
        14.500,
        12.000,
    ]
    k_exp_ik1 = length(c0ik1)
    k_pol_ik1 = length(d0ik1) - k_exp_ik1

    #methane,nitrogen,ethane
    c0ik2 = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 6, 6, 6, 6]
    d0ik2 = [1, 1, 2, 2, 4, 4, 1, 1, 1, 2, 3, 6, 2, 3, 3, 4, 4, 2, 3, 4, 5, 6, 6, 7]
    t0ik2 = [
        0.125,
        1.125,
        0.375,
        1.125,
        0.625,
        1.500,
        0.625,
        2.625,
        2.750,
        2.125,
        2.000,
        1.750,
        4.500,
        4.750,
        5.000,
        4.000,
        4.500,
        7.500,
        14.000,
        11.500,
        26.000,
        28.000,
        30.000,
        16.000,
    ]
    k_exp_ik2 = length(c0ik2)
    k_pol_ik2 = length(d0ik2) - k_exp_ik2

    #carbon_Dioxide
    c0ik3 = [1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 5, 5, 5, 6, 6]
    d0ik3 = [1, 1, 2, 3, 3, 3, 4, 5, 6, 6, 1, 4, 1, 1, 3, 3, 4, 5, 5, 5, 5, 5]
    t0ik3 = [
        0.000,
        1.250,
        1.625,
        0.375,
        0.375,
        1.375,
        1.125,
        1.375,
        0.125,
        1.625,
        3.750,
        3.500,
        7.500,
        8.000,
        6.000,
        16.000,
        11.000,
        24.000,
        26.000,
        28.000,
        24.000,
        26.000,
    ]
    k_exp_ik3 = length(c0ik3)
    k_pol_ik3 = length(d0ik3) - k_exp_ik3
    #hydrogen
    c0ik4 = [1, 1, 1, 1, 2, 2, 3, 3, 5]
    d0ik4 = [1, 1, 2, 2, 4, 1, 5, 5, 5, 1, 1, 2, 5, 1]
    t0ik4 = [0.500,0.625, 0.375, 0.625, 1.125, 2.625, 0.000, 0.250, 1.375, 4.000, 4.250, 5.000, 8.000, 8.000]
    k_exp_ik4 = length(c0ik4)
    k_pol_ik4 = length(d0ik4) - k_exp_ik4

    #water
    c0ik5 = [1, 1, 1, 2, 2, 2, 3, 5, 5]
    d0ik5 = [1, 1, 1, 2, 2, 3, 4, 1, 5, 5, 1, 2, 4, 4, 1, 1]
    t0ik5 = [
        0.500,
        1.250,
        1.875,
        0.125,
        1.500,
        1.000,
        0.750,
        1.500,
        0.625,
        2.625,
        5.000,
        4.000,
        4.500,
        3.000,
        4.000,
        6.000,
    ]
    k_exp_ik5 = length(c0ik5)
    k_pol_ik5 = length(d0ik5) - k_exp_ik5

    #helium
    c0ik6 = [1, 1, 1, 1, 1, 2, 3, 3]
    d0ik6 = [1, 1, 1, 4, 1, 3, 5, 5, 5, 2, 1, 2]
    t0ik6 =
        [0.000, 0.125, 0.750, 1.000, 0.750, 2.625, 0.125, 1.250, 2.000, 1.000, 4.500, 5.000]
    k_exp_ik6 = length(c0ik6)
    k_pol_ik6 = length(d0ik6) - k_exp_ik6

    c0ik = [
        c0ik2,
        c0ik2,
        c0ik3,
        c0ik2,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik1,
        c0ik4,
        c0ik1,
        c0ik1,
        c0ik5,
        c0ik1,
        c0ik6,
        c0ik1,
    ]

    d0ik = [
        d0ik2,
        d0ik2,
        d0ik3,
        d0ik2,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik1,
        d0ik4,
        d0ik1,
        d0ik1,
        d0ik5,
        d0ik1,
        d0ik6,
        d0ik1,
    ]

    t0ik = [
        t0ik2,
        t0ik2,
        t0ik3,
        t0ik2,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik1,
        t0ik4,
        t0ik1,
        t0ik1,
        t0ik5,
        t0ik1,
        t0ik6,
        t0ik1,
    ]

    k_pol_ik = vcat(
        k_pol_ik2,
        k_pol_ik2,
        k_pol_ik3,
        k_pol_ik2,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik4,
        k_pol_ik1,
        k_pol_ik1,
        k_pol_ik5,
        k_pol_ik1,
        k_pol_ik6,
        k_pol_ik1,
    )

    k_exp_ik = vcat(
        k_exp_ik2,
        k_exp_ik2,
        k_exp_ik3,
        k_exp_ik2,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik4,
        k_exp_ik1,
        k_exp_ik1,
        k_exp_ik5,
        k_exp_ik1,
        k_exp_ik6,
        k_exp_ik1,
    )

    vector_beta_v = [
        0.9987213770,
        0.9995180720,
        0.9975478660,
        1.0048270700,
        0.9791059720,
        1.0112403880,
        0.9483301200,
        1.0000000000,
        0.9580152940,
        0.9620508310,
        0.9947406030,
        1.0028522870,
        1.0330862920,
        1.0000000000,
        1.0000000000,
        0.9973407720,
        1.0127831690,
        1.0125990870,
        1.0000000000,
        1.0346302590,
        0.9777946340,
        0.9788801680,
        0.9744246810,
        0.9960826100,
        0.9864158300,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9725320650,
        1.0000000000,
        1.0000000000,
        0.9995217700,
        0.9103942490,
        0.9695010550,
        1.0041664120,
        1.0025257180,
        0.9968980040,
        1.1747609230,
        1.0765518820,
        1.0243114980,
        1.0607931040,
        1.0000000000,
        1.2054699760,
        1.0261693730,
        1.0000000000,
        1.0001511320,
        0.9041421590,
        1.0000000000,
        1.0000000000,
        0.9490559590,
        0.9066305640,
        0.8466475610,
        1.0083924280,
        0.9976072770,
        0.9991572050,
        1.0000000000,
        0.9938510090,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0074697260,
        1.0000000000,
        0.9956762580,
        0.9253671710,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0108179090,
        1.0000000000,
        1.0000000000,
        0.9997958680,
        0.9992431460,
        1.0449194310,
        1.0404592890,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9841042270,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9368112190,
        1.0000000000,
        1.0000000000,
        1.0008804640,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9769519680,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9081131630,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0129944310,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9846132030,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0015163710,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.7544739580,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.8289671640,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.6953583820,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9751877660,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        0.9997468470,
        1.0000000000,
        0.7956603920,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
        1.0000000000,
    ]

    beta_v = gerg_betamatrix_from_vector(vector_beta_v, inv)

    vector_beta_T = [
        0.998098830,
        1.022624490,
        0.996336508,
        0.989680305,
        0.994174910,
        0.980315756,
        0.992127525,
        1.0,
        0.981844797,
        0.977431529,
        0.957473785,
        0.947716769,
        0.937777823,
        1.0,
        1.0,
        0.987411732,
        1.063333913,
        1.011090031,
        1.0,
        0.990954281,
        1.005894529,
        1.007671428,
        1.002677329,
        0.994515234,
        0.992868130,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.956379450,
        0.957934447,
        0.946134337,
        1.0,
        1.0,
        0.997190589,
        1.004692366,
        0.692868765,
        0.999069843,
        1.013871147,
        1.033620538,
        1.018171004,
        1.023339824,
        1.027000795,
        1.019180957,
        1.0,
        1.011806317,
        1.029690780,
        1.007688620,
        1.020028790,
        0.942320195,
        1.0,
        1.0,
        0.997372205,
        1.016034583,
        0.768377630,
        0.996512863,
        0.996199694,
        0.999130554,
        1.0,
        0.998688946,
        1.0,
        1.0,
        1.0,
        0.984068272,
        1.0,
        0.970918061,
        0.932969831,
        1.0,
        1.0,
        1.0,
        0.990197354,
        1.0,
        1.0,
        1.000310289,
        0.998012298,
        0.996484021,
        0.994364425,
        1.0,
        1.0,
        1.0,
        1.0,
        0.985331233,
        1.0,
        1.0,
        1.0,
        1.0,
        0.992573556,
        1.0,
        1.0,
        1.000077547,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.993688386,
        1.0,
        1.0,
        1.0,
        1.0,
        0.985962886,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.974550548,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.962006651,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.997641010,
        1.0,
        1.0,
        1.0,
        1.0,
        0.985891113,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.988937417,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.064818089,
        1.0,
        1.049594632,
        0.897162268,
        0.973091413,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.000023103,
        1.0,
        1.025536736,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
    ]

    beta_T = gerg_betamatrix_from_vector(vector_beta_T, inv)

    gamma_ijv_vec = [
        1.013950311,
        1.002806594,
        1.006617867,
        1.038470657,
        1.045375122,
        1.054319053,
        1.124508039,
        1.343685343,
        1.052643846,
        1.156655935,
        1.116549372,
        1.141895355,
        1.146089637,
        1.018702573,
        1.0,
        1.006102927,
        1.585018334,
        1.040161207,
        0.881405683,
        1.014678542,
        1.047578256,
        1.042352891,
        1.081025408,
        1.146949309,
        1.100576129,
        1.078877166,
        1.154135439,
        1.195952177,
        1.404554090,
        1.186067025,
        1.100405929,
        1.0,
        0.970115357,
        1.008690943,
        1.094749685,
        0.997082328,
        1.256844157,
        0.932629867,
        1.002212182,
        1.032876701,
        1.047596298,
        1.222437324,
        1.081909003,
        1.068406078,
        1.116793198,
        0.851343711,
        1.164585914,
        1.104043935,
        0.973386152,
        1.183394668,
        1.152792550,
        1.0,
        1.0,
        1.542328793,
        1.024085837,
        0.864141549,
        1.029205465,
        1.003034720,
        1.006179146,
        1.006616886,
        1.026085655,
        1.045439935,
        1.169701102,
        1.057666085,
        1.071917985,
        1.143534730,
        1.098361281,
        1.106072040,
        1.0,
        1.201417898,
        1.0,
        1.030988277,
        1.0,
        1.0,
        1.003264179,
        1.001156119,
        1.019921513,
        0.999432118,
        1.057872566,
        1.079648053,
        1.102764612,
        1.199769134,
        1.053040574,
        1.074006110,
        1.0,
        1.108143673,
        1.011759763,
        1.010593999,
        1.0,
        1.0,
        1.000414440,
        1.018159650,
        1.002728434,
        1.034995284,
        1.019174227,
        1.046905515,
        1.049219137,
        1.027845529,
        1.232939523,
        1.0,
        1.084740904,
        1.223638763,
        1.033366041,
        1.0,
        1.214638734,
        1.002779804,
        1.002284353,
        1.010493989,
        1.021668316,
        1.032807063,
        1.047298475,
        1.060243344,
        1.147595688,
        1.0,
        1.087272232,
        1.0,
        0.988591117,
        1.0,
        1.0,
        1.000024335,
        1.002480637,
        1.008972412,
        1.069223964,
        1.034910633,
        1.016370338,
        1.188334783,
        1.0,
        1.119954454,
        0.956677310,
        1.076539234,
        1.0,
        1.0,
        1.002995876,
        1.009928206,
        1.017880545,
        1.028994325,
        1.039372957,
        1.184340443,
        1.0,
        1.116694577,
        1.0,
        0.835763343,
        1.0,
        1.0,
        1.001508227,
        1.006268954,
        1.020761680,
        1.013511439,
        1.243461678,
        1.0,
        1.155145836,
        1.170217596,
        1.339283552,
        1.0,
        1.0,
        1.006767176,
        1.001370076,
        1.002972346,
        1.159131722,
        1.0,
        1.190354273,
        1.0,
        1.087956749,
        1.0,
        1.0,
        1.001357085,
        1.002553544,
        1.305249405,
        1.0,
        1.219206702,
        0.599484191,
        1.0,
        1.0,
        1.0,
        1.000810520,
        1.342647661,
        1.0,
        1.252151449,
        1.0,
        1.082905109,
        1.0,
        1.0,
        1.120233729,
        1.0,
        0.870184960,
        0.551405318,
        1.171714677,
        1.0,
        1.0,
        1.0,
        1.121416201,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.143174289,
        1.0,
        1.0,
        0.993907223,
        1.0,
        1.101731308,
        1.0,
        1.159720623,
        1.014832832,
        1.0,
        1.038993495,
        1.0,
        1.0,
        1.0,
    ]

    gamma_v = gerg_betamatrix_from_vector(gamma_ijv_vec, identity)

    gamma_ijT_vec = [
        0.979273013,
        0.975665369,
        1.049707697,
        1.098655531,
        1.171607691,
        1.161117729,
        1.249173968,
        1.188899743,
        1.330570181,
        1.379850328,
        1.449245409,
        1.528532478,
        1.568231489,
        1.352643115,
        0.950000000,
        0.987473033,
        0.775810513,
        0.961155729,
        3.159776855,
        0.989843388,
        1.107654104,
        1.098650964,
        1.201264026,
        1.304886838,
        1.284462634,
        1.419029041,
        1.381770770,
        1.472607971,
        1.520975334,
        1.733280051,
        1.749119996,
        1.822157123,
        1.175696583,
        0.993425388,
        0.968808467,
        0.995157044,
        0.960174200,
        1.471831580,
        0.990034831,
        0.900949530,
        0.908772477,
        0.911498231,
        0.929982936,
        0.979217302,
        0.961218039,
        1.038675574,
        1.046169823,
        1.074455386,
        1.140671202,
        1.145512213,
        1.782924792,
        1.0,
        1.0,
        0.775453996,
        0.926018880,
        3.207456948,
        1.050971635,
        1.014730190,
        1.034832749,
        1.033283811,
        1.066665676,
        1.021150247,
        1.092177796,
        1.134532014,
        1.168636194,
        1.056033030,
        1.237191558,
        1.902008495,
        1.0,
        1.069224728,
        1.0,
        0.902736660,
        1.0,
        1.0,
        1.007392782,
        1.005250774,
        1.008344412,
        1.003269500,
        1.025657518,
        1.050044169,
        1.063694129,
        1.109973833,
        1.140905252,
        2.308215191,
        1.0,
        1.197564208,
        0.600340961,
        0.905829247,
        1.0,
        1.0,
        1.001432824,
        1.002143640,
        1.000792201,
        1.009157060,
        1.021283378,
        1.033180106,
        1.014096448,
        1.076466918,
        2.509259945,
        1.0,
        1.173916162,
        0.615512682,
        0.926156602,
        1.0,
        1.245039498,
        1.002495889,
        1.001835788,
        1.006018054,
        1.009885760,
        1.013945424,
        1.017817492,
        1.021624748,
        1.895305393,
        1.0,
        1.161390082,
        1.0,
        0.937130844,
        1.0,
        1.0,
        1.000050537,
        1.000761237,
        1.002441051,
        1.016422347,
        1.103421755,
        1.049035838,
        2.013859174,
        1.0,
        1.206043295,
        0.447666011,
        0.959065662,
        1.0,
        1.0,
        1.001204174,
        1.003194615,
        1.005647480,
        1.008191499,
        1.010825138,
        1.996386669,
        1.0,
        1.199326059,
        1.0,
        0.982651529,
        1.0,
        1.0,
        0.999762786,
        1.001633952,
        1.055369591,
        1.028939539,
        3.021197546,
        1.0,
        1.233272781,
        0.569681333,
        0.956075596,
        1.0,
        1.0,
        0.998793111,
        1.001150096,
        1.002229938,
        3.169143057,
        1.0,
        1.256123503,
        1.0,
        1.013453092,
        1.0,
        1.0,
        1.000235044,
        1.007186267,
        2.191555216,
        1.0,
        1.276565536,
        0.662072469,
        1.0,
        1.0,
        1.0,
        1.000182392,
        2.234354040,
        1.0,
        1.294070556,
        1.0,
        1.086557826,
        1.0,
        1.0,
        3.786003724,
        1.0,
        1.803567587,
        0.740416402,
        1.103693489,
        1.0,
        1.0,
        1.0,
        1.377504607,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.964767932,
        1.0,
        1.0,
        0.990430423,
        1.0,
        1.022749748,
        1.0,
        0.954215746,
        0.940587083,
        1.0,
        1.070941866,
        1.0,
        1.0,
        1.0,
    ]

    gamma_T = gerg_betamatrix_from_vector(gamma_ijT_vec, identity)

    indice1 = [1, 1, 1, 1, 1, 1, 1, 2, 2, 4, 4, 4, 5, 5, 6]
    indice2 = [2, 3, 4, 5, 6, 7, 15, 3, 4, 5, 6, 7, 6, 7, 7]

    value_eq = 1:length(indice1)
    Aij_indices_mat = zeros(Int64, 21, 21)
    for i = 1:length(value_eq)
        Aij_indices_mat[indice1[i], indice2[i]] = value_eq[i]
    end

    Aij_indices_mat = Aij_indices_mat[xsel, xsel]
    Aij_indices2 = findall(!iszero, Aij_indices_mat)
    Aij_indices = Vector{Tuple{Int64,Int64,Int64}}(undef,0)
    for i = 1:length(Aij_indices2)
        push!(Aij_indices,(Aij_indices2[i][1], Aij_indices2[i][2], Aij_indices_mat[Aij_indices2[i]]))
    end

    fij = [
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        0.771035405688,
        1.0,
        1.0,
        1.0,
        0.130424765150,
        0.281570073085,
        0.260632376098,
        0.312572600489e-1,
        -0.551609771024e-1,
        -0.551240293009e-1,
    ]

    dijk = [
        [1, 4, 1, 2, 2, 2, 2, 2, 3],
        [1, 2, 3, 1, 2, 3],
        [3, 4, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3],
        [3, 3, 4, 4, 4, 1, 1, 1, 2],
        [2, 3, 1, 1, 1, 2],
        [2, 2, 3, 1, 2, 2],
        [1, 3, 3, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
        [1, 1, 1, 2, 2, 3, 3, 4, 4, 4],
    ]

    tijk = [
        [0.000, 1.850, 7.850, 5.400, 0.000, 0.750, 2.800, 4.450, 4.250],
        [2.600, 1.950, 0.000, 3.950, 7.950, 8.000],
        [
            0.650,
            1.550,
            3.100,
            5.900,
            7.050,
            3.350,
            1.200,
            5.800,
            2.700,
            0.450,
            0.550,
            1.950,
        ],
        [1.850, 3.950, 0.000, 1.850, 3.850, 5.250, 3.850, 0.200, 6.500],
        [1.850, 1.400, 3.200, 2.500, 8.000, 3.750],
        [0.000, 0.050, 0.000, 3.650, 4.900, 4.450],
        [2.0, -1.0, 1.75, 1.4],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
        [1.000, 1.550, 1.700, 0.250, 1.350, 0.000, 1.250, 0.000, 0.700, 5.400],
    ]

    nijk = [
        [
            -0.98038985517335e-2,
            0.42487270143005e-3,
            -0.34800214576142e-1,
            -0.13333813013896,
            -0.11993694974627e-1,
            0.69243379775168e-1,
            -0.31022508148249,
            0.24495491753226,
            0.22369816716981,
        ],
        [
            -0.10859387354942,
            0.80228576727389e-1,
            -0.93303985115717e-2,
            0.40989274005848e-1,
            -0.24338019772494,
            0.23855347281124,
        ],
        [
            -0.80926050298746e-3,
            -0.75381925080059e-3,
            -0.41618768891219e-1,
            -0.23452173681569,
            0.14003840584586,
            0.63281744807738e-1,
            -0.34660425848809e-1,
            -0.23918747334251,
            0.19855255066891e-2,
            0.61777746171555e1,
            -0.69575358271105e1,
            0.10630185306388e1,
        ],
        [
            0.13746429958576e-1,
            -0.74425012129552e-2,
            -0.45516600213685e-2,
            -0.54546603350237e-2,
            0.23682016824471e-2,
            0.18007763721438,
            -0.44773942932486,
            0.19327374888200e-1,
            -0.30632197804624,
        ],
        [
            0.28661625028399,
            -0.10919833861247,
            -0.11374032082270e1,
            0.76580544237358,
            0.42638000926819e-2,
            0.17673538204534,
        ],
        [
            -0.47376518126608,
            0.48961193461001,
            -0.57011062090535e-2,
            -0.19966820041320,
            -0.69411103101723,
            0.69226192739021,
        ],
        [
            -0.25157134971934,
            -0.62203841111983e-2,
            0.88850315184396e-1,
            -0.35592212573239e-1,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
        [
            0.25574776844118e1,
            -0.79846357136353e1,
            0.47859131465806e1,
            -0.73265392369587,
            0.13805471345312e1,
            0.28349603476365,
            -0.49087385940425,
            -0.10291888921447,
            0.11836314681968,
            0.55527385721943e-4,
        ],
    ]

    etaijk = [
        [1.000, 1.000, 0.250, 0.000, 0.000, 0.000, 0.000],
        [1.000, 0.500, 0.000],
        [1.000, 1.000, 1.000, 0.875, 0.750, 0.500, 0.000, 0.000, 0.000, 0.000],
        [0.250, 0.250, 0.000, 0.000],
        [0.250, 0.250, 0.000, 0.000],
        [1.000, 1.000, 0.875],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
    ]

    epsijk = [
        [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
    ]

    betaijk = [
        [1.000, 1.000, 2.500, 3.000, 3.000, 3.000, 3.000],
        [1.000, 2.000, 3.000],
        [1.000, 1.000, 1.000, 1.250, 1.500, 2.000, 3.000, 3.000, 3.000, 3.000],
        [0.750, 1.000, 2.000, 3.000],
        [0.750, 1.000, 2.000, 3.000],
        [1.000, 1.000, 1.250],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
    ]

    gammaijk = [
        [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5, 0.5],
        [0.5, 0.5, 0.5],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
        Float64[],
    ]

    length_k_ijk = length(dijk)
    k_pol_ijk = Vector{Int64}(undef, length_k_ijk)
    k_exp_ijk = Vector{Int64}(undef, length_k_ijk)
    for i = 1:length_k_ijk
        k_exp_ijk[i] = length(etaijk[i])
        k_pol_ijk[i] = length(dijk[i]) - k_exp_ijk[i]
    end

    model = GERG2008(
        comps, #
        acentricfactor[xsel],#
        Mw[xsel], #
        rhoc[xsel], #
        1 ./ rhoc[xsel], #
        Tc[xsel], #
        pc[xsel], #
        lb_v[xsel], #
        ideal_iters, #
        nr[:, xsel], #
        zeta[:, xsel], #
        n0ik[xsel], #
        t0ik[xsel], #
        d0ik[xsel],#
        c0ik[xsel], #
        k_pol_ik[xsel],
        k_exp_ik[xsel],
        gamma_v[xsel, xsel], #
        gamma_T[xsel, xsel], #
        beta_v[xsel, xsel], #
        beta_T[xsel, xsel], #
        Aij_indices,
        fij,
        dijk,
        tijk,
        nijk,
        etaijk,
        epsijk,
        betaijk,
        gammaijk,
        k_pol_ijk,
        k_exp_ijk,
    )
    return model
end
=#


text:
@inline function signbit(x::Tuple{T,T}) where T<:AbstractFloat
    return signbit(x[1])
end

@inline function (-)(x::Tuple{T,T}) where T<:AbstractFloat
    return -x[1], -x[2]
end

@inline function abs(x::Tuple{T,T}) where T<:AbstractFloat
    return signbit(x[1]) ? -x : x
end

@inline function flipsign(x::Tuple{T,T}, y::Tuple{T,T}) where T<:AbstractFloat
    return signbit(y[1]) ? -x : x
end

@inline function copysign(x::Tuple{T,T}, y::Tuple{T,T}) where T<:AbstractFloat
    return signbit(y[1]) ? -abs(x) : abs(x)
end

function (+)(x::Tuple{T,T}, y::Tuple{T,T}) where T<:AbstractFloat
    hi, lo = two_sum(x[1], y[1])
    lohi, lolo = two_sum(x[2], y[2])
    lohi += lo
    hi, lo = fast_two_sum(hi, lohi)
    lolo += lo
    hi, lo = fast_two_sum(hi, lolo)
    return hi, lo
end

function (-)(x::Tuple{T,T}, y::Tuple{T,T}) where T<:AbstractFloat
    hi, lo = two_diff(x[1], y[1])
    lohi, lolo = two_diff(x[2], y[2])
    lohi += lo
    hi, lo = fast_two_sum(hi, lohi)
    lolo += lo
    hi, lo = fast_two_sum(hi, lolo)
    return hi, lo
end

function (*)(x::Tuple{T,T}, y::Tuple{T,T}) where T<:AbstractFloat
    lo = two_prod(x[2], y[2])
    hi = two_prod(x[1], y[2])
    hilo = lo + hi
    hi = two_prod(x[2], y[1])
    hilo = hilo + hi
    hi = two_prod(x[1], y[1])
    hilo = hilo + hi
    return hilo
end

ext: jl


text:
export gridData,
       ROIData,
       generateInputData


text:
using Oceananigans.Grids: LatitudeLongitudeGrid, ConformalCubedSphereFaceGrid
using Oceananigans.Operators: Δx_qᶜᶠᶜ, Δy_qᶠᶜᶜ, Δxᶠᶜᶜ, Δyᶜᶠᶜ, hack_sind

# Our two Coriolis schemes are energy-conserving or enstrophy-conserving
# with a "vector invariant" momentum advection scheme, but not with a "flux form"
# or "conservation form" advection scheme (which does not currently exist for
# curvilinear grids).
struct VectorInvariantEnergyConserving end
struct VectorInvariantEnstrophyConserving end

"""
    HydrostaticSphericalCoriolis{FT} <: AbstractRotation

A parameter object for constant rotation around a vertical axis.
"""
struct HydrostaticSphericalCoriolis{S, FT} <: AbstractRotation
    rotation_rate :: FT
    scheme :: S
end

"""
    HydrostaticSphericalCoriolis([FT=Float64;] rotation_rate=Ω_Earth, scheme=VectorInvariantEnergyConserving()))

Returns a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`.
By default, `rotation_rate` is assumed to be Earth's.
"""
HydrostaticSphericalCoriolis(FT::DataType=Float64; rotation_rate=Ω_Earth, scheme::S=VectorInvariantEnergyConserving()) where S =
    HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme)

@inline φᶠᶠᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.φᵃᶠᵃ[j]
@inline φᶠᶠᵃ(i, j, k, grid::ConformalCubedSphereFaceGrid) = @inbounds grid.φᶠᶠᵃ[i, j]

@inline fᶠᶠᵃ(i, j, k, grid, coriolis::HydrostaticSphericalCoriolis) =
    2 * coriolis.rotation_rate * hack_sind(φᶠᶠᵃ(i, j, k, grid))

@inline z_f_cross_U(i, j, k, grid::AbstractGrid{FT}, coriolis::HydrostaticSphericalCoriolis, U) where FT = zero(FT)

#####
##### Enstrophy-conserving scheme
#####

const VIEnstrophy = HydrostaticSphericalCoriolis{<:VectorInvariantEnstrophyConserving}

@inline x_f_cross_U(i, j, k, grid, coriolis::VIEnstrophy, U) =
    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid)

@inline y_f_cross_U(i, j, k, grid, coriolis::VIEnstrophy, U) =
    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid)

#####
##### Energy-conserving scheme
#####

const VIEnergy = HydrostaticSphericalCoriolis{<:VectorInvariantEnergyConserving}

@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)

@inline x_f_cross_U(i, j, k, grid, coriolis::VIEnergy, U) =
    @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) / Δxᶠᶜᶜ(i, j, k, grid)

@inline y_f_cross_U(i, j, k, grid, coriolis::VIEnergy, U) =
    @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) / Δyᶜᶠᶜ(i, j, k, grid)

#####
##### Show
#####

Base.show(io::IO, hydrostatic_spherical_coriolis::HydrostaticSphericalCoriolis{FT}) where FT =
    print(io, "HydrostaticSphericalCoriolis{$FT}: rotation_rate = ",
          @sprintf("%.2e",  hydrostatic_spherical_coriolis.rotation_rate))


text:
"""
    checkPrecedence(visited::BitVector,priors::Vector{Int})

Return true if the precedence constraint is satisfied

# Arguments
- `visited::BitVector`: A list of visited nodes
- `priors::Vector{Int}`: A list of nodes that must have been visited
"""
function checkPrecedence(visited::BitVector,priors::Vector{Int})
    for i in priors
        if !visited[i]
            return false
        end
    end
    return true
end

"""
    getPrecedenceNumbers(precedenceF2P::Dict{T, Vector{T}})where{T<:Int}

Return a list of pre-calculated precedence numbers

# Arguments
- `precedenceF2P::Dict{T, Vector{T}}`: The map of followers to their priors
"""
function getPrecedenceNumbers(precedenceF2P::Dict{T, Vector{T}})where{T<:Int}
    precedenceNums = Dict{Int, Int}()
    for num in keys(precedenceF2P)
        if !haskey(precedenceNums,num)
            getPrecedenceNumber!(num, precedenceF2P,precedenceNums)
        end
    end
    return precedenceNums
end

"""
    getPrecedenceNumber!(checkNum::T, precedenceF2P::Dict{T, Vector{T}},precedenceNums::Dict{Int, Int})where{T<:Int}

Return a  precedence number, but also updates the values of any missing numbers it has to find along the way

# Arguments
- `checkNum::T`: The number to calculate
- `precedenceF2P::Dict{T, Vector{T}}`: The map of followers to their priors
- `precedenceNums::Dict{Int, Int}`: The list of calculated numbers
"""
function getPrecedenceNumber!(checkNum::T, precedenceF2P::Dict{T, Vector{T}},precedenceNums::Dict{Int, Int})where{T<:Int}
    if haskey(precedenceF2P,checkNum)
        maxVal = 0
        for element in precedenceF2P[checkNum]
            if !haskey(precedenceNums,element)
                precedenceNums[element] = getPrecedenceNumber!(element, precedenceF2P,precedenceNums)
            end
            if precedenceNums[element]>maxVal
                maxVal = precedenceNums[element]
            end
        end
        precedenceNums[checkNum] = max(maxVal+1,length(precedenceF2P[checkNum]))
    else
        precedenceNums[checkNum] = 0
    end
end


text:
module PhasePlaneAnalysis
	using MATLAB
	cd(joinpath(Pkg.dir("PhasePlaneAnalysis"),"src"))
	spawn(`jupyter-notebook PhasePlaneAnalysis.ipynb`)
	mat"addpath($(pwd()))"
end # module

max_stars_repo_name: YaoLuCNS/PhasePlaneAnalysis

max_issues_repo_name: YaoLuCNS/PhasePlaneAnalysis

max_forks_repo_name: YaoLuCNS/PhasePlaneAnalysis


text:
module Courant

export advective_courant, diffusive_courant, nondiffusive_courant

using LinearAlgebra, StaticArrays
using ..VariableTemplates
using ..MoistThermodynamics
using ..PlanetParameters
import ..MoistThermodynamics: internal_energy
using ..SubgridScaleParameters
using GPUifyLoops
using CLIMA.Atmos: AtmosModel, soundspeed, vertical_unit_vector, turbulence_tensors
import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state,
                        vars_diffusive
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Grids: VerticalDirection, HorizontalDirection, EveryDirection

function advective_courant(m::BalanceLaw, state::Vars, aux::Vars,
                              diffusive::Vars, Δx, Δt, direction=VerticalDirection())
    k̂ = vertical_unit_vector(m.orientation, aux)
    if direction isa VerticalDirection
        norm_u =  abs(dot(state.ρu, k̂))/state.ρ
    elseif direction isa HorizontalDirection
        norm_u = norm( (state.ρu .- dot(state.ρu, k̂).*k̂) / state.ρ )
    else
        norm_u = norm(state.ρu/state.ρ)
    end
    return Δt * norm_u / Δx
end

function nondiffusive_courant(m::BalanceLaw, state::Vars, aux::Vars,
                              diffusive::Vars, Δx, Δt, direction=VerticalDirection())
    k̂ = vertical_unit_vector(m.orientation, aux)
    if direction isa VerticalDirection
        norm_u =  abs(dot(state.ρu, k̂))/state.ρ
    elseif direction isa HorizontalDirection
        norm_u = norm( (state.ρu .- dot(state.ρu, k̂)*k̂) / state.ρ )
    else
        norm_u = norm(state.ρu/state.ρ)
    end
    return Δt * (norm_u + soundspeed(m.moisture, m.orientation, state, aux)) / Δx
end

function diffusive_courant(m::BalanceLaw, state::Vars, aux::Vars, diffusive::Vars, Δx, Δt, direction=VerticalDirection())
    ν, τ = turbulence_tensors(m.turbulence, state, diffusive, aux, 0)
    FT = eltype(state)

    if ν isa Real
        return Δt * ν / (Δx*Δx)
    else
        k̂ = vertical_unit_vector(m.orientation, aux)
        ν_vert = dot(ν, k)

        if direction isa VerticalDirection
            return Δt * ν_vert / (Δx*Δx)
        elseif direction isa HorizontalDirection
            ν_horz = ν - ν_vert .* k
            return Δt * norm(ν_horz) / (Δx*Δx)
        else
            return Δt * norm(ν) / (Δx*Δx)
        end
    end
end

end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
export shannon_entropy, von_neumann_entropy, joint_entropy, conditional_entropy
export holevo_bound, mutual_information
export success_probability, error_probability

"""
    shannon_entropy( probabilities :: AbstractVector ) :: Float64

The classical entropy of a probability distribution:

```math
S = -\\sum_{i=1}^n p_i \\log_2(p_i)
```

A `DomainError` is thrown if input `probabilities` does not satisfy [`is_probability_distribution`](@ref).
"""
function shannon_entropy( probabilities :: AbstractVector ) :: Float64
    is_probability_distribution(probabilities) || Probabilities(probabilities)

    entropy = -1*sum(map(
        (p) -> isapprox(p, 0, atol=1e-7) ? 0 : p*log2(p),
        probabilities
    ))

    entropy
end

"""
    von_neumann_entropy( ρ :: AbstractMatrix ) :: Float64

The von neumann entropy of a density matrix:

```math
S(\\rho) = - \\sum_j  \\lambda_j \\log_2(\\lambda_j)
```

where ``\\lambda_j`` are the eigenvalues of quantum  state ``\\rho``.

A `DomainError` is thrown if `ρ` does not satisfy [`is_density_matrix`](@ref).
"""
function von_neumann_entropy(ρ :: AbstractMatrix) :: Float64
    is_density_matrix(ρ) || State(ρ)

    λs = Probabilities(eigvals(ρ))
    entropy = shannon_entropy(λs)

    entropy
end

"""
    holevo_bound(
        priors :: AbstractVector,
        ρ_states :: Vector{<:AbstractMatrix}
    ) :: Float64

The Holevo theorem places a bound on the [`mutual_information`](@ref) ``I(X : Y) \\leq \\mathcal{X}``.
It places a limit on the amount information that can be decoded from a set of quantum states.
For a mixed state ``\\rho = \\sum_i p_i \\rho_i`` the Holevo bound ``\\mathcal{X}`` is

```math
\\mathcal{X}  := S(\\rho) - \\sum_i p_i S(\\rho_i)
```

where ``S(\\rho)`` is the [`von_neumann_entropy`](@ref).

A `DomainError` is thrown  if:
* `priors` does not satisfy [`is_probability_distribution`](@ref).
* Any `ρ ∈ ρ_states` does not satisfy [`is_density_matrix`](@ref).
"""
function holevo_bound(priors::AbstractVector, ρ_states::Vector{<:AbstractMatrix}) :: Float64
    is_probability_distribution(priors) || Probabilities(priors)
    all(is_density_matrix.(ρ_states)) || State.(ρ_states)

    ρ = mixed_state(priors, ρ_states)
    ρ_ent = von_neumann_entropy(ρ)

    ρ_states_ent = von_neumann_entropy.(ρ_states)
    ρ_sum_ent = sum( priors .* ρ_states_ent )

    bound = ρ_ent - ρ_sum_ent

    bound
end

"""
    joint_entropy(priors :: AbstractVector, conditionals :: AbstractMatrix) :: Float64

Returns the entropy for the union of pdf ``P(x,y)``. The joint entropy is the [`shannon_entropy`](@ref)
of the joint probability distribution:

```math
S(X,Y) = - \\sum_{x,y} p(x,y) \\log_2(p(x,y))
```
"""
function joint_entropy(priors::AbstractVector, conditionals::AbstractMatrix) :: Float64
    is_probability_distribution(priors) || Probabilities(priors)
    is_conditional_distribution(conditionals) || Conditionals(conditionals)

    joint_probabilities = JointProbabilities(priors, conditionals)
    shannon_entropy(joint_probabilities[:])
end

"""
    conditional_entropy(priors::AbstractVector, conditionals::AbstractMatrix) :: Float64

Returns the conditional entropy for the system with specified `priors` ``p(x)`` and `conditionals` ``p(y|x)``:

```math
S(Y|X) = - \\sum_{x,y} p(x,y)\\log_2\\left(\\frac{p(y|x)}{p(x)}\\right)
```
"""
function conditional_entropy(priors::AbstractVector, conditionals::AbstractMatrix) :: Float64
    is_probability_distribution(priors) || Probabilities(priors)
    is_conditional_distribution(conditionals) || Conditionals(conditionals)

    joint_probabilities = JointProbabilities(priors, conditionals)
    (num_row, num_col) = size(joint_probabilities)

    conditional_entropy = shannon_entropy(Probabilities(joint_probabilities[:])) + sum(sum(map( row -> map( col ->
        joint_probabilities[row,col]*( isapprox(priors[col], 0, atol=1e-7) ? 0 : log2(priors[col]) ),
    1:num_col), 1:num_row)))

    conditional_entropy
end

"""
    mutual_information(
        priors :: AbstractVector,
        conditionals :: AbstractMatrix
    ) :: Float64

The entropy of the overlap between ``p(x)`` and ``p(y)``. The mutual information
is directly computed from the [`shannon_entropy`](@ref) and [`joint_entropy`](@ref):

```math
I(X : Y) = S(X) + S(Y) - S(X,Y)
```
"""
function mutual_information(priors::AbstractVector, conditionals::AbstractMatrix) :: Float64
    is_probability_distribution(priors) || Probabilities(priors)
    is_conditional_distribution(conditionals) || Conditionals(conditional)

    prior_ent = shannon_entropy(priors)

    outcome_probs = outcome_probabilities(priors,conditionals)
    outcome_ent = shannon_entropy( outcome_probs )
    joint_ent = joint_entropy(priors, conditionals)

    mutual_information = prior_ent + outcome_ent - joint_ent

    mutual_information
end

"""
    mutual_information(
        priors :: AbstractVector,
        ρ_states :: Vector{<:AbstractMatrix},
        Π :: AbstractVector
    ) :: Float64

Computes the classical mutual information for a quantum state and measurement
encoding and decoding.
The conditional probabilities are obtained from quantum states and measurements
using [`measure`](@ref).

A `DomainError` is thrown if:
* `priors` does not satisfy [`is_probability_distribution`](@ref).
* Any `ρ ∈ ρ_states` does not satisfy [`is_density_matrix`](@ref).
* `Π` does not satisfy [`is_povm`](@ref).
"""
function mutual_information(
    priors::AbstractVector,
    states::Vector{<:AbstractMatrix},
    Π::AbstractVector{<:AbstractMatrix}
) :: Float64
    is_probability_distribution(priors) || Probabilities(priors)

    if !all(ρ -> ρ isa State, states)
        states = State.(states)
    end
    if !(Π isa POVM)
        Π = POVM(Π)
    end

    conditionals = measure(Π, states)
    mutual_information(priors, conditionals)
end

"""
    success_probability(
        priors::AbstractVector,
        states::Vector{<:AbstractMatrix},
        Π::Vector{<:AbstractMatrix}
    ) :: Float64

The probability of correctly distinguishing quantum states with the specifed POVM:

```math
P_{\\text{Success}} = \\sum_{i=1}^n p_i \\text{Tr}[\\Pi_i \\rho_i].
```

The number of states must match the number POVMs.
"""
function success_probability(priors::AbstractVector,states::Vector{<:AbstractMatrix}, Π::AbstractVector{<:AbstractMatrix}) :: Float64
    is_probability_distribution(priors) || Probabilities(priors)
    all(is_density_matrix.(states)) || State.(states)
    is_povm(Π) || POVM(Π)

    sum(priors .* tr.(states .* Π))
end

"""
    error_probability(
        priors::AbstractVector,
        ρ_states::Vector{<:AbstractMatrix},
        Π::AbstractVector{<:AbstractMatrix}
    ) :: Float64

The probability of incorrectly distinguishing quantum states with the specifed POVM.
This quantity is simply obtained as ``P_{\\text{Error}} = 1 - P_{\\text{Success}}``.
"""
function error_probability(priors::AbstractVector, ρ_states::Vector{<:AbstractMatrix}, Π::AbstractVector{<:AbstractMatrix}) :: Float64
    1 .- success_probability(priors, ρ_states, Π)
end


text:
using Poisson2D
using Test
using LinearAlgebra

@testset "Direct Solve" begin
    @test include("direct_test1.jl")
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

module Materials

abstract type AbstractMaterial end
abstract type AbstractMaterialState end

export AbstractMaterial, AbstractMaterialState
export integrate_material!, update_material!, reset_material!

"""
    :+(a::T, b::T) where T <: AbstractMaterialState

Fieldwise addition for material states.
"""
@generated function Base.:+(a::T, b::T) where T <: AbstractMaterialState
   expr = [:(a.$p + b.$p) for p in fieldnames(T)]
   return :(T($(expr...)))
end

"""
    integrate_material!(material::AbstractMaterial)

Integrate one timestep. The input `material.variables` represents the old
problem state.

Abstract method. Must be implemented for each material type. When integration is
done, the method **must** write the new state into `material.variables_new`.

**Do not** write into `material.variables`; actually committing the timestep
(i.e. accepting that one step of time evolution and applying it permanently)
is the job of `update_material!`.
"""
function integrate_material!(material::M) where M <: AbstractMaterial
    error("One needs to define how to integrate material $M!")
end

"""
    update_material!(material::AbstractMaterial)

Commit the result of `integrate_material!`.

In `material`, we add `ddrivers` into `drivers`, `dparameters` into
`parameters`, and replace `variables` by `variables_new`. Then we
automatically invoke `reset_material!`.
"""
function update_material!(material::AbstractMaterial)
    material.drivers += material.ddrivers
    material.parameters += material.dparameters
    material.variables = material.variables_new
    reset_material!(material)
    return nothing
end

"""
    reset_material!(material::AbstractMaterial)

In `material`, we zero out `ddrivers`, `dparameters` and `variables_new`. This
clears out the tentative state produced when a timestep has been computed, but
has not yet been committed.

Used internally by `update_material!`.
"""
function reset_material!(material::AbstractMaterial)
    material.ddrivers = typeof(material.ddrivers)()
    material.dparameters = typeof(material.dparameters)()
    material.variables_new = typeof(material.variables_new)()
    return nothing
end

include("utilities.jl")
using .Utilities
export Symm2, Symm4
export delta, II, IT, IS, IA, IV, ID, isotropic_elasticity_tensor, isotropic_compliance_tensor
export lame, delame, debang, find_root

include("perfectplastic.jl")
using .PerfectPlasticModule
export PerfectPlastic, PerfectPlasticDriverState, PerfectPlasticParameterState, PerfectPlasticVariableState

include("chaboche.jl")
using .ChabocheModule
export Chaboche, ChabocheDriverState, ChabocheParameterState, ChabocheVariableState

include("memory.jl")
using .MemoryModule
export Memory, MemoryDriverState, MemoryParameterState, MemoryVariableState

include("DSA.jl")
using .DSAModule
export DSA, DSADriverState, DSAParameterState, DSAVariableState

include("increments.jl")
using .Increments
export uniaxial_increment!, biaxial_increment!, stress_driven_uniaxial_increment!,
       general_increment!, stress_driven_general_increment!, find_dstrain!

end


text:
module DICOM

export dcm_parse, dcm_write, lookup, lookup_vr
export @tag_str

"""
    @tag_str(s)

Return the dicom tag, corresponding to the string `s`.
```jldoctest
julia> using DICOM

julia> tag"ROI Mean"
(0x6000, 0x1302)
```
"""
macro tag_str(s)
    DICOM.fieldname_dict[s]
end

include("dcm_dict.jl")

# Create dicom dictionary - used for reading/writing DICOM files
# Keys are tuple containing hex Group and Element of DICOM entry
# Stored value is String array: [Field Name, VR, Data Length], e.g:
# Julia> DICOM.dcm_dict[(0x0008,0x0005)]
# 3-element Array{String,1}:
#  "Specific Character Set"
#  "CS"
#  "1-n"
function dcm_init()
    dcm_dict = Dict{Tuple{UInt16,UInt16},Array{String,1}}()
    for d in (_dcmdict_data_::Array{Any,1})
        dcm_dict[(UInt16(d[1][1]),UInt16(d[1][2]))] = d[2:end]
    end
    return(dcm_dict)
end

# For convenience, dictionary to get hex tag from field name, e.g:
# Julia> DICOM.fieldname_dict["Specific Character Set"]
# (0x0008, 0x0005)
function fieldname_init()
    fieldname_dict = Dict{AbstractString, Tuple{UInt16,UInt16}}()
    for d in (_dcmdict_data_::Array{Any,1})
        fieldname_dict[d[2]] = (UInt16(d[1][1]),UInt16(d[1][2]))
    end
    return(fieldname_dict)
end

const dcm_dict = dcm_init()
const fieldname_dict = fieldname_init()
_dcmdict_data_ = 0

# These "empty" values are used internally. They are returned if a search fails.
const emptyVR = "" # Can be any VR that doesn't exist
const emptyVR_lookup = ["", emptyVR, ""] # Used in lookup_vr as failure state
const emptyTag = (0x0000,0x0000)
const emptyDcmDict = Dict(DICOM.emptyTag => nothing)

const VR_names = [ "AE","AS","AT","CS","DA","DS","DT","FL","FD","IS","LO","LT","OB","OF",
       "OW","PN","SH","SL","SQ","SS","ST","TM","UI","UL","UN","US","UT" ]

# mapping UID => bigendian? explicitvr?
const meta_uids = Dict([("1.2.840.10008.1.2", (false, false)),
                  ("1.2.840.10008.1.2.1", (false, true)),
                  ("1.2.840.10008.1.2.1.99", (false, true)),
                  ("1.2.840.10008.1.2.2", (true, true))]);

"""
    lookup_vr(tag::Tuple{Integer,Integer})

Return VR value for tag from DICOM dictionary

# Example
```jldoctest
julia> lookup_vr((0x0028,0x0004))
"CS"
```
"""
function lookup_vr(gelt::Tuple{UInt16,UInt16})
    if gelt[1]&0xff00 == 0x5000
        gelt = (0x5000,gelt[2])
    elseif gelt[1]&0xff00 == 0x6000
        gelt = (0x6000,gelt[2])
    end
    r = get(dcm_dict, gelt, emptyVR_lookup)
    return(r[2])
end

function lookup(d::Dict{Tuple{UInt16,UInt16},Any}, fieldnameString::String)
    return(get(d, fieldname_dict[fieldnameString], nothing))
end

always_implicit(grp, elt) = (grp == 0xFFFE && (elt == 0xE0DD||elt == 0xE000||
                                               elt == 0xE00D))

dcm_store(st::IOStream, gelt::Tuple{UInt16,UInt16}, writef::Function) = dcm_store(st, gelt, writef, emptyVR)
function dcm_store(st::IOStream, gelt::Tuple{UInt16,UInt16}, writef::Function, vr::String)
    lentype = UInt32
    write(st, UInt16(gelt[1])) # Grp
    write(st, UInt16(gelt[2])) # Elt
    if vr !== emptyVR
        write(st, vr)
        if vr in ("OB", "OW", "OF", "SQ", "UT", "UN")
            write(st, UInt16(0))
        else
            lentype = UInt16
        end
    end
    # Write data first, then calculate length, then go back to write length
    p = position(st)
    write(st, zero(lentype)) # Placeholder for the data length
    writef(st)
    endp = position(st)
    # Remove placeholder's length, either 2 (UInt16) or 4 (UInt32) steps (1 step = 8 bits)
    if lentype == UInt32
        sz = endp-p-4
    else
        sz = endp-p-2
    end
    szWasOdd = isodd(sz) # If length is odd, round up - UInt8(0) will be written at end
    if szWasOdd
        sz+=1
    end
    seek(st, p)
    write(st, convert(lentype, max(0,sz)))
    seek(st, endp)
    if szWasOdd
        write(st, UInt8(0))
    end
end

function undefined_length(st, vr)
    data = IOBuffer()
    w1 = w2 = 0
    while true
        # read until 0xFFFE 0xE0DD
        w1 = w2
        w2 = read(st, UInt16)
        if w1 == 0xFFFE
            if w2 == 0xE0DD
                break
            end
            write(data, w1)
        end
        if w2 != 0xFFFE
            write(data, w2)
        end
    end
    skip(st, 4)
    take!(data)
end

function sequence_item(st::IOStream, evr, sz)
    item = Dict{Tuple{UInt16,UInt16},Any}()
    endpos = position(st) + sz
    while position(st) < endpos
        (gelt, data, vr) = element(st, evr)
        if isequal(gelt, (0xFFFE,0xE00D))
            break
        end
        item[gelt] = data
    end
    return item
end

function sequence_item_write(st::IOStream, evr::Bool, items::Dict{Tuple{UInt16,UInt16},Any})
    for gelt in sort(collect(keys(items)))
        element_write(st, evr, gelt, items[gelt])
    end
    write(st, UInt16[0xFFFE, 0xE00D, 0x0000, 0x0000])
end

function sequence_parse(st, evr, sz)
    sq = Array{Dict{Tuple{UInt16,UInt16},Any},1}()
    while sz > 0
        grp = read(st, UInt16)
        elt = read(st, UInt16)
        itemlen = read(st, UInt32)
        if grp==0xFFFE && elt==0xE0DD
            return sq
        end
        if grp != 0xFFFE || elt != 0xE000
            error("dicom: expected item tag in sequence")
        end
        push!(sq, sequence_item(st, evr, itemlen))
        sz -= 8 + (itemlen != 0xffffffff) * itemlen
    end
    return sq
end

function sequence_write(st::IOStream, evr::Bool, items::Array{Dict{Tuple{UInt16,UInt16},Any},1})
    for subitem in items
        if length(subitem) > 0
            dcm_store(st, (0xFFFE,0xE000), s->sequence_item_write(s, evr, subitem))
        end
    end
    write(st, UInt16[0xFFFE, 0xE0DD, 0x0000, 0x0000])
end

# always little-endian, "encapsulated" iff sz==0xffffffff
function pixeldata_parse(st::IOStream, sz, vr::String, dcm=emptyDcmDict)
    # (0x0028,0x0103) defines Pixel Representation
    isSigned = false
    f = get(dcm, (0x0028,0x0103), nothing)
    if f !== nothing
        # Data is signed if f==1
        isSigned = f == 1
    end
    # (0x0028,0x0100) defines Bits Allocated
    bitType = 16
    f = get(dcm, (0x0028,0x0100), nothing)
    if f !== nothing
        bitType = Int(f)
    else
        f = get(dcm, (0x0028,0x0101), nothing)
        bitType = f !== nothing ? Int(f) :
            vr == "OB" ? 8 : 16
    end
    if bitType == 8
        dtype = isSigned ? Int8 : UInt8
    else
        dtype = isSigned ? Int16 : UInt16
    end

    yr=1
    zr=1
    # (0028,0010) defines number of rows
    f = get(dcm, (0x0028,0x0010), nothing)
    if f !== nothing
        xr = Int(f)
    end
    # (0028,0011) defines number of columns
    f = get(dcm, (0x0028,0x0011), nothing)
    if f !== nothing
        yr = Int(f)
    end
    # (0028,0012) defines number of planes
    f = get(dcm, (0x0028,0x0012), nothing)
    if f !== nothing
        zr = Int(f)
    end
    # (0028,0008) defines number of frames
    f = get(dcm, (0x0028,0x0008), nothing)
    if f !== nothing
        zr *= Int(f)
    end
    if sz != 0xffffffff
        data =
        zr > 1 ? Array{dtype}(undef, xr, yr, zr) : Array{dtype}(undef, xr, yr)
        read!(st, data)
    else
        # start with Basic Offset Table Item
        data = Array{Any,1}(element(st, false)[2])
        while true
            grp = read(st, UInt16)
            elt = read(st, UInt16)
            xr = read(st, UInt32)
            if grp == 0xFFFE && elt == 0xE0DD
                return data
            end
            if grp != 0xFFFE || elt != 0xE000
                error("dicom: expected item tag in encapsulated pixel data")
            end
            if dtype === UInt16; xr = div(xr,2); end
            push!(data, read!(st, Array{dtype}(undef, xr)))
        end
    end
    return data
end

function pixeldata_write(st, evr, d)
    # if length(el) > 1
    #     error("dicom: compression not supported")
    # end
    nt = eltype(d)
    vr = nt === UInt8  || nt === Int8  ? "OB" :
         nt === UInt16 || nt === Int16 ? "OW" :
         nt === Float32                ? "OF" :
         error("dicom: unsupported pixel format")
    if evr !== false
        dcm_store(st, (0x7FE0,0x0010), s->write(s,d), vr)
    elseif vr != "OW"
        error("dicom: implicit VR only supports 16-bit pixels")
    else
        dcm_store(st, (0x7FE0,0x0010), s->write(s,d))
    end
end

function skip_spaces(st, endpos)
    while true
        c = read(st,Char)
        if c != ' ' || position(st) == endpos
            return c
        end
    end
end

function string_parse(st, sz, maxlen, spaces)
    endpos = position(st)+sz
    data = [ "" ]
    first = true
    while position(st) < endpos
        c = !first||spaces ? read(st,Char) : skip_spaces(st, endpos)
        if c == '\\'
            push!(data, "")
            first = true
        elseif c == '\0'
            break
        else
            data[end] = string(data[end],c)  # TODO: inefficient
            first = false
        end
    end
    if !spaces
        return map(rstrip,data)
    end
    return data
end

numeric_parse(st::IOStream, T::DataType, sz) = T[read(st, T) for i=1:div(sz,sizeof(T))]

function element(st::IOStream, evr::Bool, dcm=emptyDcmDict, dVR=Dict{Tuple{UInt16,UInt16},String}())
    lentype = UInt32
    diffvr = false
    local grp
    try
        grp = read(st, UInt16)
    catch
        return(emptyTag,0,emptyVR)
    end
    elt = read(st, UInt16)
    gelt = (grp,elt)
    if grp <= 0x0002
        evr = true
    end
    if evr && !always_implicit(grp,elt)
        vr = String(read!(st, Array{UInt8}(undef, 2)))
        if vr in ("OB", "OW", "OF", "SQ", "UT", "UN")
            skip(st, 2)
        else
            lentype = UInt16
        end
        diffvr = !isequal(vr, lookup_vr(gelt))
    else
        vr = elt == 0x0000 ? "UL" : lookup_vr(gelt)
    end
    if isodd(grp) && grp > 0x0008 && 0x0010 <= elt <+ 0x00FF
        # Private creator
        vr = "LO"
    elseif isodd(grp) && grp > 0x0008
        # Assume private
        vr = "UN"
    end
    if haskey(dVR, gelt)
        vr = dVR[gelt]
    end
    if vr === emptyVR
        if haskey(dVR, (0x0000,0x0000))
            vr = dVR[(0x0000,0x0000)]
        elseif !haskey(dVR, gelt)
            error("dicom: unknown tag ", gelt)
        end
    end

    sz = read(st,lentype)

    # Empty VR can be supplied in dVR to skip an element
    if vr == ""
        sz = isodd(sz) ? sz+1 : sz
        skip(st,sz)
        return(element(st,evr,dcm,dVR))
    end

    data =
    vr=="ST" || vr=="LT" || vr=="UT" ? String(read!(st, Array{UInt8}(undef, sz))) :

    sz==0 || vr=="XX" ? Any[] :

    vr == "SQ" ? sequence_parse(st, evr, sz) :

    gelt == (0x7FE0,0x0010) ? pixeldata_parse(st, sz, vr, dcm) :

    sz == 0xffffffff ? undefined_length(st, vr) :

    vr == "FL" ? numeric_parse(st, Float32, sz) :
    vr == "FD" ? numeric_parse(st, Float64, sz) :
    vr == "SL" ? numeric_parse(st, Int32  , sz) :
    vr == "SS" ? numeric_parse(st, Int16  , sz) :
    vr == "UL" ? numeric_parse(st, UInt32 , sz) :
    vr == "US" ? numeric_parse(st, UInt16 , sz) :

    vr == "OB" ? read!(st, Array{UInt8}(undef, sz))        :
    vr == "OF" ? read!(st, Array{Float32}(undef, div(sz,4))) :
    vr == "OW" ? read!(st, Array{UInt16}(undef, div(sz,2))) :

    vr == "AT" ? [ read!(st, Array{UInt16}(undef, 2)) for n=1:div(sz,4) ] :

    vr == "AS" ? String(read!(st, Array{UInt8}(undef, 4))) :

    vr == "DS" ? map(x->parse(Float64,x), string_parse(st, sz, 16, false)) :
    vr == "IS" ? map(x->parse(Int,x), string_parse(st, sz, 12, false)) :

    vr == "AE" ? string_parse(st, sz, 16, false) :
    vr == "CS" ? string_parse(st, sz, 16, false) :
    vr == "SH" ? string_parse(st, sz, 16, false) :
    vr == "LO" ? string_parse(st, sz, 64, false) :
    vr == "UI" ? string_parse(st, sz, 64, false) :
    vr == "PN" ? string_parse(st, sz, 64, true)  :

    vr == "DA" ? string_parse(st, sz, 10, true) :
    vr == "DT" ? string_parse(st, sz, 26, false) :
    vr == "TM" ? string_parse(st, sz, 16, false) :
    read!(st, Array{UInt8}(undef, sz))

    if isodd(sz) && sz != 0xffffffff
        skip(st, 1)
    end

    # For convenience, get rid of array if it is just acting as a container
    # Exception is "SQ", where array is part of structure
    if length(data) == 1 && vr != "SQ"
        data = data[1]
        if length(data) == 1
            data = data[1]
        end
    end

    # Return vr by default
    return(gelt, data, vr)
end

# todo: support maxlen
string_write(vals::Array{SubString{String}}, maxlen) = string_write(convert(Array{String}, vals), maxlen)
string_write(vals::SubString{String}, maxlen) = string_write(convert(String, vals), maxlen)
string_write(vals::Char, maxlen) = string_write(string(vals), maxlen)
string_write(vals::String, maxlen) = string_write([vals], maxlen)
string_write(vals::Array{String,1}, maxlen) = join(vals, '\\')

element_write(st::IOStream, evr::Bool, gelt::Tuple{UInt16,UInt16}, data::Any) = element_write(st,evr,gelt,data,lookup_vr(gelt))
function element_write(st::IOStream, evr::Bool, gelt::Tuple{UInt16,UInt16}, data::Any, vr::String)
    if vr === emptyVR
        # Element tags ending in 0x0000 are not included in dcm_dicm.jl, their vr is UL
        if gelt[2] == 0x0000
            vr = "UL"
        elseif isodd(gelt[1]) && gelt[1] > 0x0008 && 0x0010 <= gelt[2] <+ 0x00FF
                # Private creator
                vr = "LO"
        elseif isodd(gelt[1]) && gelt[1] > 0x0008
                # Assume private
                vr = "UN"
        else
            error("dicom: unknown tag ", gelt)
        end
    end
    if gelt == (0x7FE0, 0x0010)
        return pixeldata_write(st, evr, data)
    end

    if vr == "SQ"
        vr = evr ? vr : emptyVR
        return dcm_store(st, gelt,
                         s->sequence_write(s, evr, data), vr)
    end

    # Pack data into array container. This is to undo "data = data[1]" from element().
    if !isa(data,Array) && vr in ("FL","FD","SL","SS","UL","US")
        data = [data]
    end

    data =
    isempty(data) ? UInt8[] :
    vr in ("OB","OF","OW","ST","LT","UT") ? data :
    vr in ("AE", "CS", "SH", "LO", "UI", "PN", "DA", "DT", "TM") ?
        string_write(data, 0) :
    vr == "FL" ? convert(Array{Float32,1}, data) :
    vr == "FD" ? convert(Array{Float64,1}, data) :
    vr == "SL" ? convert(Array{Int32,1},   data) :
    vr == "SS" ? convert(Array{Int16,1},   data) :
    vr == "UL" ? convert(Array{UInt32,1},  data) :
    vr == "US" ? convert(Array{UInt16,1},  data) :
    vr == "AT" ? [data...] :
    vr in ("DS","IS") ? string_write(map(string,data), 0) :
    data

    if evr === false && gelt[1]>0x0002
        vr = emptyVR
    end

    dcm_store(st, gelt, s->write(s, data), vr)
end

"""
    dcm_parse(fn::AbstractString)

Reads file fn and returns a Dict
"""
function dcm_parse(fn::AbstractString, giveVR=false; header=true, maxGrp=0xffff, dVR=Dict{Tuple{UInt16,UInt16},String}())
    st = open(fn)
    if header
        # First 128 bytes are preamble - should be skipped
        skip(st, 128)
        # "DICM" identifier must be after preamble
        sig = String(read!(st,Array{UInt8}(undef, 4)))
        if sig != "DICM"
            error("dicom: invalid file header")
            # seek(st, 0)
        end
    end
    # a bit of a hack to detect explicit VR. seek past the first tag,
    # and check to see if a valid VR name is there
    skip(st, 4)
    sig = String(read!(st,Array{UInt8}(undef, 2)))
    evr = sig in VR_names
    skip(st, -6)
    dcm = Dict{Tuple{UInt16,UInt16},Any}()
    if giveVR
        dcmVR = Dict{Tuple{UInt16,UInt16},String}()
    end

    while true
        (gelt, data, vr) = element(st, evr, dcm, dVR)
        if gelt === emptyTag || gelt[1] > maxGrp
            break
        else
            dcm[gelt] = data
            if giveVR
                dcmVR[gelt] = vr
            end
        end
        # look for transfer syntax UID
        if gelt == (0x0002,0x0010)
            # Default is endian=little, explicitVR=true
            metaInfo = get(meta_uids, data, (false, true))
            evr = metaInfo[2]
            if metaInfo[1]
                # todo: set byte order to big
            else
                # todo: set byte order to little
            end
        end
    end
    close(st)
    if giveVR
        return(dcm,dcmVR)
    else
        return(dcm)
    end
end

dcm_write(fn::String, d::Dict{Tuple{UInt16,UInt16},Any}, dVR=Dict{Tuple{UInt16,UInt16},String}()) = dcm_write(open(fn,"w+"),d,dVR)
function dcm_write(st::IOStream, d::Dict{Tuple{UInt16,UInt16},Any}, dVR=Dict{Tuple{UInt16,UInt16},String}())
    write(st, zeros(UInt8, 128))
    write(st, "DICM")
    # If no dictionary containing VRs is provided, then assume implicit VR - at first
    evr = !isempty(dVR)
    if !haskey(d,(0x0002,0x0010))
        # Insert UID for our transfer syntax, if it doesn't exist
        if evr
            element_write(st, evr, (0x0002,0x0010), "1.2.840.10008.1.2.1", "UI")
        else
            element_write(st, evr, (0x0002,0x0010), "1.2.840.10008.1.2", "UI")
        end
    else
        # Otherwise, use existing transfer UID, and overwrite evr accordingly
        metaInfo = get(meta_uids, d[(0x0002,0x0010)], (false, true))
        evr = metaInfo[2]
    end
    # dVR is only used if it isn't empty and evr=true
    if evr && !isempty(dVR)
        for gelt in sort(collect(keys(d)))
            # dVR only needs to contain keys for cases where lookup_vr() fails
            haskey(dVR, gelt) ? element_write(st, evr, gelt, d[gelt], dVR[gelt]) :
                element_write(st, evr, gelt, d[gelt])
        end
    else
        for gelt in sort(collect(keys(d)))
            element_write(st, evr, gelt, d[gelt])
        end
    end
    close(st)
end

end


text:
using MicMods
using Test

@testset "MicMods.jl" begin
    # Write your tests here.
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule racket_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("racket")
JLLWrappers.@generate_main_file("racket", UUID("d1d5beed-0862-5c4e-adcd-9d7a31bb6a79"))
end  # module racket_jll


text:
# Utility to call the cache's organized by symbol
function fill_array_by_symbol_dispatch(f, c, symbols, args...)
    for (i, sym) in enumerate(symbols)
        f(c[i], Val(sym), args...)
    end
end

# The generate_perturbation function calculates the perturbation itself
# It can do used without any derivatives overhead (except, perhaps, extra memory in the cache)
function generate_perturbation(m::PerturbationModel, p_d, p_f, order::Val{1} = Val(1);
                               cache = SolverCache(m, Val(1), p_d),
                               settings = PerturbationSolverSettings())
    @assert cache.p_d_symbols == collect(Symbol.(keys(p_d)))

    p = isnothing(p_f) ? p_d : order_vector_by_symbols(merge(p_d, p_f), m.mod.p_symbols)

    # solver type provided to all callbacks
    ret = calculate_steady_state!(m, cache, settings, p)
    maybe_call_function(settings.calculate_steady_state_callback, ret, m, cache, settings,
                        p)  # before returning
    (ret == :Success) || return FirstOrderPerturbationSolution(ret, m, cache)

    ret = evaluate_first_order_functions!(m, cache, settings, p)
    maybe_call_function(settings.evaluate_functions_callback, ret, m, cache, settings, p)
    (ret == :Success) || return FirstOrderPerturbationSolution(ret, m, cache)

    ret = solve_first_order!(m, cache, settings)
    maybe_call_function(settings.solve_first_order_callback, ret, m, cache, settings)
    (ret == :Success) || return FirstOrderPerturbationSolution(ret, m, cache)
    return FirstOrderPerturbationSolution(:Success, m, cache)
end

# The generate_perturbation function calculates the perturbation itself
# It can do used without any derivatives overhead (except, perhaps, extra memory in the cache)
function generate_perturbation(m::PerturbationModel, p_d, p_f, order::Val{2};
                               cache = SolverCache(m, Val(2), p_d),
                               settings = PerturbationSolverSettings())
    @assert cache.p_d_symbols == collect(Symbol.(keys(p_d)))
    @assert cache.order == Val(2)

    p = isnothing(p_f) ? p_d : order_vector_by_symbols(merge(p_d, p_f), m.mod.p_symbols)

    # Calculate the first-order perturbation
    sol_first = generate_perturbation(m, p_d, p_f, Val(1); cache, settings)

    (sol_first.retcode == :Success) || return SecondOrderPerturbationSolution(sol_first.retcode, m, cache)

    # solver type provided to all callbacks
    ret = evaluate_second_order_functions!(m, cache, settings, p)
    (ret == :Success) || return SecondOrderPerturbationSolution(ret, m, cache)
    maybe_call_function(settings.evaluate_functions_callback, ret, m, cache, settings, p)
    ret = solve_second_order!(m, cache, settings)
    (ret == :Success) || return SecondOrderPerturbationSolution(ret, m, cache)
    return SecondOrderPerturbationSolution(:Success, m, cache)
end

function calculate_steady_state!(m::PerturbationModel, c, settings, p)
    @unpack n_y, n_x = m
    n = n_y + n_x

    (settings.print_level > 2) && println("Calculating steady state")
    try
        if !isnothing(m.mod.ȳ!) && !isnothing(m.mod.x̄!) # use closed form if possible
            m.mod.ȳ!(c.y, p)
            m.mod.x̄!(c.x, p)
            isnothing(m.mod.ȳ_p!) ||
                fill_array_by_symbol_dispatch(m.mod.ȳ_p!, c.y_p, c.p_d_symbols, p)
            isnothing(m.mod.x̄_p!) ||
                fill_array_by_symbol_dispatch(m.mod.x̄_p!, c.x_p, c.p_d_symbols, p)
        elseif !isnothing(m.mod.steady_state!) # use user-provided calculation otherwise
            m.mod.steady_state!(c.y, c.x, p)
        else # fallback is to solve system of equations from user-provided initial condition
            y_0 = zeros(n_y)
            x_0 = zeros(n_x)
            m.mod.ȳ_iv!(y_0, p)
            m.mod.x̄_iv!(x_0, p)
            w_0 = [y_0; x_0]

            if isnothing(m.mod.H̄_w!) # no jacobian
                nlsol = nlsolve((H, w) -> m.mod.H̄!(H, w, p), w_0;
                                DifferentiableStateSpaceModels.nlsolve_options(settings)...)
            else
                J_0 = zeros(n, n)
                F_0 = zeros(n)
                df = OnceDifferentiable((H, w) -> m.mod.H̄!(H, w, p),
                                        (J, w) -> m.mod.H̄_w!(J, w, p), w_0, F_0, J_0)  # TODO: the buffer to use for the w_0 is unclear?
                nlsol = nlsolve(df, w_0;
                                DifferentiableStateSpaceModels.nlsolve_options(settings)...)
            end
            if !converged(nlsol)
                if settings.print_level > 0
                    println("No steady state found\n")
                end
                return :SteadyStateFailure
            end
            settings.print_level > 1 &&
                println("Steady state found in $(nlsol.iterations) iterations\n")
            c.y .= nlsol.zero[1:n_y]
            c.x .= nlsol.zero[(n_y + 1):end]
        end
    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success
end

function evaluate_first_order_functions!(m, c, settings, p)
    (settings.print_level > 2) && println("Evaluating first-order functions into cache")
    try
        @unpack y, x = c  # Precondition: valid (y, x) steady states

        m.mod.H_yp!(c.H_yp, y, x, p)
        m.mod.H_y!(c.H_y, y, x, p)
        m.mod.H_xp!(c.H_xp, y, x, p)
        m.mod.H_x!(c.H_x, y, x, p)
        m.mod.Γ!(c.Γ, p)
        maybe_call_function(m.mod.Ω!, c.Ω, p) # supports  m.mod.Ω! = nothing
        (length(c.p_d_symbols) > 0) && m.mod.Ψ!(c.Ψ, y, x, p)
    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success  # no failing code-paths yet.
end

function evaluate_second_order_functions!(m, c, settings, p)
    (settings.print_level > 2) && println("Evaluating second-order functions into cache")
    try
        @unpack y, x = c  # Precondition: valid (y, x) steady states
        (length(c.p_d_symbols) == 0) && m.mod.Ψ!(c.Ψ, y, x, p)  # would have been called otherwise in first_order_functions
        m.mod.Ψ!(c.Ψ, y, x, p)
        m.mod.Ψ_yp!(c.Ψ_yp, y, x, p)
        m.mod.Ψ_y!(c.Ψ_y, y, x, p)
        m.mod.Ψ_xp!(c.Ψ_xp, y, x, p)
        m.mod.Ψ_x!(c.Ψ_x, y, x, p)
    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success  # no failing code-paths yet.
end

function solve_first_order!(m, c, settings)
    @unpack ϵ_BK, print_level = settings
    @unpack n_x, n_y, n_p, n_ϵ = m
    n = n_x + n_y

    (settings.print_level > 2) && println("Solving first order perturbation")
    try
        buff = c.first_order_solver_buffer
        buff.A .= [c.H_xp c.H_yp]
        buff.B .= [c.H_x c.H_y]
        s = schur!(buff.A, buff.B) # Generalized Schur decomposition, inplace using buffers
        # The generalized eigenvalues λ_i are S_ii / T_ii
        # Following Blanchard-Kahn condition, we reorder the Schur so that
        # S_22 ./ T_22 < 1, ie, the eigenvalues < 1 come last
        # inds = [s.α[i] / s.β[i] >= 1 for i in 1:n]
        inds = abs.(s.α) .>= (1 - ϵ_BK) .* abs.(s.β)
        if sum(inds) != n_x
            if print_level > 0
                # @show (n_x, sum(inds), inds, abs.(s.α), abs.(s.β)) # move to print_level > 1?
                println("Blanchard-Kahn condition not satisfied\n")
            end
            return :BlanchardKahnFailure
        end

        ordschur!(s, inds)
        # In Julia A = QSZ' and B = QTZ'

        b = 1:n_x
        l = (n_x + 1):n

        # Extract the Schur components to real matrices, for inplace factorizations/etc.
        Z = s.Z'

        c.g_x .= real(-Z[l, l] \ Z[l, b])
        blob = Z[b, b] .+ Z[b, l] * c.g_x
        c.h_x .= real(-blob \ (s.S[b, b] \ (s.T[b, b] * blob)))
        # buff.Z .= real(s.Z')
        # buff.Z_ll .= buff.Z[l,l] #preallocate for buffer for inplace LU.  No known structure?

        # Both of these are upper-triangular, helpful for fast linsolve
        # buff.S_bb .= UpperTriangular(real(s.S[b,b]))
        # buff.T_bb .= UpperTriangular(real(s.T[b,b]))

        # TODO: Check if RecursiveFactorization.jl is faster or slower than using MKL/etc. Add toggle
        #Z_ll = lu!(buff.Z_ll)
        # Z_ll = RecursiveFactorization.lu!(buff.Z_ll)
        # c.g_x .= ldiv!(Z_ll, buff.Z[l, b])
        # c.g_x .*= -1

        # The following is an as-inplace-as-possible version of
        # blob = buff.Z[b, b] .+ buff.Z[b, l] * c.g_x
        # c.h_x .= -blob \ (buff.S_bb \ (buff.T_bb * blob))
        # temp = buff.Z[b, b] .+ buff.Z[b, l] * c.g_x
        # mul!(c.h_x, buff.T_bb, temp)  # doing everything in place
        # ldiv!(buff.S_bb, c.h_x)  # no factorization required since triangular
        # temp_lu = lu!(temp)
        # ldiv!(temp_lu, c.h_x)
        # c.h_x .*= -1

        # fill in Σ, Ω.
        c.Σ .= Symmetric(c.Γ * c.Γ')

        # Q transforms
        mul!(c.C_1, c.Q, vcat(c.g_x, I))

        # Stationary Distribution
        c.η_Σ_sq .= Symmetric(c.η * c.Σ * c.η')  # used in 2nd order as well

        V = cholesky(lyapd(c.h_x, c.η_Σ_sq); check = false) #inplace wouldn't help since allocating for lyapd.  Can use lyapds! perhaps, but would need work and have low payoffs

        # no inplace assignment or copy for cholesky, so reach inside internals for now.  Assumes same uplo flag
        c.V.factors .= V.factors

        # eta * Gamma
        mul!(c.B, c.η, c.Γ)

        # @exfiltrate  # flip on to see intermediate calculations.  TURN OFF BEFORE PROFILING
    catch e
        if !is_linear_algebra_exception(e)
            (settings.print_level > 2) && println("Rethrowing exception")
            rethrow(e)
        else
            settings.print_level == 0 || display(e)
            return :Failure # generic failure
        end
    end
    return :Success
end

function solve_second_order!(m, c, settings)
    @unpack n_x, n_y, n_p, n_ϵ, n_z, η = m
    n = n_x + n_y

    buff = c.second_order_solver_buffer

    try
        # "Sylvester prep for _xx"
        buff.A .= [c.H_y c.H_xp + c.H_yp * c.g_x]
        buff.C[:, 1:n_y] .= c.H_yp
        ws = GeneralizedSylvesterWs(n, n, n_x, 2)
        # TODO: Tullio/etc. quadratic form trickier for any of this?
        buff.R .= vcat(c.g_x * c.h_x, c.g_x, c.h_x, c.I_x)
        for i in 1:n
            buff.E[i, :] = vec(buff.R' * c.Ψ[i] * buff.R) # (24), flip the sign for gsylv
        end
        buff.E .*= -1

        # Sylvester
        generalized_sylvester_solver!(buff.A, buff.C, c.h_x, buff.E, 2, ws)
        c.g_xx .= reshape(buff.E[1:n_y, :], n_y, n_x, n_x)
        c.h_xx .= reshape(buff.E[(n_y + 1):end, :], n_x, n_x, n_x)

        # Linear equations for _σσ
        buff.A_σ .= [c.H_yp + c.H_y c.H_xp + c.H_yp * c.g_x]
        C_σ = zeros(n)
        H_yp_g = c.H_yp * buff.E[1:n_y, :]
        buff.R_σ .= vcat(c.g_x, zeros(n_y, n_x), c.I_x, zeros(n_x, n_x))
        for i in 1:n # (29), flip the sign for (34)
            C_σ[i] -= dot(buff.R_σ' * c.Ψ[i] * buff.R_σ, c.η_Σ_sq)
            C_σ[i] -= dot(H_yp_g[i, :], c.η_Σ_sq)
        end
        #A_σ_lu = lu!(buff.A_σ) # modifes the buff.A_σ
        A_σ_lu = RecursiveFactorization.lu!(buff.A_σ)
        ldiv!(A_σ_lu, C_σ) # solve (34) inplace. X_σ = C_σ after modification
        c.g_σσ .= C_σ[1:n_y]
        c.h_σσ .= C_σ[(n_y + 1):end]

        c.C_0 .= 0.5 * c.Q * vcat(c.g_σσ, zeros(n_x))

        # TODO: This looks like a tullio thing
        fill!(c.C_2, zero(eltype(c.C_2)))  # reset as we need to use `+=`
        for i in 1:n_z
            for j in 1:n_y
                c.C_2[i, :, :] += 0.5 * c.Q[i, j] * c.g_xx[j, :, :]
            end
        end
    catch e
        settings.print_level == 0 || display(e)
        return :Failure # generic failure
    end
    #    @exfiltrate  # flip on to see intermediate calculations.  TURN OFF BEFORE PROFILING
    return :Success
end