In your formula for lat:
lat = Atan(y_i / (dist_axis * (1 - f) ^ 2)) * 180 / PI
I think (1-f)^2 should be just (1-f) because sqrt(1-epsilon^2)=1-f since 1-epsilon^2=b^2/a^2=(1-f)^2 See for instance Chauvenet, A Manual of Spherical and Practical Astronomy, formula 539, p. 491. Do you agree or am I wrong?
Poul Riis Denmark