; PLOT PROJECTIONS OF FIELD THROUGH SIMULATION VOLUME
function project, PROJECTION=proj, FIELD=field, VERBOSE=verbose, CACHE=cache
;@common_blocks.inc
   verbose = KEYWORD_SET(verbose)
   slice = proj.slice
   xyi = slice.xyi & xi = slice.xi & yi = slice.yi & zi = slice.zi & axes = slice.axes
   grids = *proj.grids & ii_l = *proj.ii_l & ii_r = *proj.ii_r & ii_d = *proj.ii_d & delta_distance = *proj.delta_distance & dp = *proj.dp
   depth_values = *proj.depth_values & grid_index = *proj.grid_index & grid_values = *proj.grid_values & unique_depths = *proj.unique_depths & n_depths = proj.n_depths
   grid_xy = *proj.grid_xy & grid_wh = *proj.grid_wh
   border = proj.border

   cache_good = 0
   if N_ELEMENTS(cache) ne 0 then begin
       cache_good = 1
       if TOTAL(cache.axes eq slice.axes) eq 3 then cache_good = 2
   endif
   

   if cache_good eq 0 then begin
       get_projection_data, PROJECTION=proj, FIELDS=[field]
       grid_data = *proj.grid_data

       if verbose then print, 'Fetching field...'
       field_index = (where(strtrim(grid_data[0].data_fields,2) eq field))[0]
       emissivity_data = PTRARR(N_ELEMENTS(grid_data), /ALLOCATE_HEAP)
       for j=0L,n_elements(grid_data)-1 do begin
           *emissivity_data[j] = fn(grid_data[j], field)
       endfor
       cache = {emissivity_data:PTR_NEW(emissivity_data), axes:slice.axes}
   endif else begin
       emissivity_data = *cache.emissivity_data
       if cache_good eq 1 then begin
           swap = axes_diff(OLD=cache.axes, NEW=slice.axes)
           for i=0L,N_ELEMENTS(*cache.emissivity_data)-1 do begin
               transpose_grid, GRID=emissivity_data[i], COUNT=ii_d[i,cache.axes mod 3]+2*proj.border, AXES=swap
           endfor
           cache.axes = slice.axes
       endif
   endelse

   slice_data = DBLARR(proj.slice.image_size,proj.slice.image_size)

   if verbose then print, 'Projecting...'
   ; ITERATE OVER DEPTHS AND GRIDS
   ds = slice_data
   emissivity = slice_data
   emission = slice_data
   per_steradian = 1 / (4.D * 3.14159265358979D)

   w = proj.slice.image_size

   interp = border > 0 < 2
   for j=0L, n_depths-1 do begin
       indices            = where(depth_values eq depth_values[unique_depths[j]])
       layer_grids        = grid_values[indices]
       layer_grid_indices = grid_index[indices]
       layer_order = sort(grids[layer_grids].level)
       ds[*] = 0.D
       emissivity[*] = 0.D
       for jj=0L, N_ELEMENTS(layer_grids)-1 do begin
           i = layer_grids[layer_order[jj]] & xy = grid_xy[*,i] & wh = grid_wh[*,i]
;print, jj, i
;help, *emissivity_data[i], layer_grid_indices[jj]
sss = size(*emissivity_data[i])
if layer_grid_indices[jj]+border le sss[3]-1 then begin
           splat_data, TARGET=emissivity, SOURCE=emissivity_data[i], COORDS=grid_xy[*,i], DIMS=grid_wh[*,i], PIXEL=dp[*,i], LAYER=layer_grid_indices[jj], SHAPE=ii_d[i,xyi]+2*border, INTERP=interp

           splat_const_data, TARGET=ds, DATA=delta_distance[zi,i] * proj.distance_unit * proj.slice.z_rate, COORDS=grid_xy[*,i]
        endif else begin
           print, 'layer value and grid dimensions do not match'
           print, layer_grid_indices[jj]+border 
           print, sss
           print, j, jj
        endelse
       endfor
       ; CALCULATE EMISSION
       ; emission calculation - I1 = I0 + j * ds

       z = depth_values[unique_depths[j]]-.5D
       su = (LONG(z*(w/2.)*proj.slice.su) + w) mod w
       sv = (LONG(z*(w/2.)*proj.slice.sv) + w) mod w
       need_shift = (su ne 0) or (sv ne 0)

       ; APPLY RESULTS TO SLICE
       emission[0,0]   = per_steradian * emissivity * ds
       if need_shift then emission[0,0] = SHIFT(emission, su, sv)
       slice_data[0,0] = slice_data + emission

       if (j mod 40) eq 0 then begin
           check_cancel, stopit=stopit
           ;if verbose then $
           ;    show_intermediate_image, slice_data, 'Intensity'
       endif
       if stopit then begin
           stopit = 0
           break
       endif
   endfor
   slice_data = warp_high_quality(slice_data, proj.slice) 
   if verbose then show_intermediate_image, slice_data, 'Intensity'

   slice_data = (slice_data) > 0
   if shouldbe_logged(Field) then slice_data = alog10(slice_data)
   if verbose then print, 'Projection done.'
   print, mean(slice_data)
;breaki
   return, slice_data
end

function construct_projection, FIELD=field, INTERPOLATION=interp
;@common_blocks.inc
   if N_ELEMENTS(interp) eq 0 then interp = 1

   proj = build_projection(INTERPOLATION=interp)
   return, project(PROJECTION=proj, FIELD=field, verbose=1)
end