derivs.hpp

        auto Ar02 = tdx::template get_Ar02<be>(model, T, rhotot, molefrac);
        auto Z = 1.0 + Ar01;
        auto dZdrho = (Ar01 + Ar02)/rhotot; // (dZ/rhotot)_{T,x}
        return std::make_tuple(log(Z), Z, dZdrho);
    }
    
    /***
    * \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. molar density at constant temperature and mole fractions
    *
    * \f[
    * \deriv{ \ln\vec\phi}{\rho}{T,\vec x} = \frac{1}{RT}H(\Psi_r)\vec x - \frac{1}{Z}\deriv{Z}{\rho}{T,\vec x}
    * \f]
    */
    template<ADBackends be = ADBackends::autodiff>
    static auto get_d_ln_fugacity_coefficients_drho_constTmolefracs(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = forceeval(rhovec.sum());
        auto molefrac = (rhovec / rhotot).eval();
        auto R = model.R(molefrac);
        auto [lnZ, Z, dZdrho] = get_lnZ_Z_dZdrho(model, T, rhovec);
        auto hessian = build_Psir_Hessian_autodiff(model, T, rhovec);
        return forceeval((1/(R*T)*(hessian*molefrac.matrix()).array() - dZdrho/Z).eval());
    }
    
    /***
    * \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. molar volume at constant temperature and mole fractions
    *
    * \f[
    * \deriv{ \ln\vec\phi}{v}{T,\vec x} = \deriv{ \ln\vec\phi}{\rho}{T,\vec x}\deriv{ \rho}{v}{}
    * \f]
    */
    template<ADBackends be = ADBackends::autodiff>
    static auto get_d_ln_fugacity_coefficients_dv_constTmolefracs(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = forceeval(rhovec.sum());
        auto drhodv = -rhotot*rhotot; //  rho = 1/v; drho/dv = -1/v^2 = -rho^2
        return get_d_ln_fugacity_coefficients_drho_constTmolefracs(model, T, rhovec)*drhodv;
    }
    
    /***
    * \brief Calculate the derivative of the natural logarithm of fugacity coefficient of each component w.r.t. mole fraction of each component, at constant temperature and molar density
    *
    * \f[
    * \deriv{ \ln\vec\phi}{\vec x}{T, \rho}
    * \f]
    */
    template<ADBackends be = ADBackends::autodiff>
    static auto get_d_ln_fugacity_coefficients_dmolefracs_constTrho(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = forceeval(rhovec.sum());
        auto molefrac = (rhovec / rhotot).eval();
        auto R = model.R(molefrac);
        
        using tdx = TDXDerivatives<Model, Scalar, VectorType>;
        auto Ar01 = tdx::template get_Ar01<be>(model, T, rhotot, molefrac);
        auto Z = 1.0 + Ar01;
        VectorType dZdx = rhotot*build_d2alphardrhodxi_constT(model, T, rhotot, molefrac);
        Eigen::RowVector<decltype(rhotot), Eigen::Dynamic> dZdx_Z = dZdx/Z;
        
        // Starting matrix is from the first term
        auto hessian = build_Psir_Hessian_autodiff(model, T, rhovec);
        Eigen::ArrayXXd out = rhotot/(R*T)*hessian;
        
        // Then each row gets the second part
        out.rowwise() -= dZdx_Z.array();
        return out;
    }
    
    template<ADBackends be = ADBackends::autodiff>
    static auto get_d_ln_fugacity_coefficients_dmolefracs_constTrho1(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = forceeval(rhovec.sum());
        auto molefrac = (rhovec / rhotot).eval();
        auto R = model.R(molefrac);
        
        auto hessian = build_Psir_Hessian_autodiff(model, T, rhovec);
        // Starting matrix is from the first term
        Eigen::ArrayXXd out = 1/(R*T)*rhotot*hessian;
        return out;
    }

    /***
    * \brief Calculate the temperature derivative of the chemical potential of each component
    * \note: Some contributions to the ideal gas part are missing (reference state and cp0), but are not relevant to phase equilibria
    */
    static auto get_dchempotdT_autodiff(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = rhovec.sum();
        auto molefrac = (rhovec / rhotot).eval();
        auto rhorefideal = 1.0;
        return (build_d2PsirdTdrhoi_autodiff(model, T, rhovec) + model.R(molefrac)*(rhorefideal + log(rhovec/rhorefideal))).eval();
    }

    /***
    * \brief Calculate the temperature derivative of the pressure at constant molar concentrations
    */
    static auto get_dpdT_constrhovec(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = rhovec.sum();
        auto molefrac = (rhovec / rhotot).eval();
        auto dPsirdT = get_dPsirdT_constrhovec(model, T, rhovec);
        return rhotot*model.R(molefrac) - dPsirdT + rhovec.matrix().dot(build_d2PsirdTdrhoi_autodiff(model, T, rhovec).matrix());
    }

    /***
    * \brief Calculate the molar concentration derivatives of the pressure at constant temperature
    */
    static auto get_dpdrhovec_constT(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = rhovec.sum();
        auto molefrac = (rhovec / rhotot).eval();
        auto RT = model.R(molefrac)*T;
        auto [func, grad, hessian] = build_Psir_fgradHessian_autodiff(model, T, rhovec); // The hessian matrix
        return (RT + (hessian*rhovec.matrix()).array()).eval(); // at constant temperature
    }

    /***
    * \brief Calculate the partial molar volumes of each component
    * 
    * \f[
    * \hat v_i = \left(\frac{\partial V}{\partial n_i}\right)_{T,V,n_{j \neq i}}
    * \f]
    */
    static auto get_partial_molar_volumes(const Model& model, const Scalar& T, const VectorType& rhovec) {
        auto rhotot = rhovec.sum();
        auto molefrac = (rhovec / rhotot).eval();
        using tdx = TDXDerivatives<Model, Scalar, VectorType>;
        auto RT = model.R(molefrac)*T;

        auto dpdrhovec = get_dpdrhovec_constT(model, T, rhovec);
        auto numerator = -rhotot*dpdrhovec;
        
        auto ders = tdx::template get_Ar0n<2>(model, T, rhotot, molefrac);
        auto denominator = -pow2(rhotot)*RT*(1 + 2*ders[1] + ders[2]);
        return (numerator/denominator).eval();
    }
    
    static VectorType get_Psir_sigma_derivs(const Model& model, const Scalar& T, const VectorType& rhovec, const VectorType& v) {
        autodiff::Real<4, double> sigma = 0.0;
        auto rhovecad = rhovec.template cast<decltype(sigma)>(), vad = v.template cast<decltype(sigma)>();
        auto wrapper = [&rhovecad, &vad, &T, &model](const auto& sigma_1) {
            auto rhovecused = (rhovecad + sigma_1 * vad).eval();
            auto rhotot = rhovecused.sum();
            auto molefrac = (rhovecused / rhotot).eval();
            return forceeval(model.alphar(T, rhotot, molefrac) * model.R(molefrac) * T * rhotot);
        };
        auto der = derivatives(wrapper, along(1), at(sigma));
        VectorType ret(der.size());
        for (auto i = 0; i < ret.size(); ++i){ ret[i] = der[i];}
        return ret;
    }
};

template<int Nderivsmax, AlphaWrapperOption opt>
class DerivativeHolderSquare{
    
public:
    Eigen::Array<double, Nderivsmax+1, Nderivsmax+1> derivs;
    
    template<typename Model, typename Scalar, typename VecType>
    DerivativeHolderSquare(const Model& model, const Scalar& T, const Scalar& rho, const VecType& z) {
        using tdx = TDXDerivatives<decltype(model), Scalar, VecType>;
        static_assert(Nderivsmax == 2, "It's gotta be 2 for now");
        AlphaCallWrapper<opt, Model> wrapper(model);
        
        auto AX02 = tdx::template get_Agen0n<2>(wrapper, T, rho, z);
        derivs(0, 0) = AX02[0];
        derivs(0, 1) = AX02[1];
        derivs(0, 2) = AX02[2];
        
        auto AX20 = tdx::template get_Agenn0<2>(wrapper, T, rho, z);
        derivs(0, 0) = AX20[0];
        derivs(1, 0) = AX20[1];
        derivs(2, 0) = AX20[2];
        
        derivs(1, 1) = tdx::template get_Agenxy<1,1>(wrapper, T, rho, z);
    }
};



}; // namespace teqp