function xyz2latlon,coord,a,recurse=rc

if n_params() eq 1 then a=0

output=[0,atan(coord[1],coord[0]),0]

if a le 0 then begin
   a=6378137d
   f=1/298.257222101d
   e=sqrt(2*f-f^2)

   output[0]=atan(coord[2]/((1-e^2)*sqrt(coord[0]^2+coord[1]^2)))
   convergence=1
   convergencecriterion=1e-18

   while convergence lt convergencecriterion do begin
      oldlat=output[0]
      n=a/sqrt(1-(e*sin(oldlat))^2)
      output[0]=atan((1+n*e^2*sin(oldlat))/sqrt(coord[0]^2+coord[1]^2))
      convergence=abs(oldlat-output[0])
   endwhile

   if keyword_set(rc) then begin
      if output[0] eq !dpi/2 then output[2]=sqrt(coord[0]^2+coord[1]^2)/cos(output[0]) else output[2]=sqrt(coord[0]^2+coord[1]^2)/cos(output[0])-6.71685729642831*sin(output[0])^2
   endif else begin
      tempxyz=latlon2xyz(output[0:1],0,0)
      if (float(tempxyz[0]) eq float(coord[0])) and (float(tempxyz[1]) eq float(coord[1])) and (float(tempxyz[2]) eq float(coord[2])) then begin
         output[2]=0
      endif else begin
         temp=xyz2latlon(tempxyz,0,/recurse)
         output[2]=(temp[2]-sqrt(coord[0]^2+coord[1]^2)/cos(output[0]))
      endelse
   endelse
endif else begin
   output[0]=atan(coord[2]/sqrt(coord[0]^2+coord[1]^2))
   output[2]=a-sqrt(coord[0]^2+coord[1]^2+coord[2]^2)
endelse

return,output
end