function is_derived_field_number, name, dfields
  is = -1
  for j=0, N_elements(dfields)-1 do begin
     if STRMATCH((*dfields[j])[0], name) then is = j
  endfor
return, is
end

pro pack_as_list, use_g, a
; pack all data into 1D arrays 
   zero = -1E30
   zero_solution_under_subgrid, use_g, a, ZERO=zero

  nvar = a[0].n_dfields
  npoints = total(a.np_total)
  da = dblarr(nvar,npoints)
  for j=0,Nvar-1 do begin
      oruni = 0L
      for i=0L,N_elements(a)-1 do begin
          ind = where(*a[i].data[j] ne zero, count)
          runi = oruni+count
          if count gt 0 then da[j,oruni:runi-1] = (*a[i].data[j])[ind]
          oruni = runi
      endfor
  endfor
  for j=0,Nvar-1 do a[0].data[j] = ptr_new(da[j,0:runi-1])
  a[0].np_total = runi
  for i=1L,N_elements(a)-1 do begin
      for j=0,Nvar do ptr_free, a[i].data[j]
  endfor
  a = a[0]
end


function identify_unique_field_names, VariableNames, dfields
  Nvar = N_elements(VariableNames)
  Nderived = N_elements(dfields)
  field_names = ['']
  match = -1
  for i=0,Nvar-1 do begin
     match = is_derived_field_number(VariableNames[i], dfields)
     if match ge 0 then begin
         tfields = (*dfields[match])[1:*]
         for j=0,N_elements(tfields)-1 do begin
             matchagain = is_derived_field_number(tfields[j], dfields)
             if matchagain ge 0 then $
               field_names = [field_names, (*dfields[matchagain])[1:*]]  $
             else field_names = [field_names, tfields[j]]
         endfor
     endif else field_names = [field_names, VariableNames[i]]
  endfor
  field_names = field_names[1:*]
  nfn = N_elements(field_names)
; remove duplicate entries
  field_names = remove_additional_occurance(field_names)
  return,field_names
end

function determine_data_format_from_file_name, name
  data_format = 5
  if strmatch(name, '*.grid*') then data_format = 4    ; hdf4 
  if strmatch(name, '*.cpu*') then data_format = 5     ; packed format
  if strmatch(name, '*_P*.hdf5*') then data_format = 6  ; movie format
  if strmatch(name, '*.speck*') then data_format = 1  ; partiview particle format
  if strmatch(name, '*snapshot_*') then data_format = 2  ; gadget particle format
  if (data_format eq 4) then begin
     spawn, 'h5ls '+name, res
print, name
     if (N_elements(res) gt 5) then data_format = 3
  endif

  if (data_format le 0) then begin
     print, 'determine_data_format_from_file_name: unrecognized data format! ', name
  endif

print, 'data_format:',data_format     
  return, data_format
end

function radial_velocity, cen, x,y,z,vx,vy,vz, cons_vel=cons_vel
x = x-cen[0]
y = y-cen[1]
z = z-cen[2]
r=sqrt(x^2+y^2+z^2)
if not keyword_set(cons_vel) then begin
    ind = (where(r eq min(r)))[0] ; subtract central velocity
    vc = [vx[ind], vy[ind], vz[ind]]
    if 1 then print, 'subtract central velocity:', vc
endif else vc =[0.,0.,0.]

vr = ((vx-vc[0])*x+(vy-vc[1])*y+(vz-vc[2])*z)/r
return,vr
end
function tangential_velocity, cen, x,y,z,vx,vy,vz
x = x-cen[0]
y = y-cen[1]
z = z-cen[2]
r=sqrt(x^2+y^2+z^2)
vt = sqrt( (vy*z-vz*y)^2+(vz*x-vx*z)^2+(vx*y-vy*x)^2 )/r
return,vt
end

function add_dim, a
  dim = (size(a,/dimensions))
  if N_elements(dim) eq 1 then a = reform(a,dim[0],1,1)
  if N_elements(dim) eq 2 then a = reform(a,dim[0],dim[1],1)
  return, a
end

