;+
;
;  Name: CALCULATE_LSHELL
;  
;  Purpose: Calculates the l shell value given spacecraft position data in GSM coordinates
;  
;  Inputs: Spacecraft position vector, an array of [time, x, y, z]
;          where time is double precision, and x, y, z, are in GSM coordinates and have units in RE
;
;  Outputs: The l shell value
;
;  Usage: lshell = calculate_lshell(pos_gsm_re)
;  
;  See Also: tkm2re.pro
;
;  History:
;
;  Notes: This routine requires IDL geopack routines
;
; $LastChangedBy:  $
; $LastChangedDate:  $
; $LastChangedRevision:  $
; $URL:  $
;-

FUNCTION calculate_lshell, gsm_re

    compile_opt idl2


lshell=make_array(n_elements(gsm_re[0,*]), /double)

FOR i=0,n_elements(gsm_re[0,*])-1 do begin

  time=time_struct(gsm_re[0,i])
   
  ;recalculate geomagnetic dipole
  geopack_recalc, time.year, time.doy, time.hour, time.min, time.sec, tilt=tilt
  
  ;trace field line from position to equator
  geopack_trace, gsm_re[1,i], gsm_re[2,i], gsm_re[3,i], 1, 0, out_foot_x, $
      out_foot_y, out_foot_z,/igrf,/equator,/refine,rlim=60.0D

  ;calculate radial distance at equator
  lshell[i] = sqrt(out_foot_x^2+out_foot_y^2+out_foot_z^2)

ENDFOR

RETURN, lshell

END