1function tau = tauf(taup, e2)
2%TAUF   tan(phi)
3%
4%   TAUF(taup, e2) returns tangent of phi in terms of taup the tangent of
5%   chi.  e2, the square of the eccentricity, is a scalar; taup can be any
6%   shape.
7
8  numit = 5;
9  e2m = 1 - e2;
10  tau = taup / e2m;
11  stol = 0.1 * sqrt(eps) * max(1, abs(taup));
12  g = isfinite(tau);
13  for i = 1 : numit
14    if ~any(g), break, end
15    tau1 = hypot(1, tau);
16    sig = sinh( eatanhe( tau ./ tau1, e2 ) );
17    taupa = hypot(1, sig) .* tau - sig .* tau1;
18    dtau = (taup - taupa) .* (1 + e2m .* tau.^2) ./ ...
19           (e2m * tau1 .* hypot(1, taupa));
20    tau(g) = tau(g) + dtau(g);
21    g = g & abs(dtau) >= stol;
22  end
23end
24