function render_image
@common_blocks.inc
   interpolate =  interpolate_i
; interpolate all values
  logit = shouldbe_logged(list_str[var_index])

; get triangles:
   secmin = min(z_stuetz, ind, max=zmax)
   IF verbose THEN  print, 'min, max before griding:', secmin, zmax
;   print, x_stuetz[ind], y_stuetz[ind], dx_stuetz[ind]
;   z_stuetz[ind] = 10.
;  print, 'b:',b
   result = 0.
   result = DBLARR(xy_sl_size,xy_sl_size)
   result[*,*] = secmin

   limits =  FLOAT([0., 0., 1., 1.])
   os = where(slice_ori lt DomainLeftEdge *1.0001)
   
   IF (interpolate eq 1) THEN BEGIN
      gs =  FLOAT([1./FLOAT(xy_sl_size), 1./FLOAT(xy_sl_size)])
      tx_stuetz = FLOAT((x_stuetz-slice_size(0))/(slice_size(2)-slice_size(0))) > 0 < 1
      ty_stuetz = FLOAT((y_stuetz-slice_size(1))/(slice_size(3)-slice_size(1))) > 0 < 1
      minx = min(tx_stuetz, max=maxx)
      miny = min(ty_stuetz, max=maxy)
 ;     tx_stuetz = (tx_stuetz-minx)/(maxx-minx)
 ;     ty_stuetz = (ty_stuetz-miny)/(maxy-miny)
 ;     print, 'minx,miny,maxx,maxy:', minx,miny,maxx,maxy
;      ind =  where((tx_stuetz GE 0 ) AND (tx_stuetz LE 1 ) AND $
;                   (ty_stuetz GE 0 ) AND (ty_stuetz LE 1 ))
;      ind =  where((tx_stuetz GE -0.05 ) AND (tx_stuetz LE 1.05 ) AND $
;                   (ty_stuetz GE -0.05 ) AND (ty_stuetz LE 1.05 ))
     
      tx_stuetz =    tx_stuetz[*] > 0 < 1
      ty_stuetz =    ty_stuetz[*] > 0 < 1
      tz_stuetz  =    z_stuetz[*]

      TRIANGULATE, tx_stuetz, ty_stuetz, triangles, b
   ENDIF
   if n_elements(result) eq 0 then result = fltarr(xy_sl_size,xy_sl_size)
   if (interpolate eq 1) then result = $
;    TRIGRID(tx_stuetz,ty_stuetz, FLOAT(tz_stuetz), triangles, INPUT=result, $
;            gs, limits, min_value=secmin)  <  max(z_stuetz)
       griddata(tx_stuetz,ty_stuetz,tz_stuetz, dimension=xy_sl_size,/linear, triangles=triangles, $
   missing=secmin,/grid) > secmin

;   if (interpolate eq 2) then result = $
;      griddata(tx_stuetz,ty_stuetz,tz_stuetz, method='Kriging',dimension=xy_sl_size,/linear, triangles=triangles)  > secmin
;    TRIGRID(tx_stuetz,ty_stuetz, FLOAT(tz_stuetz), triangles, INPUT=result, $
;            gs, limits, /quintic, min_value=secmin) <  max(z_stuetz)
;            gs, limits, /quintic, min_value=secmin, extrapolate=b) <  max(z_stuetz)

   if (interpolate eq 0 or interpolate ge 2) then BEGIN 
; no interpolation:  
      dxts =  DOUBLE(xy_sl_size)/(DOUBLE(slice_size[2])-DOUBLE(slice_size[0]))
      dyts =  DOUBLE(xy_sl_size)/(DOUBLE(slice_size[3])-DOUBLE(slice_size[1]))

      if interpolate eq 3 then tdx_stuetz = dx_stuetz-2.0001/dxts > 1./dxts else tdx_stuetz=dx_stuetz

      ind = sort(1.D/dx_stuetz)
      if interpolate eq 2 then begin
         xmi =  LONG64((x_stuetz[ind]-DOUBLE(slice_size(0)))*dxts) >  0  < (xy_sl_size-1)
         ymi =  LONG64((y_stuetz[ind]-DOUBLE(slice_size(1)))*dyts) > 0 < (xy_sl_size-1)
      endif else begin
         xmi =  LONG64((x_stuetz[ind]-0.5D*tdx_stuetz[ind]-DOUBLE(slice_size(0)))*dxts) >  0  < (xy_sl_size-1)
         ymi =  LONG64((y_stuetz[ind]-0.5D*tdx_stuetz[ind]-DOUBLE(slice_size(1)))*dyts) > 0 < (xy_sl_size-1)
      endelse
     
      xma =  LONG64((x_stuetz[ind]+0.5D*tdx_stuetz[ind]-DOUBLE(slice_size(0)))*dxts) < (xy_sl_size-1) > 0
      yma =  LONG64((y_stuetz[ind]+0.5D*tdx_stuetz[ind]-DOUBLE(slice_size(1)))*dyts) < (xy_sl_size-1) > 0
