pro setup_ic,DMONLY=DMONLY,HYDROGRID=HYDROGRID,split=split @setconst nx=0L ny=0L nz=0L dx=0.0 a0=0.0 o0=0.0 l0=0.0 h0=0.0 openr,1,'ic_xc.dat',/f77 readu,1,nx,ny,nz,dx,a0,o0,l0,h0 close,1 pos=fltarr(3,nx,ny,nz) vel=fltarr(3,nx,ny,nz) dummy=fltarr(nx,ny) sub=['x','y','z'] FOR i=0,2 DO BEGIN openr,1,'ic_'+sub(i)+'c.dat',/f77 readu,1,nx,ny,nz,dx,a0,o0,l0,h0 openr,2,'ic_v'+sub(i)+'c.dat',/f77 readu,2,nx,ny,nz,dx,a0,o0,l0,h0 FOR j=0,nz-1 DO BEGIN readu,1,dummy pos(i,*,*,j)=dummy(*,*) readu,2,dummy vel(i,*,*,j)=dummy(*,*) END close,1 close,2 END Omega0=double(o0) OmegaLambda=double(l0) HubbleParam=double(h0/100) BoxSize=double(ROUND(dx*nx*HubbleParam*1000)) ; setuniverse,o=Omega0,l=OmegaLambda,h=HubbleParam ; mtot=HubbleParam*(rocrit(0.0)*Omega0)/solmas/1e10/nx/ny/nz*mparsck^3*(BoxSize/1000/HubbleParam)^3 mtot=2.5703901/nx/ny/nz*(BoxSize/1000/HubbleParam)^3 npart=[0L,nx*ny*nz,0L,0L,0L,0L] massarr=[0.0d0,mtot,0.0d0,0.0d0,0.0d0,0.0d0] time=double(a0) redshift=double(1./a0)-1.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} print,'n1..3:',nx,ny,nz print,'Cos:',dx,a0,o0,l0,h0 ; we need to use Int (u/a dt) which at high z i ; approx u/a 3/2 dt ; @setconst ; dt=z2t(h.redshift) ; disp=vel*1e5*dt/mparsck ; npos=(pos*a0-1.5*disp)/a0 ; plot,pos(0,*,*,0),pos(1,*,*,0),psym=4,xst=1,yst=1,xr=[0,10],yr=[0,10] ; FOR i=0,10 DO BEGIN ; oplot,[0,100],[1,1]*i*dx+dx/2,color=mycolor(1) ; oplot,[1,1]*i*dx+dx/2,[0,100],color=mycolor(1) ; END ; oplot,npos(0,*,*,0),npos(1,*,*,0),psym=4,col=mycolor(4) ob0=0.043 IF NOT(IS_DEF(DMONLY)) THEN BEGIN h.npart(0)=nx*ny*nz mgas=mtot*ob0/o0 mdm=mtot-mgas h.massarr[0]=mgas h.massarr[1]=mdm dxdm=mgas*dx/2/mtot dxgas=mdm*dx/2/mtot IF IS_DEF(HYDROGRID) THEN dxdm=0.0 END ELSE BEGIN mdm=mtot mgas=0.0d0 dxdm=0.0 dxgas=0.0 END ntot=h.npart(0)+h.npart(1) nsing=h.npart(0) pos_all=fltarr(3,ntot) vel_all=fltarr(3,ntot) id=lindgen(ntot)+1L IF NOT(IS_DEF(DMONLY)) THEN dummy=fltarr(nsing) count=0L FOR ix=0,nx-1 DO BEGIN FOR iy=0,ny-1 DO BEGIN FOR iz=0,nz-1 DO BEGIN FOR k=0,2 DO BEGIN IF NOT(IS_DEF(DMONLY)) THEN BEGIN IF IS_DEF(HYDROGRID) THEN BEGIN p=[ix,iy,iz] pp=(float(p(k))+0.5)*dx pos_all(k,count)=pp vel_all(k,count)=0.0 END ELSE BEGIN pos_all(k,count)=pos(k,ix,iy,iz)-dxgas vel_all(k,count)=vel(k,ix,iy,iz) END END pos_all(k,nsing+count)=pos(k,ix,iy,iz)+dxdm vel_all(k,nsing+count)=vel(k,ix,iy,iz) END count=count+1L END END END vel_all=vel_all/sqrt(a0) pos_all=pos_all*h0/100*1000. plot,pos_all(0,nsing:ntot-1),pos_all(1,nsing:ntot-1),psym=3,xr=[0,2000],yr=[0,2000] IF NOT(IS_DEF(DMONLY)) THEN oplot,pos_all(0,0:nsing-1),pos_all(1,0:nsing-1),psym=4,col=mycolor(1) jj=where(pos_all LT 0) IF jj(0) GE 0 THEN pos_all(jj)=pos_all(jj)+100000. jj=where(pos_all GT 100000.) IF jj(0) GE 0 THEN pos_all(jj)=pos_all(jj)-100000. IF NOT(IS_DEF(DMONLY)) THEN name='IC.dat' ELSE name='IC_dm.dat' 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 !!! 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,dummy(ngaswritten:ngaswritten+ngastowrite-1),'U ' add_block,fname,dummy(ngaswritten:ngaswritten+ngastowrite-1),'RHO ' 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 ' IF NOT(IS_DEF(DMONLY)) THEN BEGIN add_block,name,dummy,'U ' add_block,name,dummy,'RHO ' END END stop end