; very crude gradients avoiding using ghostzones
function div_x, vx, dx
  dim = (size(vx,/dimensions))
  if N_elements(dim) ne 3 then vx = add_dim(vx)
  dim = dim[0]
  div_x = (shift(vx,[-1,0,0])-shift(vx,[1,0,0]))/2./dx
; forward and backward differences where we do not have ghostzones
  if dim gt 1 then  div_x[0,*,*] = (vx[1,*,*] - vx[0,*,*])/dx
  if dim gt 1 then  div_x[dim-1,*,*] = (vx[dim-1,*,*] - vx[dim-2,*,*])/dx
  return, div_x
end
function div_y, vy, dy
  dim = (size(vy,/dimensions))
  if N_elements(dim) ne 3 then vy = add_dim(vy)
  dim = dim[1]
  div_y = (shift(vy,[0,-1,0])-shift(vy,[0,1,0]))/2./dy
  if dim gt 1 then  div_y[*,0,*] = (vy[*,1,*] - vy[*,0,*])/dy
  if dim gt 1 then  div_y[*,dim-1,*] = (vy[*,dim-1,*] - vy[*,dim-2,*])/dy
  return, div_y
end
function div_z, vz, dz
  dim = (size(vz,/dimensions))
  if N_elements(dim) ne 3 then begin
     vz = add_dim(vz)
     dim = (size(vz,/dimensions))
  endif
  dim = dim[2]
  div_z = (shift(vz,[0,0,-1])-shift(vz,[0,0,1]))/2./dz
  if dim gt 1 then  div_z[*,*,0] = (vz[*,*,1] - vz[*,*,0])/dz
  if dim gt 1 then  div_z[*,*,dim-1] = (vz[*,*,dim-1] - vz[*,*,dim-2])/dz
  return, div_z
end
function div_v, vx,vy,vz, dx,dy,dz
  div = div_x(vx,dx)+div_y(vy,dy)+div_z(vz,dz)
return, div
end


function abs_div_v, vx,vy,vz, dx,dy,dz
  div = sqrt(div_x(vx,dx)^2+div_y(vy,dy)^2+div_z(vz,dz)^2)
return, div
end

function enclosed_gas_mass, cen,x,y,z, mass
x = x-cen[0]
y = y-cen[1]
z = z-cen[2]
r=sqrt(x^2+y^2+z^2)
ind = sort(r)
Menc = (TOTAL((mass)[ind],/double,/cumulative))[sort(ind)]
return,Menc
end

function rho_slope, rho, x, y, z, dx, dy, dz, center=c
if N_elements(c) ne 3 then c = [0,0,0]
; gradient
  gx = div_x(rho, dx)
  gy = div_y(rho, dy)
  gz = div_z(rho, dz)
;radius
;  r = sqrt((x-c[0])^2+(y-c[1])^2+(z-c[2])^2)
; slope
  alpha = (gx*(x-c[0])+gy*(y-c[1])+gz*(z-c[2]))/rho
  return, alpha
end


function fn, aa, name
@common_blocks.inc
ind = where(strtrim(aa[0].data_fields,2) eq strtrim(name,2))
if ind[0] eq -1 then begin
   pind = where(name eq *parti.names)
   if pind[0] ne -1 then return, parti.d[pind]
   print, 'no field with name >'+name+'< found in ', aa
   return, 0
endif

return, (*aa[0].data[ind[0]])
end

function get_data, gi, VariableNames, viewdim, DERIVED_FIELDS=derived_fields,$
                   LEFT_INDEX=start, RIGHT_INDEX=righti, redshift=zred, $
                   LEFT_EDGE=ledge, RIGHT_EDGE=redge, interpolate=interp, $
                   return_as_list=ral
@common_blocks.inc
; for every new derived field give name and the original fields they
; are derived from.  
; viewdim is the view direction. 
; to add your own: 
;     1)  increment the 'Nderived =' statement
;     2)  define the name and the fields it gets derived from in an
;         array of strings
;     3) look to the end of this file to add the function
; if you set  /Derived_Fields this routines just returns the
; names of the derived fields defined here [enables to only have one
; place where they are mentioned]
; vvvvvvvv modify below vvvvvvvvvvvvv

