function readdat, filepath, datestr, timestr, datacount

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

;read the GPS NMEA header
nmea_string = string('')
readf, dat_fh, nmea_string
point_lun, dat_fh, 0

;parse date and time 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)

;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(200) ;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 channel count and data count
channelcount = 2
datacount = (filesize - (last_gps + 2)) / (channelcount * 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 bytes
;from the beginning of the file. Skipping 2 bytes of newline 
bdata     = uintarr(datacount * channelcount)
data      = assoc(dat_fh, bdata, last_gps+2)

;the following assigns the first data array to bx, swapping byte ordering in
;the process. MIGHT NOT BE NECESSARY SINCE LINUX AND MAC ARE BOTH LITTLE ENDIAN
dbdt      = data(0)

;=========================================================
;CHECKING THE GPS BIT:
;gps_bit  = fix((dbdt and (2^15)))/2^15
;indices = indgen(864000) *2
;ccc= gps_bit(indices)   ; these are the gps bit indices - 0's and 1's
;aaa=where(ccc eq 1)
;ddd = (aaa-shift(aaa,1))-10
;eee = where(abs(ddd) gt 0) *10

;kkk=intarr(432000)
;kkk(eee) = 1

;=========================================================

; the 32-bit time stamp data and their indices
time_arr = ulonarr(datacount*2/22)
time_inds1 = ulonarr(datacount*2/22)
time_inds2 = ulonarr(datacount*2/22)
k=0UL
for i = 0, datacount*2-1, 22 do begin
    time1 = ulong(dbdt[0:1], 0)
    time_arr(k) = ulong(dbdt[i:i+1], 0)
    time_inds1(k) = i;
    time_inds2(k) = i+1;
    k = k+1
endfor
time_inds = [time_inds1, time_inds2]

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

;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-n_elements(time_arr)) * 2 ;to get an array of even numbered indices

;then
data_arr  = fltarr(3, datacount-n_elements(time_arr))
data_arr(0, *) = dbdt(x_indices)
data_arr(1, *) = dbdt(x_indices+1)

; Create fractional time stamps for each set of x-y data points
k=0UL
for i = 0, n_elements(data_arr(0,*))-1, 10 do begin
     data_arr(2, i:i+9) = time_arr(k)+indgen(10)/10.
     k = k+1
endfor

;========================================================================
;plot,data_arr(0,000:4300),yrange=[-1,1]
;oplot,ccc(2000:9000)+2
;oplot,kkk(2000:9000)-2
;======================================================================

;to be done:
;- I have not done anything to determine array sizes, which is something that
;does need to be completed before implemeneting the assoc() function. The
;file information needs to be queried for this, as well as the size of the
;NMEA string.

; Also, I have no function to query the number of channels that are sampled.
;I believe this is contained in the 3 bytes following the NMEA string
return, data_arr
end