pro make_periodic, x, dle, ds ind = where(x lt dle) if ind[0] ge 0 then x[ind] = x[ind]+ds ind = where(x gt dle+ds) if ind[0] ge 0 then x[ind] = x[ind]-ds end pro read_gadget_file, filename @common_blocks.inc vel.read = 0 mhd = 0 npart=lonarr(6) massarr=dblarr(6) gadgettime=0.0D redshift=0.0D flag_sfr=0L flag_feedback=0L npartTotal=lonarr(6) flag_use_entropy_instead_u = 0L bytesleft=256-6*4 - 6*8 - 8 - 8 - 2*4-6*4 - 4 la=intarr(bytesleft/2) openr,lun,filename,/f77_unformatted, /get_lun readu,lun,npart,massarr,gadgettime,redshift,flag_sfr,flag_feedback,npartTotal,flag_use_entropy_instead_u,la print, 'flag_use_entropy_instead_u:', flag_use_entropy_instead_u N=total(npart) pos=fltarr(3,N) gadgetvel=fltarr(3,N) id=lonarr(N) hubble = 0.5 lengthunit = 3.086d24;/0.5 velunit = 1.e5 massunit = 1.989d43;/0.5 massunit_in_ms = massunit/1.989d33 timeunit = lengthunit/velunit densityunit = massunit/lengthunit^3 energyunit = (lengthunit/timeunit)^2 ind=where((npart gt 0) and (massarr eq 0)) if ind(0) ne -1 then begin Nwithmass= total(npart(ind)) mass=fltarr(Nwithmass) endif else begin Nwithmass= 0 mass = fltarr(total(npart)) if npart[0] gt 0 then mass[0:npart[0]-1] = massarr[0] if npart[1] gt 0 then mass[npart[0]:npart[0]+npart[1]-1] = massarr[1] endelse readu,lun,pos readu,lun,gadgetvel readu,lun,id if Nwithmass gt 0 then begin readu,lun,mass endif NGas=npart[0] names = [''] if Ngas gt 0 then begin u=fltarr(Ngas) readu,lun,u if N_elements(gamma) eq 0 then gamma = 5./3. ; if (flag_use_entropy_instead_u ne 0) then u=(u*rho^gamma)/((gamma-1.)*rho) if ((mhd eq 1) and (not EOF(lun))) then begin Bfield=fltarr(3,Ngas) readu,lun,Bfield endif rho=fltarr(Ngas) readu,lun,rho hsml=fltarr(Ngas) if not EOF(lun) then begin readu,lun,hsml names = ['hsml'] endif endif close,lun free_lun, lun print, ngas, ' gas particles read' particle_type = intarr(total(npart)) particle_type[*] = 1 particle_type[0:max([Ngas-1,0])] = 2 dle = [min(pos[0,*], max=dre0),min(pos[1,*], max=dre1),min(pos[2,*], max=dre2)] dre = [dre0,dre1,dre2] DomainLeftEdge = dle DomainRightEdge = dre print, 'Domain:', dle,dre ;find maximum density particle if Ngas gt 0 then begin densind = where(*parti.names eq 'Density') if densind[0] ne -1 then begin dummy = max(rho, ind) print, dummy center = [pos[0,ind[0]],pos[1,ind[0]],pos[2,ind[0]]] endif endif center = 0.5*[1.,1.,1.] center = [1.,1.,1.] ; center = 0.5*[1.,1.,0.] ; center = (dle+dre)/2 c = center x = pos[0,*]-c[0] y = pos[1,*]-c[1] ; z = 0 z = pos[2,*]-c[2] ; ; ds = dre-dle ; make_periodic, x, dle, ds ; make_periodic, y, dle, ds ; make_periodic, z, dle, ds names = [names[*],'Radius',' Radial Velocity'] radius = reform(sqrt(x^2+y^2+z^2)) ; *lengthunit ;breaki x = 0 y=0 z=0 ; get fieldnames print, 'Linear Momentum' print, total(gadgetvel[0,*]), total(gadgetvel[1,*]), total(gadgetvel[2,*]) ;print, 'Angular Momentum:' if Ngas gt 0 then print, 'fraction of volume:', total(1./rho)/Ngas posvel = ['particle_position_x','particle_position_y','particle_position_z','particle mass', 'particle id', $ 'particle_velocity_x', 'particle_velocity_y', 'particle_velocity_z'] if names[0] eq '' then names = posvel else names = [posvel[*], names[*]] if (Ngas gt 0) then begin names = [names[*], 'Density', 'u', 'p', 'entropy', 'Temperature','particle_type'] if (mhd eq 1) then names = [names[*], 'Bx', 'By', 'Bz'] print, 'assuming gamma=', gamma,' to derive temperature' print, 'the units have not been debugged yet for this ..' print, 'if that is not your gamma change this in read_gadget_file.pro' endif Nfields = N_elements(names) print, Nfields, names list_str = names count=0L parti.d[count++] = ptr_new(pos[0,*],/no_copy) parti.d[count++] = ptr_new(pos[1,*],/no_copy) parti.d[count++] = ptr_new(pos[2,*],/no_copy) parti.d[count++] = ptr_new(mass[*],/no_copy) parti.d[count++] = ptr_new(id[*],/no_copy) parti.d[count++] = ptr_new(gadgetvel[0,*],/no_copy) parti.d[count++] = ptr_new(gadgetvel[1,*],/no_copy) parti.d[count++] = ptr_new(gadgetvel[2,*],/no_copy) if N_elements(hsml) gt 0 then parti.d[count++] = ptr_new(hsml[*],/no_copy) if N_elements(radius) gt 0 then begin parti.d[count++] = ptr_new(radius[*],/no_copy) parti.d[count++] = ptr_new(((pos[0,*]-c[0])*gadgetvel[0,*]+ $ (pos[1,*]-c[1])*gadgetvel[1,*]+ $ (pos[2,*]-c[2])*gadgetvel[2,*])/radius,/no_copy) endif ;if N_elements(rho) gt 0 then parti.d[count++] =ptr_new(rho[*]/2.7766217 ,/no_copy) ; units of mean density ... if N_elements(rho) gt 0 then parti.d[count++] = ptr_new(rho[*] ,/no_copy) ; if N_elements(u) gt 0 then begin parti.d[count++]= ptr_new(u[*],/no_copy) parti.d[count++]= ptr_new(u[*]*(gamma-1.)*rho,/no_copy) ; pressure mu_gadget = 0.588235 ; parti.d[count++]= ; ptr_new(alog(1.67e-24*mu_gadget*energyunit*u*(gamma-1)/1.38e-16/(rho*massunit_in_ms*hubble*hubble/gadgettime^3)^(gamma-1)),/no_copy) ; entropy parti.d[count++]= ptr_new(u[*]*(gamma-1.)/rho^(gamma-1),/no_copy) ; entropy parti.d[count++]= ptr_new(1.67e-24*mu_gadget*energyunit*u*(gamma-1)/1.38e-16,/no_copy) ;temperature print, 'multiplying u with:', 1.67e-24*mu_gadget*energyunit/1.38e-16*(gamma-1), ' to get temperature' endif if N_elements(particle_type) gt 0 then parti.d[count++]= ptr_new(particle_type,/no_copy) if mhd eq 1 then begin parti.d[count++]= ptr_new(Bfield[0,*],/no_copy) parti.d[count++]= ptr_new(Bfield[1,*],/no_copy) parti.d[count++]= ptr_new(Bfield[2,*],/no_copy) endif ;breaki goto, skip_voronoi ;compute voronoi densities and volumes vv = fltarr(total(npart)) vd = fltarr(total(npart)) print, "computing voronoi volume and density" if Ngas gt 0 then begin print, "gas particles" vv[0:Ngas-1] = voronoi_volume(pos[0,0:Ngas-1], pos[1,0:Ngas-1],pos[2,0:Ngas-1], fractional=1) vd[0:Ngas-1] = 1.D/vv[0:Ngas-1]/Ngas ; normalized so mean density = 1 (assumed equal mass) endif if total(npart) gt Ngas then begin print, "other particles" vv[Ngas:*] = voronoi_volume(pos[0,Ngas:*], pos[1,Ngas:*],pos[2,Ngas:*], fractional=1) vd[Ngas:*] = 1.D/vv[Ngas:*]/(total(npart)-Ngas) ; normalized so mean density = 1 (assumed equal mass) endif parti.d[count++]= ptr_new(vv[*],/no_copy) parti.d[count++]= ptr_new(vd[*],/no_copy) print, "done with voronoi computation" names = [names[*], 'voronoi_volume', 'voronoi_density'] skip_voronoi: parti.names = ptr_new(names) for i=0,Nfields-1 do parti.read[i] = 1 print, 'total number of particles:', n 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 $ } print, 'Time:', gadgettime info.time = gadgettime info.num_particle = total(Npart) info.hier_file = filename info.num_baryon_fields = Nfields info.baryon_file = filename info.particle_file = filename info.left_edge = dle info.right_edge = dre info.dim = [100,100,100] all_grid_info = info grid_info = info times = [gadgettime] redshifts = [redshift] end