Nderived = 50
id = 0L
dfields = PTRARR(Nderived) 
dfields[id++] = PTR_NEW(["Density^2", "Density"])
dfields[id++] = PTR_NEW(["Balmer emission", "HII_Density", "Electron_Density", "Temperature"])
dfields[id++] = PTR_NEW(["Approximate Balmer", "Density", "Temperature"])
dfields[id++] = PTR_NEW(["abs(B)", "Bx", "By", "Bz"])
dfields[id++] = PTR_NEW(["H2 Fraction", "Density", "H2I_Density"])
dfields[id++] = PTR_NEW(["HI Fraction", "Density", "HI_Density"])
dfields[id++] = PTR_NEW(["HII Fraction", "Density", "HII_Density"])
dfields[id++] = PTR_NEW(["Electron Fraction", "Density", "Electron_Density"])
dfields[id++] = PTR_NEW(["Metallicity (Solar)", "Density", "SN_Colour"])
dfields[id++] = PTR_NEW(["Metallicity (solar)", "Density", "Metal_Density"])
dfields[id++] = PTR_NEW(["Relativistic Gamma", "x-velocity", "y-velocity", "z-velocity"])
dfields[id++] = PTR_NEW(["Plasma beta", "Density", "GasEnergy", "Bx", "By", "Bz"])
dfields[id++] = PTR_NEW(["Alfven speed", "Density",  "Bx", "By", "Bz"])
dfields[id++] = PTR_NEW(["rms 3D velocity", "x-velocity", "y-velocity", "z-velocity"])
dfields[id++] = PTR_NEW(["Tb (kinetic)", "HI_Density", "Temperature"])
dfields[id++] = PTR_NEW(["Tb (fully coupled)", "HI_Density", "Temperature"])
dfields[id++] = PTR_NEW(["Rotation Measure", "B||", "Electron_Density"])
dfields[id++] = PTR_NEW(["Approx. Rotation Meas.", "B||", "Density"])
dfields[id++] = PTR_NEW(["Dynamical time (gas)", "Density"])
dfields[id++] = PTR_NEW(["Radius", "x", "y", "z"])
dfields[id++] = PTR_NEW(["Cylindrical Radius", "x", "y", "z"])
dfields[id++] = PTR_NEW(["radial velocity", "x-velocity", "y-velocity", "z-velocity", "x", "y", "z"])
dfields[id++] = PTR_NEW(["tangential velocity", "x-velocity", "y-velocity", "z-velocity", "x", "y", "z"])
dfields[id++] = PTR_NEW(["enclosed gas mass", "mass", "x","y","z"])
dfields[id++] = PTR_NEW(["H number density", "HI_Density", "HII_Density"])
dfields[id++] = PTR_NEW(["Pressure", "Density", "Gas_Energy"])
dfields[id++] = PTR_NEW(["Gas Surface Density", "Density", "x", "y", "z"])
dfields[id++] = PTR_NEW(["radial dM/dt", "Density", "x-velocity", "y-velocity", "z-velocity", "x", "y", "z"])
dfields[id++] = PTR_NEW(["Sound speed", "Temperature"])
dfields[id++] = PTR_NEW(["Jeans mass", "Density", "Temperature"])
dfields[id++] = PTR_NEW(["|v_particles|^10", "Particle_x-velocity", $
                  "Particle_y-velocity", "Particle_z-velocity"])
dfields[id++] = PTR_NEW(["abs(ParticleVelocity)", "Particle_x-velocity", $
                  "Particle_y-velocity", "Particle_z-velocity"])
