Since this is more or less a tehnical issue, please open an issue on my local Math fork: GitHub - rok-cesnovar/math: The Stan Math Library is a C++ template library for automatic differentiation of any order using forward, reverse, and mixed modes. It includes a range of built-in functions for probabilistic modeling, linear algebra, and equation solving. and we can go into details there.
We need to make the function templated and work with stan::math::var. That is all. It currently uses Eigen::MatrixXd everywhere and we need to make sure those accept vars.
We dont need to write a /rev version if it.
Its much eaasier for me to add this to a fork of Math and make stan call it from there as you propose. The current process is a bit clumsy for me, but obviously much easier to plug and play with off-the-shelf Rstan.