function readlcdata, infile, outunit=outunit, double=double, noextcorr=noextcorr

; Smart ASCII reader allowing varying column order and variable delimitation.
; Example files:
;
; #DEFAULT (GRBlog - no header required):
; &| ref  |   t   | filt | mag | emag | lim | refsys | psym
; GCN1494 | 0.013 | V | 15.51  | 0.14 | no  | USNO   | 8 
; GCN1497 | 0.105 | R | 17.70  | 0.09 | no
; GCN1498 | 0.109 | R | 16.6   |      | yes
;
; #EXAMPLE 1
; @inunit=days
; @outunit=days
; @expunit=secs
; & tstart exp mag  emag lim
;  0.105  180 15.6 0.5  no
;  0.110  300 17.7 0.7  no

; Count how many lines in the file
openr, 1, infile
iline = ''
nlines = 0
while not eof(1) do begin
  readf, 1, iline
  nlines = nlines + 1
endwhile
close, 1

; Set up structure
if keyword_set(double) eq 0 then begin
  obs = replicate({source:'', telescope:'', reference:'', t:-1.0, filter:'', mag:0.0, emag:0.0, $
                 uplim:0, refsys:'', findex:-1, flux:0.0, fluxerr:0.0, $
                 tstart:-1.0, tend:-1.0, exp:0.0, cts:0.0, ects:0.0, psym:-1},$
                 nlines)
endif else begin
  obs = replicate({source:'', telescope:'', reference:'', t:-1.0D, filter:'', mag:0.0D, emag:0.0D, $
                 uplim:0, refsys:'', findex:-1, flux:0.0D, fluxerr:0.0D, $
                 tstart:-1.0D, tend:-1.0D, exp:0.0D, cts:0.0D, ects:0.0D, psym:-1},$
                 nlines)
endelse
; source is now a catch-all that can be interpreted as either telescope, instrument, reference/citation, or various
; random uses.  telescope refers to the telescope and reference to the citation exclusively



; Read in the data from the table
npts = 0
openr, 1, infile

; defaults
splitchar = '|'
dsrc = 0
dtel = -1
dref = -1
dt = 1
dtstart = -1
dtend = -1
dexp = -1
ddt = -1
ddtplus = -1
ddtminus = -1
dfilt = 2
dmag = 3
demag = 4
dflux = -1
deflux = -1
dcts = -1
dects = -1
dlim = 5
ddet = -1
dsys = 6
dpsym = 7
dband = -1
dutstart = -1
dutmid = -1
dutend = -1
inunit = 'd'
ebv = 0
skip = 0
if n_elements(outunit) gt 0 then begin
  outunit = strmid(outunit,0,1)
  fixoutunit = 1   ;you can only set the output (storage) unit if it's 
endif else begin   ;    not set by calling function
  outunit = inunit
  fixoutunit = 0
