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 PRO set_parents, gi parents = where(gi.ngnl NE 0) if (parents[0] eq -1) then return np = N_ELEMENTS(parents)-1 FOR i = 0L, np DO BEGIN cp = parents[i] cg = long(gi[cp].ngnl-1) stopit = 0 WHILE (NOT stopit) DO BEGIN gi[cg].parent = cp+1 ; while we're at it, also set the start and end indices tp = gi[cp] pd = tp.dim pdx = (tp.right_edge-tp.left_edge)/DOUBLE(pd) for k=0,2 do begin gi[cg].parent_s_index[k] = ROUND((gi[cg].left_edge[k] - tp.left_edge[k])/pdx[k]) < (pd[k]-1) gi[cg].parent_e_index[k] = ROUND((gi[cg].right_edge[k]- tp.left_edge[k])/pdx[k])-1 < (pd[k]-1) endfor IF (gi[cg].ngtl NE 0) THEN cg = gi[cg].ngtl-1 ELSE stopit=1 IF (vec_GT_vec(gi[cp].Left_edge, gi[cg].Left_edge) OR $ vec_LT_vec(gi[cp].Right_edge, gi[cg].Right_edge)) THEN BEGIN print, 'ok. thats wrong', gi[cp].Right_edge, gi[cg].Right_edge ENDIF ENDWHILE ENDFOR ind = where(gi.parent eq gi.num, count) if count gt 0 then begin print, 'found ', count,' grids that think they are their own parents' gi[ind].parent = 0 endif RETURN END function READ_GRID_INFO, file_name, REDSHIFT=red, DATA_DIR=data_dir, NOORDER=noorder, Gamma=gamma ; extract interesting information from hierarchy file ; 3d data: data_dim = 3 info = {grid, $ num: 0L, $ time: 0.D, $ timestep: 0L, $ parent: -1L, $ parent_s_index: INTARR(data_dim), $ parent_e_index: INTARR(data_dim), $ hier_file: '', $ level: 0, $ dim: INTARR(data_dim), $ Start_index: INTARR(data_dim), $ End_index: INTARR(data_dim), $ Left_edge: DBLARR(data_dim), $ Right_edge: DBLARR(data_dim), $ num_baryon_fields: 0 , $ num_particle: LONG(0) , $ baryon_file: '', $ particle_file:'', $ ngnl: 0L, $ ngtl: 0L $ } gamma = 5./3. ; unless told otherwise help_string = '' print, 'reading grid info from ', file_name i = 0L h_particle_file = "" h_baryon_file = "" h_paticle_file = "" temp_grid_info = info if strmatch(file_name, '*.hierarchy') then BEGIN ; hdf4 or packed hdf5 format fstring = read_ascii_file_return_list_of_strings(file_name[0], /verbose, nlines=nlines) dim_lines = fstring[where((strmatch(fstring, 'GridDim*') eq 1), ngrids)] temp_grid_info = Replicate(info, ngrids) wherestart = where(strmatch(fstring, 'GridStart*') eq 1) whereleft = where(strmatch(fstring, 'GridLeft*') eq 1) startindex = fstring[wherestart] endindex = fstring[wherestart+1] leftedge = fstring[whereleft] rightedge = fstring[whereleft+1] nfields = fstring[where(strmatch(fstring, 'NumberOfBaryon*') eq 1)] equalpos = strpos(nfields[0],"=") nfields = fix(strmid(nfields, equalpos+1)) time = fstring[where(strmatch(fstring, 'Time*') eq 1)] equalpos = strpos(time[0],"=") time = FLOAT(strmid(nfields, equalpos+1)) ngnl = LONARR(ngrids) hngnl = fstring[where(strmatch(fstring, '*NextGridNextLeve*') eq 1)] ind = LONARR(N_elements(hngnl)) equalpos = strpos(hngnl,"=") for ii=0L, ngrids-1 do BEGIN mi = LONG(STRMID(hngnl[ii], 14, strpos(hngnl[ii], ']')-14))-1 mngnl = LONG(strmid(hngnl[ii], equalpos[ii]+1)) ngnl[mi] = mngnl ; print, mi, mngnl, ' ', hngnl[ii] endfor ngtl = LONARR(ngrids) hngtl = fstring[where(strmatch(fstring, '*NextGridThisLeve*') eq 1)] ind = LONARR(ngrids) equalpos = strpos(hngtl,"=") for ii=0L, ngrids-1 do $ ngtl(LONG(STRMID(hngtl[ii], 14, strpos(hngtl[ii], ']')-14))-1) = long(strmid(hngtl[ii], equalpos[ii]+1)) baryonfile = where(strmatch(fstring, 'BaryonFile*') eq 1, count) if (count gt 0) then begin baryonfile = fstring[baryonfile] equalpos = strpos(baryonfile[0],"=") baryonfile = strtrim(strmid(baryonfile, equalpos+1),2) ; this last slash in practice is always at the same position ... (we hope) lastslash = strpos(baryonfile[0],"/", /REVERSE_SEARCH) baryonfile = strmid(baryonfile, lastslash+1) endif else begin baryonfile = strarr(ngrids) endelse npart = fstring[where(strmatch(fstring, 'NumberOfPart*') eq 1)] npartvalue = get_int_value_arr(npart) nonzeros = where(npartvalue ne 0, nonzero_count) if (nonzero_count eq ngrids) then begin partfile = fstring[where(strmatch(fstring, 'ParticleFileName*') eq 1)] equalpos = strpos(partfile[0],"=") partfile = strtrim(strmid(partfile, equalpos+1),2) lastslash = strpos(partfile[0],"/", /REVERSE_SEARCH) partfile = strmid(partfile, lastslash+1) endif else begin partfile = strarr(ngrids) if (nonzero_count ne 0) then begin linenums = where(strmatch(fstring, 'ParticleFileNam*') eq 1) if N_elements(linenums) eq N_elements(nonzeros) then $ partfile[nonzeros] = fstring[linenums] $ else partfile = fstring[linenums] equalpos = strpos(partfile[nonzeros[0]],"=") partfile = strtrim(strmid(partfile, equalpos[0]+1),2) lastslash = strpos(partfile[0],"/", /REVERSE_SEARCH) partfile = strmid(partfile, lastslash+1) endif taskstr = fstring[where(strmatch(fstring, 'Task*') eq 1)] tasks = get_int_value_arr(taskstr) if (N_elements(tasks) gt 0) then begin dot_is_at = StrPos(file_name, '.',/reverse_search) slash_is_at = StrPos(file_name, '/',/reverse_search) l_file_name = STRMID(file_name, slash_is_at+1, $ dot_is_at-slash_is_at-1) zeros = where(npartvalue eq 0, zero_count) if (zero_count gt 0) then $ partfile[zeros] = l_file_name + $ STRING(tasks[zeros], FORMAT='(".cpu",I4.4)') endif endelse dot_is_at = StrPos(file_name, '.',/reverse_search) l_file_name = STRMID(file_name, 0, dot_is_at) slash_is_at = StrPos(l_file_name, '/',/reverse_search) data_dir = STRMID(file_name, 0, slash_is_at+1) l_file_name = STRMID(l_file_name, slash_is_at+1,300) temp_grid_info.num = lindgen(ngrids)+1 temp_grid_info.hier_file = file_name[0] temp_grid_info.start_index = get_intarr_arr(startindex) temp_grid_info.end_index = get_intarr_arr(endindex) temp_grid_info.dim = temp_grid_info.end_index-temp_grid_info.start_index+1 temp_grid_info.left_edge = get_floatarr_arr(leftedge) temp_grid_info.right_edge = get_floatarr_arr(rightedge) temp_grid_info.num_baryon_fields = nfields temp_grid_info.num_particle = npartvalue ;; Hack to read nBody runs if (nfields[0] gt 0) then $ temp_grid_info.baryon_file = data_dir+baryonfile $ else $ temp_grid_info.baryon_file = data_dir+partfile temp_grid_info.particle_file = data_dir+partfile temp_grid_info.ngnl = ngnl temp_grid_info.ngtl = ngtl temp_grid_info.time = time if NOT KEYWORD_SET(noorder) then begin temp_grid_info = set_level_info(temp_grid_info) ; temp_grid_info = order_grids(temp_grid_info) ENDIF ENDIF $ ; if hdf4 or hdf5 packed format else BEGIN ; movie format temp_grid_info = $ make_hierarchy_from_movie_format_file(file_name, REDSHIFT=red) ; red = -temp_grid_info[0].time ; print, 'REDSHIFT:', red, ' from file:', file_name end if (N_elements(ngnl) gt 1) then set_parents, temp_grid_info ; return grid parameters if N_elements(red) eq 0 then red = -1 RETURN, temp_grid_info END ;.compile TOOLS/read_grid_info