Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85290766
addgas
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 28, 01:37
Size
9 KB
Mime Type
text/x-python
Expires
Mon, Sep 30, 01:37 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21155991
Attached To
rGTOOLS Gtools
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
Log In to Comment