PRO ryu2snap,dm_only=dm_only,split=split @setconst ; ng=64L ; zri=34.63170 ; ng=128L ; zri = 44.77303 ; ng=256L ; zri = 55.92095 ng=512L zri=67.99017 split=4 rcube = 100 fh=0.76 fhe=0.24 Omega0=0.27D OmegaLambda=0.73D HubbleParam=0.7D BoxSize=rcube*1000.0D urho = 1.879e-29*HubbleParam^2*Omega0 uvel = 1.224e+7*rcube*(1.+zri)^0.5*Omega0^0.5 utmp = 1.814e+6*rcube^2*(1.+zri)*Omega0/(2.*fh+0.75*fhe) n=ng/2 dummy=0L openr,1,'dm_pos_512-256',/f77 readu,1,dummy readu,1,dummy pos=fltarr(3,n,n,n) dum=fltarr(n,n) FOR iz=0,n-1 DO BEGIN FOR k=0,2 DO BEGIN readu,1,dum pos(k,*,*,iz)=dum(*,*) END END close,1 vel=fltarr(3,n,n,n) vel(*,*,*,*)=0.0 pos=pos/(2*n)*100000. dum=fltarr(ng,ng) openr,1,'fluid08a',/f77 readu,1,dummy readu,1,dummy rho=fltarr(ng,ng,ng) FOR iz=0,ng-1 DO BEGIN readu,1,dum rho(*,*,iz)=dum(*,*) END close,1 openr,1,'fluid08e',/f77 readu,1,dummy readu,1,dummy temp=fltarr(ng,ng,ng) FOR iz=0,ng-1 DO BEGIN readu,1,dum temp(*,*,iz)=dum(*,*) END close,1 openr,1,'fluid08b',/f77 readu,1,dummy readu,1,dummy openr,2,'fluid08c',/f77 readu,2,dummy readu,2,dummy openr,3,'fluid08d',/f77 readu,3,dummy readu,3,dummy velg=fltarr(3,ng,ng,ng) FOR iz=0,ng-1 DO BEGIN FOR k=0,2 DO BEGIN readu,1+k,dum velg(k,*,*,iz)=dum(*,*) END END close,1 close,2 close,3 temp=temp/rho*utmp FOR k=0,2 DO velg(k,*,*,*)=velg(k,*,*,*)/rho(*,*,*)*uvel rho=rho*urho ; setuniverse,o=Omega0,l=OmegaLambda,h=HubbleParam ; mdm=(1-0.159)*HubbleParam*(rocrit(0.0)*Omega0)/solmas/1e10/n/n/n*mparsck^3*(BoxSize/1000/HubbleParam)^3 mdm=7493849.7/n/n/n IF NOT(IS_DEF(dmonly)) THEN BEGIN mdm=(1.-0.159)*mdm npart=[ng*ng*ng,n*n*n,0L,0L,0L,0L] END ELSE BEGIN npart=[0L,n*n*n,0L,0L,0L,0L] END massarr=[0.0d0,mdm,0.0d0,0.0d0,0.0d0,0.0d0] time=1.0D redshift=0.0D flag_sfr=0L flag_feedback=0L partTotal=npart flag_cooling=0L num_files=1L bytesleft=256-6*4 - 6*8 - 8 - 8 - 2*4 - 6*4 - 2*4 - 4*8 la=bytarr(bytesleft) h = { head , npart:npart,$ massarr:massarr,$ time:time,$ redshift:redshift,$ flag_sfr:flag_sfr,$ flag_feedback:flag_feedback,$ partTotal:partTotal,$ flag_cooling:flag_cooling,$ num_files:num_files,$ BoxSize:BoxSize,$ Omega0:Omega0,$ OmegaLambda:OmegaLambda,$ HubbleParam:HubbleParam,$ la:la} ngas=ng*ng*ng ntot=n*n*n+ngas pos_all=fltarr(3,ntot) vel_all=fltarr(3,ntot) id=lindgen(ntot)+1L dens=fltarr(ngas) mass=fltarr(ngas) u=fltarr(ngas) hsml=fltarr(ngas) count=0L g1 = 5.0 / 3.0 - 1. xNe=1.0 yhelium = (1. - fh) / (4. * fh) mu = (1. + 4. * yhelium) / (1. + 3. * yhelium + xNe) dVphys=(rcube*mparsck/HubbleParam/ng)^3 dV=(BoxSize/ng)^3 FOR ix=0,ng-1 DO BEGIN FOR iy=0,ng-1 DO BEGIN FOR iz=0,ng-1 DO BEGIN FOR k=0,2 DO BEGIN vel_all(k,count)=velg(k,ix,iy,iz)/v_unit END pos_all(0,count)=(ix+0.5)/float(ng)*BoxSize pos_all(1,count)=(iy+0.5)/float(ng)*BoxSize pos_all(2,count)=(iz+0.5)/float(ng)*BoxSize mass(count)=rho(ix,iy,iz)*dVphys/1e10/solmas*HubbleParam dens(count)=mass(count)/dV hsml(count)=1.8*BoxSize/ng u(count)=temp(ix,iy,iz)*kcgs*m_unit/mu/e_unit/mp/g1 count=count+1L END END END FOR ix=0,n-1 DO BEGIN FOR iy=0,n-1 DO BEGIN FOR iz=0,n-1 DO BEGIN FOR k=0,2 DO BEGIN pos_all(k,count)=pos(k,ix,iy,iz) vel_all(k,count)=vel(k,ix,iy,iz) END count=count+1L END END END mtot=TOTAL(mass) print,'mean_rho from mass (gas):',mtot*1e10*solmas/0.7/(100*mparsck/0.7)^3,rocrit(0.0)*0.27*0.159 name='snap_010' IF IS_DEF(split) THEN BEGIN ngaswritten=0L ndmwritten=0L nperfile=(h.npart(0)+h.npart(1))/split h.num_files=split FOR i=0,split-1 DO BEGIN fname=name+'.'+string(i,FORMAT='(i1)') IF i LT split-1 THEN ntowrite=nperfile ELSE ntowrite=h.npart(0)+h.npart(1)-nperfile*(split-1) IF ngaswritten LT h.npart(0) THEN BEGIN IF ngaswritten+ntowrite LT h.npart(0) THEN ngastowrite=ntowrite ELSE ngastowrite=h.npart(0)-ngaswritten END ELSE BEGIN ngastowrite=0L END IF ngastowrite LT ntowrite THEN ndmtowrite=ntowrite-ngastowrite ELSE ndmtowrite=0L IF ndmwritten+ndmtowrite GT h.npart(1) THEN BEGIN print,'PANICK !!! Should not happen (dm) !!! stop END IF ngaswritten+ngastowrite GT h.npart(0) THEN BEGIN print,'PANICK !!! Should not happen (gas) !!! stop END ht=h ht.npart(0)=ngastowrite ht.npart(1)=ndmtowrite write_head,fname,ht add_block,fname,pos_all(*,ngaswritten+ndmwritten:ngaswritten+ndmwritten+ngastowrite+ndmtowrite-1),'POS ' add_block,fname,vel_all(*,ngaswritten+ndmwritten:ngaswritten+ndmwritten+ngastowrite+ndmtowrite-1),'VEL ' add_block,fname,id(ngaswritten+ndmwritten:ngaswritten+ndmwritten+ngastowrite+ndmtowrite-1),'ID ' IF ngastowrite GT 0 THEN BEGIN add_block,fname,mass(ngaswritten:ngaswritten+ngastowrite-1),'MASS' add_block,fname,u(ngaswritten:ngaswritten+ngastowrite-1),'U ' add_block,fname,dens(ngaswritten:ngaswritten+ngastowrite-1),'RHO ' add_block,fname,hsml(ngaswritten:ngaswritten+ngastowrite-1),'HSML' END ngaswritten=ngaswritten+ngastowrite ndmwritten=ndmwritten+ndmtowrite END print,'DM: ',ndmwritten,' == ',h.npart(1),' ??' print,'GAS: ',ngaswritten,' == ',h.npart(0),' ??' END ELSE BEGIN write_head,name,h add_block,name,pos_all,'POS ' add_block,name,vel_all,'VEL ' add_block,name,id,'ID ' add_block,name,mass,'MASS' add_block,name,dens,'RHO ' add_block,name,u,'U ' add_block,name,hsml,'HSML' END stop end pro test @setconst readnew,'snap_010',h,'HEAD' readnew,'snap_010',rho,'RHO' IF h.massarr(0) EQ 0 THEN BEGIN IF h.num_files GT 1 THEN BEGIN m=fltarr(h.npart(0)) offs=0L FOR i=0,h.num_files-1 DO BEGIN readnew,'snap_010.'+string(i,FORMAT='(i1)'),ht,'HEAD' readnew,'snap_010.'+string(i,FORMAT='(i1)'),mt,'MASS' m(offs:offs+ht.npart(0)-1)=mt(0:ht.npart(0)-1) offs=offs+ht.npart(0) END END ELSE BEGIN readnew,'snap_010',m,'MASS' m=m(0:h.npart(0)-1) END mtot=TOTAL(m) END ELSE BEGIN mtot=h.npart(0)*h.massarr(0) m=h.massarr(0) END dens=density(rho) dVphys=m*1e10*solmas/0.7/dens print,'mean_rho from mass (dm): ',h.npart(1)*h.massarr(1)*1e10*solmas/0.7/(100*mparsck/0.7)^3,rocrit(0.0)*0.27*(1.-0.159) print,'mean_rho from mass (gas):',mtot*1e10*solmas/0.7/(100*mparsck/0.7)^3,rocrit(0.0)*0.27*0.159 print,'mean_rho from rho:',TOTAL(dens)/h.npart(0),TOTAL(dens*dVphys)/TOTAL(dVphys) print,'vol from rho: ',total(m/rho),1e5^3 print,'vol from dens:',total(dVphys),(100*mparsck/0.7)^3 end pro test_ryu @setconst rcube = 100 fh=0.76 fhe=0.24 Omega0=0.27D OmegaLambda=0.73D HubbleParam=0.7D BoxSize=rcube*1000.0D nam_tab=['64-32','128-64','256-128','512-256'] ; nam_tab=['64-32','128-64','256-128'] ng_tab=[64L,128L,256L,512L] zri_tab=[34.63170,44.77303,55.92095,67.99017] sst=SIZE(nam_tab) FOR sim=0,sst(1)-1 DO BEGIN ng=ng_tab(sim) zri=zri_tab(sim) urho = 1.879e-29*HubbleParam^2*Omega0 dummy=0L dum=fltarr(ng,ng) openr,1,nam_tab(sim)+'/fluid08a',/f77 readu,1,dummy readu,1,dummy rho=fltarr(ng,ng,ng) FOR iz=0,ng-1 DO BEGIN readu,1,dum rho(*,*,iz)=dum(*,*)*urho END close,1 rhomean=TOTAL(rho)/ng^3 print,nam_tab(sim),rhomean,rocrit(0.0)*0.27*0.159 END end