dfields[id++] = PTR_NEW(["DarkMatter_mass", "volume", "Dark_Matter_Density"])
dfields[id++] = PTR_NEW(["Current", "B_Vorticity1", "B_Vorticity2", "B_Vorticity3"])
dfields[id++] = PTR_NEW(["Pressure (from Etot)", "TotalEnergy", "Density", "x-velocity", "y-velocity", "z-velocity"])
dfields[id++] = PTR_NEW(["Entropy", "Density", "Temperature"])
dfields[id++] = PTR_NEW(["div V", "x-velocity", "y-velocity", "z-velocity", "dx", "dy", "dz"])
dfields[id++] = PTR_NEW(["density slope", "Density", "x", "y", "z", "dx", "dy", "dz"])
dfields[id++] = PTR_NEW(["DM density slope", "Dark_Matter_Density", "x", "y", "z", "dx", "dy", "dz"])
dfields[id++] = PTR_NEW(["C^18 O (2-1)", "Density"])
dfields[id++] = PTR_NEW(["N_2 H^+ (1-0)", "Density"])
dfields[id++] = PTR_NEW(["N H_3 (1-1)", "Density"])
dfields[id++] = PTR_NEW(["Gas Fraction", "Density", "Dark_Matter_Density"])

dfields = dfields[0:id-1]

; ^^^^^^^^ modify above ^^^^^^^^^^
bs = ["Bx", "By", "Bz"]
rcd = (slice_ori eq 0.) ; is 0 for the slice coordinate 

if KEYWORD_SET(derived_fields) then return, dfields ; return field_names & stop

if N_elements(interp) eq 0 then interp = 0
num_of_grids = N_ELEMENTS(gi) 
   index_range = gi.End_index-gi.Start_index+1
   count = transpose(index_range)          ; read all active zones as default
;   if num_of_grids eq 1 then count = transpose(count)
   if N_elements(start) eq 0 then start = INTARR(num_of_grids, 3)
;   start = count
;   start[*] = 0

IF N_elements(ledge) eq 3 then begin ; compute start and count from left_edge and right_edge
   if verbose then print, 'get_data: compute indices from left and right_edge'
   delta_distance = (gi.Right_edge - gi.Left_edge)/FLOAT(index_range)
   start = intarr(num_of_grids,3) & stride = intarr(num_of_grids,3) & righti = intarr(num_of_grids,3)
   for i=0,2 do begin
      start[*,i]  =(FLOOR((ledge[i] - gi.Left_edge[i])/delta_distance[i,*]))  > 0 < (index_range[i,*]-1)
      righti[*,i] =(FLOOR((redge[i] - gi.Left_edge[i])/delta_distance[i,*]))  > 0 < (index_range[i,*]-1)
  endfor

endif

if N_ELEMENTS(righti) gt 2 then begin
    count = righti-start+1
endif


if N_elements(count) gt 0 then count = count > 1
Nvar = N_ELEMENTS(VariableNames)
data_format = determine_data_format_from_file_name(gi[0].baryon_file)
if data_format lt 0 then return,0
;if data_format eq 6 then dataisvertexc = 1 else dataisvertexc = 0 ;
;assume movie data is vertex centered
dataisvertexc = 0
if dataisvertexc then count += 1
; determine all data fields to read and read data
field_names = identify_unique_field_names(VariableNames, dfields)

ind = where(field_names eq 'B||', vcount)
if vcount gt 0 then field_names[ind] = bs[viewdim]


if field_names[0] ne 'Dummy' then begin
      if data_format eq 4 then read_all_grid_hdf,  $
         gi, field_names, a, VERBOSE=verbose, START=start, COUNT=count
      if ((data_format eq 5) OR (data_format eq 6)) then $
         read_all_grid_hdf5, gi, field_names, a, $
                             /VERBOSE, START=start, COUNT=count
      if (data_format eq 3) then $
         read_all_grid_hdf5, gi, field_names, a, $
                             /VERBOSE, START=start, COUNT=count

   endif else begin ; create dummy field 
        readall        = {readall,                         $
                    file_name: "",                   $
                    Np_total: 0L,                    $
                    dims: [0,0,0],                   $
                    li: [0,0,0],                     $  ; left index
                    ri: [0,0,0],                     $  ; right index
                    N_dfields: 0B,                   $
                    DATA_FIELDS: STRARR(100),        $
                    Data:        PTRARR(100)   }
        a = replicate(readall, N_elements(gi))
        for j=0L,NUM_OF_GRIDS-1 do begin 
           ra = gi[j].end_index-gi[j].start_index + 1
           tmp = fltarr(ra[0], ra[1], ra[2])
           if verbose then print, 'debug:', ra, N_elements(tmp)
           tmp[*] = 1.
           a[j].li = [0, 0, 0]
           a[j].ri = ra[*]-1
           a[j].data[0] = ptr_new(tmp)
           a[j].dims[*] = ra[*]
           a[j].np_total = N_elements(tmp)
        endfor
        a[*].n_dfields = 1
        return, a
     endelse

   dataisvertexc = 0
