pro read_pds

; Paul Withers, 2004.09.08
; Center for Space Physics, Boston University

; Process the PDS merimu_1001 data for mer1 = merb

common marsblk
; Properties of Mars and physical constants
; Defined in marsprops

common sctblk
; Properties of spacecraft and datasets
; Defined in sctprops

pdsfile = sctpropsstruct.pdsfile
datafile = sctpropsstruct.datafile
nraw = sctpropsstruct.nraw
sa = sctpropsstruct.sa
nsplit = sctpropsstruct.nsplit
t0 = sctpropsstruct.entrystate.t
i0 = sctpropsstruct.i0
i1 = sctpropsstruct.i1

junk=''
data_array = strarr(nraw) 
t_array = dblarr(nraw)
accx_array = t_array
accy_array = t_array
accz_array = t_array

; Read in raw data file
openr, lun1, pdsfile, /get_lun
i=0
while i lt nraw do begin
readf, lun1, junk
data_array(i) = junk
xxx = str_sep(strcompress(junk), ' ')
accx_array(i) = xxx(12)
accy_array(i) = xxx(13)
accz_array(i) = xxx(14)
t_array(i) = double(strmid(xxx(0),5,8)) ; fix the start time...
i=i+1
endwhile
free_lun, lun1

t = t_array

; No data from IMU-B after backshell separation
t = t(0:nsplit)
accx_array = accx_array(0:nsplit)
accy_array = accy_array(0:nsplit)
accz_array = accz_array(0:nsplit)

; reject duplicate data points
xx = uniq(t)
t = t(xx)
accx_array = accx_array(xx)
accy_array = accy_array(xx)
accz_array = accz_array(xx)

; reject data points at t=xxx.0000 that don't fit in sequence
xx = where((t - floor(t) ne 0.0) and (t - floor(t) ne 0.25) and $
 (t - floor(t) ne 0.5) and (t - floor(t) ne 0.75) )
;xx = where(t ne floor(t))
t = t(xx)
accx_array = accx_array(xx)
accy_array = accy_array(xx)
accz_array = accz_array(xx)

; fill in a blank
t = [t(0:i0), t(i0) + (t(i1)-t(i0))* ([1,2]/3.), t(i1:*)]
accx_array = [accx_array(0:i0), accx_array(i0) + (accx_array(i1)-accx_array(i0))* ([1,2]/3.), accx_array(i1:*)]
accy_array = [accy_array(0:i0), accy_array(i0) + (accy_array(i1)-accy_array(i0))* ([1,2]/3.), accy_array(i1:*)]
accz_array = [accz_array(0:i0), accz_array(i0) + (accz_array(i1)-accz_array(i0))* ([1,2]/3.), accz_array(i1:*)]


; subtract pre-entry mean and see what happens...
;aa = where(t gt t0-700. and t le t0-100.)
;accx_array = accx_array - mean(accx_array(aa))
;accy_array = accy_array - mean(accy_array(aa))
;accz_array = accz_array - mean(accz_array(aa))


; reject points before entry state, accept possible 0.25s offset here...
aa = where(t ge t0)
t = t(aa)
accx_array = accx_array(aa)
accy_array = accy_array(aa)
accz_array = accz_array(aa)

; To determine start time, work with (ax^2+ay^2+az^2)^0.5 only...
a = sqrt( accx_array^2 + accy_array^2 + accz_array^2 )

; Flip sign, do not always use |a| in rest of code...
accz_array = (-1) * accz_array

saccx_array = accx_array*0.
saccy_array = saccx_array
saccz_array = saccx_array
saccx_array(*) = sa
saccy_array(*) = sa
saccz_array(*) = sa


openw, lun1, datafile, /get_lun
printf, lun1, t
printf, lun1, accx_array ; a*0.
printf, lun1, accy_array ; a*0.
printf, lun1, accz_array ; a
printf, lun1, saccx_array 
printf, lun1, saccy_array 
printf, lun1, saccz_array 
free_lun, lun1

print, 'n = ', n_elements(t)
print, 'min(t) = ', min(t)

n = n_elements(t)
sctpropsstruct.n = n
print, '(-) n = ', n
; Number of timesteps/measurements

end