All those functions are built into Stan, so you can just write your own _lpdf as a Stan function. It’d be more efficient to do it in C++ with vectorization and analytical gradients, but it shouldn’t be bad encoded directly in Stan. You might be able to code it up in terms of sufficient statistics to make things faster.
No, I don’t have experience with this distribution.
Nope, I haven’t tried something like this so my only suggestion is to make sure to log and simplify on paper before you try to code it (sort of a duh, but I thought I’d say it, never write a density that has a tgamma in it) but it looks like fun stuff. What kind of data are you using this for?
Good point about Mathematica, I recently did some similar ones so you’re likely to end up with some incomplete gamma functions. If you do, those are gamma_p in Stan and you likely will need a bugfix that’s in-progress.