interp=0
   if interp and not dataisvertexc then begin ; map to vertex centered data
      for j=0L,NUM_OF_GRIDS-1 do begin 
         for i=0,Nvar-1 do begin
            d = a[j].dims ; pad data by one cell each direction
            data = fltarr(d[0]+1,d[1]+1,d[2]+1)
            if min(d) gt 2 then $
               data[*,*,*] = interpolate((*a[j].data[i])[*,*,*], $ ; trilinear interpolation
                                         [0,findgen(d[0])+0.5], [0,findgen(d[1])+0.5], [0,findgen(d[2])+0.5], /grid) $
            
            else data[0:(d[0]-1),0:(d[1]-1),0:(d[2]-1)]=(*a[j].data[i])[*,*,*]
            a[j].dims= size(data, /DIMENSIONS)
            a[j].ri  = a[j].li + a[j].dims-1
            a[j].np_total = N_elements(data)
            ptr_free, a[j].data[i]
            a[j].data[i] = ptr_new(data)
         endfor
      endfor
      ; now the first and last row/columns can be set from parent
;      if NUM_OF_GRIDS gt 1 then np = histogram(gi.parent, reverse_indices=r,min=0)
      for j=0L,NUM_OF_GRIDS-1 do begin 
         pin = where(gi[j].parent eq gi.num, count)
         d = a[j].dims
         for i=0,Nvar-1 do begin
            if count eq 1 then begin
               li = (((gi[j].parent_s_index[*]) - a[pin].li[*]) > 0 ) < (a[pin].dims -1)
               ri = li
               for k=0L,N_elements(li)-1 do  $
                  ri[k] = ((min([gi[j].parent_e_index[k]+1, a[pin].ri[k]])-a[pin].li[k]) $
                           < (a[pin].dims[k])-1) > (li[k] + 1)
               print, 'li,ri:', li, ri
               print, 'd', d
               print, 'a[pin].dims:',a[pin].dims
               print, 'a[pin].li/ri:',a[pin].li, a[pin].ri
               if min(d) gt 2  then begin
; should be doing something like the following ... however it does not
; work this specific way ...
                   
    (*a[j].data[i])[0,*,*] =congrid((*a[pin].data[i])[li[0],li[1]:ri[1],li[2]:ri[2]],1, d[1], d[2],/minus_one)
    (*a[j].data[i])[*,0,*] =congrid((*a[pin].data[i])[li[0]:ri[0],li[1],li[2]:ri[2]],d[0], 1, d[2],/minus_one)
    (*a[j].data[i])[*,*,0] =congrid((*a[pin].data[i])[li[0]:ri[0],li[1]:ri[1],li[2]],d[0], d[1], 1,/minus_one)
;    (*a[j].data[i])[0,*,*] = 15.
;    (*a[j].data[i])[*,0,*] = 15.
;    (*a[j].data[i])[*,*,0] = 15.

