NavList:

Tony,

What makes the formula look messy is that altitude, or rather zenith distance, is "hidden" behind the archav function. If you allow zenith distance to be a separate variable, z, it will look better:

hav a = sec φ · csc z · sqrt{ hav[ 180° - (φ+δ+z) ] · hav(φ-δ+z) }
where
hav z = hav(φ-δ) + [ 1 - hav(φ-δ) - hav(φ+δ) ] · hav t

Then I utilized the fact that hav(180°-x) = 1 - hav x, which adds to the mess.

Further, 
sec φ = 1 / cos φ = 1 / (1-2·hav φ)
and
csc z = 1 / sin z = 1 / sqrt[ hav(2·z) ]

But, as you wrote, "no matter how long or complicated" ...

If you try to avoid using z at all, then you'll end up with tan and cot functions which do not, as far as I am aware, easily translate to haversines. And finding their inverses ...

The simplest method is, I believe, 

hav a = [ hav p - hav(φ-h) ] / [ 1 - hav(φ-h) - hav(φ+h) ]
where
p is polar distance, p = 90° - δ, and h is altitude, h = 90° - z.

Lars