; OJ__define.pro v1.0 2000/05/23- ; 2007/12 integrate in jaques ; AUTHOR: Tom Abel ; ; NAME: OJ ; ; PURPOSE: ; This object class provides the data structure and methods to ; handle structure adaptive mesh refinement data of the ; enzo code. ; ; Methods so far allow to extract polygon and volume data objects ; ; CALLING SEQUENCE: ; To inititally create ; oOJ = OBJ_NEW('OJ') ; this will spawn GUI file selection to get filename of ; .hierarchy file ; you can supply a n exisiting hierarchy ; oOJ = OBJ_NEW('OJ', grid_info) ; or give a filename ; oOJ = OBJ_NEW('OJ', '~/data/RedshiftOutput0001.hierarchy') ; ; ; METHODS: ; ; I OJ::Get_GridStructure() ; ; grid_info = oOJ -> Get_GridStructure() ; returns the grid_info structure ; ; II OJ::Get_Polygons() ; ; an object array of IDLgrModels defined analogous to ; ; oGridPolys = OBJARR(maximum_number_OF_levels) ; FOR i = 0,maximum_number_OF_levels-1 DO BEGIN ; oGridPolys[i] = OBJ_NEW('IDLgrModel') ; ENDFOR ; ; will be returned by ; oGridPolys[0:*] = oOJ -> Get_Polygons() ; ; to request polygons for level 2,3,6,8,and 10 with the ; subvolume with Left_edge = [0.3, 0.4, 0.4] and ; Right_edge[0.9,0.8,0.9] use ; ; oGridPolys[0:*] = oOJ -> $ ; Get_Polygons(LEVELS=[2,3,6,8,10], ; LEFT_EDGE = [0.3, 0.4, 0.4], RIGHT_EDGE[0.9,0.8,0.9]) ; ; III OJ::Get_Slice_Data ; ; the heart of Jacques are these slices: ; ; returns a thin slice described by slice_ccord = [0, 0.45, 0] ; (in this example slice of constant y coordinate = 0.45) which ; is interpolated at the local resolution. ; ; slice = oOJ -> Get_Slice_Data(slice_coord) ; returns the slice using all available levels and volume ; (i.e all defined in the initial call to OBJ_NEW('OJ') ; ; Additional control keywords exist: ; slice = oOJ -> Get_Slice_Data(slice_coord, LEVELS=[2,3,6,8,10], $ ; LEFT_EDGE = [0.3D, 0.4D, 0.4D], RIGHT_EDGE[0.9D,0.8D,0.9D],$ ; DATA_FIELDS=[0,1,2,3,4,10,11,12] ) ; ; slice is a structure defined as: ; slice_struct = {slice, $ ; min_level: 0B, $ ; max_level: 0B, $ ; Np_total: 0L, $ ; N_d_dields: 0B, $ ; slice_coord:[0D,0D,0D], $ ; LEFT_EDGE: [0.D, 0.D, 0.D], $ ; RIGHT_EDGE: [1.D, 1.D, 1.D], $ ; DATA_FIELDS: BYTARR(100), $ ; X_PTR: PTRARR(1), $ ; Y_PTR: PTRARR(1), $ ; Z_PTR: PTRARR(100), $ ; LEVEL_PTR: PTRARR(1), $ ; dx_LEVEL: DBLARR(100) ) ; ; III OJ::Get_3D_Data ; ; returns the requested data fields from all files ; like particle data ; ; the file test_3D.pro in this directory shows how to use this method ; ; threeD = oOJ -> Get_3D_Data() ; returns 3d particle-like data using all available levels and volume ; (i.e all defined in the initial call to OBJ_NEW('OJ') ; ; Additional control keywords exist: ; slice = oOJ -> Get_Slice_Data(slice_coord, LEVELS=[2,3,6,8,10], $ ; LEFT_EDGE = [0.3D, 0.4D, 0.4D], RIGHT_EDGE[0.9D,0.8D,0.9D],$ ; DATA_FIELDS=[0,1,2,3,4,10,11,12] ) ; ; threeD is a structure defined as: ; threeD = {threeD, $ ; min_level: 0B, $ ; max_level: 0B, $ ; Np_total: 0L, $ ; N_dfields: 0B, $ ; LEFT_EDGE: [0.D, 0.D, 0.D], $ ; RIGHT_EDGE: [1.D, 1.D, 1.D], $ ; DATA_FIELDS: BYTARR(100), $ ; X_PTR: PTR_NEW(), $ ; Y_PTR: PTR_NEW(), $ ; Z_PTR: PTR_NEW(), $ ; D_PTR: PTRARR(100), $ ; LEVEL_PTR: PTR_NEW(), $ ; dx_LEVEL: DBLARR(100) } ; ; ; III OJ::Get_Orthogonal_Projection ; ; ; returns aprojection along one axis described by oproj_ccord = [0, 0.45, 0] ; (in this example a projection along y centered at = 0.45) which ; is interpolated at the local resolution and has a width given ; by oproj_width ; ; slice = oOJ -> Get_Orthogonal_Projection(slice_coord) ; returns the slice using all available levels and volume ; (i.e all defined in the initial call to OBJ_NEW('OJ') ; ; Additional control keywords exist: ; proj = oOJ -> Get_Orthogonal_Projection(oproj_coord, $ ; oproj_width , LEVELS=[2,3,6,8,10], $ ; LEFT_EDGE = [0.3D, 0.4D, 0.4D], RIGHT_EDGE[0.9D,0.8D,0.9D],$ ; DATA_FIELDS=[0,1,2,3,4,10,11,12] ) ; ; alternatively to the array of levels you can specify min_level ; and max_level to only project between min_level and max_level ; ; proj is a structure defined as: ; slice_struct = {proj, $ ; min_level: 0B, $ ; max_level: 0B, $ ; Np_total: 0L, $ ; N_dfields: 0B, $ ; oproj_coord:[0D,0D,0D], $ ; oproj_width:0D, $ ; LEFT_EDGE: [0.D, 0.D, 0.D], $ ; RIGHT_EDGE: [1.D, 1.D, 1.D], $ ; DATA_FIELDS: BYTARR(100), $ ; X_PTR: PTRARR(1), $ ; Y_PTR: PTRARR(1), $ ; Z_PTR: PTRARR(100), $ ; LEVEL_PTR: PTRARR(1), $ ; dx_LEVEL: DBLARR(100) ) ; FUNCTION vec_LT_vec, u, v h = u LT v RETURN, MIN([(h),1]) END FUNCTION vec_LE_vec, u, v h = u LE v RETURN, MIN([(h),1]) END FUNCTION vec_GE_vec, u, v h = u GE v RETURN, MIN([(h),1]) END FUNCTION vec_GT_vec, u, v h = u GT v RETURN, MIN([(h),1]) END FUNCTION OJ::INIT, gi_orfname @common_blocks.inc IF N_ELEMENTS(gi_orfname) le 0 then begin print, 'OJ:INIT called with invalid grid_info of filename' return, 0 endif ; intialize OJ ; if it is a string assume it is the filename and read the hierarchy if size(gi_orfname, /TYPE) eq 7 then begin grid_info = read_grid_info(gi_orfname) end else grid_info = gi_orfname left_edge=[1e50,1e50,1e50] right_edge= -1.*[1e50,1e50,1e50] for i=0,N_elements(grid_info[0].Left_edge)-1 do begin Left_edge[i] = min(grid_info.left_edge[i]) Right_edge[i] = max(grid_info.left_edge[i]) endfor ; II) Initialize Object Structure from what we know so far self.file_base = grid_info[0].hier_file self.Left_edge = Left_edge self.Right_edge = Right_edge ; III) read .hierarchy file self.gi = grid_info[*] self.num_of_grids = N_elements(grid_info) ci = self.num_of_grids print, 'OJ: ',ci, ' grids included.' self.max_level = max(self.gi[0:(ci-1)].level) self.min_level = min(self.gi[0:(ci-1)].level) print, 'amr levels: max, min',self.max_level, self.min_level RETURN, 1 END PRO OJ::CLEANUP ptr_free, self.data[*] self = 0 END FUNCTION OJ::Get_GridStructure RETURN, self.gi[0:(self.num_OF_grids-1)] END FUNCTION OJ::GET_Polygons, LEVELS = levels, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge ; Construct IDLgrModel array for the Polygons represnting a grid ; I) Rough Check of INPUT data: ; IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF (MAX(Left_edge) GT 1.D) THEN BEGIN print, 'OJ: Left_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Left_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Polygons() : Left_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_Polygons() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF (MAX(Right_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_Polygons() : Right_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Right_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Polygons() : Right_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_Polygons() : Warning: Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_Polygons() : Right_edge is not larger than Left_edge' RETURN,0 ENDIF ENDIF ELSE Right_edge = self.Right_edge IF (N_ELEMENTS(LEVELS) GT 0) THEN BEGIN IF (max(LEVELS) GT self.max_level) THEN BEGIN print, 'OJ::Get_Polygons() : requested level larger than ', self.max_level RETURN, 0 ENDIF IF (min(LEVELS) GT self.min_level) THEN BEGIN print, 'OJ::Get_Polygons() : requested level smaller than ', self.min_level RETURN, 0 ENDIF ENDIF ELSE $ levels = indgen(self.max_level-self.min_level+1)+self.min_level ; II) Define a color scheme for polygons gout_colors = bytarr(40,3) gout_colors = $ [[255,255,255],[ 0,255,255],[ 0, 0,255],[255,255, 0], $ [255, 0, 0],[255, 0,255],[ 0, 0,125],[ 0,125,125], $ [125,125, 0],[125, 0, 0],[125, 0,125],[125,125,125], $ [255,255,255],[ 0,255,255],[ 0, 0,255],[255,255, 0], $ [255, 0, 0],[255, 0,255],[ 0, 0,125],[ 0,125,125], $ [125,125, 0],[125, 0, 0],[125, 0,125],[125,125,125], $ [255,255,255],[ 0,255,255],[ 0, 0,255],[255,255, 0], $ [255, 0, 0],[255, 0,255],[ 0, 0,125],[ 0,125,125], $ [125,125, 0],[125, 0, 0],[125, 0,125],[125,125,125], $ [255,255,255],[ 0,255,255],[ 0, 0,255],[255,255, 0]] ; III) Create Polygons nl = N_ELEMENTS(levels) oGridPolys = OBJARR(nl) FOR i = 0,nl-1 DO BEGIN oGridPolys[i] = OBJ_NEW('IDLgrModel') ENDFOR ng = self.num_of_grids grid_polys = OBJARR(ng,4) tgi = self.gi[0:ng] ci = 0 FOR i = 0,ng-1 DO BEGIN IF (vec_GE_vec(tgi[i].Left_edge, Left_edge) AND $ vec_LE_vec(tgi[i].Right_edge,Right_edge) AND $ (total(levels EQ tgi[i].level) GE 0)) THEN BEGIN lx = tgi[i].Left_edge[0] ly = tgi[i].Left_edge[1] lz = tgi[i].Left_edge[2] rx = tgi[i].Right_edge[0] ry = tgi[i].Right_edge[1] rz = tgi[i].Right_edge[2] goutcol = gout_colors(*,tgi[i].level) grid_polys[i,0] = OBJ_NEW('IDLgrPolygon', $ [lx,lx,lx,lx], $ [ly,ry,ry,ly], $ [lz,lz,rz,rz], $ COLOR = goutcol, SHADING = 1, style = 1) grid_polys[i,1] = OBJ_NEW('IDLgrPolygon', $ [rx,rx,rx,rx], $ [ly,ry,ry,ly], $ [lz,lz,rz,rz], $ COLOR = goutcol, SHADING = 1, style = 1) grid_polys[i,2] = OBJ_NEW('IDLgrPolygon', $ [rx,lx,lx,rx], $ [ry,ry,ly,ly], $ [lz,lz,lz,lz], $ COLOR = goutcol, SHADING = 1, style = 1) grid_polys[i,3] = OBJ_NEW('IDLgrPolygon', $ [rx,lx,lx,rx], $ [ry,ry,ly,ly], $ [rz,rz,rz,rz], $ COLOR = goutcol, SHADING = 1, style = 1) level_ind = where(levels EQ tgi[i].level) FOR j = 0,3 DO $ oGridPolys[level_ind[0]] -> Add, grid_polys[i,j] ci = ci + 1 ENDIF ENDFOR RETURN, oGridPolys END FUNCTION OJ::GET_Volumes, MIN_LEVEL = min_level, MAX_LEVEL = max_level,$ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ RGB_COLORS = rgb ; Construct IDLgrModel including 3D Volumes on specified (or all) levels ; I) Rough Check of INPUT data: ; IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF (MAX(Left_edge) GT 1.D) THEN BEGIN print, 'OJ: Left_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Left_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Volumes() : Left_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_Volumes() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF (MAX(Right_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_Volumes() : Right_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Right_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Volumes() : Right_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_Volumes() : Warning: Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_Volumes() : Right_edge is not larger than Left_edge' RETURN,0 ENDIF ENDIF ELSE Right_edge = self.Right_edge IF (N_ELEMENTS(MIN_LEVEL) GT 0) THEN BEGIN IF (MIN_LEVEL GT self.max_level) THEN BEGIN print, 'OJ::Get_Volumes() : requested minimum level larger than ', self.max_level print, ' will use', self.max_level, ' as minimum level.' min_level = self.max_level ENDIF IF (MIN_LEVEL LT self.min_level) THEN BEGIN print, 'OJ::Get_Volumes() : requested minimum level smaller than ', self.min_level print, ' will use', self.min_level, ' as minimum level.' min_level = self.min_level ENDIF ENDIF ELSE min_level = self.min_level IF (N_ELEMENTS(MAX_LEVEL) GT 0) THEN BEGIN IF (MAX_LEVEL GT self.max_level) THEN BEGIN print, 'OJ::Get_Volumes() : requested maximum level larger than ', self.max_level print, ' will use', self.max_level, ' as maximum level.' max_level = self.max_level ENDIF IF (MAX_LEVEL LT self.min_level) THEN BEGIN print, 'OJ::Get_Volumes() : requested minimum level smaller than ', self.min_level print, ' will use', self.min_level, ' as minimum level.' max_level = self.min_level ENDIF ENDIF ELSE max_level = self.max_level set_default_rgb = 1 IF (N_ELEMENTS(rgb) GT 0) THEN BEGIN hs = size(rgb) IF total(hs[0:2] NE [2,256,3]) THEN BEGIN print, 'OJ::Get_Volumes() : something wrong with supplied rgb table', $ ' Using default rgb table instead' ENDIF ELSE set_default_rgb = 0 ENDIF IF set_default_rgb THEN BEGIN rgb = bytarr(256,3) rgb[*,0]=BYTSCL(255-FIX(indgen(256)^2)) rgb[*,1]=90 rgb[128:*,2]=2*indgen(128) ; rgb[*,0]=0 ; rgb[*,1]=0 ; rgb[*,2]=indgen(256) ENDIF ; III) Create Volumes nl = (max_level-min_level+1) oGridVolumes = OBJ_NEW('IDLgrModel') ng = self.num_of_grids grid_dens = OBJARR(ng) tgi = self.gi[0:ng] ci = 0 parents = tgi.parent i = 0 WHILE ((tgi[i].level LE max_level) AND (i LT ng) ) DO BEGIN IF ((tgi[i].level GE min_level) AND $ vec_GE_vec(tgi[i].Left_edge, Left_edge) AND $ vec_LE_vec(tgi[i].Right_edge,Right_edge)) THEN BEGIN read_hdf, tgi[i].baryon_file, 1, data pi = where(parents EQ tgi[i].num) IF (pi[0] GE 0) THEN $ FOR j = 0,N_ELEMENTS(pi)-1 DO $ IF ((tgi[pi[j]].level LE max_level) AND $ (tgi[pi[j]].level GE min_level)) THEN BEGIN ts = tgi[pi[j]].parent_s_index te = tgi[pi[j]].parent_e_index data[ts[0]:te[0],ts[1]:te[1],ts[2]:te[2]] = 1.e-10 ; print, 'zeroing:',ts,te,' of grid:',tgi[i].num, ' because of ', tgi[pi[j]].num ENDIF grid_dens[i] = OBJ_NEW('IDLgrVolume',$ bytscl((data), min = 3.e-2,max = 1.),$ /NO_COPY) ; bytscl((data/2.^(tgi[i].level-min_level)),$ ; min = .05,max = 30.), /NO_COPY) ; bytscl(alog10(data/2.^(tgi[i].level-min_level)),$ ; min = -1.,max = 1.5), /NO_COPY) cc = (tgi[i].Right_edge-tgi[i].Left_edge)/$ DOUBLE(tgi[i].End_index-tgi[i].Start_index+1) ; print, i, tgi[i].level, cc op = indgen(256) op = FIX(FLOAT(op)/2.^(tgi[i].level-min_level)) grid_dens[i] -> SetProperty, RENDER_STEP = [1,1,1], $ XCOORD_CONV=[tgi[i].Left_edge[0],cc[0]],$ YCOORD_CONV=[tgi[i].Left_edge[1],cc[1]],$ ZCOORD_CONV=[tgi[i].Left_edge[2],cc[2]], $ COMPOSITE_FUNCTION = 0, ZERO_OPACITY_SKIP = 1, $ RGB_TABLE0 = rgb , OPACITY_TABLE0 = op , $ /INTERPOLATE oGridVolumes -> Add, grid_dens[i] ci = ci+1 ENDIF i = i+1 ENDWHILE print, 'OJ::Get_Volumes() : included volume data from ', ci, 'grids.' RETURN, oGridVolumes END FUNCTION OJ::GET_Slice_data, slice_coord, $ MIN_LEVEL = min_level, MAX_LEVEL = max_level, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ LEVELS=Levels,DATA_FIELDS=data_fields ; Construct slice data on specified (or all) levels slice = {slice, $ min_level: 0B, $ max_level: 0B, $ Np_total: 0L, $ N_dfields: 0B, $ slice_coord:[0D,0D,0D], $ LEFT_EDGE: [0.D, 0.D, 0.D], $ RIGHT_EDGE: [1.D, 1.D, 1.D], $ DATA_FIELDS: BYTARR(100), $ X_PTR: PTR_NEW(), $ Y_PTR: PTR_NEW(), $ Z_PTR: PTRARR(100), $ LEVEL_PTR: PTR_NEW(), $ dx_LEVEL: DBLARR(100) } ; I) Rough Check of INPUT data: ; ; check slice_coord to be 3 element array: IF (N_ELEMENTS(slice_coord) NE 3) THEN BEGIN print, 'OJ: call to Get_Slice_Data with incorect slice_coord:',$ slice_coord RETURN,0 ENDIF this_ind = WHERE(slice_coord GT 1.D-30) ; check slice_coord to have only one non-zero component IF (N_ELEMENTS(this_ind) NE 1) THEN BEGIN print, 'OJ: call to Get_Slice_Data with incorect slice_coord:',$ slice_coord RETURN,0 ENDIF IF (this_ind[0] EQ -1) THEN BEGIN print, 'OJ: call to Get_Slice_Data with incorect slice_coord:',$ slice_coord RETURN,0 ENDIF ; check region definitions: IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF (MAX(Left_edge) GT 1.D) THEN BEGIN print, 'OJ: Left_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Left_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Slice_Data() : Left_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_Slice_Data() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF (MAX(Right_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_Slice_Data() : Right_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Right_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Slice_Data() : Right_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_Slice_Data() : Warning:\n', $ 'Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_Slice_Data() : Right_edge is not larger than Left_edge' RETURN,0 ENDIF ENDIF ELSE Right_edge = self.Right_edge ; check Level selection N_levels_requested = N_ELEMENTS(Levels) IF (N_levels_requested GT 0) THEN BEGIN IF (MAX(Levels) GT self.max_level) THEN BEGIN print, 'OJ::Get_Slice_Data() : Levels :', Levels, $ 'contain entry exceeding self.max_level:', self.max_level RETURN, 0 ENDIF IF (MIN(Levels) LT self.min_level) THEN BEGIN print, 'OJ::Get_Slice_Data() : Levels :', Levels, $ 'contain entry less than self.min_level:', self.min_level RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Levels)) NE N_levels_requested) THEN BEGIN print, 'OJ::Get_Slice_Data() : Levels :', Levels, $ 'are not unique!' RETURN, 0 ENDIF ENDIF ELSE BEGIN levels = INDGEN(self.max_level-self.min_level+1)+self.min_level ENDELSE Levels = Levels[sort(Levels)] ; sort them, just in case N_df_requested = N_ELEMENTS(data_fields) IF (N_df_requested GT 0) THEN BEGIN IF (MAX(data_fields) GT self.gi[0].num_baryon_fields) THEN BEGIN print, 'OJ::Get_Slice_Data() : Data_fields :', data_fields, $ 'warning: contain entry exceeding self.num_baryon_fields:', self.gi[0].num_baryon_fields ; RETURN, 0 ENDIF IF (MIN(data_fields) LT 0) THEN BEGIN print, 'OJ::Get_Slice_Data() : Data_fields :', data_fields, $ 'contain entry less than 0' RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Data_fields)) NE N_df_requested) THEN BEGIN print, 'OJ::Get_Slice_Data() : Data_fields :', data_fields, $ 'were not unique!'+'... fixed!' Data_fields = Data_fields(UNIQ(Data_fields)) ENDIF ENDIF ELSE BEGIN ; get everything then Data_fields = INDGEN(self.gi[0].num_baryon_fields) ENDELSE axes = (slice_coord gt 1.D-30) const_sub = where(axes eq max(axes)) const_sub = const_sub(0) other_subs = where(axes eq 0) xi = other_subs[0] & yi = other_subs[1] & zi = const_sub x_stuetz = DBLARR(1000L*1000L) y_stuetz = x_stuetz z_stuetz = DBLARR(N_df_requested,1000L*1000L) c_level = BYTARR(1000L*1000L) dx_level = DBLARR(100) ; III) Create Slice max_level = MAX(Levels, MIN=min_level) ; max and minimum levels nl = N_levels_requested ; number of requested levels to be included ng = self.num_of_grids ; total number of grids available tgi = self.gi[0:(ng-1)] ; grid info ci = 0L i = 0L run_i = -1L ; running index of data points print, "OJ::Get_Slice_Data(): reading data files and interpolating slice data values ..." ; select grids to use ; first the ones that are in the level list: use_g = tgi[where_is_in(tgi.level,Levels)] ; now the ones that contain the slice index_range_a = use_g.End_index-use_g.Start_index+1 d_dis = (use_g.Right_edge - use_g.Left_edge)/ $ (use_g.End_index-use_g.Start_index+1) c_dis = d_dis[zi,*] ; ind = where(((use_g[*].Left_edge[zi]+0.5D*c_dis[*]) LT slice_coord[zi]) AND $ ; ((use_g[*].Right_edge[zi]-0.5D*c_dis[*]) GT slice_coord[zi])) ind = lindgen(N_elements(use_g)) if (N_elements(d_dis[*,0]) gt 2) then begin ind = where(((use_g[*].Left_edge[zi]) LT slice_coord[zi]) AND $ ((use_g[*].Right_edge[zi]) GT slice_coord[zi])) use_g =use_g(ind) endif ung = N_ELEMENTS(use_g) x_dis = d_dis(xi,ind[*]) y_dis = d_dis(yi,ind[*]) c_dis = c_dis(ind[*]) index_range_a[*,0:(ung-1)] = index_range_a(*,ind[*]) print, ' using', ung,' grids from Levels ', Levels WHILE (i LT ung) DO BEGIN ; now what is its dx cg = use_g[i] index_range = index_range_a[*,i] dx = x_dis[i] dy = y_dis[i] ; <- yes these are necessary because of unequal round off errors dz = c_dis[i] lind = where(cg.level EQ Levels) lco = (slice_coord[zi] - cg.Left_edge[zi] - 0.5D * dx)/dx > 0.D rco = (slice_coord[zi] - cg.Left_edge[zi] + 0.5D * dy)/dy ; IF ((lco GE 0.D) AND (rco LE (DOUBLE(index_range[zi])-0.0D))) $ ; two indeces for interpolation? ; THEN BEGIN ; ok this grid will be used ci = ci+1 ; which indeces of the constant component are nearest: ic_l = FLOOR(lco) ic_r = FLOOR(rco) < (cg.end_index[zi] - cg.start_index[zi]) data = fltarr(N_df_requested, index_range[0],index_range[1],index_range[2]) for df=0,N_df_requested-1 do begin hdat = 0. read_hdf, cg.baryon_file, Data_fields[df], hdat ;, /verbose data[df, *, *, *] = hdat[*,*,*] endfor ; list of grids that have the current grid as parent: dx_level[cg.level] = dx pi = where(use_g.parent EQ cg.num, nchild) ; Zero the data where children exist for ic = 0, nchild-1 do begin ii = [use_g[pi[ic]].parent_s_index[0], $ use_g[pi[ic]].parent_e_index[0]] jj = [use_g[pi[ic]].parent_s_index[1], $ use_g[pi[ic]].parent_e_index[1]] kk = [use_g[pi[ic]].parent_s_index[2], $ use_g[pi[ic]].parent_e_index[2]] data[*, ii[0]:ii[1], jj[0]:jj[1], kk[0]:kk[1]] = 0.0 endfor ;; END zero cells "ic" ; now loop over both other indeces, compute x_stuetz, y_stuetz and interpolate z_stuetz FOR j=0,index_range[xi]-1 DO $ FOR k=0,index_range[yi]-1 DO BEGIN ; scaled coords in [0.,1] xco = cg.Left_edge[xi] + (0.5D + DOUBLE(j))*dx yco = cg.Left_edge[yi] + (0.5D + DOUBLE(k))*dy ; now make sure there are no finer grids covering this region ; (only from the selected grids -- where data ne 0) cp = DBLARR(3) cp[other_subs] = [xco, yco] cp[zi] = slice_coord[zi] ldi=INTARR(3) ldi[other_subs] = [j,k] ldi[zi ] = [ic_l] rdi=ldi rdi[zi ] = [ic_r] if (data[0, ldi[0], ldi[1], ldi[2]] ne 0 and $ data[0, rdi[0], rdi[1], rdi[2]] ne 0) then begin run_i = run_i + 1 x_stuetz[run_i] = cp[xi] y_stuetz[run_i] = cp[yi] for df=0,N_df_requested-1 do begin ld = data[df, ldi[0],ldi[1],ldi[2]] rd = data[df, rdi[0],rdi[1],rdi[2]] ; linear interpolation to get correct data value z_stuetz[df, run_i] = ld + (rd-ld)/(rco-lco)* dz * $ (cp[zi]-lco*dz) endfor c_level[run_i] = BYTE(cg.level) ENDIF ENDFOR ; ENDIF i = i+1 ENDWHILE ; reduce arrays: ; x_stuetz = TEMPORARY(x_stuetz[0:run_i]) ; y_stuetz = TEMPORARY(y_stuetz[0:run_i]) ; c_level = TEMPORARY(c_level [0:run_i]) nonzero = where(z_stuetz[0,*] ne 0, nvalid) if (nvalid gt 0) then begin x_stuetz = (x_stuetz[nonzero]) y_stuetz = (y_stuetz[nonzero]) c_level = (c_level [nonzero]) endif slice.min_level = min_level slice.max_level = max_level slice.Np_total = nvalid slice.N_dfields = N_df_requested slice.slice_coord = slice_coord slice.Left_edge = Left_edge slice.Right_edge = Right_edge slice.Data_fields= Data_fields slice.x_ptr = PTR_NEW(x_stuetz) slice.y_ptr = PTR_NEW(y_stuetz) for df=0,N_df_requested-1 do begin z_s = reform(z_stuetz[df, nonzero]) slice.z_ptr[df] = PTR_NEW(z_s) endfor slice.Level_ptr= PTR_NEW(c_Level) slice.dx_level = dx_level print, 'OJ::Get_Slice_Data() : included ', nvalid,' (', $ sqrt(FLOAT(nvalid)), '^2) data points from ', ci, ' grids.' RETURN, slice END ;-------------------------------------------------------------------------------- ; Written by John H. Wise ; Reads particle data ; Created: October 23, 2003 function OJ::GET_Particle_Data, $ MIN_LEVEL = min_level, MAX_LEVEL = max_level, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ LEVELS=Levels,DATA_FIELDS=data_fields, $ STARS = stars struct = {part, $ min_level: 0B, $ max_level: 0B, $ Np_total: PTR_NEW(), $ N_dfields: 0B, $ LEFT_EDGE: [0.D, 0.D, 0.D], $ RIGHT_EDGE: [1.D, 1.D, 1.D], $ DATA_FIELDS: BYTARR(100), $ level: 0, $ X_PTR: PTR_NEW(), $ Y_PTR: PTR_NEW(), $ Z_PTR: PTR_NEW(), $ D_PTR: PTRARR(100) } ; I) Rough Check of INPUT data: ; ; If there's star formation, have 3 extra particle attributes if keyword_set(stars) then num_pattr = 11 else num_pattr = 8 ; check region definitions: IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF (MAX(Left_edge) GT 1.D) THEN BEGIN print, 'OJ: Left_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Left_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Particle_Data() : Left_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_Particle_Data() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF (MAX(Right_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_Particle_Data() : Right_edge contains component greater 1.' RETURN,0 ENDIF IF (MIN(Right_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Particle_Data() : Right_edge contains component smaller than 0..' RETURN,0 ENDIF IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_Particle_Data() : Warning:\n', $ 'Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_Particle_Data() : Right_edge is not larger than Left_edge' RETURN,0 ENDIF ENDIF ELSE Right_edge = self.Right_edge ; check Level selection N_levels_requested = N_ELEMENTS(Levels) IF (N_levels_requested GT 0) THEN BEGIN IF (MAX(Levels) GT self.max_level) THEN BEGIN print, 'OJ::Get_Particle_Data() : Levels :', Levels, $ 'contain entry exceeding self.max_level:', self.max_level RETURN, 0 ENDIF IF (MIN(Levels) LT self.min_level) THEN BEGIN print, 'OJ::Get_Particle_Data() : Levels :', Levels, $ 'contain entry less than self.min_level:', self.min_level RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Levels)) NE N_levels_requested) THEN BEGIN print, 'OJ::Get_Particle_Data() : Levels :', Levels, $ 'are not unique!' RETURN, 0 ENDIF ENDIF ELSE BEGIN levels = INDGEN(self.max_level-self.min_level+1)+self.min_level ENDELSE Levels = Levels[sort(Levels)] ; sort them, just in case N_df_requested = N_ELEMENTS(data_fields) IF (N_df_requested GT 0) THEN BEGIN IF (MAX(data_fields) GT self.gi[0].num_baryon_fields) THEN BEGIN print, 'OJ::Get_Particle_Data() : Data_fields :', data_fields, $ 'warning: contain entry exceeding self.num_particle_attr:', self.gi[0].num_particle_attr ; RETURN, 0 ENDIF IF (MIN(data_fields) LT 0) THEN BEGIN print, 'OJ::Get_Particle_Data() : Data_fields :', data_fields, $ 'contain entry less than 0' RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Data_fields)) NE N_df_requested) THEN BEGIN print, 'OJ::Get_Particle_Data() : Data_fields :', data_fields, $ 'were not unique!'+'... fixed!' Data_fields = Data_fields(UNIQ(Data_fields)) ENDIF ENDIF ELSE BEGIN ; get everything then Data_fields = INDGEN(num_pattr) + 1 ENDELSE max_level = max(levels, min=min_level) ; Max and min levels nl = N_levels_requested ; Number of requested levels to be included ng = self.num_of_grids ; Total number of grids available tgi = self.gi[0:(ng-1)] ; Grid Info part = replicate(struct, ng) print, "OJ:Get_Particle_Data(): reading data files..." ; Look for the ones that are in the level list use_g = tgi[where_is_in(tgi.level,levels)] d_dis = (use_g.Right_edge - use_g.Left_edge)/ $ (use_g.End_index-use_g.Start_index+1) ; Number of grids ung = n_elements(use_g) print, ' using', ung, ' grids from levels ', levels ; Read the Particle Data pdat = get_particle_data( use_g[*].baryon_file, Data_fields, $ self.gi[0].num_baryon_fields+2);, /verbose ii = 0 ; grid counter while (ii lt ung) do begin part[ii].Np_total = pdat[ii].Np_total part[ii].x_ptr = pdat[ii].d_ptr[0] part[ii].y_ptr = pdat[ii].d_ptr[1] part[ii].z_ptr = pdat[ii].d_ptr[2] part[ii].level = use_g[ii].level for jj = 3, num_pattr do begin part[ii].d_ptr[jj-3] = pdat[ii].d_ptr[jj] endfor ii = ii + 1 endwhile part.min_level = min_level part.max_level = max_level part.N_dfields = N_df_requested part.left_edge = left_edge part.right_edge = right_edge part.data_fields = data_fields return, part end ;-------------------------------------------------------------------------------- function OJ::nextavailabledatacontainer return, min(where(self.data_fields eq '')) end function OJ::dataindexfromname, field return, (where(self.data_fields eq field))[0] end pro OJ::r3D, $ MIN_LEVEL = min_level, MAX_LEVEL = max_level, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ LEVELS=Levels,DATA_FIELDS=data_fields, $ CENTER=center, VIEWDIM=viewdim, CUBE=cube ; Retrieve all grid data on specified (or all) levels threeD = {threeD, $ min_level: 0B, $ max_level: 0B, $ Np_total: 0L, $ N_dfields: 0B, $ LEFT_EDGE: [0.D, 0.D, 0.D], $ RIGHT_EDGE: [1.D, 1.D, 1.D], $ DATA_FIELDS: STRARR(100), $ XP: PTR_NEW(), $ ; x coord pointer YP: PTR_NEW(), $ ZP: PTR_NEW(), $ DP: PTRARR(100), $ ; data pointer LP: PTR_NEW(), $ ; level pointer dxl: DBLARR(100) } ; dx of level ; I) Rough Check of INPUT data: ; IF (N_ELEMENTS(DATA_DIM) EQ 0) THEN data_dim=3 ; check region definitions: IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF (MAX(Left_edge) GT 1.D) THEN BEGIN print, 'OJ: Left_edge contains component greater 1.' RETURN ENDIF IF (MIN(Left_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_3D_Data() : Left_edge contains component smaller than 0..' RETURN ENDIF IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_3D_Data() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF (MAX(Right_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_3D_Data() : Right_edge contains component greater 1.' RETURN ENDIF IF (MIN(Right_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_3D_Data() : Right_edge contains component smaller than 0..' RETURN ENDIF IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_3D_Data() : Warning:\n', $ 'Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_3D_Data() : Right_edge is not larger than Left_edge' RETURN ENDIF ENDIF ELSE Right_edge = self.Right_edge ; check Level selection N_levels_requested = N_ELEMENTS(Levels) IF (N_levels_requested GT 0) THEN BEGIN IF (MAX(Levels) GT self.max_level) THEN BEGIN print, 'OJ::Get_3D_Data() : Levels :', Levels, $ 'contain entry exceeding self.max_level:', self.max_level RETURN ENDIF IF (MIN(Levels) LT self.min_level) THEN BEGIN print, 'OJ::Get_3D_Data() : Levels :', Levels, $ 'contain entry less than self.min_level:', self.min_level RETURN ENDIF IF (N_ELEMENTS(UNIQ(Levels)) NE N_levels_requested) THEN BEGIN print, 'OJ::Get_3D_Data() : Levels :', Levels, $ 'are not unique!' RETURN ENDIF ENDIF ELSE BEGIN levels = INDGEN(self.max_level-self.min_level+1)+self.min_level ENDELSE Levels = Levels[sort(Levels)] ; sort them, just in case N_df_requested = N_ELEMENTS(data_fields) IF (N_df_requested GT 0) THEN BEGIN IF (N_ELEMENTS(UNIQ(Data_fields)) NE N_df_requested) THEN BEGIN print, 'OJ::Get_3D_Data() : Data_fields :', data_fields, $ 'were not unique!'+'... fixed!' Data_fields = Data_fields(UNIQ(Data_fields)) N_df_requested = N_ELEMENTS(data_fields) ENDIF ENDIF ELSE BEGIN print, 'OJ::g3D called without data_fields to read.' ENDELSE if N_ELEMENTS(viewdim) ne 1 then viewdim=0 if N_elements(center) ne 3 then center = [0.5,0.5,0.5] ; III) Create 3D data max_level = MAX(Levels, MIN=min_level) ; max and minimum levels nl = N_levels_requested ; number of requested levels to be included ng = self.num_of_grids ; total number of grids available tgi = self.gi[0:(ng-1)] ; grid info ii = 0 run_i = 0L ; running index of data points print, "OJ::Get_3D_Data(): reading data files ..." ; select grids to use ; first the ones that are in the level list: use_g = tgi[where_is_in(tgi.level,Levels)] if N_elements(cube) ge 1 then begin print, 'selecting grids only in cube' ocube = cube cube = fltarr(2) cube[0] = ocube[*] if N_elements(ocube) lt 2 then cube[1] = 256. if total(use_g[0].dim gt 1) eq 3 then $ ; select if 3D otherwise use_g = grid_in_cube(use_g, center[0:2], cube[0], FIX(cube[1]), /noresolution) endif d_dis = (use_g.Right_edge - use_g.Left_edge)/ $ (use_g.End_index-use_g.Start_index+1) ; number of grids ung = N_ELEMENTS(use_g) print, ' using', ung,' grids from Levels ', use_g[uniq(use_g.Level)].level ; read all data dummy = data_fields if n_elements(cube) ge 1 then begin lef = center-cube[0]/2. rig = center+cube[0]/2. hdat = get_data(use_g, data_fields, viewdim, LEFT_EDGE=lef, RIGHT_EDGE=rig,/return_as_list) endif else hdat = get_data(use_g, data_fields, viewdim,/return_as_list) n_df_requested = N_elements(data_fields) ; how many data points in total? tot_cells = hdat[0].np_total ; allocate some memory if n_df_requested ge 1 then d_stuetz = DBLARR(N_df_requested,tot_cells) data_dim = N_ELEMENTS(where(self.gi[0].dim gt 0)) x_dis = d_dis[0,*] for df=0,n_df_requested-1 do d_stuetz[df,*] = (*hdat[0].data[df])[*] ; reduce arrays: for df=0,N_df_requested-1 do begin exi = (where(self.data_fields eq data_fields[df]))[0] if exi eq -1 then begin in = self->OJ::nextavailabledatacontainer() self.data_fields[in] = data_fields[df] self.data[in] = ptr_new(d_stuetz[df, *]) endif else begin ptr_free, self.data[exi] self.data[exi] = ptr_new(d_stuetz[df, *]) ENDELSE endfor print, 'OJ::Get_3D_Data() : included ', tot_cells,' (',(FLOAT(tot_cells))^(0.3334),$ '^3) data points from ', ung, ' grids.' RETURN END function OJ::data3D, field if N_elements(field) ne 1 then begin print, 'OJ::data called with a request for more than one field' print, 'returning 0.' return, 0 endif in = self->OJ::dataindexfromname(field) if in ne -1 then begin if ptr_valid(self.data[in]) then return, (*self.data[in]) ; if data already known then return pointer to it endif else self->OJ::r3D, DATA_FIELDS=field ; else read it first return, (*self.data[self->OJ::dataindexfromname(field)]) end pro OJ::delete3Dcache ind = where(self.data_fields ne '', count) if count gt 0 then begin for i=0,count-1 do begin ptr_free, self.data[i] self.data[i] = ptr_new() self.data_fields[i] = '' endfor endif return end function OJ::Get_Property, DATA_FIELDS=data_fields if KEYWORD_SET(data_fields) then return, self.data_fields end FUNCTION OJ::GET_2D_data, $ MIN_LEVEL = min_level, MAX_LEVEL = max_level, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ LEVELS=Levels,DATA_FIELDS=data_fields ; Retrieve all grid data on specified (or all) levels twoD = {twoD, $ min_level: 0B, $ max_level: 0B, $ Np_total: 0L, $ N_dfields: 0B, $ LEFT_EDGE: [0.D, 0.D], $ RIGHT_EDGE: [1.D, 1.D], $ DATA_FIELDS: BYTARR(100), $ X_PTR: PTR_NEW(), $ Y_PTR: PTR_NEW(), $ D_PTR: PTRARR(100), $ LEVEL_PTR: PTR_NEW(), $ dx_LEVEL: DBLARR(100) } ; I) Rough Check of INPUT data: ; IF (N_ELEMENTS(DATA_DIM) EQ 0) THEN data_dim=2 ; check region definitions: IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_2D_Data() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_2D_Data() : Warning:\n', $ 'Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_2D_Data() : Right_edge is not larger than Left_edge' RETURN,0 ENDIF ENDIF ELSE Right_edge = self.Right_edge ; check Level selection N_levels_requested = N_ELEMENTS(Levels) IF (N_levels_requested GT 0) THEN BEGIN IF (MAX(Levels) GT self.max_level) THEN BEGIN print, 'OJ::Get_2D_Data() : Levels :', Levels, $ 'contain entry exceeding self.max_level:', self.max_level RETURN, 0 ENDIF IF (MIN(Levels) LT self.min_level) THEN BEGIN print, 'OJ::Get_2D_Data() : Levels :', Levels, $ 'contain entry less than self.min_level:', self.min_level RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Levels)) NE N_levels_requested) THEN BEGIN print, 'OJ::Get_2D_Data() : Levels :', Levels, $ 'are not unique!' RETURN, 0 ENDIF ENDIF ELSE BEGIN levels = INDGEN(self.max_level-self.min_level+1)+self.min_level ENDELSE Levels = Levels[sort(Levels)] ; sort them, just in case N_df_requested = N_ELEMENTS(data_fields) IF (N_df_requested GT 0) THEN BEGIN IF (MAX(data_fields) GT self.gi[0].num_baryon_fields) THEN BEGIN print, 'OJ::Get_2D_Data() : Data_fields :', data_fields, $ 'warning: contain entry exceeding self.num_baryon_fields:', self.gi[0].num_baryon_fields ; RETURN, 0 ENDIF IF (MIN(data_fields) LT 0) THEN BEGIN print, 'OJ::Get_2D_Data() : Data_fields :', data_fields, $ 'contain entry less than 0' RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Data_fields)) NE N_df_requested) THEN BEGIN print, 'OJ::Get_2D_Data() : Data_fields :', data_fields, $ 'were not unique!'+'... fixed!' Data_fields = Data_fields(UNIQ(Data_fields)) ENDIF ENDIF ELSE BEGIN ; get everything then Data_fields = INDGEN(self.gi[0].num_baryon_fields) ENDELSE ; III) Create 2D data max_level = MAX(Levels, MIN=min_level) ; max and minimum levels nl = N_levels_requested ; number of requested levels to be included ng = self.num_of_grids ; total number of grids available tgi = self.gi[0:(ng-1)] ; grid info ci = 0L ii = 0L run_i = 0L ; running index of data points print, "OJ::Get_2D_Data(): reading data files ..." ; select grids to use ; first the ones that are in the level list: use_g = tgi[where_is_in(tgi.level,Levels)] d_dis = (use_g.Right_edge - use_g.Left_edge)/(use_g.End_index-use_g.Start_index+1) ; number of grids ung = N_ELEMENTS(use_g) print, ' using', ung,' grids from Levels ', Levels ; read all data read_all_grid_hdf, use_g[*].baryon_file, Data_fields, hdat, /verbose ; how many dta points in total? tot_cells = TOTAL(hdat.np_total) ; allocate some memory x_stuetz = DBLARR(tot_cells) y_stuetz = x_stuetz d_stuetz = DBLARR(N_df_requested,tot_cells) c_level = BYTARR(tot_cells) dx_level = DBLARR(100) data_dim = N_ELEMENTS(where(self.gi[0].dim gt 0)) if data_dim ne 2 then begin print, 'AMRData::get_2d_data(): data_dim not 2D! ' ; breaki endif x_dis = d_dis[0,*] y_dis = d_dis[1,*] WHILE (ii LT ung) DO BEGIN ; now what is its dx cg = use_g[ii] dx = x_dis[ii] dy = y_dis[ii] ci = ci+1 dx_level[cg.level] = dx ; list of grids that have the current grid as parent: pi = where(use_g.parent EQ cg.num) npi = N_ELEMENTS(pi) ; mark regions covered by finer grids if (pi[0] ge 0) then $ for cnt=0,npi-1 DO (*hdat[ii].d_ptr[0])[ $ use_g(pi[cnt]).parent_s_index[0]:use_g(pi[cnt]).parent_e_index[0], $ use_g(pi[cnt]).parent_s_index[1]:use_g(pi[cnt]).parent_e_index[1]] = -1.d50 ; derive x, y coordinates of cells d = hdat[ii].dims in = indgen(hdat[ii].np_total, /long) i = in mod d[0] j = (in - i)/d[0] xco = cg.Left_edge[0] + (0.5D + DOUBLE(i))*dx yco = cg.Left_edge[1] + (0.5D + DOUBLE(j))*dy nump = hdat[ii].np_total x_stuetz[run_i:(run_i+nump-1)] = xco[*] y_stuetz[run_i:(run_i+nump-1)] = yco[*] ;breaki for df=0,n_df_requested-1 do $ d_stuetz[df,run_i:(run_i+nump-1)] = (*hdat[ii].d_ptr[df])[*] c_level[run_i:(run_i+nump-1)] = BYTE(cg.level) run_i = run_i + nump ii = ii+1L ; increment grid ENDWHILE ind = where(d_stuetz[0,0L:run_i-1L] gt -1.d50) overlap = run_i-N_elements(ind) print, overlap ,' points dismissed because of finer grids info available' ; reduce arrays: x_stuetz = TEMPORARY(x_stuetz[ind]) y_stuetz = TEMPORARY(y_stuetz[ind]) twoD.min_level = min_level twoD.max_level = max_level twoD.Np_total = N_elements(ind) twoD.N_dfields = N_df_requested twoD.Left_edge = Left_edge[0:1] twoD.Right_edge = Right_edge[0:1] twoD.Data_fields= Data_fields twoD.x_ptr = PTR_NEW(x_stuetz) twoD.y_ptr = PTR_NEW(y_stuetz) twoD.dx_level = dx_level z_s = fltarr(twoD.Np_total) for df=0,N_df_requested-1 do begin z_s[*] = d_stuetz[df, ind] twoD.d_ptr[df] = PTR_NEW(z_s) endfor twoD.Level_ptr= PTR_NEW(c_Level[ind]) print, 'OJ::Get_2D_Data() : included ', run_i,' (',(FLOAT(run_i))^(0.5),$ '^2) data points from ', ci, ' grids.' RETURN, twoD END FUNCTION OJ::Image_from_2D_data, $ coords, xysize, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ LEVELS=Levels,DATA_FIELDS=data_fields ; Retrieve all grid data on specified (or all) levels ; ; OJ->Image_from_2D_data([0.0, 0.0, 1., 1.], [512,512]) ; makes and 512x512 image of the domain images = {images, $ XSIZE: 0, $ YSIZE: 0, $ N_dfields: 0B, $ LEFT_EDGE: [0.D, 0.D], $ RIGHT_EDGE: [1.D, 1.D], $ DATA_FIELDS: BYTARR(100), $ D_PTR: PTRARR(100) $ } ; I) Rough Check of INPUT data: ; IF (N_ELEMENTS(DATA_DIM) EQ 0) THEN data_dim=2 ; check region definitions: IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF NOT vec_GE_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_2D_Data() : Warning: Left_edge is not larger or equal than Left_edge of OJ' ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF NOT vec_LE_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_2D_Data() : Warning:\n', $ 'Right_edge is not smaller or equal than Right_edge of OJ' ENDIF IF NOT vec_LT_vec(Left_edge,Right_edge) THEN BEGIN print, 'OJ::Get_2D_Data() : Right_edge is not larger than Left_edge' RETURN,0 ENDIF ENDIF ELSE Right_edge = self.Right_edge N_df_requested = N_ELEMENTS(data_fields) IF (N_df_requested GT 0) THEN BEGIN IF (MAX(data_fields) GT self.gi[0].num_baryon_fields) THEN BEGIN print, 'OJ::Get_2D_Data() : Data_fields :', data_fields, $ 'warning: contain entry exceeding self.num_baryon_fields:', self.gi[0].num_baryon_fields ; RETURN, 0 ENDIF IF (MIN(data_fields) LT 0) THEN BEGIN print, 'OJ::Get_2D_Data() : Data_fields :', data_fields, $ 'contain entry less than 0' RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Data_fields)) NE N_df_requested) THEN BEGIN print, 'OJ::Get_2D_Data() : Data_fields :', data_fields, $ 'were not unique!'+'... fixed!' Data_fields = Data_fields(UNIQ(Data_fields)) ENDIF ENDIF ELSE BEGIN ; get everything then Data_fields = INDGEN(self.gi[0].num_baryon_fields) ENDELSE coord=coords ; III) Create 2D data ; check Level selection npixels = xysize[0]*xysize[1] pixeldx = (coords[2]-coords[0])/float(xysize[0]) pixeldy = (coords[3]-coords[1])/float(xysize[1]) minpixeldx = min([coords[2]-coords[0], coords[3]-coords[1]])/float(max(xysize)) ng = self.num_of_grids ; total number of grids available tgi = self.gi[0:(ng-1)] ; grid info d_dis = (tgi.Right_edge - tgi.Left_edge)/(tgi.End_index-tgi.Start_index+1) ci = 0L ii = 0L run_i = 0L ; running index of data points print, "OJ::Image_from_2D_Data(): reading data files ..." ; select grids to use ; in image and large enough pixels to be visible ; assume grids are order from coarse to finest already ; ind = where((d_dis[0,*] ge minpixeldx) AND $ ((tgi[*].Left_edge[0]) LT coords[2]) AND $ ((tgi[*].Left_edge[1]) LT coords[3]) AND $ ((tgi[*].Right_edge[0]) GT coords[0]) AND $ ((tgi[*].Right_edge[1]) GT coords[1]) ) use_g = tgi[ind] d_dis = d_dis[*,ind] ; number of grids ung = N_ELEMENTS(use_g) print, ' using', ung,' grids with dx larger than ', minpixeldx ; read all data read_all_grid_hdf, use_g[*].baryon_file, Data_fields, hdat, /verbose ; how many data points in total? tot_cells = TOTAL(hdat.np_total) ; allocate some memory d_stuetz = DBLARR(N_df_requested,xysize[0],xysize[1]) x_dis = d_dis[0,*] y_dis = d_dis[1,*] WHILE (ii LT ung) DO BEGIN ; now what is its dx cg = use_g[ii] ci = ci+1 dx = x_dis[ii] dy = y_dis[ii] ; what is the first coordinate on image xl = coord[2]-coord[0] yl = coord[3]-coord[1] cgxl = cg.Right_edge[0]-cg.Left_edge[0] cgyl = cg.Right_edge[1]-cg.Left_edge[1] cgdxpixel = dx cgdypixel = dy scalex = cgdxpixel/pixeldx scaley = cgdypixel/pixeldy ; find cells of current grid that cintribute to image cimi = max([FLOOR((coord[0]-cg.Left_edge[0])/dx), 0]) cima = min([CEIL((coord[2]-cg.Left_edge[0])/dx), (cg.end_index[0]-cg.start_index[0])]) cjmi = max([FLOOR((coord[1]-cg.Left_edge[1])/dy), 0]) cjma = min([CEIL((coord[3]-cg.Left_edge[1])/dy), (cg.end_index[1]-cg.start_index[1])]) imi = max([ROUND((cg.Left_edge[0]-coord[0])/pixeldx), 0]) ima = min([ROUND((cg.Right_edge[0]-coord[0])/pixeldx), xysize[0]])-1 jmi = max([ROUND((cg.left_edge[1]-coord[1])/pixeldy), 0]) jma = min([ROUND((cg.right_edge[1]-coord[1])/pixeldy), xysize[1]])-1 ndx = cima-cimi+1 ndy = cjma-cjmi+1 slex = cg.left_edge[0] + float(cimi)*dx sley = cg.left_edge[1] + float(cjmi)*dy simi = max([ROUND((coord[0]-slex)/pixeldx), 0]) sima = min([ROUND((coord[2]-slex)/pixeldx), ROUND(ndx*scalex)])-1 sjmi = max([ROUND((coord[1]-sley)/pixeldy), 0]) sjma = min([ROUND((coord[3]-sley)/pixeldy), ROUND(ndy*scaley)])-1 for df=0,n_df_requested-1 do BEGIN stamp = (*hdat[ii].d_ptr[df])[cimi:cima,cjmi:cjma] image = congrid(stamp,ndx*scalex,ndy*scaley) d_stuetz[df,imi:ima,jmi:jma] = image[simi:sima,sjmi:sjma] ENDFOR ii = ii+1L ; increment grid ENDWHILE images.xsize = xysize[0] images.ysize = xysize[1] images.N_dfields = N_df_requested images.Left_edge = [coord[0],coord[2]] images.Right_edge = [coord[1],coord[3]] images.Data_fields= Data_fields image = fltarr(xysize[0], xysize[1]) for df=0,N_df_requested-1 do begin image[*,*] = d_stuetz[df, *,*] images.d_ptr[df] = PTR_NEW(image) endfor RETURN, images END FUNCTION OJ::GET_Orthogonal_Projection, oproj_coord, oproj_width, $ MIN_LEVEL = min_level, MAX_LEVEL = max_level, $ LEFT_EDGE = Left_edge, RIGHT_EDGE = Right_edge, $ LEVELS=Levels,DATA_FIELDS=data_fields ; Construct slice data on specified (or all) levels slice = {oproject, $ min_level: 0B, $ max_level: 0B, $ Np_total: 0L, $ N_dfields: 0B, $ oproj_coord:[0D,0D,0D], $ ; which axis midplane coord oproj_width: 0D, $ ; width of the projected slice in [0.,1.] LEFT_EDGE: [0.D, 0.D, 0.D], $ RIGHT_EDGE: [1.D, 1.D, 1.D], $ DATA_FIELDS: BYTARR(100), $ X_PTR: PTR_NEW(), $ Y_PTR: PTR_NEW(), $ Z_PTR: PTRARR(100), $ LEVEL_PTR: PTR_NEW(), $ dx_LEVEL: DBLARR(100) } ; I) Rough Check of INPUT data: ; ; check oproj_coord to be 3 element array: IF (N_ELEMENTS(oproj_coord) NE 3) THEN BEGIN print, 'OJ: call to Get_Orthogonal_Projection with incorect oproj_coord:',$ oproj_coord oproj_coord = [0.D, 0.5D, 0.D] print, 'using ', oproj_coord ENDIF this_ind = WHERE(oproj_coord GT 1.D-30) axes = (oproj_coord gt 1.D-30) const_sub = where(axes eq max(axes)) const_sub = const_sub(0) other_subs = where(axes eq 0) xi = other_subs[0] & yi = other_subs[1] & zi = const_sub ; check slice_coord to have only one non-zero component IF (N_ELEMENTS(this_ind) NE 1) THEN BEGIN print, 'OJ: call to Get_Orthogonal_Projection with incorect oproj_coord:',$ oproj_coord RETURN,0 ENDIF IF (this_ind[0] EQ -1) THEN BEGIN print, 'OJ: call to Get_Orthogonal_Projection with incorect oproj_coord:',$ oproj_coord RETURN,0 ENDIF IF (N_ELEMENTS(oproj_width) NE 1) THEN BEGIN print, 'OJ: call to Get_Orthogonal_Projection with incorect oproj_width:',$ oproj_coord oproj_width=1.D print, ' using', oproj_width, ' instead' ENDIF ; check region definitions: IF (N_ELEMENTS(Left_edge) NE 0) THEN BEGIN IF (MAX(Left_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection: Left_edge contains component greater 1.' Left_edge = self.Left_edge print, ' using:', self.Left_edge ENDIF IF (MIN(Left_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Left_edge contains component smaller than 0..' Left_edge = self.Left_edge print, ' using:', self.Left_edge ENDIF IF NOT vec_LT_vec(Left_edge,self.Left_edge) THEN BEGIN print, 'OJ::Get_Slice_Data() : Warning: Left_edge is smaller than Left_edge of OJ' print, Left_edge, self.left_edge ENDIF ENDIF ELSE Left_edge = self.Left_edge IF (N_ELEMENTS(Right_edge) NE 0) THEN BEGIN IF (MAX(Right_edge) GT 1.D) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection: Right_edge contains component greater 1.' Right_edge = self.Right_edge print, ' using:', self.Right_edge ENDIF IF (MIN(Right_edge) LT 0.D) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Right_edge contains component smaller than 0..' Right_edge = self.Right_edge print, ' using:', self.Right_edge ENDIF IF NOT vec_LT_vec(Right_edge,self.Right_edge) THEN BEGIN print, 'OJ::Get_Slice_Data() : Warning: Right_edge is larger than Right_edge of OJ' ENDIF ENDIF ELSE Right_edge = self.Right_edge ; check Level selection N_levels_requested = N_ELEMENTS(Levels) IF (N_levels_requested GT 0) THEN BEGIN IF (MAX(Levels) GT self.max_level) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Levels :', Levels, $ 'contain entry exceeding self.max_level:', self.max_level RETURN, 0 ENDIF IF (MIN(Levels) LT self.min_level) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Levels :', Levels, $ 'contain entry less than self.min_level:', self.min_level RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Levels)) NE N_levels_requested) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Levels :', Levels, $ 'are not unique!' RETURN, 0 ENDIF ENDIF ELSE BEGIN IF ((N_ELEMENTS(min_level) GE 0) AND (N_ELEMENTS(max_level) GE 0)) THEN BEGIN min_level = max([self.min_level, min_level]) max_level = min([self.max_level, max_level]) levels = INDGEN(max_level-min_level+1)+min_level ENDIF ELSE levels = INDGEN(self.max_level-self.min_level+1)+self.min_level ENDELSE Levels = Levels[sort(Levels)] ; sort them, just in case N_levels_requested = N_ELEMENTS(Levels) ; recount N_df_requested = N_ELEMENTS(data_fields) IF (N_df_requested GT 0) THEN BEGIN IF (MAX(data_fields) GT self.gi[0].num_baryon_fields) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Data_fields :', data_fields, $ 'contain entry exceeding self.num_baryon_fields:', self.gi[0].num_baryon_fields ; RETURN, 0 ENDIF IF (MIN(data_fields) LT 0) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Data_fields :', data_fields, $ 'contain entry less than 0' RETURN, 0 ENDIF IF (N_ELEMENTS(UNIQ(Data_fields)) NE N_df_requested) THEN BEGIN print, 'OJ::Get_Orthogonal_Projection() : Data_fields :', data_fields, $ 'were not unique!'+'... fixed!' Data_fields = Data_fields(UNIQ(Data_fields)) ENDIF ENDIF ELSE BEGIN ; get everything then Data_fields = INDGEN(self.gi[0].num_baryon_fields) ENDELSE Data_fields = Data_fields[sort(Levels)] ; sort them for faster access x_stuetz = DBLARR(1000L*1000L) ; x coordinates y_stuetz = x_stuetz ; y coordinates z_stuetz = x_stuetz ; data value w_stuetz = x_stuetz ; weight c_level = BYTARR(1000L*1000L) ; point from which level? dx_level = DBLARR(100) ; cell size on this level ; III) Create Projection max_level = MAX(Levels, MIN=min_level) ; max and minimum levels nl = N_levels_requested ; number of requested levels to be included ng = self.num_of_grids ; total number of grids available tgi = self.gi[0:(ng-1)] ; grid info print, "OJ::Get_Orthogonal_Projection(): reading data files and interpolating slice data values ..." ; select grids to use ; first the ones that are in the level list: use_g = tgi[where_is_in(tgi.level,Levels)] ; now the ones that contain the slice index_range_a = use_g.End_index-use_g.Start_index+1 d_dis = (use_g.Right_edge - use_g.Left_edge)/ $ (use_g.End_index-use_g.Start_index+1) c_dis = d_dis[zi,*] lep = oproj_coord[zi]-0.5D*oproj_width > 0.D rip = lep + oproj_width < 1.D ind = where(((use_g[*].Left_edge[zi]) LT rip) AND $ ((use_g[*].Right_edge[zi]) GT lep)) use_g =use_g(ind) ung = N_ELEMENTS(use_g) x_dis = d_dis(xi,ind[*]) y_dis = d_dis(yi,ind[*]) c_dis = c_dis(ind[*]) index_range_a[*,0:(ung-1)] = index_range_a(*,ind[*]) print, ' using', ung,' grids from Levels ', Levels ; create projections of indiviual grids gd = PTRARR(ung) print, ' reading and projecting of individual grids' print, ' and creating full projection' plida=PTRARR(ng+1) ; pointer array for use with indeces FOR i=0,ng DO plida[i] = ptr_new([0L]) ci = 0 i = ung - 1L run_i = -1L ; running index of data points ; now start at lowest level and construct list of points WHILE (i GE 0) DO BEGIN ; now what is its dx cg = use_g[i] index_range = index_range_a[*,i] dx = x_dis[i] dy = y_dis[i] ; <- necessary because of unequal round off errors dz = c_dis[i] weighted = 0 ;read_hdf, cg.baryon_file, 1, weight ; density ;read_hdf, cg.baryon_file, 16, data ; temperature data = 0. read_hdf, cg.baryon_file, 1, data ; density ; read_hdf, cg.baryon_file, 14, sdata ; H2 ; data = sdata/temporary(data) ; H2 ; list of grids that have the current grid as parent: dx_level[cg.level] = dx pi = where(use_g.parent EQ cg.num) s = size(data) bs = s[1:3]-1 IF (pi[0] GE 0) THEN BEGIN FOR jj=0,N_ELEMENTS(pi)-1 DO BEGIN ; zeroing all data values where childs are: CASE zi OF 2: BEGIN data[(use_g(pi[jj]).parent_s_index[xi]):(use_g(pi[jj]).parent_e_index[xi]), $ (use_g(pi[jj]).parent_s_index[yi]):(use_g(pi[jj]).parent_e_index[yi]), $ (use_g(pi[jj]).parent_s_index[zi]):(use_g(pi[jj]).parent_e_index[zi])] = 0.D IF weighted THEN $ weight[(use_g(pi[jj]).parent_s_index[xi]):(use_g(pi[jj]).parent_e_index[xi]), $ (use_g(pi[jj]).parent_s_index[yi]):(use_g(pi[jj]).parent_e_index[yi]), $ (use_g(pi[jj]).parent_s_index[zi]):(use_g(pi[jj]).parent_e_index[zi])] = 0.D END 1: BEGIN data[(use_g(pi[jj]).parent_s_index[xi]):(use_g(pi[jj]).parent_e_index[xi]), $ (use_g(pi[jj]).parent_s_index[zi]):(use_g(pi[jj]).parent_e_index[zi]), $ (use_g(pi[jj]).parent_s_index[yi]):(use_g(pi[jj]).parent_e_index[yi])] = 0.D IF weighted THEN $ weight[(use_g(pi[jj]).parent_s_index[xi]):(use_g(pi[jj]).parent_e_index[xi]), $ (use_g(pi[jj]).parent_s_index[zi]):(use_g(pi[jj]).parent_e_index[zi]), $ (use_g(pi[jj]).parent_s_index[yi]):(use_g(pi[jj]).parent_e_index[yi])] = 0.D END 0: BEGIN data[(use_g(pi[jj]).parent_s_index[zi]):(use_g(pi[jj]).parent_e_index[zi]), $ (use_g(pi[jj]).parent_s_index[xi]):(use_g(pi[jj]).parent_e_index[xi]), $ (use_g(pi[jj]).parent_s_index[yi]):(use_g(pi[jj]).parent_e_index[yi])] = 0.D IF weighted THEN $ weight[(use_g(pi[jj]).parent_s_index[zi]):(use_g(pi[jj]).parent_e_index[zi]), $ (use_g(pi[jj]).parent_s_index[xi]):(use_g(pi[jj]).parent_e_index[xi]), $ (use_g(pi[jj]).parent_s_index[yi]):(use_g(pi[jj]).parent_e_index[yi])] = 0.D END ELSE: PRINT, 'warning : constant index is not in [0,1,2]:',zi ENDCASE ENDFOR ENDIF s=size(data) lco = (lep - (cg.Left_edge[zi] + 0.5D * dz))/dz > 0.D rco = (rip - (cg.Left_edge[zi] + 0.5D * dz))/dz < DOUBLE(s[zi+1]-1) totdz = min([cg.Right_edge[zi],rip]) - max([cg.Left_edge[zi],lep]) lc = ROUND(lco) rc = ROUND(rco) IF weighted THEN BEGIN IF (zi EQ 0) THEN weight2d=weight[lc:rc,*,*] IF (zi EQ 1) THEN weight2d=weight[*,lc:rc,*] IF (zi EQ 2) THEN weight2d=weight[*,*,lc:rc] ENDIF IF (zi EQ 0) THEN data2d=data[lc:rc,*,*] IF (zi EQ 1) THEN data2d=data[*,lc:rc,*] IF (zi EQ 2) THEN data2d=data[*,*,lc:rc] ; print, size(data2d) ; integrator = totdz/DOUBLE(rc-lc+1) integrator = totdz/DOUBLE(rc-lc+1) IF weighted THEN BEGIN data2d = TOTAL(data2d*weight2d, zi+1)*integrator weight2d = TOTAL(weight2d, zi+1)*integrator ENDIF ELSE data2d = TOTAL(data2d, zi+1)*integrator s = size(data2d) ci = ci+1 dx_level[cg.level] = dx ; now loop over both other indeces, compute x_stuetz, y_stuetz and interpolate z_stuetz ilev = LINDGEN(ung-1-i+1)+i ; ilev = where(use_g.level GE cg.level) Imparentof = where(use_g.parent EQ cg.num) FOR j=0,index_range[xi]-1 DO $ FOR k=0,index_range[yi]-1 DO BEGIN ; scaled coords in [0.,1] xco = cg.Left_edge[xi] + (0.5D + DOUBLE(j))*dx yco = cg.Left_edge[yi] + (0.5D + DOUBLE(k))*dy cp = DBLARR(3) cp[other_subs] = [xco, yco] cp[zi] = oproj_coord[zi] forget = 0 add_new = 1 IF (NOT forget) THEN BEGIN ; of the grids we have already read are there any on my level or higher ? iogr = where((use_g[ilev].left_edge[xi] lt (cp[xi]+0.5D*dx)) AND $ (use_g[ilev].right_edge[xi] gt (cp[xi]-0.5D*dx)) AND $ (use_g[ilev].left_edge[yi] lt (cp[yi]+0.5D*dy)) AND $ (use_g[ilev].right_edge[yi] gt (cp[yi]-0.5D*dy))) IF (iogr[0] GE 0) THEN BEGIN nbggr = N_ELEMENTS(iogr) FOR kk=0,(nbggr-1) DO BEGIN ; add to existing point: runco = (*plida[use_g[ilev[iogr[kk]]].num]) IF (((runco[0]) ne -1) AND (ilev[iogr[kk]] GT i)) THEN BEGIN ; fully contained ? Then don't do search for points. IF ((use_g[ilev[iogr[kk]]].left_edge[xi] GT (cp[xi]-0.5D*dx)) AND $ (use_g[ilev[iogr[kk]]].right_edge[xi] LT (cp[xi]+0.5D*dx)) AND $ (use_g[ilev[iogr[kk]]].left_edge[yi] GT (cp[yi]-0.5D*dy)) AND $ (use_g[ilev[iogr[kk]]].right_edge[yi] LT (cp[yi]+0.5D*dy))) THEN $ theones=runco $ ELSE $ ; we have to ... theones= where((x_stuetz(runco) GT (cp[xi]-0.5D*dx)) AND $ (x_stuetz(runco) LT (cp[xi]+0.5D*dx)) AND $ (y_stuetz(runco) GT (cp[yi]-0.5D*dy)) AND $ (y_stuetz(runco) LT (cp[yi]+0.5D*dy))) IF (theones[0] GE 0) THEN BEGIN z_stuetz(runco(theones)) = z_stuetz(runco(theones)) + $ data2d[j,k] IF weighted THEN $ w_stuetz(runco(theones)) = w_stuetz(runco(theones)) + $ weight2d[j,k] add_new = 0 ; print, 'grid:',use_g[iogr[kk]].num, cg.num, theones, 'runco:',runco ; print, theones ENDIF ENDIF ENDFOR ENDIF IF (add_new) THEN BEGIN run_i = run_i + 1 x_stuetz[run_i] = cp[xi] y_stuetz[run_i] = cp[yi] z_stuetz[run_i] = data2d[j,k] IF weighted THEN w_stuetz[run_i] = weight2d[j,k] (*plida[cg.num])= [(*plida[cg.num]), run_i] ; print, N_ELEMENTS(*plida[cg.level]) c_level[run_i] = BYTE(cg.level) rg = cg ENDIF ENDIF ENDFOR print, i, cg.num, cg.level, run_i i = i-1 ENDWHILE ; reduce arrays: ind = REVERSE(where(z_stuetz[0:run_i] NE 0.D)) x_stuetz = x_stuetz[ind] y_stuetz = y_stuetz[ind] IF weighted THEN z_stuetz = z_stuetz[ind]/w_stuetz[ind] $ ELSE z_stuetz = z_stuetz[ind] c_level = c_level[ind] run_i = N_ELEMENTS(ind)-1 slice.min_level = min_level slice.max_level = max_level slice.Np_total = run_i+1 slice.N_dfields = 1 slice.oproj_coord = oproj_coord slice.oproj_width = oproj_width slice.Left_edge = Left_edge slice.Right_edge = Right_edge slice.Data_fields= Data_fields slice.x_ptr = PTR_NEW(x_stuetz) slice.y_ptr = PTR_NEW(y_stuetz) slice.z_ptr[0] = PTR_NEW(z_stuetz) slice.Level_ptr= PTR_NEW(c_Level) slice.dx_level = dx_level print, 'OJ::Get_Orthogonal_Projection() : included ', run_i+1,' (',sqrt(FLOAT(run_i+1)),$ '^2) data points from ', ci, ' grids.' RETURN, slice END PRO OJ__define data_max_dim = 3 num_pattr = 8 ; # of Particle Attributes: ; Star formation = 11, None = 8 max_num_of_grids = 150000L max_num_of_levels = 45 info = {grid, $ num: 0L, $ time: 0.D, $ timestep: 0L, $ parent: -1L, $ parent_s_index: INTARR(data_max_dim), $ parent_e_index: INTARR(data_max_dim), $ hier_file: '', $ level: 0, $ dim: INTARR(data_max_dim), $ Start_index: INTARR(data_max_dim), $ End_index: INTARR(data_max_dim), $ Left_edge: DBLARR(data_max_dim), $ Right_edge: DBLARR(data_max_dim), $ num_baryon_fields: 0, $ num_particle: LONG(0), $ baryon_file: '', $ particle_file:'' , $ ; num_particle_attr: num_pattr, $ ngnl: 0L, $ ngtl: 0L $ } grid_info = Replicate(info, max_num_of_grids) n_datafields = 50 struct = {OJ, $ file_base:'', $ Left_edge:DBLARR(data_max_dim), $ Right_edge:DBLARR(data_max_dim), $ Redshift:0.D, $ list_str:strarr(n_datafields), $ min_level:0, $ max_level:0, $ refined_by:0, $ num_of_grids:LONG(0), $ gi:grid_info, $ data:PTRARR(n_datafields), $ data_fields:strarr(n_datafields)$ } END