; Returns the eccentric anomaly given the mean anomaly and eccentricity
; 
function e_anomaly,ma1,ecc
  n = floor(ma1/2./!dpi)
  ma = ma1-n*2d*!dpi
  ea = ma
  i=0
  mx=1
  err_accept = (size(/type,ma1) eq 5) ? 1e-14  : 1e-7
  while mx gt err_accept do begin
     ea = ma + ecc*sin(ea)
     ea = ea - (ea-ecc*sin(ea)-ma)/(1-ecc*cos(ea))
     err = ea-ecc*sin(ea) - ma
     mx = max(abs(err))
;     plot,ma,abs(err),yrange=[1e-20,1e3],/ylog,ps=-3 & wait,.1
     i=i+1
     if i gt 30 then begin
       print,'Warning, eccentric anomaly error: ',mx
       mx = 0.
     endif
  endwhile
;  print,i
  return,ea+n*2d*!dpi
end