function eq=elevation(z0,Rp,thetap,thetas)


a=-asinh(cot(thetas));

eq=(cos(thetap)+((Rp-z0)/sqrt(Rp^2-(Rp-z0)^2))*sin(thetap))*(asinh((((Rp-z0)/sqrt(Rp^2-(Rp-z0)^2))*cos(thetap)-sin(thetap))/(((Rp-z0)/sqrt(Rp^2-(Rp-z0)^2))*sin(thetap)+cos(thetap)))-a)-Rp/(2*Rp-z0);