pro construct_interpolated_slice_data, grid_info, slice_ori, slice_coord,$ data_dir,sds_num,slice_size, slice_data, SMOOTH=smooth, SECMIN=SECMIN, $ INTERPOLATE = interpolate, LOG=logit if (not keyword_set(INTERPOLATE)) then begin construct_slice_data, grid_info, slice_ori, slice_coord,$ data_dir,sds_num,slice_size, slice_data, SMOOTH=smooth, SECMIN=SECMIN return endif ; 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_coord has 4 components specifying the left and right ; edge of the slice ; slice_size: is a two component vector giving the x and y size in pixels ; note that slice_ori has to have two elements set to zero since we ; only do orthogonal slices. num_dim = 3 ; rough check of input data slice_data = fltarr(slice_size(0), slice_size(1)) ; slice_data(*,*) = 0. 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 (Fix(TOTAL(axes)) gt 1) then BEGIN print, 'construct_slice_data: slice_ori:', $ slice_ori,' not allowed !' STOP ENDIF IF (N_ELEMENTS(grid_info) lt 1) THEN BEGIN print, 'construct_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 = [slice_coord(0),slice_coord(1)] max_right = ([slice_coord(2),slice_coord(3)]) print, min_left scale_up = min(1./(max_right-min_left)) ; use only grids that will show up on the image: b = 0.7*(max_right-min_left)/FLOAT(slice_size) ; determine grids that are at least partially in the slice grid_flag = grid_in_slice(grid_info, slice_ori, slice_coord) index_range = grid_info.End_index-grid_info.Start_index+1 delta_distance = (grid_info.Right_edge - $ grid_info.Left_edge)/(index_range) u_g = grid_info(FIX(where((delta_distance(other_subs,*) ge max(b)) $ and (grid_flag gt 0) ) )) dim_grid_info = size(u_g) num_of_grids = dim_grid_info(1) ; initialize minimum of data min_data = 1.e30 FOR i =0, num_of_grids-1 DO BEGIN if ((u_g(i).Left_edge(other_subs(0)) le min_left(0))$ and (u_g(i).Left_edge(other_subs(1)) le min_left(1)) $ and (u_g(i).Right_edge(other_subs(0)) ge max_right(0)) $ and (u_g(i).Right_edge(other_subs(1)) ge max_right(1))) THEN $ fully_contained_above = i ENDFOR 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)/FLOAT(index_range) ; print, size(u_g) ; print, 'max:', min_left, max_right, scale_up x_stuetz = FLTARR(100000) y_stuetz = FLTARR(100000) z_stuetz = FLTARR(100000) n_stuetz = 0 this_level=0 le_o = u_g(0:num_of_grids-1).left_edge(other_subs(0)) le_t = u_g(0:num_of_grids-1).left_edge(other_subs(1)) ri_o = u_g(0:num_of_grids-1).right_edge(other_subs(0)) ri_t = u_g(0:num_of_grids-1).right_edge(other_subs(1)) dd_co = max(delta_distance(*,fully_contained_above)) dd_pix = FIX(max(dd_co*scale_up*slice_size)) dd_pix_r= dd_pix slice_ori_l = fltarr(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 ; if not on top grid modify bounds to get rid of edge effects Num_pix_x = slice_size(0) Num_pix_y = slice_size(1) IF (min(slice_ori_l(other_subs)-dd_co) ge 0.) THEN BEGIN slice_ori_l(other_subs)=slice_ori_l(other_subs)-dd_co Num_pix_x = Num_pix_x + dd_pix Num_pix_y = Num_pix_y + dd_pix ENDIF ELSE BEGIN dd_pix = 0 ENDELSE if (max(slice_ori_r(other_subs)+dd_co) le 1.) then begin slice_ori_r(other_subs)= (slice_ori_r(other_subs) +dd_co) Num_pix_x = Num_pix_x + dd_pix_r Num_pix_y = Num_pix_y + dd_pix_r ENDIF ELSE BEGIN dd_pix_r = 0 ENDELSE FOR i=fully_contained_above, num_of_grids-1 DO BEGIN ; check this delta_dist = delta_distance(*,i) data_points = index_range(*,i) ; yep check again !! i_s = ROUND((slice_ori_l(const_sub) - 0.5*delta_dist(const_sub) - $ u_g(i).Left_edge(const_sub))/delta_dist(const_sub) >0 $ < (data_points(const_sub)-1)) i_l = CEIL( (slice_ori_l(other_subs) - $ u_g(i).Left_edge(other_subs)) $ /delta_dist(other_subs) ) > 0 i_r = (FLOOR((slice_ori_r(other_subs) - $ u_g(i).Left_edge(other_subs))/ $ delta_dist(other_subs) - 1) < $ (data_points(other_subs)-1) ) > 0 if ((min(i_r-i_l) ge 2) ) THEN BEGIN read_hdf, data_dir+u_g(i).baryon_file, sds_num, temp s = size(temp) if (min([s(2),s(1)]) gt 3) THEN BEGIN CASE const_sub OF 0: temp = REFORM(temp(i_s, i_l(0):i_r(0),i_l(1):i_r(1))) 1: temp = REFORM(temp(i_l(0):i_r(0),i_s,i_l(1):i_r(1))) 2: temp = REFORM(temp(i_l(0):i_r(0), i_l(1):i_r(1),i_s)) ELSE: ENDCASE Leftedgeh = [u_g(i).Left_edge(other_subs(0)),$ u_g(i).Left_edge(other_subs(1))] ; print, 'size(temp):',size(temp), i_l, i_r FOR is=i_l(0),i_r(0) DO BEGIN FOR js=i_l(1),i_r(1) DO BEGIN x_stuetz_h = Leftedgeh(0) + $ float(is)*delta_dist(other_subs(0))+ $ 0.5*delta_dist(other_subs(0)) y_stuetz_h = Leftedgeh(1) + $ float(js)*delta_dist(other_subs(1))+ $ 0.5*delta_dist(other_subs(1)) ; make sure there won't be a finer grid close to this location: ; but don't worry for last grid: con = 0 if (i lt num_of_grids-1) then begin con_o = (le_o((i+1):(num_of_grids-1)) lt x_stuetz_h) and $ (le_t((i+1):(num_of_grids-1)) lt y_stuetz_h) con_t = (ri_o((i+1):(num_of_grids-1)) gt x_stuetz_h) and $ (ri_t((i+1):(num_of_grids-1)) gt y_stuetz_h) ; stop endif if ((max(con_o and con_t) lt 1) and $ ((temp(is-i_l(0),js-i_l(1)) ge 1.e-7) or $ (temp(is-i_l(0),js-i_l(1)) le (-1.e-7)))) then begin x_stuetz(n_stuetz) = x_stuetz_h y_stuetz(n_stuetz) = y_stuetz_h z_stuetz(n_stuetz) = temp(is-i_l(0),js-i_l(1)) n_stuetz = n_stuetz+1 endif ENDFOR ENDFOR ; print, size(con_o),size(con_t),size(le_o),size(x_stuetz_h) ENDIF ENDIF ENDFOR print, size(con), con print, 'n_stuetz:', n_stuetz ; interpolate all values ; get triangles: x_stuetz = x_stuetz(0:n_stuetz-1) y_stuetz = y_stuetz(0:n_stuetz-1) z_stuetz = z_stuetz(0:n_stuetz-1) SECMIN = MIN(z_stuetz(where(z_stuetz gt min(z_stuetz)))) if (keyword_set(logit)) THEN z_stuetz = $ ALOG10(abs(z_stuetz) > secmin)*z_stuetz/(abs(z_stuetz) > 1.e-30) TRIANGULATE, x_stuetz, y_stuetz, triangles, b ; print, 'b:',b result = FLTARR(Num_pix_x,Num_pix_y) if (keyword_set(INTERPOLATE)) then begin if (interpolate eq 2) then result = $ TRIGRID(x_stuetz,y_stuetz, z_stuetz, triangles, INPUT=result, $ NX=Num_pix_x,Ny=Num_pix_y, /quintic, min_value=-10. ) if (interpolate eq 1) then result = $ TRIGRID(x_stuetz,y_stuetz, z_stuetz, triangles, INPUT=result, $ NX=Num_pix_x,Ny=Num_pix_y, min_value=-10. ) endif print, 'size(result)',size(result), 'dd_pix:', dd_pix, slice_size(0) ; result = tri_surf(z_stuetz(0:n_stuetz-1), $ ; x_stuetz(0:n_stuetz-1), $ ; y_stuetz(0:n_stuetz-1), $ ; NX=Num_pix_x,Ny=Num_pix_y $ ; ) > min(z_stuetz(0:n_stuetz-1)) slice_data = result(dd_pix:(dd_pix+slice_size(0)-1), $ dd_pix:(dd_pix+slice_size(1)-1)) return END ; .compile construct_interpolated_slice_data.pro