function fields_to_indices, DATA=a, FIELDS=fields
    data_names = strtrim(a[0].data_fields,2)
    if N_ELEMENTS(fields) eq 0 then $
        indices = LINDGEN(a[0].n_dfields) $
    else begin
        indices = LONARR(N_ELEMENTS(fields))
        for i=0L,N_ELEMENTS(fields) do $
            indices[i] = (where(data_names eq fields[i]))[0]
        indices = indices[uniq(indices)]
        indices = indices[where(indices ne -1)]
    endelse
    return, indices
end

pro ghost_grid, PROJECTION=proj, DATA=a, THICKNESS=cells, FIELDS=fields
    info = *proj.grids
    ii_l = TRANSPOSE(*proj.ii_l) & ii_r = TRANSPOSE(*proj.ii_r) & ii_d = TRANSPOSE(*proj.ii_d)
    indices = fields_to_indices(DATA=a, FIELDS=fields)
    level_order = sort([info.level])

    ; get the parent indices
    pp = LONARR(n_elements(info))
    for i=0L,n_elements(info)-1 do pp[i] = (where(info.num eq info[i].parent))[0]

    if not KEYWORD_SET(cells) then cells = 0
    iW = ii_d + 2*cells
    ; the size ratio of child cells to parent cells
    scale = 2^(info.level - info[pp].level)
    ; calculate how many layers of parent cells we need
    parent_cells = 1L + (cells + scale/2 - 1) / scale
    ; calculate the integer coordinates of all the grids
    dims = grid_int_coords(INFO=info)
    ghosted = INTARR(N_ELEMENTS(info))
    ; determine which grids are neighbors
    neighbors = grid_neighbors(INFO=info, DIMS=dims, CELLS=cells)
    ; now, for each grid,
    for jj=0L,n_elements(info)-1 do begin
        j = level_order[jj]
        p = pp[j]
        ; if p < 0 then use this grid instead of a parent
        if p lt 0 then begin
            p = 0
            parent_start = [0,0,0]
            parent_end   = ii_d[*,j] - 1
            ; find the indices of the points of the new grid
            X = LINDGEN(iW[0,j]) - cells
            Y = LINDGEN(iW[1,j]) - cells
            Z = LINDGEN(iW[2,j]) - cells
        endif else begin
            ; note that the parent has already been ghosted!
            parent_start = ((info[j].parent_s_index - parent_cells[j]) > (ii_l[*,p]-cells) < (ii_r[*,p]+cells)) - ii_l[*,p] + cells
            parent_end   = ((info[j].parent_e_index + parent_cells[j]) > (ii_l[*,p]-cells) < (ii_r[*,p]+cells)) - ii_l[*,p] + cells
            ; find the coordinates of the points of the new grid
            ; in the parent grid's index space
            X = (DINDGEN(iW[0,j]) - cells + .5) / scale[j] - .5 + parent_cells[j] - ii_l[0,p]
            Y = (DINDGEN(iW[1,j]) - cells + .5) / scale[j] - .5 + parent_cells[j] - ii_l[1,p]
            Z = (DINDGEN(iW[2,j]) - cells + .5) / scale[j] - .5 + parent_cells[j] - ii_l[2,p]
        endelse
        ; decide where to splat the neighbors
        n_ind = where(neighbors[*,j] ge 0)
        if n_ind[0] ge 0 then begin
            n = neighbors[n_ind, j]
            ; how much can we get from this neighbor?
            ; start with the box dims[0:2,n]+ii_l[0:2,n] : dims[0:2,n]+ii_r[0:2,n]
            ; clip to dims[0:2,j]+ii_l[0:2,j]-cells : dims[0:2,j]+ii_r[0:2,j]+cells to get absolute index coords
            ; subtract dims[0:2,n]+ii_l[0:2,n] to get neighbor index coords
            ; subtract dims[0:2,j]+ii_l[0:2,j]-cells to get local splat coords
            left    = CMREPLICATE(dims[0:2,j]+ii_l[0:2,j]-cells, N_ELEMENTS(n))
            right   = CMREPLICATE(dims[0:2,j]+ii_r[0:2,j]+cells, N_ELEMENTS(n))
            n_left  = (dims[0:2,n]+ii_l[0:2,n]) > left < right
            n_right = (dims[0:2,n]+ii_r[0:2,n]) > left < right
            ; nn_s: neighbor's neighbor-index coords, start
            ; ns_e: neighbor's splat-index coords, end
            nn_s = n_left  - (dims[0:2,n]+ii_l[0:2,n])
            nn_e = n_right - (dims[0:2,n]+ii_l[0:2,n])
            ns_s = n_left  - left
            ns_e = n_right - left
            ; some of the neighbors may already have been ghosted
            n_ghosted = where(ghosted[n])
            if n_ghosted[0] ge 0 then begin
                nn_s[*, n_ghosted] += cells
                nn_e[*, n_ghosted] += cells
            endif
        endif else n = -1

        ; for each field on grid j
        for kk=0L,N_ELEMENTS(indices)-1 do begin
            if cells gt 0 then begin
                ; get the region in question
                pk = (*a[p].data[indices[kk]])[parent_start[0]:parent_end[0], parent_start[1]:parent_end[1], parent_start[2]:parent_end[2]]
                ; interpolate
                k = interpolate(pk, X, Y, Z, /GRID)
                ; splat neighbors
                if n[0] ge 0 then begin
                    for ni=0L,N_ELEMENTS(n)-1 do begin
                        k[ns_s[0,ni]:ns_e[0,ni], ns_s[1,ni]:ns_e[1,ni], ns_s[2,ni]:ns_e[2,ni]] = $
                            (*a[n[ni]].data[indices[kk]])[nn_s[0,ni]:nn_e[0,ni], nn_s[1,ni]:nn_e[1,ni], nn_s[2,ni]:nn_e[2,ni]]
                    endfor
                endif
            endif else $
                k = FLTARR(iW[0, j], iW[1,j], iW[2, j])
            ; splat
            k[cells,cells,cells] = *a[j].data[indices[kk]]
            ; store
            *a[j].data[indices[kk]] = k
        endfor
        ghosted[j] = 1
    endfor
end