FUNCTION get_vel, hts
;
; $Id: get_vel.pro,v 1.1 2007/11/08 18:30:33 nathan Exp $
;
; $Log: get_vel.pro,v $
; Revision 1.1  2007/11/08 18:30:33  nathan
; moved from adam/jprogs
;
   x = DOUBLE(hts.time)
   y = hts.height
   n = hts.plot_num
   
   ;** Remove adjacent points taken at the same time
   q = UNIQ((hts.time - hts(0).time) / (hts(0).time))
   x = x(q)
   y = y(q)
   n = n(q)
   
   ;** Take derivatives along each trajectory
   z = 0
   indeces = 0
   print, 'line 19: n-elements z = ', N_ELEMENTS(z)
   FOR i=0, MAX(n) DO BEGIN
      q = WHERE(n eq i)
      nm = n_elements(q)
         IF (nm ge 5) THEN BEGIN
            p = poly_fit(x(q) - x(MIN(q)), y(q), 2)
            vp = [p(1), 2.* p(2)]
            z = [z, (poly(x(q) - x(MIN(q)), vp))]
            indeces = [indeces, q]
         ;** Do a linear fit if there are less than five data points.  
         ENDIF ELSE BEGIN
            POPUP_MSG,'A track contained too few data points for a quadratic fit:' $
                      +' using linear fit.', TITLE = 'Warning'
            p = SVDFIT(x(q) - x(MIN(q)), y(q), 2)
            vp = p(1)
            z = [z, (vp*x(q))]
            indeces = [indeces, q]
         ENDELSE
   ENDFOR  
   
   ;** Remove dummy first elements
   z = z(1:*)
   indeces = indeces(1:*)
   x = x(indeces)
   y = y(indeces)
   
   ;** Convert the velocities to km/s (from R/s)
;   z = z*1390000. / 2
   
   
   
   q = WHERE((z ge 0) AND (z le 3000))
   IF (q(0) ne -1) THEN BEGIN
      x = x(q)
      y = y(q)
      z = z(q)
   ENDIF
   
   vstrct = {vel:0.0, time:0.0, height:0.0}
   vels = REPLICATE(vstrct, n_elements(x))
   
   vels.vel = z
   vels.time = x
   vels.height = y
   
   RETURN, vels
END