Third-order high-temperature Birch-Murnaghan equation of state

#pragma version=1.3

// set bracketing values and tolerance for volume
// absolute brackets will be (1 +/- kBracketFactor) * V0
// should be okay for EOS for solids
static constant kBracketFactor=0.1
static constant kTol=1e-5

// V = volume at P, T
// Tref = reference temperature
// V0 = volume at Tref
// thermal expansion coefficient = c1 + c2*T + c3/T + c4/T^2 + c5/T^0.5
// K0 is the bulk modulus at Tref
// Kprime is the pressure deriviative of bulk modulus
// dKdT is temperature derivative of bulk modulus

// For room temperature EOS, set a, b, c, and dKdT to zero.

// 3rd order high-temperature Birch-Murnaghan equation of state, returns P
threadsafe function HTBM(V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT
   
    variable KT=K0+dKdT*(T-Tref)
    variable V0T=V0*exp(c1*(T-Tref)+0.5*c2*(T^2-Tref^2)+c3*ln(T/Tref)-c4*(T^-1-Tref^-1)+2*c5*(T^0.5-Tref^0.5))
    return 3/2*KT*((V0T/V)^(7/3) - (V0T/V)^(5/3))*(1 + 3/4*(Kprime-4)*((V0T/V)^(2/3) - 1))
end

// wrapper function for FindRoots
threadsafe function HTBM_wrapper(wave w, variable V)
    // w contains V0, T, Tref, a, b, c, dKdT, K0, Kprime
    return HTBM(V, w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9], w[10])
end

// solve for V that satisfies P=pressure
threadsafe function solveHTBM(pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT, [limH,limL])
    variable pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT, limL, limH
    limL = ParamIsDefault(limL) ? (1-kBracketFactor)*V0 : limL
    limH = ParamIsDefault(limH) ? (1+kBracketFactor)*V0 : limH
   
    Make /free/D w={temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT}
    FindRoots /H=(limH)/L=(limL) /Q /T=(kTol)/Z=(pressure) HTBM_wrapper, w
    return V_flag==0 ? V_Root : NaN
end

// calculate integral of V dP from P=P1 to P=P2 at temperature = T
threadsafe function integrateVolume(P1, P2, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable P1, P2, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT
   
    variable V2 = solveHTBM(P2, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable V1 = solveHTBM(P1, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    return P2*V2 - P1*V1 - integrateP(V2, T, Tref, V1, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
end

// isothermal integral of 3rd order Birch-Murnaghan wrt V from V0 to V at T=T
threadsafe function integrateP(V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT
   
    variable KT=K0+dKdT*(T-Tref)
    variable V0T=V0*exp(c1*(T-Tref)+0.5*c2*(T^2-Tref^2)+c3*ln(T/Tref)-c4*(T^-1-Tref^-1)+2*c5*(T^0.5-Tref^0.5))
   
    variable integral = 0
    integral = -(Kprime-4)*V0T^3*V^-2 - (2-3*(Kprime-4))*V0T^(7/3)*V^(-4/3) - (3*(Kprime-4)-4)*V0T^(5/3)*V^(-2/3)
    integral += V0T*((Kprime-4) + (2-3*(Kprime-4)) + (3*(Kprime-4)-4))
    integral *= 9/16*KT
    return integral
end