; only the following is what seems to be wrong ... track this down and
; interpolation should work ... 
;    (*a[j].data[i])[d[0]-1,*,*] =congrid((*a[pin].data[i])[ri[0],li[1]:ri[1],li[2]:ri[2]],1, d[1], d[2]);,/minus_one)
;    (*a[j].data[i])[*,d[1]-1,*] =congrid((*a[pin].data[i])[li[0]:ri[0],ri[1],li[2]:ri[2]],d[0], 1, d[2]);,/minus_one)
;    (*a[j].data[i])[*,*,d[2]-1] =congrid((*a[pin].data[i])[li[0]:ri[0],li[1]:ri[1],ri[2]],d[0], d[1], 1);,/minus_one)
               endif
            endif
            if count eq 0 then begin ; no parent found, just extend flat
               (*a[j].data[i])[0,*,*] = (*a[j].data[i])[1,*,*]
               (*a[j].data[i])[*,0,*] = (*a[j].data[i])[*,1,*]
               (*a[j].data[i])[*,*,0] = (*a[j].data[i])[*,*,1]
               (*a[j].data[i])[d[0]-1,*,*] = (*a[j].data[i])[d[0]-2,*,*]
               (*a[j].data[i])[*,d[1]-1,*] = (*a[j].data[i])[*,d[1]-2,*]
               (*a[j].data[i])[*,*,d[2]-1] = (*a[j].data[i])[*,*,d[2]-2]
            endif 
         endfor
     endfor
   endif


;      breaki
; Sometimes we need the redshift
   z = 0.0
   if N_elements(redshifts) eq N_elements(times) then $
      if redshifts[time_index] ne -1 then z=redshifts[time_index]

   ini = [0]
   ins = [0]
   if verbose then print, field_names
   for j=0L,N_elements(a)-1 DO BEGIN
      a[j].data = shift(a[j].data, nvar)
      a[j].data_fields = shift(a[j].data_fields, nvar)
      if a[j].np_total gt 0 then begin
         for i=0,Nvar-1 do begin ; allocate memory
            d = (size(*a[j].data[Nvar]))
            tmp = fltarr(d[1],d[2],d[3])
            a[j].data[i] = ptr_new(tmp)
         endfor
      endif
   endfor

   for i=0, Nvar-1 do begin     ; fill in fields and data
      ii = N_elements(a[0].data)
      match = is_derived_field_number(VariableNames[i], dfields)
      if match ge 0 then begin
         for j=0L,N_elements(a)-1 DO BEGIN
            case match of
; vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv defined derived functions here vvvvvvvvvvvvvvvvvvvvvvvvv
;          0: *a[j].data[i] = fn(a[j],"HI_Density")*fn(a[j],"Electron_Density") / $
;                             fn(a[j],"Temperature")^0.9 
;          1: *a[j].data[i] =
;          fn(a[j],"Density")^2/(fn(a[j],"Temperature")^0.9 )
          0: *a[j].data[i] = fn(a[j],"Density")^2
          1: begin
            urho = get_unit("Electron_Density")
            *a[j].data[i] = balmer_emission(fn(a[j], "Electron_Density")*urho, $
                                            fn(a[j], "Temperature")) * $
                            fn(a[j], "HII_Density") * fn(a[j], "Electron_Density") * $
                            urho^2
          end
          2: begin ;; Assume everything below 0.8 eV is neutral (x_e = 1e-6)
            urho = get_unit("Density")
            ;; temperature in eV
            neutral = WHERE(fn(a[j], "Temperature")/11605.0 lt 0.8, nneutral)
            nelec = fn(a[j], "Density")
            if (nneutral gt 0) then nelec[neutral] *= 1e-6
            *a[j].data[i] = balmer_emission(nelec*urho, fn(a[j], "Temperature")) * $
                            (nelec*urho)^2
          end
          3: *a[j].data[i] = sqrt(fn(a[j],"Bx")^2+fn(a[j],"By")^2+fn(a[j],"Bz")^2)
          4: *a[j].data[i] = fn(a[j],"H2I_Density")/fn(a[j],"Density")
          5: *a[j].data[i] = fn(a[j],"HI_Density")/fn(a[j],"Density")
          6: *a[j].data[i] = fn(a[j],"HII_Density")/fn(a[j],"Density")
          7: *a[j].data[i] = fn(a[j], "Electron_Density")/fn(a[j],"Density")
          8: *a[j].data[i] = fn(a[j], "SN_Colour")/fn(a[j],"Density")/0.0204
          9: *a[j].data[i] = fn(a[j], "Metal_Density")/fn(a[j],"Density")/0.0204
         10: *a[j].data[i] = 1./(1-sqrt(fn(a[j], "x-velocity")^2 + $
                             fn(a[j], "y-velocity")^2 + fn(a[j], "z-velocity")^2 ))
         11: *a[j].data[i] = fn(a[j],"Density")*fn(a[j],"GasEnergy")*(gamma-1.) / $
                             (fn(a[j],"Bx")^2+fn(a[j],"By")^2+fn(a[j],"Bz")^2)/2
         12: *a[j].data[i] = ((fn(a[j],"Density")) / $
                             (fn(a[j],"Bx")^2+fn(a[j],"By")^2+fn(a[j],"Bz")^2))^(-0.5)
         13: *a[j].data[i] = sqrt( $
                 fn(a[j], "x-velocity")^2 + fn(a[j], "y-velocity")^2 + fn(a[j], "z-velocity")^2 )
         14: begin
            urho = get_unit("HI_Density")
            *a[j].data[i] = delta_tb(fn(a[j], "HI_Density")*urho, $
                                     fn(a[j], "Temperature"), z)
         end
         15: begin
            urho = get_unit("HI_Density")
            *a[j].data[i] = delta_tb(fn(a[j], "HI_Density")*urho, $
                                     fn(a[j], "Temperature"), $
                                     z, /FULLYCOUPLED)
         end
         16: *a[j].data[i] = fn(a[j], field_names[0]) * fn(a[j], "Electron_Density")
         17: *a[j].data[i] = fn(a[j], field_names[0]) * fn(a[j], "Density")
         18: *a[j].data[i] = sqrt(!DPI*3./(32.*fn(a[j], "Density")))
         19: *a[j].data[i] = sqrt((center[0]-fn(a[j], "x"))^2+(center[1]-fn(a[j], "y"))^2 $  ; radius
                                  +(center[2]-fn(a[j], "z"))^2)
         20: *a[j].data[i] = sqrt(rcd[0]*(center[0]-fn(a[j], "x"))^2 + $
                                  rcd[1]*(center[1]-fn(a[j], "y"))^2 + $ 
                                  rcd[2]*(center[2]-fn(a[j], "z"))^2)    ; cylindrical radius
         21: *a[j].data[i] = radial_velocity(center, fn(a[j], "x"),fn(a[j], "y"),fn(a[j], "z"),fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"), cons_vel=[0.,0.,0.])
         22: *a[j].data[i] = tangential_velocity(center, fn(a[j], "x"),fn(a[j], "y"),fn(a[j], "z"),fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"))
;         22: *a[j].data[i] = abs_div_v(fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"), $
;                                   fn(a[j], "dx"),fn(a[j], "dy"),fn(a[j], "dz"))
         23: *a[j].data[i] = enclosed_gas_mass(center,  fn(a[j], "x"),fn(a[j], "y"),fn(a[j], "z"),fn(a[j], "mass"))
         24: *a[j].data[i] = (fn(a[j], "HI_Density") + fn(a[j], "HII_Density")) 
         25: *a[j].data[i] = fn(a[j],"Density")*fn(a[j],"Gas_Energy")*(gamma-1.)
;         26: *a[j].data[i] = vector-divergence(fn(a[j], "x-velocity"),
;                             fn(a[j], "y-velocity"), fn(a[j], "z-velocity"))
         27: *a[j].data[i] = fn(a[j],"Density")*4.*!Pi*$
                radial_velocity(center, fn(a[j], "x"),fn(a[j], "y"),$ ; rho 4Pi vr
     fn(a[j], "z"),fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"),cons_vel=[0,0,0]) * $
     (center[0]-fn(a[j], "x"))^2+(center[1]-fn(a[j], "y"))^2+(center[2]-fn(a[j], "z"))^2 ;R^2

         28: begin
            uvel = get_unit("Velocity", /force_use_units)
            *a[j].data[i] = 1.1725e-1 * sqrt(fn(a[j], "Temperature")) / uvel
         end
         29: begin
            umass = get_unit("mass", /force_use_units)
            *a[j].data[i] = 19.0315 * (fn(a[j], "Temperature"))^1.5 / $
                            sqrt(fn(a[j], "Density")) / umass
         end

         30: *a[j].data[i] = (fn(a[j],"Particle_x-velocity")^2+fn(a[j],"Particle_y-velocity")^2 $
                              +fn(a[j],"Particle_z-velocity")^2)^5
         31: *a[j].data[i] = sqrt(fn(a[j],"Particle_x-velocity")^2+fn(a[j],"Particle_y-velocity")^2 $
                              +fn(a[j],"Particle_z-velocity")^2)
         32: *a[j].data[i] = fn(a[j],"Dark_Matter_Density")*fn(a[j],"volume")
         33: *a[j].data[i] = sqrt(fn(a[j],"B_Vorticity1")^2+fn(a[j],"B_Vorticity2")^2 $
                                  +fn(a[j],"B_Vorticity3")^2)
         34: *a[j].data[i] =  1.D*(Gamma-1.)*(fn(a[j], "TotalEnergy") - 0.5D*($
                              DOUBLE(fn(a[j], "x-velocity"))^2 + $
                              DOUBLE(fn(a[j], "y-velocity"))^2 + $
                              DOUBLE(fn(a[j], "z-velocity")^2 )))*fn(a[j], "Density")
         35: *a[j].data[i] =  ALOG(1.D*(0.58/0.6*fn(a[j], "Temperature")) /(fn(a[j], "Density")*6.9416d10)^(gamma-1.)) ; entropy in funny units ... assuming was mu=0.6 in enzo sim

         36: *a[j].data[i] = div_v(fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"), $
                                   fn(a[j], "dx"),fn(a[j], "dy"),fn(a[j], "dz"))
          
         37: *a[j].data[i] = rho_slope(fn(a[j], "Density"),fn(a[j], "x"), $
                                       fn(a[j], "y"), fn(a[j], "z"), $
                                       fn(a[j], "dx"),fn(a[j], "dy"),fn(a[j], "dz"), $
                                      center=center)
         38: *a[j].data[i] = rho_slope(fn(a[j], "Dark_Matter_Density"),fn(a[j], "x"), $
                                       fn(a[j], "y"), fn(a[j], "z"), $
                                       fn(a[j], "dx"),fn(a[j], "dy"),fn(a[j], "dz"), $
                                      center=center)
          

         39: *a[j].data[i] = molecule_thin_emission(fn(a[j], "Density")*$
                                                  get_unit("Density", /force_use_units), 0)
         40: *a[j].data[i] = molecule_thin_emission(fn(a[j], "Density")*$
                                                  get_unit("Density", /force_use_units), 1)
         41: *a[j].data[i] = molecule_thin_emission(fn(a[j], "Density")*$
                                                  get_unit("Density", /force_use_units), 2)

;         39: begin
;            dmrho = fn(a[j], "Dark_Matter_Density")
;            zero = WHERE(dmrho eq 0.0, nzero)
;            if (nzero gt 0) then dmrho[zero] = 1e20  ; force the fraction to zero
;            *a[j].data[i] = fn(a[j], "Density") / (fn(a[j], "Density")+dmrho)
;         end



; ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ define derived functions above ^^^^^^^^^^^^^^^^^^^^^^^^^
      endcase
            a[j].data_fields[i] = (*dfields[match])[0]
         endfor                 ; loop over grids
      endif else begin              ; is derived variable
         for j=0L,N_elements(a)-1 DO  BEGIN
            if a[j].np_total gt 0 then *a[j].data[i] = fn(a[j], VariableNames[i])
            a[j].data_fields[i] = VariableNames[i]
         endfor         
      endelse
   endfor                       ; loop over variables

     if keyword_set(ral) then pack_as_list, gi, a

   for j=0L,N_elements(a)-1 DO BEGIN
      heap_free, a[j].data[Nvar:*] 
      a[j].data_fields[Nvar:*] = ''
   endfor
   a.n_dfields = Nvar

return, a
end