;      print, xmi[ind], xma[ind], ymi[ind], yma[ind]
      if interpolate eq 0 or interpolate eq 3 then begin
         if (data_format eq 20) then $
            FOR i=0L,n_stuetz-1L DO result[xmi[i]:xma[i],ymi[i]:yma[i]] += z_stuetz[ind[i]] $
         else $
            FOR i=0L,n_stuetz-1L DO result[xmi[i]:xma[i],ymi[i]:yma[i]] = z_stuetz[ind[i]]
      endif
      if interpolate eq 2 then FOR i=0L,n_stuetz-1L DO result[xmi[i],ymi[i]] = z_stuetz[ind[i]]
   ENDIF
   
   IF verbose THEN print, 'max after griding:',max(result), min(result)

   RETURN, result
END
;.compile construct_interpolated_slice_data


function construct_interpolated_slice_data
; PLOT DERIVED QUANTITIES such as angular momentum, Entropy and Pressure ...
@common_blocks.inc
; using grid_info read in data (found in data_dir directory) of all
; necessary grids for a slice ; specified by the 3 component vector
; slice_ori. slice_size has 4 components specifying the left and right
; edge of the slice
; note that slice_ori has to have two elements set to zero since we
; only do orthogonal slices.
; 
  interpolate =  interpolate_i

; Don't render if we don't have any baryon data loaded.
  slice_data = DBLARR(xy_sl_size, xy_sl_size)
  if (list_str[0] eq "No baryon data" or N_ELEMENTS(grid_info) eq 0) then begin
     print, "Warning: no baryon data to render."
     return, slice_data
  endif


; should I log it ? Yes unless
  logit = shouldbe_logged(list_str[var_index])
 
  num_dim = 3
; rough check of input data
  axes = (slice_ori gt 1.e-30)
  const_sub =  where(axes eq max(axes))
  const_sub = const_sub(0)
  other_subs = where(axes eq 0)

  if (data_format eq 2) then begin ; gadget data
     depthp = DOUBLE(slice_value+depth_value/2.) 
     depthm = DOUBLE(slice_value-depth_value/2.) 
     ind = where(((*parti.d[const_sub]) le depthp) and ((*parti.d[const_sub]) ge depthm))
     if ind[0] eq -1 then return, 0
     x_stuetz = (*parti.d[other_subs[0]])[ind]
     y_stuetz = (*parti.d[other_subs[1]])[ind]
     z_stuetz = (*parti.d[var_index])[ind]
     hsmlind = (where((*parti.names) eq 'hsml'))[0]
     dx_stuetz= (*parti.d[hsmlind])[ind]/(parti.size)
     n_stuetz = N_elements(ind)
     goto, enough
  endif

  
  IF (Fix(TOTAL(axes)) gt 1) then BEGIN
      print, 'construct_interpolated_slice_data: slice_ori:', $
        slice_ori,' not allowed !'
      STOP
  ENDIF
  IF (N_ELEMENTS(grid_info) lt 1) THEN BEGIN
      print, 'construct_interpolated_slice_data: no grid information in grid_info:',$
        grid_info
      STOP
  ENDIF
  dim_grid_info = size(grid_info)
  num_of_grids  = dim_grid_info[1]

; min left edge and scaling
  min_left  = DOUBLE([slice_size(0),slice_size(1)])
  max_right = DOUBLE([slice_size(2),slice_size(3)]) 
  scale_up    = min(1.D/(max_right-min_left))

  index_range    = grid_info.End_index-grid_info.Start_index+1
  delta_distance = DOUBLE(grid_info.Right_edge - grid_info.Left_edge)/DOUBLE(index_range)

