; PLOT PROJECTIONS OF EMISSION LINES THROUGH SIMULATION VOLUME function project_with_absorption, PROJECTION=proj, LINE_INDEX=line_index, OPACITY=opacity, STARS_ONLY=stars_only, VERBOSE=verbose, CACHE=cache 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 compute_emission = (stars_only eq 0) and not opacity 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=[dust_cross_section(/NEEDED_FIELDS), $ *(line_emissivity(/AVAILABLE_LINES))[line_index].needed_fields] grid_data = *proj.grid_data if verbose then print, 'Computing additional derived fields...' lambda = (line_emissivity(/AVAILABLE_LINES))[line_index].lambda * 1d-7 ; convert nm to cm absorption_data = dust_cross_section(LAMBDA=lambda, DATA=grid_data) emissivity_data = -1 if compute_emission then $ emissivity_data = line_emissivity(LINE_INDEX=line_index, DATA=grid_data) cache = {absorption_data:PTR_NEW(absorption_data), emissivity_data:PTR_NEW(emissivity_data), axes:slice.axes} endif else begin absorption_data = *cache.absorption_data if compute_emission then $ 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.absorption_data)-1 do begin transpose_grid, GRID=absorption_data[i], COUNT=ii_d[i,cache.axes mod 3]+2*proj.border, AXES=swap if compute_emission then $ 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 alpha_v = slice_data ds = slice_data dtau = slice_data extinction = slice_data efficiency = slice_data if compute_emission then begin emissivity = slice_data emission = slice_data per_steradian = 1 / (4.D * 3.14159265358979D) endif w = proj.slice.image_size if (stars_only gt 0) and PTR_VALID(proj.part_data) then begin star_data = *proj.part_data n_stars = (size(star_data))[1] star_xyzw = [[star_data[*,0:2]], [DBLARR(n_stars) + 1d]] b = proj.slice.bounds_proj l = b[*,0] d = b[*,1]-b[*,0] ; form the shear permutation matrix P = DBLARR(4,4) P[axes mod 3,[0,1,2]] = (axes lt 3) * 2 - 1 & P[3,3] = 1 SP = [[1,0,proj.slice.su,0],[0,1,proj.slice.sv,0],[0,0,1,0],[0,0,0,1]] ## P ; adjust for origin offset SP_off = [[1,0,0,.5],[0,1,0,.5],[0,0,1,.5],[0,0,0,1]] ## SP ## [[1,0,0,-.5],[0,1,0,-.5],[0,0,1,-.5],[0,0,0,1]] ; scale x and y axes by image dimensions and discard w coordinate SP_scl = [[w,0,0,1],[0,w,0,1],[0,0,1,0]] ## [[1d/d[0],0,0,-l[0]/d[0]],[0,1d/d[1],0,-l[1]/d[1]],[0,0,1,0],[0,0,0,1]] ## SP_off ; now calculate our final star positions star_xyz = SP_scl ## star_xyzw ind = reverse(sort(star_xyz[*,2])) star_x = star_xyz[ind,0] > 1. < w star_y = star_xyz[ind,1] > 1. < w star_z = star_xyz[ind,2] star_M = star_data[ind,3] star_L = star_data[ind,4] star_T = star_data[ind,5] star_radiates = star_data[ind,6] ; decide to which depths the stars belong and build index lists star_depth_indices = LONG(INTERPOL(DINDGEN(n_depths), depth_values[unique_depths], star_z)) < (n_depths-1) > 0 hi = HISTOGRAM(star_depth_indices, REVERSE_INDICES=star_ri, MIN=0, MAX=n_depths-1) endif 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) alpha_v[*] = 0.D ds[*] = 0.D if compute_emission then 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] splat_data, TARGET=alpha_v, SOURCE=absorption_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 if compute_emission then $ 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] endfor ; CALCULATE EXTINCTION AND EMISSION ; extinction calculation - I1 = I0 * e^(-alpha ds) + j/alpha * (1-e^(-alpha ds)) dtau[0,0] = alpha_v*ds extinction[0,0] = exp(-dtau) 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) if not opacity then begin ; APPLY RESULTS TO SLICE if stars_only eq 0 then begin ix = where(dtau gt 0.01) ; substitute Taylor polynomial where calculation is poorly conditioned ;efficiency[0,0] = ds-alpha_v*ds^2/2.+alpha_v^2*ds^3/6.-alpha_v^3*ds^4/24. efficiency[0,0] = ds*(-.5D*dtau*(.33333D*dtau*(.25D*dtau-1.D)+1.D)+1.D) ; faster if ix[0] ge 0 then efficiency[ix] = (1.D - extinction[ix]) / alpha_v[ix] emission[0,0] = per_steradian * emissivity * efficiency if need_shift then extinction[0,0] = SHIFT(extinction, su, sv) if need_shift then emission[0,0] = SHIFT(emission, su, sv) slice_data[0,0] = slice_data * extinction + emission endif else if N_ELEMENTS(star_x) gt 0 then begin ind_start = star_ri[j] ind_end = star_ri[j+1]-1 if ind_end ge ind_start then ind = star_ri[ind_start:ind_end] else ind = -1 if ind[0] ge 0 then begin if stars_only eq 2 then ind2 = where(star_radiates[ind]) $ else ind2 = where(not star_radiates[ind]) if ind2[0] ge 0 then ind = ind[ind2] else ind = -1 endif if need_shift then extinction[0,0] = SHIFT(extinction, su, sv) if ind[0] ge 0 then begin starlight = hist2d(star_x[ind], star_y[ind], star_L[ind], MIN1=1.,MIN2=1.,MAX1=w,MAX2=w) slice_data[0,0] = slice_data * extinction + starlight endif else slice_data[0,0] = slice_data * extinction endif endif else begin if need_shift then extinction[0,0] = SHIFT(extinction, su, sv) slice_data[0,0] = 1.D - (1.D - slice_data) * extinction endelse 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 if stars_only eq 0 then $ slice_data = warp_high_quality(slice_data, proj.slice) $ else $ slice_data = warp_points(slice_data, proj.slice) if verbose then show_intermediate_image, slice_data, 'Intensity' if not opacity then begin slice_data = (slice_data/proj.intensity_unit) > 0 if shouldbe_logged('Intensity') then slice_data = alog10(slice_data) endif if verbose then print, 'Projection with absorption done.' return, slice_data end function construct_projection_with_absorption, OPACITY=opacity, STARS_ONLY=stars_only, INTERPOLATION=interp @common_blocks.inc opacity = KEYWORD_SET(opacity) if not KEYWORD_SET(stars_only) then stars_only = 0 if N_ELEMENTS(interp) eq 0 then interp = 1 proj = build_projection(INTERPOLATION=interp) return, project_with_absorption(PROJECTION=proj, LINE_INDEX=line_index, OPACITY=opacity, STARS_ONLY=stars_only) end