function jac=xyzJacobian(x,y,z)

% algorithm from /usr/local/postprocessing/src/xyztogeo.c
% use jac to rotate vectors as
%  
%   n         x             x          n
%   e = jac * y       or    y = jac' * e
%   u         z             z          u
%
%   or 
%
%   neu = xyz * jac'  or   xyz = neu * jac

earth_rad = 6378137.0;
f = 1.0 / 298.257222101;
deg2rad = pi / 180.0;
twopi = pi * 2.0;
tolerance = 0.0001;

eccsq = 2.0 * f - f * f;
eq_radius = sqrt(x*x+y*y);
lat_p = atan2(z, eq_radius);
lon_i = atan2(y,x);
if lon_i < 0.0 
   lon_i = lon_i + twopi;
end

h_p = 0.0;
niter=0;
converged = 0;
while (converged == 0)
   rad_curve = earth_rad / sqrt(1.0 - eccsq * sin(lat_p) * sin(lat_p));
   rad_lat = eq_radius * ( 1.0 - eccsq * rad_curve/ (rad_curve + h_p) );
   lat_i = atan2(z, rad_lat);

   if abs(lat_i) < (pi / 4.0) 
        h_i = eq_radius / cos(lat_i) - rad_curve;
   else
        h_i = z / sin(lat_i) - (1.0 - eccsq)*rad_curve;
   end

   if abs(h_i-h_p) < tolerance &  abs(lat_i-lat_p)*rad_curve < tolerance 
        converged = 1;
   end

   niter=niter+1;
   if niter > 50
        'xyztogeo error : failure to converge'
        converged = 1;
   end

   h_p = h_i;
   lat_p = lat_i;
end

jac = [-sin(lat_i)*cos(lon_i)   -sin(lat_i)*sin(lon_i)    cos(lat_i);
       -sin(lon_i)               cos(lon_i)               0.0;
        cos(lat_i)*cos(lon_i)    cos(lat_i)*sin(lon_i)    sin(lat_i)];

