function get_variable_names, gi, OLD_FIELD=old_field, var_index=var_index, particle_names=pnames
  dfields = get_data(0,0,0,DERIVED_FIELDS=1)
  pdfields = get_particle_data(0,0,DERIVED_FIELDS=1)

  if N_elements(gi) eq 0 then return, 'none'
  if max(gi.num_particle) gt 0 then fname = gi[(where(gi.num_particle gt 0))[0]].baryon_file $
  else fname = gi[0].baryon_file

  list = list_field_names(fname)
  data_format = determine_data_format_from_file_name(fname)
  if data_format eq 1 then begin
     pnames = list
     return, ['none']
  endif 
  if data_format eq 2 then begin ; gadget data
     pnames = list
     return, list               ; also allow in variables to allow interpolated
                                ; images of gadget data 
  endif 


  list = [list[*],  'x', 'y', 'z', 'dx', 'dy', 'dz', 'volume','mass']
  pnames = ''
  ind = where(strmatch(list, 'particle*') or strmatch(list, 'dynamical_time') $
              or strmatch(list, 'creation_time') $
              or strmatch(list, 'metallicity_fraction') $
             )

  print, 'getVAR', ind
  if ind[0] ne -1 then begin 
     pnames = list[ind]
     ii = bytarr(N_elements(list))
     ii[ind] = 1
     opposite = where(ii eq 0)
     if opposite[0] eq -1 then list = "No baryon data" else list = list[opposite]
     ind = sort(pnames)
     rev = sort(ind)
     pnames = (pnames[uniq(pnames, ind)])[rev]
; now make sure positions are first:
     partpos = where(strmatch(pnames, 'particle_position*'))
     for i=0, N_elements(partpos)-1 do begin
        he = pnames[i]
        pnames[i] = pnames[partpos[i]]
        pnames[partpos[i]] = he
     endfor
     if total(strmatch(pnames, '*mass')) and where(pnames eq 'particle massdensity') eq -1 then pnames=[pnames[*], 'particle massdensity']
  endif 

; select the derived fields that can be computed with the list of available data 
  df = ['']
  dfp = df
  for i=0,N_elements(dfields)-1 do begin
     found = 0
     for j=1,N_elements(*dfields[i])-1 do BEGIN
        if (where((*dfields[i])[j] eq 'Radius'))[0] ge 0 then found += 1
        for k=0,N_elements(list)-1 DO BEGIN
           if (list[k] eq (*dfields[i])[j]) then found += 1
;           print, (*dfields[i])[j], (where((*dfields[i])[j] eq ['x', 'y', 'z','volume','mass']))[0]
; to see whether B|| exist just check that there is at least one B component
           if (((*dfields[i])[j] eq 'B||') and list[k] eq "Bx") then found += 1
        endfor
     ENDFOR
     if found eq j-1 then df = [df, (*dfields[i])[0]]     
  endfor
; now the same for particles
  for i=0,N_elements(pdfields)-1 do begin
     found = 0
     for j=1,N_elements(*pdfields[i])-1 do $
        for k=0,N_elements(pnames)-1 DO BEGIN
           if (pnames[k] eq (*pdfields[i])[j]) then found += 1
        endfor
     if found eq j-1 then dfp = [dfp, (*pdfields[i])[0]]     
  endfor
  if N_elements(df) gt 1 then list =  [list, df[1:*]]
  if N_elements(dfp) gt 1 then pnames =  [pnames, dfp[1:*]]
  print, 'available variables:', list, pnames
  if KEYWORD_SET(oldfield) then begin
     ind = where(list eq old_field, count)
     var_index = (count gt 0) ? ind[0] : 0
  endif
  return, list
end

function get_all_variable_names, flist, $
                                 OLD_FIELD=old_field, var_index=var_index, $
                                 particle_names=pnames

  nf = N_elements(flist)
  for i=0,nf-1 do begin
     list = get_variable_names(flist[i],particle_names=plist)
     list = remove_additional_occurance(flist)
;     plist = remove_additional_occurance(plist)
  endfor
  if KEYWORD_SET(oldfield) then begin
     ind = where(list eq old_field, count)
     var_index = (count gt 0) ? ind[0] : 0
  endif
return, list
end