pro setup_cloud

rho_mean_today = 9.2070437e-30            ; in g/cm^3

m_cluster_today = 1d15 * 1.9890000d+33    ; 1e15 Msol in g

r_cloud = (m_cluster_today / rho_mean_today / 4 / !Pi * 3) ^ (1. / 3.)      ; in cm

r_cloud_code = r_cloud /  3.085678e21        ; in code units

N=  21L

fout="box.ic"

Np_max=   N^3

pos= fltarr(3,Np_max)
vel= fltarr(3,Np_max)
mass= fltarr(Np_max)
id= lonarr(Np_max)

ip = 0L
for i=0L, N-1 do begin
  x = (i/(N-1.0) - 0.5) * 2 * r_cloud_code
  for j=0L, N-1 do begin
    y = (j/(N-1.0) - 0.5) * 2 * r_cloud_code
    for k=0L, N-1 do begin
      z = (k/(N-1.0) - 0.5) * 2 * r_cloud_code
      r = sqrt(x^2+y^2+z^2)
      if(r LE r_cloud_code) then begin
         pos(0,ip) = x
         pos(1,ip) = y
         pos(2,ip) = z
         id(ip) = ip +1
         ip = ip + 1
      end
    end
  end
end

vel(*,*) = 0
mass(*) = m_cluster_today / 1.989d43 / ip      ; in code units

; reduce arrays
pos = pos(*,0:ip-1)
vel = vel(*,0:ip-1)
id = id(0:ip-1)
mass = mass(0:ip-1)


npart=lonarr(6)
massarr=dblarr(6)
time=0.0D
redshift=0.0D
flag_sfr=0L
flag_feedback=0L
partTotal=lonarr(6)
flag_cooling=0L
num_files=1L
BoxSize=1.0D
Omega0=0.0D
OmegaLambda=0.0D
HubbleParam=0.7D
flag_stellarage=0L;          /*!< flags whether the file contains formation times of star particles */
flag_metals=0L;              /*!< flags whether the file contains metallicity values for gas and star particles */
npartTotalHighWord=lonarr(6);
Labels=bytarr(2,15)
bytesleft=256-6*4 - 6*8 - 8 - 8 - 2*4 - 6*4 - 2*4 - 4*8 - 2*4 - 6*4 - 2*15
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,$
             flag_stellarage:flag_stellarage,$
             flag_metals:flag_metals,$
             npartTotalHighWord:npartTotalHighWord,$
	     Labels:Labels,$
             la:la}

h.npart(1)=ip
h.partTotal(1)=ip

write_head,fout,h
add_block,fout,pos,'POS '
add_block,fout,vel,'VEL '
add_block,fout,id,'ID  '
add_block,fout,mass,'MASS'

stop

end