;Calculate osculating orbital elements for a cartesian state vector
;Input
;  pv - Six element array, ECI position and velocity of spacecraft. First three
;       elements are position and second three are velocity
;  gm - optional gravity constant. Default constant is Earth WGS-84 km and seconds
;  /verbose - debugging output
;return
;  an array of six primary orbital elements and some secondary elements, in this order
;   [a,e,i,raan,ap,ta,alat]
;   a=semimajor axis (same length unit as position vector)
;   e=eccentricity
;   i=inclination (radians)
;   raan=right ascension of ascending node (radians)
;   ap=argument of periapse (radians)
;   ta=true anomaly (radians)
;   alat=argument of latitude (radians)
;
;note:
; 1) Use ECI coordinates. ECEF will silently give you the wrong answer.
; 2) Use kilometers and km/s as units for input vector. This should be the natural units of scpv()
; 3) To calculate the argument of latitude:
;      pv=scpv(time); Not scpv(time,/ecef)!
;      ele=pv2ele(pv)
;      arg_lat=ele[4]+ele[5]; argument of latitude is argument of periapse plus true anomaly
;      if arg_lat lt 0 then arg_lat+=2*!DPI
;      if arg_lat gt 2*!DPI then arg_lat-=2*!DPI
function pv2ele,pv,GM,verbose=verbose
  if n_elements(gm) eq 0 then gm=398600.4418
  z=[0,0,1]
  rv=pv[0:2]
  vv=pv[3:5]
  r=vlength(rv)
  v=vlength(vv)
  if keyword_set(verbose) then print,string(rv,r,format='(%"Position: <%f, %f, %f> (%f) km")')
  if keyword_set(verbose) then print,string(vv,v,format='(%"Velocity: <%f, %f, %f> (%f) km/s")')
  hv=crossp_grid(rv,vv)
  h=vlength(hv)
  if keyword_set(verbose) then print,string(hv,h,format='(%"AngMmntm: <%e, %e, %e> (%e) km^2/s")')
  nv=crossp_grid(z,hv)
  if keyword_set(verbose) then print,string(nv,  format='(%"ANodeVec: <%e, %e, %e>")')
  GMEVCoeff1=V*V-GM/R;
  GMEVCoeff2=-dotp(RV,VV);
  GMEV=((RV*GMEVCoeff1)+(VV*GMEVCoeff2));
  if keyword_set(verbose) then print,string(GMEV,format='(%"GM*EcVec: <%e, %e, %e>")')
  EV=GMEV/GM;
  if keyword_set(verbose) then print,string(EV,  format='(%"EcVec:    <%e, %e, %e>")')
  E=norm(EV);
  P=(H^2)/GM;
  A=P/(1-E*E);
  if e lt 1 then begin
    Apoapse=A*(1+e)
    Period=2*!DPI/sqrt(A*A*A/GM)
  end else begin
    Apoapse=!values.d_infinity
    Period=Apoapse
  end
  Periapse=A*(1-E)
  I=vangle(HV,z)
  if(I eq 0) or I eq !DPI then begin
    LAN=0
    if(E eq 0) then begin
      AP=0
    end else begin
      AP=atan(GMEV[1],GMEV[0])
    end
  end else begin
    LAN=atan(nv[1],nv[0])
    if lan lt 0 then lan+=2*!DPI
    if e eq 0 then begin
      ap=0
    end else begin
      ap=vangle(nv,gmev)
      if gmev[2]<0 then begin
        ap=2*!DPI-ap
      end
    end
  end
  if E eq 0 then begin
    TA=vangle(nv,rv)
    if rv[2]<0 then ta=2*!DPI-TA
  end else begin
    TA=vangle(rv,gmev)
    if dotp(rv,vv) lt 0 then TA=2*!DPI-TA
  end
  alat=ap+ta; argument of latitude is argument of periapse plus true anomaly
  if alat lt 0 then alat+=2*!DPI
  if alat gt 2*!DPI then alat-=2*!DPI
  return,[a,e,i,lan,ap,ta,alat]
end