pro sort_time, day, time, subs
;
; Input day and time, program return sorted day and time and array of
; pointers to the elements of the original array(subs)
; Original subscript order is maintained when values are equal.
;
;
ms_day = 86400000            ; 1 day = 86400000 msec
n = n_elements(day)          ; number of elements in day array
;
;  Check the time array, to make sure there are no more than 1 day in
;  this array.  If there are more than 1 day in time array, modify the 
;  day array by increment 1 day, and modify time array by subtract 
;  864000000 msec(1 day).
;  Also generate subs array, containing the subscript of original array.
;
nd = time/ms_day                   ; number of days in time array
day = day + nd                     ; increment day array by nd
time = time - nd*ms_day            ; new time array
subs = long(indgen(n))             ; array index
;
;  Original subscript order is maintained when values are equal.
;
day_sort_indx = bsort(day, day_sort)  ; day_sort is day in sorted order
;  
day  = day_sort                       ; day, time, and subs are in
time = time(day_sort_indx)            ; ascending order
subs = subs(day_sort_indx)
;
;  Find the "same day", sort in order of the time within the "same day"
;  group
;
ptr_test = 0                          ; pointer of testing location
done = 0
while ( done eq 0) do begin
   day_test = day(ptr_test)
   sameday = where(day(ptr_test:n-1) eq day_test)
   n_sd = n_elements(sameday)
;
   if (n_sd eq 1) then begin       ; no other day are same as day_test
      ptr_test = ptr_test+1        ; no same day, continue the test
   endif else begin
      ptr_end = ptr_test + n_sd - 1      ; ending pointer of same day
;
;     sort from ptr_test to ptr_end in time. 
;
      time_sort_indx = bsort(time(ptr_test:ptr_end), time_sort)
;    
;     time_sort_indx is from 0 to number of "same day", not an actual index.
;     Actual index is time_sort_indx plus ptr_test.
;
      time_sort_indx = time_sort_indx + ptr_test     ; actual subscript
      day(ptr_test:ptr_end) = day(time_sort_indx)
      time(ptr_test:ptr_end) = time(time_sort_indx)
      subs(ptr_test:ptr_end) = subs(time_sort_indx)
      ptr_test = ptr_end + 1             ; new testing position
   endelse           
   if (ptr_test gt n-2) then done = 1
endwhile
end