; determine grids that are at least partially in the slice
  grid_flag      = grid_in_slice(grid_info, DOUBLE(slice_ori),  DOUBLE(slice_size))
  del_pix =   DOUBLE(max([xy_sl_size,xy_sl_size]))*delta_distance/(max_right(0)-min_left(0))

  min_data = 1.e30
  select = where((del_pix[other_subs[0],*] GE FLOAT(fastmode)) AND $
                 (del_pix[other_subs[0],*] LE FLOAT(xy_sl_size)/3) AND (grid_flag gt 0), count)
  if count gt 0 then u_g = grid_info[select] else u_g = grid_info
  
  if N_ELEMENTS(where(grid_info[0].dim[*] gt 0)) eq 2 then begin
     print, 'make sure to press the z-button in the GUI!!!'
     u_g = grid_info
  endif

; print,  del_pix
  dim_grid_info = size(u_g)
  num_of_grids  = dim_grid_info[1]

  fully_contained_above = min(where((u_g(*).Left_edge(other_subs(0)) le (min_left(0)))$
                                    and  (u_g(*).Left_edge(other_subs(1)) le (min_left(1))) $
                                    and (u_g(*).Right_edge(other_subs(0)) ge (max_right(0))) $
                                    and (u_g(*).Right_edge(other_subs(1)) ge (max_right(1))))) > 0

;  print, 'slice is fully contained in grids after index:', fully_contained_above
  
  index_range = u_g.End_index-u_g.Start_index+1
  delta_distance = (u_g.Right_edge - u_g.Left_edge)/DOUBLE(index_range)

  slice_ori_l                = dblarr(3)
  slice_ori_l(const_sub)     = slice_ori(const_sub)
  slice_ori_l(other_subs) = min_left
  slice_ori_r                = slice_ori_l 
  slice_ori_r(other_subs) = max_right
; 
  ii_s = FLOOR( 0. + $
               (slice_ori_l[const_sub] - u_g.Left_edge(const_sub))/delta_distance[const_sub,*] >0 ) $
              < (index_range[const_sub,*]-1)

  read_two = 0
  read_two = (total(["div V", "density slope", "DM density slope"] eq list_str[var_index]) gt 0) < 1
  print, "read_two", read_two
  if read_two then begin
     ind = where((ii_s eq (index_range[const_sub,*]-1)) eq 1, count)
     if count gt 0 then  ii_s[ind] -= 1
  endif

  ii_sp1 = (ii_s+1) ; plus one to have two indices to interpolate from

  ii_l1 =   (FLOOR((slice_ori_l[other_subs[0]] - u_g.Left_edge[other_subs[0]]) $
                 /delta_distance[other_subs[0],*]) -1) > 0
  ii_l2 =   (FLOOR((slice_ori_l[other_subs[1]] - u_g.Left_edge[other_subs[1]]) $
                 /delta_distance[other_subs[1],*]) -1) > 0
  ii_r1   = ((FLOOR((slice_ori_r[other_subs[0]] - u_g.Left_edge(other_subs[0]))/ $
                 delta_distance[other_subs[0],*]) ) < (index_range[other_subs[0],*]-1))  > 0
  ii_r2   = ((FLOOR((slice_ori_r[other_subs[1]] - u_g.Left_edge(other_subs[1]))/ $
                 delta_distance[other_subs[1],*]) ) < (index_range[other_subs[1],*]-1))  > 0

start = intarr(num_of_grids,3) & stride = intarr(num_of_grids,3) & righti = intarr(num_of_grids,3)
start[*,const_sub]     = ii_s[*]
start[*,other_subs[0]] = ii_l1
start[*,other_subs[1]] = ii_l2
if read_two then righti[*,const_sub]   = ii_sp1[*] $
else righti[*,const_sub]     = ii_s[*]

righti[*,other_subs[0]] = ii_r1
righti[*,other_subs[1]] = ii_r2


; read all selected grid data
  a = get_data(u_g, list_str[var_index], const_sub, LEFT_INDEX=start, RIGHT_INDEX=righti)
  if total(a.np_total) lt 1 then return, 0
  if verbose then print, 'constructing image'
  zero = -1e-30
  zero_solution_under_subgrid, u_g, a, ZERO=zero

  npoints =  TOTAL(float(a.dims[0])*a.dims[1]*a.dims[2])
  if verbose then  print, 'at most will have ', npoints, ' in slice'
  x_stuetz =  DBLARR(npoints) &  y_stuetz =  DBLARR(npoints) &  z_stuetz =  DBLARR(npoints)
  dx_stuetz = DBLARR(npoints)


  n_stuetz = 0L
  FOR i=0L,N_elements(a)-1 DO BEGIN
