pro select_stars, BOUNDS=bounds, STAR_DATA=data
@common_blocks.inc
    xyz = [[*parti.d[0]], [*parti.d[1]], [*parti.d[2]]]
    ind = where((xyz[*,0] ge bounds[0,0]) and (xyz[*,0] le bounds[0,1]) and $
                (xyz[*,1] ge bounds[1,0]) and (xyz[*,1] le bounds[1,1]) and $
                (xyz[*,2] ge bounds[2,0]) and (xyz[*,2] le bounds[2,1]))
    if ind[0] lt 0 then begin
        xyz = -1
        return
    endif

    xyz = xyz[ind,*]

    ; L/L0 = (M/M0) ^ 3.5
    ; R/R0 = (M/M0) ^ 0.8
    ; L/L0 = (R/R0)^2 * (T/T0)^4
    ; (M/M0) ^ 3.5 = (M/M0) ^ 1.6 * (T/T0)^4
    ; T/T0 = (M/M0)^0.475
    M = fp('particle mass', ind) * get_unit('mass', /FORCE_USE_UNITS) ; solar masses
    L = M^3.5   * 3.846e33 ; mass-luminosity relation - ergs/sec
    T = M^0.475 * 5778     ; Kelvin
    ; the flag for UV luminous stars is located (temporarily?) in 'dynamical_time'
    radiates = M gt 30. ;fp('dynamical_time', ind) ne 0.
    ; particles with masses less than 0.012 solar masses are too small to undergo fusion
    ; but those stars will be trillions of times dimmer than the brightest stars
    ; we'll cut off stars below 1000 K as well since they aren't very luminous in the visible
    ; spectrum

    ind = where((M gt 0.012) and (T gt 1000))

    if ind[0] ge 0 then begin
        xyz = xyz[ind,*] & M = M[ind] & L = L[ind] & T = T[ind] & radiates = radiates[ind]
        data = [[xyz],[M],[L],[T],[radiates]]
    endif
end

pro get_projection_data, PROJECTION=proj, FIELDS=fields_needed, NOPARTS=noparts
@common_blocks.inc
    xi = proj.slice.xi & yi = proj.slice.yi & zi = proj.slice.zi
    if not PTR_VALID(proj.part_data) and not KEYWORD_SET(noparts) then begin
        if N_ELEMENTS(*parti.names) ne 0 and (*parti.names)[0] ne '' then begin
            ;per_area_per_steradian = proj.slice.image_size^2 / $
            ;    ((bounds[xi,1]-bounds[xi,0]) * (bounds[yi,1]-bounds[yi,0]) * $
            ;    get_unit('Length')^2 * 4 * !pi)
            read_particle_data, FIELDS=[0,1,2,(where(*parti.names eq 'particle mass'))[0],(where(*parti.names eq 'dynamical_time'))[0]]
            select_stars, BOUNDS=proj.slice.bounds_data, STAR_DATA=star_data
            if N_ELEMENTS(star_data) gt 0 then $
                proj.part_data = PTR_NEW(TEMPORARY(star_data))
        endif
    endif

    fields_needed = fields_needed[uniq(fields_needed)]
    do_read = 1

    if PTR_VALID(proj.grid_data) then begin
        fields_existing = (*proj.grid_data)[0].data_fields[0:(*proj.grid_data)[0].n_dfields-1]
        keep = INTARR(N_ELEMENTS(fields_needed))
        for i=0L, N_ELEMENTS(fields_needed)-1 do $
            keep[i] = (where(fields_needed[i] eq fields_existing))[0] lt 0

        if TOTAL(keep) gt 0 then $
            fields_needed = fields_needed[where(keep)] $
        else $
            do_read = 0
    endif

    if do_read then begin
        print, 'Reading ', string(FORMAT='(I0)', N_ELEMENTS(fields_needed)), ' data fields from ', string(FORMAT="(I0)", N_ELEMENTS(*proj.grids)), ' grids...'
        ;print, fields_needed
        new_grid_data = get_data(*proj.grids, fields_needed, zi, LEFT_INDEX=*proj.ii_l, RIGHT_INDEX=*proj.ii_r)

        if proj.border gt 0 then begin
            print, 'Computing ghost values...'
            ghost_grid, PROJECTION=proj, DATA=new_grid_data, THICKNESS=proj.border
        endif
        print, 'Transposing data...'
        transpose_data, GRID_DATA=new_grid_data, COUNT=*proj.ii_d+2*proj.border, INDICES=LINDGEN(new_grid_data[0].n_dfields), AXES=proj.slice.axes
            
        if not PTR_VALID(proj.grid_data) then $
            grid_data = new_grid_data $
        else begin
            grid_data = *proj.grid_data
            k = grid_data[0].n_dfields
            n = new_grid_data[0].n_dfields
            grid_data.n_dfields += new_grid_data.n_dfields
            ; IDL doesn't support the obvious direct array indexing assignments here, so we use a temporary variable
            b = grid_data.data_fields   & b[k:k+n-1,*] = new_grid_data.data_fields[0:n-1] & grid_data.data_fields = b
            b = grid_data.data          & b[k:k+n-1,*] = new_grid_data.data[0:n-1]        & grid_data.data = b
        endelse
        proj.grid_data = PTR_NEW(grid_data)
    endif
end

function build_projection, INTERPOLATION=interp, IMAGE_SIZE=image_size, NOSWAPPING=noswapping
@common_blocks.inc
    if not KEYWORD_SET(image_size) then image_size = xy_sl_size
    if not KEYWORD_SET(interp) then interp = interpolate_i
    slice = define_slice(SLICE_ORI=slice_ori, SLICE_SIZE=slice_size, SLICE_VALUE=slice_value, DEPTH_VALUE=depth_value, IMAGE_SIZE=image_size)
    if KEYWORD_SET(noswapping) then begin
        slice.axes = [0,1,2]
        slice.xi = 0 & slice.yi = 1 & slice.zi = 2 & slice.xyi = [0,1]
        slice.bounds_proj = [0,0,0,1,1,1]
        slice.bounds_data = [0,0,0,1,1,1]
    endif
    proj = define_projection(GRIDS=grid_info, SLICE=slice, FASTMODE=fastmode)
    proj.border = interp
    return, proj
end