Page MenuHomec4science
Files F101348214

addgas
No OneTemporary
Actions

File Metadata

Created
Sat, Feb 8, 07:24

addgas
View Options

#!/usr/bin/env python
from numpy import *
import os,sys,string,copy,types
from optparse import OptionParser
#################################
def readblock(f,data_type,shape=None,skip=None,byteorder=sys.byteorder):
#################################
'''
data_type = Int,Float32,Float
or
data_type = array
shape = tuple
'''
import types
try:
nb1 = fromstring(f.read(4),int32)
if sys.byteorder != byteorder:
nb1.byteswap()
nb1 = nb1[0]
except IndexError:
raise "ReadBlockError"
if type(data_type) == types.TupleType:
data = []
for tpe in data_type:
if type(tpe) == types.IntType:
val = f.read(tpe)
elif tpe == int32:
val = fromstring(f.read(4),tpe)
val = val
elif tpe == float64:
val = fromstring(f.read(8),tpe)
val = val[0]
else: # ancienne methode
val = fromstring(f.read(tpe.bytes),tpe)
if sys.byteorder != byteorder:
val.byteswap()
val = val[0]
data.append(val)
else:
data = fromstring(f.read(nb1),data_type)
if sys.byteorder != byteorder:
data.byteswap()
nb2 = fromstring(f.read(4),int32)
if sys.byteorder != byteorder:
nb2.byteswap()
nb2 = nb2[0]
if nb1 != nb2:
raise "ReadBlockError"
# reshape if needed
if shape != None:
data.shape=shape
return data
#################################
def writeblock(f,data,byteorder=sys.byteorder,save_memory=0):
#################################
'''
data = array
or
data = ((x,Float32),(y,Int),(z,Float32),(label,40))
shape = tuple
'''
if type(data) == types.TupleType:
# first, compute nbytes
nbytes = 0
for dat in data:
if type(dat[0])==types.StringType:
nbytes = nbytes + dat[1]
else:
#nbytes = nbytes + array([dat[0]],dat[1]).type().bytes* array([dat[0]],dat[1]).size()
pass
nbytes = 256
nbytes = array([nbytes],int32)
# write block
if sys.byteorder != byteorder:
nbytes.byteswap()
f.write(nbytes.tostring())
for dat in data:
if type(dat[0])==types.StringType:
f.write(string.ljust(dat[0],dat[1])[:dat[1]])
elif dat[1] == int32:
ar = array([dat[0]],dat[1])
f.write(ar.tostring())
elif dat[1] == float64:
ar = array([dat[0]],dat[1])
f.write(ar.tostring())
else:
ar = array([dat[0]],dat[1])
if sys.byteorder != byteorder:
ar.byteswap()
f.write(ar.tostring())
f.write(nbytes.tostring())
else:
# write block
nbytes = array([data.nbytes],int32)
if sys.byteorder != byteorder:
nbytes.byteswap()
data.byteswap()
f.write(nbytes.tostring())
if save_memory:
n = len(data)
f.write(data[:n/2].tostring())
f.write(data[n/2:].tostring())
f.write(nbytes.tostring())
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser.add_option("--fb",
action="store",
dest="fb",
type="float",
default = 0.167587,
help="baryonic fraction",
metavar=" FLOAT")
parser.add_option("-n",
action="store",
dest="n",
type="int",
default = 512,
help="grid size",
metavar=" INT")
parser.add_option("-L",
action="store",
dest="L",
type="float",
default = 2000,
help="box size",
metavar=" FLOAT")
(options, args) = parser.parse_args()
files = args
return files,options
######################################################################################3
#
# M A I N
#
######################################################################################
files,opt = parse_options()
file1 = files[0]
file2 = files[1]
n = opt.n
L = opt.L
f = open(file1)
isize = os.path.getsize(file1)
# read header
tpl = (24,48,float64,float64,int32,int32,24,int32,int32,float64,float64,float64,float64,int32,int32,24,int32,float64,52)
npart,massarr,atime,redshift,flag_sfr,flag_feedback,nall,flag_cooling,num_files,boxsize,omega0,omegalambda,hubbleparam,flag_age,flag_metals,nallhw,flag_entr_ic,critical_energy_spec,empty = readblock(f,tpl,byteorder='little')
npart = fromstring(npart,int32)
massarr = fromstring(massarr,float64)
nall = fromstring(nall,int32)
nallhw = fromstring(nallhw,float64)
nbody = sum(npart)
# read pos
pos = readblock(f,float32,shape=(nbody,3))
vel = readblock(f,float32,shape=(nbody,3))
num = readblock(f,int32)
fsize = f.tell()
if fsize != isize:
raise "ReadError","file %s not red completely"%(file1)
f.close()
######################################################################################3
# now, we can work
######################################################################################
f = opt.fb # baryon fraction
mg = massarr[1]*f
md = massarr[1]*(1-f)
##################
# init
##################
iz,iy,ix = indices((n,n,n))
ix = ravel(ix)
iy = ravel(iy)
iz = ravel(iz)
gas_pos = zeros((n**3,3),float32)
gas_vel = zeros((n**3,3),float32)
# 0,0,0
print "0,0,0"
idx = fmod(ix+0,n)
idy = fmod(iy+0,n)
idz = fmod(iz+0,n)
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 1,0,0
print "1,0,0"
idx = fmod(ix+1,n)
idy = fmod(iy+0,n)
idz = fmod(iz+0,n)
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 1,1,0
print "1,1,0"
idx = fmod(ix+1,n)
idy = fmod(iy+1,n)
idz = fmod(iz+0,n)
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 0,1,0
print "0,1,0"
idx = fmod(ix+0,n)
idy = fmod(iy+1,n)
idz = fmod(iz+0,n)
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 0,0,1
print "0,0,1"
idx = fmod(ix+0,n)
idy = fmod(iy+0,n)
idz = fmod(iz+1,n)
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 1,0,1
print "1,0,1"
idx = fmod(ix+1,n)
idy = fmod(iy+0,n)
idz = fmod(iz+1,n)
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 1,1,1
print "1,1,1"
idx = fmod(ix+1,n)
idy = fmod(iy+1,n)
idz = fmod(iz+1,n)
j = idz*n*n + idy*n + idx
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# 0,1,1
print "0,1,1"
idx = fmod(ix+0,n)
idy = fmod(iy+1,n)
idz = fmod(iz+1,n)
j = idz*n*n + idy*n + idx
lx = (idx - fmod(ix+0,n)<0)*L
ly = (idy - fmod(iy+0,n)<0)*L
lz = (idz - fmod(iz+0,n)<0)*L
dpos = transpose(array([lx,ly,lz],float32))
j = idz*n*n + idy*n + idx
gas_pos = gas_pos + take(pos,j,axis=0)+dpos
gas_vel = gas_vel + take(vel,j,axis=0)
# do the mean
gas_pos = gas_pos/8.
gas_vel = gas_vel/8.
print "concatenate"
pos = concatenate((gas_pos,pos))
vel = concatenate((gas_vel,vel))
print "done"
u = zeros(nbody).astype(float32)
num = arange(nbody+nbody).astype(int32)
npart = array([nbody,0,nbody,0,0,0],int32)
nall = array([nbody,0,nbody,0,0,0],int32)
massarr = array([mg,0,md,0,0,0],float64)
#massarr[0] = 0.035574646727636612
#massarr[1] = 0.16770904885885832
nbody = nbody + nbody
# write final model
npart = npart.tostring()
massarr = massarr.tostring()
nall = nall.tostring()
nallhw = nallhw.tostring()
#empty = 52*" "
f = open(file2,'w')
tpl = ((npart,24),(massarr,48),(atime,float64),(redshift,float64),(flag_sfr,int32),(flag_feedback,int32),(nall,24),(flag_cooling,int32),(num_files,int32),(boxsize,float64),(omega0,float64),(omegalambda,float64),(hubbleparam,float64),(flag_age,int32),(flag_metals,int32),(nallhw,24),(flag_entr_ic,int32),(critical_energy_spec,float64),(empty,52))
writeblock(f,tpl)
#writeblock(f,pos.astype(float32))
#writeblock(f,vel.astype(float32))
#writeblock(f,num.astype(int32))
#writeblock(f,u.astype(float32))
writeblock(f,pos.astype(float32),save_memory=1)
writeblock(f,vel.astype(float32),save_memory=1)
writeblock(f,num.astype(int32),save_memory=1)
writeblock(f,u.astype(float32),save_memory=1)
f.close()

Event Timeline

c4science · Help