; check this
     i_s   = ii_s[i]
     i_sp1 = ii_sp1[i]
     i_l = [ii_l1[i], ii_l2[i]]
     i_r = [ii_r1[i], ii_r2[i]]
     if a[i].np_total eq 0 then continue

     b = (*a[i].data[0])

     if read_two then begin
        ind = where(total(b ne zero, const_sub+1) gt 1, Np)
        xdiff = (slice_ori[const_sub] - $
                (u_g[i].left_edge[const_sub] + (0.+i_s)*delta_distance[const_sub, i]) $
                < 2.*delta_distance[const_sub, i]) $
                > (- 0.*delta_distance[const_sub, i])  
        
        if const_sub eq 0 then slope = double(b[1, *, *] - b[0, *, *])/delta_distance[const_sub, i]
        if const_sub eq 1 then slope = double(b[*, 1, *] - b[*, 0, *])/delta_distance[const_sub, i]
        if const_sub eq 2 then slope = double(b[*, *, 1] - b[*, *, 0])/delta_distance[const_sub, i]
        
        if const_sub eq 0 then b = b[0, *, *] + slope*xdiff
        if const_sub eq 1 then b = b[*, 0, *] + slope*xdiff
        if const_sub eq 2 then b = b[*, *, 0] + slope*xdiff

        print, delta_distance[const_sub, i], xdiff/delta_distance[const_sub, i], i_s, i_sp1, index_range[const_sub,i]-1

        if (Np gt 0) then begin
           c = b
           c[*] = zero
           c[ind] = b[ind]
           b = c
        endif

;        b = REFORM(b)
     endif

     ind = where(b ne zero, Np)

     s = size(b, /dimensions)
     id = i_r[0]-i_l[0]+1
     jd = i_r[1]-i_l[1]+1
     if id lt 0 or jd lt 0 then begin
        print, 'construct_interpolated_slcie_data: does this data have at least 2 dimensions?'
        return, slice_data
     endif

     xy_st = array_indices([id,jd], LINDGEN(LONG(id)*jd), /DIMENSIONS)

     if (Np gt 0) then begin
         x_stuetz_h = u_g[i].Left_edge[other_subs[0]] + $
           DOUBLE(0.5+i_l[0]+xy_st[0,ind])*delta_distance[other_subs[0],i]
         y_stuetz_h =  u_g[i].Left_edge[other_subs[1]] + $
           DOUBLE(0.5+i_l[1]+xy_st[1,ind])*delta_distance[other_subs[1],i]
         nne = N_stuetz+np-1
         x_stuetz(n_stuetz:nne) = x_stuetz_h
         y_stuetz(n_stuetz:nne) = y_stuetz_h
         z_stuetz(n_stuetz:nne) = b[ind]
         dx_stuetz(n_stuetz:nne) =  delta_distance[0,i]
         N_stuetz += np
     endif
  ENDFOR                        ; over grids ...
  
enough:
  if verbose then print, 'n_stuetz:', n_stuetz
  zmax = 0.
  if (n_stuetz gt 0) then begin
     ind = lindgen(n_stuetz)
;     ind = ind[sort(z_stuetz[0:n_stuetz-1])]
     x_stuetz = x_stuetz[ind]
     y_stuetz = y_stuetz[ind]
     z_stuetz = z_stuetz[ind]
     dx_stuetz= dx_stuetz[ind]
     secmin =  min(z_stuetz, MAX=zmax)    
     if verbose then print, 'min/max data:', secmin,zmax
     IF logit THEN BEGIN 
;        z_stuetz =  ALOG10(abs(z_stuetz) >
;        abs(secmin))*z_stuetz/(abs(z_stuetz) > 1.e-40)
        z_stuetz =  ALOG10(z_stuetz) 
        secmin = ALOG10(secmin)
        zmax   = ALOG10(zmax)
        if verbose then print, 'construct_interpolated_slice_data: Took log of data value',secmin,zmax
     ENDIF
     
     if (secmin gt zmax) then slice_data[*] = secmin $
     else slice_data = render_image()

     schlieren = 0
     if ((schlieren eq 1) and (interpolate_i eq 1)) then $
        slice_data = alog10(sqrt((shift(slice_data, [-1,0])-shift(slice_data, [0,0]))^2 +  $
                          (shift(slice_data, [0,-1])-shift(slice_data, [0,0]))^2)/sqrt(2.))
        
  endif 
  
  if verbose then print, 'used ', n_stuetz, '=', SQRT(DOUBLE(n_stuetz)),'^2 points.'



  RETURN, slice_data
END
;.compile construct_interpolated_slice_data