endelse
if keyword_set(outunit) eq 0 then outunit = 'd'
expunit = 's'
fluxconv = 1.0 ; from counts
fluxunitconv = 1.0
    unitconv=((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+ $
               (outunit eq 's')*24.*60.*60. + (outunit eq 'k')*24*3.6 ) / $
              ((inunit eq 'd')+(inunit eq 'h')*24.+(inunit eq 'm')*24.*60.+$
               (inunit eq 's')*24.*60.*60. + (inunit eq 'k')*24*3.6 )
    expconv = ((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+$
               (outunit eq 's')*24.*60.*60. + (outunit eq 'k')*24*3.6 ) / $
              ((expunit eq 'd')+(expunit eq 'h')*24.+(expunit eq 'm')*24.*60.+$
               (expunit eq 's')*24.*60.*60. + (expunit eq 'k')*24*3.6 )
    sec2outunit = ((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+ $
               (outunit eq 's')*24.*60.*60. + (outunit eq 'k')*24*3.6 ) / (24.*60.*60.)

defaultfilt = 'unk'
defaultsrc = ''
defaultref = ''
defaulttel = ''
defaultsym = -1
defaultlim = 0

i = 0
for l = 0, nlines-1 do begin
  readf, 1, iline

  ; ignore leading/trailing spaces
  iline = strtrim(iline,2)

  ; ignore blank lines and comments
  if strlen(iline) eq 0 then continue
  if strmid(iline,0,1) eq '#' then continue
  iline = (strsplit(iline,'#',/extract)) [0]
  if strlen(iline) le 1 then continue

  ; special commands start with @
  if strmid(iline,0,1) eq '@' then begin
    vstatement = clip(strsplit(strmid(iline,1,1000), '=', /extract))
    if n_elements(vstatement) gt 1 then v = vstatement[1]
    if vstatement[0] eq 'inunit' then inunit = strmid(vstatement[1],0,1)
    if vstatement[0] eq 'expunit' then expunit = strmid(vstatement[1],0,1)
    if vstatement[0] eq 'outunit' and fixoutunit eq 0 then outunit = strmid(vstatement[1],0,1)
    unitconv=((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+ $
               (outunit eq 's')*24.*60.*60. + (outunit eq 'k')*24*3.6) / $
              ((inunit eq 'd')+(inunit eq 'h')*24.+(inunit eq 'm')*24.*60.+$
               (inunit eq 's')*24.*60.*60.+ (inunit eq 'k')*24*3.6 )
    if vstatement[0] eq 'expunit' then expunit = strmid(vstatement[1],0,1)
    expconv = ((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+$
               (outunit eq 's')*24.*60.*60. + (outunit eq 'k')*24*3.6) / $
              ((expunit eq 'd')+(expunit eq 'h')*24.+(expunit eq 'm')*24.*60.+$
               (expunit eq 's')*24.*60.*60. + (expunit eq 'k')*24*3.6 )
    if vstatement[0] eq 'fluxconv' then fluxconv = vstatement[1]
    if vstatement[0] eq 'fluxunitconv' then fluxunitconv = vstatement[1]
    if strmid(vstatement[0],0,4) eq 'filt' then defaultfilt = vstatement[1]
    if strmid(vstatement[0],0,4) eq 'band' then defaultfilt = vstatement[1]
    if        vstatement[0]    eq 'source' then defaultsrc = vstatement[1]
    if strmid(vstatement[0],0,4) eq 'inst' then defaultsrc = vstatement[1]
    if strmid(vstatement[0],0,3) eq 'ref' then  defaultref = vstatement[1]
    if strmid(vstatement[0],0,3) eq 'tel' then  defaulttel = vstatement[1]
    if strmid(vstatement[0],0,3) eq 'sym' then  defaultsym = vstatement[1]
    if strmid(vstatement[0],0,3) eq 'utb' then utburst = vstatement[1]
    if strmid(vstatement[0],0,3) eq 'ebv' then ebv = float(vstatement[1])
    if strmid(vstatement[0],0,2) eq 'av' then ebv = float(vstatement[1])/3.1
    if strmid(vstatement[0],0,4) eq 'skip' then skip = 1
    if strmid(vstatement[0],0,6) eq 'unskip' then skip = 0
    if strmid(vstatement[0],0,3) eq 'lim' then defaultlim = (v eq 'yes' or v eq 'y')-(v eq 'X' or v eq 'x')
    if strmid(vstatement[0],0,4) eq 'flag' then defaultlim = 255
    if strmid(vstatement[0],0,6) eq 'unflag' then defaultlim = 0
    continue
  endif

  if skip gt 0 then continue
   
  ; change the delimiter or the order/meaning of columns with %
  if strmid(iline,0,1) eq '%' then begin
    splitchar = strmid(iline,1,1)
    if splitchar eq ' ' then begin
      idata = strtrim(strsplit(strmid(iline,2),/extract),2)
    endif else begin
      idata = strtrim(strsplit(strmid(iline,2),splitchar,/extract),2)
    endelse
    ;idata = idata[1:*]

    nd = n_elements(idata)
    
    dsrc = -1
    dtel = -1
    dref = -1
    dband = -1
    dfilt = -1
    dt = -1
    dtstart = -1
    dtend = -1
    ddt = -1
    ddtplus = -1
    ddtminus = -1
    dexp = -1
    dmag = -1
    demag = -1
    dflux = -1
    deflux = -1
    dcts = -1
    dects = -1
    dlim = -1
    ddet = -1
    dsys = -1
    dpsym = -1
    dband = -1
    dutstart = -1
    dutmid = -1
    dutend = -1

    for d = 0, nd-1 do begin
      en = clip(idata[d])
      if strpos(en, 'src') ge 0 or strpos(en, 'sou') ge 0 then dsrc = d
      if (strpos(en, 'ref') ge 0 and strpos(en,'refsys') ne 0) then dref = d
      if strpos(en, 'tel') ge 0 then dtel = d
      if en eq 't' or en eq 'tburst' or en eq 'tmid' then dt = d
      if strpos(en, 'tst') ge 0 or strpos(en, 't_st') ge 0 then dtstart = d
      if strpos(en, 'tend') ge 0 or strpos(en, 't_end') ge 0 then dtend = d
      if strpos(en, 'dt') ge 0 then ddt = d
      if strpos(en, 'exp') ge 0 then dexp = d
      if strpos(en, 'filt') ge 0 then dfilt = d
      if strpos(en, 'band') ge 0 then dband = d
      if strmid(en,0,3) eq 'mag' then dmag = d
      if strmid(en,0,1) eq 'e' and strpos(en, 'mag') gt 0 then demag = d
      if strmid(en,0,1) eq 'm' and strpos(en, 'err') gt 0 then demag = d
      if (dmag eq d-1 or dmag eq d-2) and strpos(en, 'unc') eq 0 then demag = d
      if strpos(en, 'lim') ge 0 then dlim = d
      if strpos(en, 'det') ge 0 then ddet = d
      if strpos(en, 'sym') ge 0 then dpsym = d
      if en eq 'flux' then dflux = d
      if strmid(en,0,1) eq 'e' and strpos(en, 'flux') gt 0 then deflux = d
      if strmid(en,0,1) eq 'f' and strpos(en, 'err') gt 0 then deflux = d
      if strpos(en, 'ct') eq 0 then dcts = d
      if strpos(en, 'ect') ge 0 then dects = d
      if strpos(en, 'utstart') ge 0 then dutstart = d
      if strpos(en, 'utmid') ge 0 then dutmid = d
      if strpos(en, 'utend') ge 0 then dutend = d
    endfor

    continue
  endif

  if strmid(iline,0,1) eq '!' then begin ;special XRT file interpretation
    splitchar = string(9B)
    dt = 0
    ;ddt = 1
    ddtplus = 1
    ddtminus = 2
    dfilt = -1
    dsys = -1
    dmag = -1
    dlim = -1
    dpsym = -1
    dcts = 3
    dects = 4
    inunit = 's'
    expunit = 's'
    ;if n_elements(isanalyser) eq 0 then isanalyser = ((strpos(iline,'ANAL')) [0]) gt 0
    fluxconv = 1.33e+11
    ;if isanalyser then fluxconv *= 1e-3
    unitconv=((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+ $
               (outunit eq 's')*24.*60.*60.) / $
              ((inunit eq 'd')+(inunit eq 'h')*24.+(inunit eq 'm')*24.*60.+$
               (inunit eq 's')*24.*60.*60.)
    expconv = ((outunit eq 'd')+(outunit eq 'h')*24.+(outunit eq 'm')*24.*60.+$
               (outunit eq 's')*24.*60.*60.) / $
              ((expunit eq 'd')+(expunit eq 'h')*24.+(expunit eq 'm')*24.*60.+$
               (expunit eq 's')*24.*60.*60.)

    continue
  endif
  if strmid(iline,0,4) eq 'READ' then continue
  if strmid(iline,0,2) eq 'OR' then continue
  if strmid(iline,0,2) eq 'NO' then continue

  ; properties to read directly from file
  
    if splitchar ne ' ' then $
      idata = strtrim(strsplit(iline,splitchar,/extract),2) $
    else $
      idata = strtrim(strsplit(iline,/extract),2)

    ned = n_elements(idata)

    if dsrc ge 0 and ned gt dsrc    then obs[i].source = idata[dsrc]
       if obs[i].source eq '' then obs[i].source = defaultsrc
      if strpos(obs[i].source, 'GCN') eq 0 and $
         strpos(obs[i].source, 'GCNR') ne 0 then obs[i].source = 'GCN'
    if dref ge 0 and ned gt dref    then obs[i].reference = idata[dref]
       if obs[i].reference eq '' then obs[i].reference = defaultref
    if dtel ge 0 and ned gt dtel then obs[i].telescope = idata[dtel]
       if obs[i].telescope eq '' then obs[i].telescope = defaulttel

    if dt ge 0 and   ned gt dt      then obs[i].t = double(idata[dt])*unitconv
    if dtstart ge 0 and ned gt dtstart then $
                                obs[i].tstart = double(idata[dtstart])*unitconv
    if dtend ge 0 and ned gt dtend  then $
                                    obs[i].tend = double(idata[dtend])*unitconv
    if dexp ge 0 and ned gt dexp    then obs[i].exp=double(idata[dexp])*expconv
    if ddt ge 0 and ned gt ddt      then obs[i].exp=double(idata[ddt])*expconv
    if (ddtplus ge 0 and ddtminus ge 0 and ned gt ddtminus and ned gt ddtplus) then begin
      obs[i].exp=     ( abs(double(idata[ddtplus])) + abs(double(idata[ddtminus])) )*expconv
      obs[i].tstart = obs[i].t - abs(double(idata[ddtminus]))*expconv
      obs[i].tend   = obs[i].t + abs(double(idata[ddtplus]))*expconv
      obs[i].t = (obs[i].tstart + obs[i].tend)/2 ; may not always want to do this...?
    endif
    if dfilt ge 0 and ned gt dfilt  then obs[i].filter = idata[dfilt]
      if obs[i].filter eq 'Rc' then obs[i].filter = 'R'
      if obs[i].filter eq 'Ic' then obs[i].filter = 'I'
      if obs[i].filter eq '' then obs[i].filter = defaultfilt
    if dband ge 0 and ned gt dband  then obs[i].filter = idata[dband]
    ;print, '|'+idata[dmag]+'|'
    if dmag ge 0 and ned gt dmag    then obs[i].mag = double(idata[dmag])
    if demag ge 0 and ned gt demag  then obs[i].emag = double(idata[demag])
    if dflux ge 0 and ned gt dflux  then obs[i].flux = double(idata[dflux]) * fluxunitconv
    if deflux ge 0 and ned gt deflux then obs[i].fluxerr =double(idata[deflux]) * fluxunitconv
    if dcts ge 0 and ned gt dcts    then obs[i].cts = double(idata[dcts])
    if dects ge 0 and ned gt dects  then obs[i].ects = double(idata[dects])
      obs[i].uplim = defaultlim
    if dlim ge 0 and ned gt dlim    then obs[i].uplim = $
      (idata[dlim] eq 'yes' or idata[dlim] eq 'y' or idata[dlim] eq '1')-(idata[dlim] eq 'X' or idata[dlim] eq 'x')
    if ddet ge 0 and ned gt ddet    then obs[i].uplim = $
   (strmid(idata[dlim],0,1) eq 'n')-(idata[dlim] eq 'X' or idata[dlim] eq 'x')
    if dsys ge 0 and ned gt dsys    then obs[i].refsys = idata[dsys]
    if dpsym ge 0 and ned gt dpsym  then obs[i].psym = idata[dpsym]
      if dpsym lt 0 or dpsym ge ned then obs[i].psym = defaultsym
    if dutstart ge 0 and ned gt dutstart then obs[i].tstart = utdiff(utburst,idata[dutstart])*sec2outunit
    if dutmid ge 0 and ned gt dutmid then obs[i].t = utdiff(utburst,idata[dutmid])*sec2outunit
    if dutend ge 0 and ned gt dutend then obs[i].tend = utdiff(utburst,idata[dutend])*sec2outunit

  ; if UVOT, set to a UVOT filter
  if obs[i].source eq 'UVOT' then begin
    if obs[i].filter eq 'u' or obs[i].filter eq 'U' then obs[i].filter = 'UVU'
    if obs[i].filter eq 'b' then obs[i].filter = 'UVB'
    if obs[i].filter eq 'v' then obs[i].filter = 'UVV'
    if obs[i].filter eq 'uvm1' then obs[i].filter = 'UVM1'
    if obs[i].filter eq 'uvw1' then obs[i].filter = 'UVW1'
    if obs[i].filter eq 'uvm2' then obs[i].filter = 'UVM2'
  endif

  ; properties to calculate from other properties

    ;flux (if given as cts)
    if obs[i].flux le 0.0 and obs[i].cts gt 0.0 then begin
      obs.flux = obs.cts * fluxconv
      obs.fluxerr = obs.ects * fluxconv
    endif

    ;flux (if given as mag)
    filt = (strsplit(obs[i].filter,'_',/extract)) [0]
    k = zeropt(filt)
    if obs[i].flux le 0. and obs[i].mag gt 0. then begin 
      obs[i].flux = 10.0^(k-obs[i].mag/2.5)
      ;print, obs[i].filter, obs[i].mag, obs[i].flux
      obs[i].fluxerr = obs[i].flux - 10.0^(k-(obs[i].mag+obs[i].emag)/2.5)
    endif

    ;mag (if given as flux)
    if obs[i].mag le 0. and obs[i].flux gt 0. then begin
      obs[i].mag = 2.5*(k - alog10(obs[i].flux))
      obs[i].emag = obs[i].mag - $
            (2.5*(k - alog10(obs[i].flux+obs[i].fluxerr)))
    endif

    ;print, obs[i].filter, ' ', filt, '  ', obs[i].mag, obs[i].flux

    ;print, obs[i].tstart, obs[i].tend, obs[i].t

    ; find exposure time and midtime if not given but t1, t2 are given
    if obs[i].tstart ne -1.0 and obs[i].tend ne -1.0 then begin
      if obs[i].t eq -1.0 then $
        obs[i].t = obs[i].tstart + (obs[i].tend - obs[i].tstart)/2
      if obs[i].exp le 0.0 then $
        obs[i].exp = obs[i].tend - obs[i].tstart 
      ;note:  sometimes exp != tend-tstart due to dead time
    endif

    ; find midpoint/start time if not given
    ;   this is a bug if units are not the same!
    if obs[i].exp gt 0.0 then begin
      if obs[i].t eq -1.0 and obs[i].tstart gt 0.0 then obs[i].t = obs[i].tstart + obs[i].exp/2
      if obs[i].tstart eq -1.0 and obs[i].t gt 0.0 then obs[i].tstart = obs[i].t - obs[i].exp/2
      if obs[i].tend eq -1.0 and obs[i].t gt 0.0 then obs[i].tend = obs[i].t + obs[i].exp/2
    endif else begin  ;exposure time is unknown - assume to be zero
      if obs[i].t eq -1.0 and obs[i].tstart gt 0.0 then obs[i].t = obs[i].tstart 
      if obs[i].tstart eq -1.0 and obs[i].t gt 0.0 then obs[i].tstart = obs[i].t
      if obs[i].tend eq -1.0 and obs[i].t gt 0.0 then obs[i].tend = obs[i].t
    endelse

    ;print, obs[i].tstart, obs[i].tend, obs[i].t

    ;print, obs[i].t, obs[i].tstart, obs[i].tend, obs[i].exp

    ;print, obs[i].tstart, obs[i].tend, obs[i].t
    ;print, ''

    ;print, obs[i].filter, obs[i].tstart, obs[i].mag, obs[i].emag

  i = i + 1

endfor


close, 1

obs = obs[0:i-1]

if ebv gt 0 and keyword_set(noextcorr) eq 0 then begin
  for j = 0, n_elements(obs)-1 do begin
    fm_unred, filtwv(obs[j].filter),1.0,ebv,extfactor
    obs[j].flux = obs[j].flux * extfactor
    obs[j].fluxerr = obs[j].fluxerr * extfactor
  endfor
endif

return, obs

end