pro ulf_read_dat, $
   filepath, $
   datestr, $
   timestr, $
   latstr, $
   latdstr, $
   lonstr, $
   londstr, $
   datacount, $
   data_arr


;open the data file:
dat_fh = 1
openr, dat_fh, filepath, /get_lun ;open the data file

;read the GPS NMEA header
;example:
;$GPRMC,084222,A,7803.6598,N,01337.2191,E,0.000,0.0,310813,2.9,E*78
;RMC          recommended mimimum sentence C
;084222       fix taken at 08:24:22
;A            status: a active; v void
;7803.6598,N  lat - 78 deg. 03.6598' N
;01337.2191,E lon - 13 deg. 37.2191' E
;0.000        speed
;0.0          track angle
;310813       date - 31 of August, 2013
;2.9,E        magnetic variation
;*78          checksum
nmea_string = string('')
readf, dat_fh, nmea_string
point_lun, dat_fh, 0
print, nmea_string
;parse date, time, latitude, longitude from the NMEA string
nmea_chars = strarr(strlen(nmea_string))
for i=0, strlen(nmea_string) - 1 do begin
    nmea_chars(i) = strmid(nmea_string, i, 1)
endfor
delim = where(nmea_chars eq ',')
timestr = strmid(nmea_string, delim(0)+1, 6)
datestr = strmid(nmea_string, delim(8)+1, 6)
year = strmid(datestr, 4, 2)

;subtract 2^10 (7168) day rollover uncertainty offset
if year eq '35' then begin
   month = strcompress(fix(strmid(datestr,2,2)), /remove_all)
   day = strcompress(fix(strmid(datestr,0,2)), /remove_all)
   py_str = '/usr/bin/python -c "from datetime import timedelta; from datetime import datetime; d = timedelta(days=-7168); print datetime(2035,'+month+', '+day+') + d"'
   spawn, py_str, date_new
   year_new = strmid(date_new, 2, 2)
   month_new = strmid(date_new, 5, 2)
   day_new = strmid(date_new, 8, 2)
   datestr = day_new+month_new+year_new
endif
if year eq '36' then begin
   month = strcompress(fix(strmid(datestr,2,2)), /remove_all)
   day = strcompress(fix(strmid(datestr,0,2)), /remove_all)
   py_str = '/usr/bin/python -c "from datetime import timedelta; from datetime import datetime; d = timedelta(days=-7168); print datetime(2036,'+month+', '+day+') + d"'
   spawn, py_str, date_new
   year_new = strmid(date_new, 2, 2)
   month_new = strmid(date_new, 5, 2)
   day_new = strmid(date_new, 8, 2)
   datestr = day_new+month_new+year_new
endif
;3 * 2^10 rollover
if year eq '75' then begin
   month = strcompress(fix(strmid(datestr,2,2)), /remove_all)
   day = strcompress(fix(strmid(datestr,0,2)), /remove_all)
   py_str = '/usr/bin/python -c "from datetime import timedelta; from datetime import datetime; d = timedelta(days=-21504); print datetime(2035,'+month+', '+day+') + d"'
   spawn, py_str, date_new
   year_new = strmid(date_new, 2, 2)
   month_new = strmid(date_new, 5, 2)
   day_new = strmid(date_new, 8, 2)
   datestr = day_new+month_new+year_new
endif

print, datestr + ' ' + timestr
latstr = strmid(nmea_string, delim(2)+1, 9)
latdstr = strmid(nmea_string, delim(3)+1, 1)

lonstr = strmid(nmea_string, delim(4)+1, 10)
londstr = strmid(nmea_string, delim(5)+1, 1)

;get the size of the file
status = FSTAT(dat_fh)
filesize = status.size

;the next line associates the variable gps_end as a byte array - in other
;words, ALL of the data are interpreted as bytes at this point. This is done
;in order to locate the CR/LF after the NMEA string.
cr        = bytarr(100) ;the maximum length of NMEA string is 80
gps_end   = assoc(dat_fh, cr)

;next line searches the gps_end array for occurrences of a line feed.
last_gps  = where(gps_end(0) eq 10)

;next line gets the index number of the first occurrence of a line feed.
last_gps  = last_gps(0)

;get the channels and clock divider from the 3 bytes after the end of
;the nmea string
bheader = bytarr(3)
assoc_header = assoc(dat_fh, bheader, last_gps+1)
;find number of channel bits set
header = assoc_header[0]
print, header, format='(3z5)'
num_chan = bit_population(header[0])
if num_chan ne 2 then begin
   print, 'more than 2 channels'
   stop
endif
clock_divider = header[1] * 256 + header[2]
;sample rate in Hz is 10KHz clock / number of channels / clock_divider
sample_rate = 10000 / num_chan / clock_divider

;get the channel count and data count
datacount = (filesize - (last_gps + 4)) / (num_chan * 2)

;The next line associates the variable 'data' as a uintarr. in addition, the
;association is made such that the first integer is located at last_gps+4 bytes
;from the beginning of the file. This comes from Paul's information
bdata     = uintarr(datacount * num_chan)
data      = assoc(dat_fh, bdata, last_gps+4)

;the following assigns the first data array to bx, swapping byte ordering in
;the process.
dbdt      = swap_endian(data[0])

;finally, the 16-bit data are masked to 12 bits, converted to signed integers,
; have the midrange point (2048) subtracted and are scaled.
;dbdt      = (fix((dbdt and (2^12-1))) - 2048)/204.8

;don't forget to close the data file!!!!
free_lun, dat_fh

;at this point, all of the data have been read in, although the various
;channels still need to be extracted. This would be done something like:

x_indices = lindgen(datacount) * 2 ;to get an array of even numbered indices

;then
data_arr  = uintarr(num_chan, datacount)
print, 'datacount: ', datacount
data_arr[0, *] = dbdt[x_indices]
data_arr[1, *] = dbdt[x_indices+1]

gps_sync = ishft((reform(data_arr[0, *]) AND '8000'x), -15)
w_sync = where(gps_sync eq 1, sync_count)

if sync_count ge 2 then begin

   sync_delta = w_sync[1:sync_count-1] - w_sync[0:sync_count-2]
   w_no_sync = where(sync_delta ne 10, no_sync_count)
   if no_sync_count gt 0 then print, max(sync_delta)

endif else begin

   print, 'no sync?'

endelse

end