;+
;NAME:
; clean_time
;CALL: clean_time,datanames [,SORT=SORT]
;PURPOSE:
; Check and fix the time arrays of the input variables.
; If the time array cannot be fixed, remove the corresponding
; time and data values.
;WARNING: This program is not an artificial intelligence. It will not
; work on all problems. Particularly if there are lots of errors
; piled up together or seperated by only a few time intervals.
; It should be able to handle up to 40 consecuous time errors if
; there is a bit of clean time (~40 points) before the next errors.
;INPUTS: Array of strings. Each string should be associated with a data
; quantity. (see the store_data and get_data routines)
; alternately, datanames could be one data structure.
; eg: get_data,'Ve_3dp',dat=dat & clean_time,dat
;OPTIONAL INPUTS: seed: an input to the filter routine. Default is 5.
; seed is a filtering width factor. If the default fails, choose
; another seed that is not a multiple of 5. Seed will affect
; run time.
;KEYWORDS:
; SORT: this is a kludge. assume all the time data is good.
; this will sort the data monotonically.
; PLOT: plot the new time array to see if it is acceptable.
; this is mostly a debugging tool
; LOUD: Print lots of messages as the program chugs along.
;OUTPUTS: Prints a strange report: eg: no errors found, data
; cleaned, xx% data loss, data unrecoverable, too many errors.
;SIDE EFFECTS:If sucessful in finding errors, will alter the data array
; and time array.
;LAST MODIFICATION: @(#)clean_time.pro 1.5 95/12/04
;AUTHOR: Frank V. Marcoline
;-
FUNCTION nglitches,array ;find places where array is not monotonic
l = n_elements(array)
delta = where( (array(1:l-1l)-array(0:l-2l)) LT 0, count )
return,count
END
FUNCTION rm_ele,arr,ele ;remove element ele from array arr
n = n_elements(arr)
IF ((ele GE n) OR (ele LT 0)) THEN BEGIN
print,'Element ',strcompress(ele),$
' out of range: 0 < # < ',strcompress(n-1)
return,-1
ENDIF
IF ele EQ 0 THEN arr = arr(1:n-1) $
ELSE IF ele EQ n-1 THEN arr = arr(0:n-2) $
ELSE arr = [arr(0:ele-1),arr(ele+1:n-1)]
return,arr
END
PRO clean_time,datanames,seed,PLOT=PLOT,LOUD=LOUD,SORT=SORT
on_error,0
IF n_params() NE 2 THEN seed = 5 ELSE seed = long(seed)
type = data_type(datanames)
CASE type OF
8: nparams = 1
7: nparams = (reverse(size(datanames)))(0)
0: BEGIN
print,'No defined parameters entered'
return
ENDCASE
ELSE: nparams = 1
ENDCASE
IF NOT keyword_set(PLOT) THEN PLOT = 0
IF NOT keyword_set(LOUD) THEN LOUD = 0
IF NOT keyword_set(SORT) THEN SORT = 0
FOR i=0l,nparams-1l DO BEGIN
CASE type OF
8: time = datanames.x
7: begin
print,'Dataname: ',datanames(i)
get_data,datanames(i),dat=data,index=index
IF index EQ 0 THEN BEGIN
print,'Unsuccessful attempt to get_data: ',datanames(i)
return
ENDIF ELSE time = data.x
ENDCASE
ELSE: time = datanames
ENDCASE
IF PLOT THEN plot,time MOD time(0),title='Time arary.'
length = n_elements(time)
;;the median is more likely to be representive of the average time step
;;than the average is, because time spikes are often large shifts
IF length GT 1 THEN delta_t = median(-ts_diff(time,1)) $
ELSE BEGIN
print,'Time array contains only one data point'
return
ENDELSE
IF LOUD THEN print,'Median time step:',delta_t
IF NOT SORT THEN BEGIN
;;simplest test, is the time array monotonic and increasing?
glitches = nglitches(time)
IF glitches EQ 0 THEN BEGIN
IF LOUD THEN print,'Time array monotonically increasing.'
ENDIF ELSE BEGIN
IF LOUD THEN print,'Time array not monotonically increasing'
print,'Number of faults in time array: ',glitches
dt = -ts_diff(time,1)
dt = dt(0:length-2)
spike = where((dt GT 5*delta_t) OR (dt LT -5*delta_t),spikes)+1
IF LOUD THEN print,'Spikes and gaps: ',strcompress(spikes)
IF LOUD THEN print,spike
IF spikes GT 1 THEN gaps = lonarr(spikes) ELSE gaps = 0
FOR j=0l,spikes-2 DO BEGIN
IF time(spike(j)) GT time(spike(j)-1) THEN BEGIN
;;test1:see if a time jump has no dips within the next 40 points
test1 = where(time(spike(j)+1:((spike(j)+41)<(length-1))) $
LT time(spike(j)),n1)
;;test2:see if a time jump has no jumps within the previous 40 points
test2 = where(time(((spike(j)-41)>0):(spike(j)-1)) $
GT time(spike(j)-1),n2)
; print,test1,test2
IF n1+n2 EQ 0 THEN BEGIN ;this is a gap, not a spike
gaps(j) = 1
; print,j,spike(j)
ENDIF
ENDIF
ENDFOR
spike = spike(where(gaps eq 0,spikes)) ;remove gaps, leave just spikes
IF LOUD THEN print,'Spikes:',strcompress(spikes),' Gaps: ',strcompress(total(gaps))
; IF LOUD THEN print,'Spike:',spike
;;try to get rid of spikes
k = 0
FOR j=0l,spikes-1 DO BEGIN
IF k EQ 0 THEN BEGIN
IF dt(spike(j)-1) GT 0 THEN $ ;spike, not a dip
REPEAT k = k+1 $ ;find return
UNTIL ((dt(spike(j+k)-1) LT 0) OR (k+j GE spikes-1)) $
ELSE $ ;this is a dip
REPEAT k = k+1 $
UNTIL ((dt(spike(j+k)-1) GT 0) OR (k+j GE spikes-1)) ;find return
;this statement represents a failure of accounting, please fix me
IF (k+j) LT spikes THEN BEGIN
span = spike(j+k)-(spike(j)-1)
tspan = time(spike(j+k))-time(spike(j)-1)
avg = tspan/span
ENDIF
FOR l=1l,span-1 DO time(spike(j)+l-1) = time(spike(j)-1)+(l*avg)
ENDIF ELSE k = k-1
ENDFOR
ENDELSE
glitches = nglitches(time)
print,'Number of glitches left: ',glitches
; med = median(time,70) ;construct a replacement time series
; ;70 wil handle up to 32 contiguous time errors
; trouble = where(ts_diff(med,1) eq 0,spots)+1
; IF spots GT 1 THEN BEGIN
; trouble = trouble(0:spots-2)
; spots = spots-1
; ENDIF
; inarow = lonarr(spots)+1
; k = 0
; FOR j=0,spots-2 DO BEGIN
; IF trouble(j+1)-trouble(j) EQ 1 THEN inarow(k) = inarow(k)+1 $
; ELSE k = k+1
; ENDFOR
; inarow = inarow(where(inarow GT 1),spikes)
; FOR j=seed,length/seed,seed DO $
; IF nglitches(med) GT 0 THEN med = median(med,j*seed)
; print,'glitches after median filtering:',nglitches(med)
; ;;we knocked out the salt and pepper noise. now replace the points
; ;;with reasonable values.
;
; span = 0
; print,'spots',spots
; FOR j=0,spots-1 DO BEGIN
; k = 0
; REPEAT BEGIN
; k = k+1
; IF trouble(j)+k LT length THEN $
; span = time(trouble(j)+k)-time(trouble(j)-1)
; ENDREP UNTIL ((k GE 34) OR (span GT 0.5*moments(0)))
; print,k
; IF span GT 0 THEN BEGIN ;fix the time elements
; FOR l=0,k-1 DO BEGIN
; time(trouble(j)+l) = time(trouble(j)-1)+span*(l+1.0)/(k+1)
; ENDFOR ;loop variable l
; ENDIF ELSE BEGIN
; ;;throw out a bunch of junk...
; ENDELSE
; ENDFOR ;loop variable j
ENDIF ;end of: if not sort
IF PLOT THEN plot,time MOD time(0),title='Repaired time array.'
;prepair for departure...
IF SORT THEN BEGIN
time_index = sort(time)
time = time(time_index)
ENDIF
CASE type OF
8: BEGIN
IF SORT THEN BEGIN
datanames.y = datanames.y(time_index,*)
ENDIF
datanames.x = time
ENDCASE
7: BEGIN
data.x = time
IF SORT THEN BEGIN
data.y = data.y(time_index,*)
ENDIF
store_data,datanames(i),dat=data
ENDCASE
ELSE: datanames = time
ENDCASE
ENDFOR ;loop variable i
RETURN
END