I also thought I'd try and use your query as well as the calculation I used previously and that also works, with almost the same results (minor differences that I am happy to overlook, probably due to the 3961 vs 3963.1676)
3963.1676 * (2 * atn2(sqrt((sin(radians(cast(po.Lat as float) - cast(ps.Lat as float))/2) * sin(radians(cast(po.Lat as float) - cast(ps.Lat as float))/2) + cos(radians(ps.Lat))
* cos(radians(po.Lat)) * sin(radians(cast(po.Long as float) - cast(ps.Long as float))/2) * sin(radians(cast(po.Long as float) - cast(ps.Long as float))/2))), sqrt(1 - (sin(radians(cast(po.Lat as float) - cast(ps.Lat as float))/2) * sin(radians(cast(po.Lat as float) - cast(ps.Lat as float))/2) + cos(radians(ps.Lat))
* cos(radians(po.Lat)) * sin(radians(cast(po.Long as float) - cast(ps.Long as float))/2) * sin(radians(cast(po.Long as float) - cast(ps.Long as float))/2)))))
So if anyone else is interested, either one works pretty much.