diff --git a/Gtools/io.py b/Gtools/io.py index 7c59fbb..a690311 100644 --- a/Gtools/io.py +++ b/Gtools/io.py @@ -1,1481 +1,1534 @@ # standard modules import os,sys,string,copy,types # array module from numpy import * #import numarray.strings as str import pNbody as nbd from pNbody import param +import pickle + ################################# def read_params(file): ################################# ''' Read params Gadget file and return the content in a dictionary ''' f = open(file) lines = f.readlines() f.close() # remove empty lines lines = filter(lambda l:l!='\n', lines) # remove trailing lines = map(string.strip, lines) # remove comments lines = filter(lambda x:x[0]!='%', lines) # split lines elts = map(string.split,lines) # make dictionary params = {} for e in elts: try : params[e[0]]=float(e[1]) except ValueError: params[e[0]]= e[1] return params ################################# def read_info(file): ################################# ''' Read info Gadget file ''' f = open(file,'r') lines = f.readlines() f.close() lines = filter(lambda l:l!='\n',lines) lines = map(string.strip,lines) elts = map(string.split,lines) try: Step = array(map(lambda x:float(x[2][:-1]),elts)) Time = array(map(lambda x:float(x[4][:-1]),elts)) Systemstep = array(map(lambda x:float(x[6] ),elts)) #Numactive = array(map(lambda x:float(x[8] ),elts)) except: Step = array(map(lambda x:float(x[2][:-1]),elts)) Time = array(map(lambda x:float(x[4][:-1]),elts)) Systemstep = array(map(lambda x:float(x[8][:-1]),elts)) #return Step,Time,Systemstep,Numactive return Step,Time,Systemstep ################################# def read_cpu(file): ################################# ''' Read cpu Gadget file ''' Time = array([]) Step = array([]) CPUs = array([]) CPU_Total = array([]) CPU_Gravity = array([]) CPU_Hydro = array([]) CPU_Domain = array([]) CPU_Potential = array([]) CPU_Predict = array([]) CPU_TimeLine = array([]) CPU_Snapshot = array([]) CPU_TreeWalk = array([]) CPU_TreeConstruction = array([]) CPU_CommSum = array([]) CPU_Imbalance = array([]) CPU_HydCompWalk = array([]) CPU_HydCommSumm = array([]) CPU_HydImbalance = array([]) CPU_EnsureNgb = array([]) CPU_PM = array([]) CPU_Peano = array([]) f = open(file,'r') lines = f.readlines() f.close() lines = map(string.strip,lines) lines1 = filter(lambda l:l[:4]=='Step',lines) lines2 = filter(lambda l:l[:4]!='Step',lines) elts1 = map(string.split,lines1) Step = array(map(lambda x:float(x[1][:-1]),elts1))[1:] Time = array(map(lambda x:float(x[3][:-1]),elts1))[1:] CPUs = array(map(lambda x:float(x[5] ),elts1))[1:] elts2 = array(map(lambda l:map(float,string.split(l)),lines2)) CPU_Total = elts2[:,0][1:] CPU_Gravity = elts2[:,1][1:] CPU_Hydro = elts2[:,2][1:] CPU_Domain = elts2[:,3][1:] CPU_Potential = elts2[:,4][1:] CPU_Predict = elts2[:,5][1:] CPU_TimeLine = elts2[:,6][1:] CPU_Snapshot = elts2[:,7][1:] CPU_TreeWalk = elts2[:,8][1:] CPU_TreeConstruction = elts2[:,9][1:] CPU_CommSum = elts2[:,10][1:] CPU_Imbalance = elts2[:,11][1:] CPU_HydCompWalk = elts2[:,12][1:] CPU_HydCommSumm = elts2[:,13][1:] CPU_HydImbalance = elts2[:,14][1:] CPU_EnsureNgb = elts2[:,15][1:] CPU_PM = elts2[:,16][1:] CPU_Peano = elts2[:,17][1:] return Step,Time,CPUs,CPU_Total,CPU_Gravity,CPU_Hydro,CPU_Domain,CPU_Potential,CPU_Predict,CPU_TimeLine,CPU_Snapshot,CPU_TreeWalk,CPU_TreeConstruction,CPU_CommSum,CPU_Imbalance,CPU_HydCompWalk,CPU_HydCommSumm,CPU_HydImbalance,CPU_EnsureNgb,CPU_PM,CPU_Peano ################################# def read_new_energy(file): ################################# ''' Read energy energy Gadget file ''' def toFloatList(l): return map(float,l) # read the file f = open(file,'r') header = f.readline() if header[0] != '#': raise "NotNewEgyFileError","file %s not a new energy file."%(file) # create dict from header header = string.strip(header[2:]) elts = string.split(header) iobs = {} i = 0 for elt in elts: iobs[elt]=i i = i + 1 # read end of file lines = f.readlines() f.close() # remove trailing lines = map(string.strip, lines) # split lines = map(string.split, lines) # convert into float lines = map(toFloatList, lines) # convert into array lines = array(map(array, lines)) # convert into array lines = transpose(lines) vals = {} for key in iobs.keys(): vals[key] = lines[iobs[key]] return vals ################################# def read_energy(file,iobs=None): ################################# ''' Read energy Gadget file ''' def toFloatList(l): return map(float,l) def GetNelts(self): f = open(file,'r') line = f.readline() nelts = len(string.split(line)) return nelts nlts = GetNelts(file) if nlts==28: iobs = {'Time':0,'EnergyInt' : 1,'EnergyPot' : 2,'EnergyKin' : 3, 'EnergyInt1': 4,'EnergyPot1': 5,'EnergyKin1': 6, 'EnergyInt2': 7,'EnergyPot2': 8,'EnergyKin2': 9, 'EnergyInt3':10,'EnergyPot3':11,'EnergyKin3':12, 'EnergyInt4':13,'EnergyPot4':14,'EnergyKin4':15, 'EnergyInt5':16,'EnergyPot5':17,'EnergyKin5':18, 'EnergyInt6':19,'EnergyPot6':20,'EnergyKin6':21, 'MassComp1' :22,'MassComp2' :23,'MassComp3':24,'MassComp4':25,'MassComp5':26,'MassComp6':27} elif nlts==35: iobs = {'Time':0,'EnergyInt' : 1,'EnergyRadSph' : 2,'EnergyPot' : 3,'EnergyKin' : 4, 'EnergyInt1': 5,'EnergyRadSph1': 6,'EnergyPot1': 7,'EnergyKin1': 8, 'EnergyInt2': 9,'EnergyRadSph2':10,'EnergyPot2':11,'EnergyKin2':12, 'EnergyInt3':13,'EnergyRadSph3':14,'EnergyPot3':15,'EnergyKin3':16, 'EnergyInt4':17,'EnergyRadSph4':18,'EnergyPot4':19,'EnergyKin4':20, 'EnergyInt5':21,'EnergyRadSph5':22,'EnergyPot5':23,'EnergyKin5':24, 'EnergyInt6':25,'EnergyRadSph6':26,'EnergyPot6':27,'EnergyKin6':28, 'MassComp1' :29,'MassComp2':30,'MassComp3':31,'MassComp4':32,'MassComp5':33,'MassComp6':34} elif nlts==42: iobs = {'Time':0,'EnergyInt' : 1,'EnergyRadSph' : 2,'EnergyRadSticky' : 3,'EnergyPot' : 4,'EnergyKin' : 5, 'EnergyInt1': 6,'EnergyRadSph1': 7,'EnergyRadSticky1': 8,'EnergyPot1': 9,'EnergyKin1':10, 'EnergyInt2':11,'EnergyRadSph2':12,'EnergyRadSticky2':13,'EnergyPot2':14,'EnergyKin2':15, 'EnergyInt3':16,'EnergyRadSph3':17,'EnergyRadSticky3':18,'EnergyPot3':19,'EnergyKin3':20, 'EnergyInt4':21,'EnergyRadSph4':22,'EnergyRadSticky4':23,'EnergyPot4':24,'EnergyKin4':25, 'EnergyInt5':26,'EnergyRadSph5':27,'EnergyRadSticky5':28,'EnergyPot5':29,'EnergyKin5':30, 'EnergyInt6':31,'EnergyRadSph6':32,'EnergyRadSticky6':33,'EnergyPot6':34,'EnergyKin6':35, 'MassComp1' :36,'MassComp2':37,'MassComp3':38,'MassComp4':39,'MassComp5':40,'MassComp6':41} f = open(file) lines = f.readlines() f.close() # remove trailing lines = map(string.strip, lines) # split lines = map(string.split, lines) # convert into float lines = map(toFloatList, lines) # convert into array lines = array(map(array, lines)) # convert into array lines = transpose(lines) vals = {} for key in iobs.keys(): vals[key] = lines[iobs[key]] return vals ################################# def read_phase(file): ################################# ''' Read phase Gadget file ''' f = open(file,'r') f.readline() lines = f.readlines() f.close() lines = filter(lambda l:l!='\n',lines) lines = map(string.strip,lines) elts = map(string.split,lines) Step = array(map(lambda x:float(x[1][:-1]),elts)) Time = array(map(lambda x:float(x[3] ),elts)) GasPart = array(map(lambda x:float(x[5] ),elts)) SphPart = array(map(lambda x:float(x[7] ),elts)) StickyPart = array(map(lambda x:float(x[9] ),elts)) DarkPart = array(map(lambda x:float(x[11] ),elts)) return Step,Time,GasPart,SphPart,StickyPart,DarkPart ################################# def read_accretion(file): ################################# ''' Read accretion Gadget file ''' f = open(file,'r') lines = f.readlines() f.close() lines = filter(lambda l:l!='\n',lines) lines = map(string.strip,lines) elts = map(string.split,lines) elts = str.array(elts) vals = {} for i in range(elts.shape[1]/2): name = elts[:,2*i][0] if name[-1]==':': name = name[:-1] vect = elts[:,2*i+1].eval() if name != 'Step': vals[name] = vect.astype(float) else: vals[name] = vect return vals ################################# def read_bondi(file): ################################# ''' Read bondi Gadget file ''' f = open(file,'r') lines = f.readlines() f.close() lines = filter(lambda l:l!='\n',lines) lines = map(string.strip,lines) elts = map(string.split,lines) elts = str.array(elts) vals = {} for i in range(elts.shape[1]/2): name = elts[:,2*i][0] if name[-1]==':': name = name[:-1] vect = elts[:,2*i+1].eval() if name != 'Step': vals[name] = vect.astype(float) else: vals[name] = vect return vals ################################# def read_bubble(file): ################################# ''' Read bubble Gadget file ''' f = open(file,'r') lines = f.readlines() f.close() lines = filter(lambda l:l!='\n',lines) lines = map(string.strip,lines) elts = map(string.split,lines) elts = str.array(elts) vals = {} for i in range(elts.shape[1]/2): name = elts[:,2*i][0] if name[-1]==':': name = name[:-1] vect = elts[:,2*i+1].eval() if name != 'Step': vals[name] = vect.astype(float) else: vals[name] = vect return vals ################################# def read_cooling(file): ################################# ''' Read cooling file ''' f = open(file,'r') f.readline() f.readline() lines = f.readlines() f.close() lines = map(string.strip,lines) elts = map(string.split,lines) logT = array(map(lambda x:float(x[0]),elts)) logL0 = array(map(lambda x:float(x[1]),elts)) logL1 = array(map(lambda x:float(x[2]),elts)) logL2 = array(map(lambda x:float(x[3]),elts)) logL3 = array(map(lambda x:float(x[4]),elts)) logL4 = array(map(lambda x:float(x[5]),elts)) logL5 = array(map(lambda x:float(x[6]),elts)) logL6 = array(map(lambda x:float(x[7]),elts)) return logT,logL0,logL1,logL2,logL3,logL4,logL5,logL6 ################################# def read_den(file): ################################# ''' read .den file (hop) ''' f = open(file) N = fromstring(f.read(4),int32)[0] # number of particles dens = fromstring(f.read(4*N),float32) # density of all particles f.close() return dens ################################# def read_tag(file): ################################# ''' read .tag file (hop) ''' f = open(file) N = fromstring(f.read(4),int32)[0] # number of particles Ng = fromstring(f.read(4),int32)[0] # number of groups tags = fromstring(f.read(4*N),int32) # group of all particles f.close() return Ng,tags ################################# def read_size(file): ################################# ''' read .size file (hop) ''' f = open(file) N = int(f.readline()) # number of particles Npg = int(f.readline()) # number of particles in groups Ng = int(f.readline()) # number of groups lines = f.readlines() f.close() lines = map(string.strip,lines) elts = map(string.split,lines) gid = array(map(lambda x:int(x[0]),elts)) gnp = array(map(lambda x:int(x[1]),elts)) return N,Npg,Ng,gid,gnp ################################# def read_gbound(file): ################################# ''' read .gbound file (hop) ''' f = open(file) Ng = int(f.readline()) # number of groups f.readline() f.readline() f.readline() f.readline() f.readline() f.readline() f.readline() f.readline() f.readline() f.readline() lines = f.readlines() f.close() lines = lines[:Ng] lines = map(string.strip,lines) elts = map(string.split,lines) gid = array(map(lambda x:int(x[0]),elts)) gnp = array(map(lambda x:int(x[1]),elts)) gx = array(map(lambda x:float(x[3]),elts)) gy = array(map(lambda x:float(x[4]),elts)) gz = array(map(lambda x:float(x[5]),elts)) return gid,gnp,gx,gy,gz ################################# def read_track(file): ################################# ''' read .trk file return a dictionary of elements ''' f = open(file,'r') lines = f.readlines() f.close() elts = {} for line in lines: line = string.strip(line) name = string.split(line)[0] val = string.split(line)[2] try: elts['%s'%(name)] = float(val) except: elts['%s'%(name)] = val return elts ################################# def read_tracks(files): ################################# ''' read a set of .trk file return a global dictionary ''' ################################################ # read first file and create main dic ################################################ elts = read_track(files[0]) mainDic = {} for key in elts.keys(): if type(elts[key]) == types.FloatType: mainDic[key] = array([],float) else: mainDic[key] = [] ################################################ # read all files ################################################ for file in files: elts = read_track(file) for key in elts.keys(): if mainDic.has_key(key): if type(elts[key]) == float: mainDic[key] = concatenate( ( mainDic[key], array([elts[key]]) ) ) else: mainDic[key].append(elts[key]) return mainDic ################################# def write_nbodyparam(file,dic,eye,size): ################################# """ write an nbody parameters file """ def RotateAround(angle,axis,point,pos): x = pos # center point x = x-point # construction of the rotation matrix norm = sqrt(axis[0]**2 + axis[1]**2 + axis[2]**2) if norm == 0: return x sn = sin(-angle/2.) e0 = cos(-angle/2.) e1 = axis[0]*sn/norm e2 = axis[1]*sn/norm e3 = axis[2]*sn/norm a = zeros((3,3),float) a[0,0] = e0**2 + e1**2 - e2**2 - e3**2 a[1,0] = 2.*(e1*e2 + e0*e3) a[2,0] = 2.*(e1*e3 - e0*e2) a[0,1] = 2.*(e1*e2 - e0*e3) a[1,1] = e0**2 - e1**2 + e2**2 - e3**2 a[2,1] = 2.*(e2*e3 + e0*e1) a[0,2] = 2.*(e1*e3 + e0*e2) a[1,2] = 2.*(e2*e3 - e0*e1) a[2,2] = e0**2 - e1**2 - e2**2 + e3**2 a = a.astype(float) # multiply x and v x = dot(x,a) # decenter point x = x+point pos = x return pos params = param.Params(nbd.PARAMETERFILE,None) gwinShapeX = size[0] gwinShapeY = size[1] ############################# # projection param ############################# if (dic['Observer_0:ProjectionMode']==1): # frutum gwinPerspectiveTop = dic['Observer_0:PerspectiveNear']*tan(dic['Observer_0:PerspectiveFov']*0.5*pi/180.) gwinPerspectiveRight = gwinPerspectiveTop*float(gwinShapeX)/float(gwinShapeY) else: # ortho gwinPerspectiveRight = dic['Observer_0:PerspectiveRight'] gwinPerspectiveTop = gwinShapeY/gwinShapeX*gwinPerspectiveRight gwinClip1 = dic['Observer_0:PerspectiveNear'] gwinClip2 = dic['Observer_0:PerspectiveFar'] M = zeros((4,3),float) axis = zeros((3,),float) point = zeros((3,),float) M[0,0] = dic['Observer_0:M0'] M[0,1] = dic['Observer_0:M1'] M[0,2] = dic['Observer_0:M2'] M[1,0] = dic['Observer_0:M4'] M[1,1] = dic['Observer_0:M5'] M[1,2] = dic['Observer_0:M6'] M[2,0] = dic['Observer_0:M8'] M[2,1] = dic['Observer_0:M9'] M[2,2] = dic['Observer_0:M10'] M[3,0] = dic['Observer_0:M12'] M[3,1] = dic['Observer_0:M13'] M[3,2] = dic['Observer_0:M14'] # compute stereo axis = M[3]-M[0]; r = sqrt(axis[0]*axis[0]+axis[1]*axis[1]+axis[2]*axis[2]); axis[0] = axis[0]/r; axis[1] = axis[1]/r; axis[2] = axis[2]/r; if (eye=='right'): M = M + axis* dic['Observer_0:EyeDist'] elif (eye=='left'): M = M - axis* dic['Observer_0:EyeDist'] else: pass # rotate, to be compatible with gwin EYE = 0 PTS = 4 axis = M[0] - M[1] point = M[0] M = RotateAround(pi/2.,axis,point,M); N = M M = zeros((16,),float) M[0] = N[0,0] M[1] = N[0,1] M[2] = N[0,2] M[3] = 0 M[4] = N[1,0] M[5] = N[1,1] M[6] = N[1,2] M[7] = 0 M[8] = N[2,0] M[9] = N[2,1] M[10] = N[2,2] M[11] = 0 M[12] = N[3,0] M[13] = N[3,1] M[14] = N[3,2] M[15] = 0 dist = sqrt( pow(M[0]-M[8],2) + pow(M[1]-M[9],2) + pow(M[2]-M[10],2)); # head obs1 = M[0]; obs2 = M[1]; obs3 = M[2]; # lookat point */ norm = sqrt( pow(M[0]-M[4],2) + pow(M[1]-M[5],2) + pow(M[2]-M[6],2)); obs4 = obs1 + (M[4] -obs1)/norm*dist; obs5 = obs2 + (M[5] -obs2)/norm*dist; obs6 = obs3 + (M[6] -obs3)/norm*dist; # arm norm = sqrt( pow(M[0]-M[8],2) + pow(M[1]-M[9],2) + pow(M[2]-M[10],2)); obs7 = obs1 + (M[8] - obs1)/norm; obs8 = obs2 + (M[9] - obs2)/norm; obs9 = obs3 + (M[10]- obs3)/norm; # head norm = sqrt( pow(M[0]-M[12],2) + pow(M[1]-M[13],2) + pow(M[2]-M[14],2)); obs10 = obs1 + (M[12] - obs1)/norm; obs11 = obs2 + (M[13] - obs2)/norm; obs12 = obs3 + (M[14] - obs3)/norm; obs = array([obs1,obs2,obs3,obs4,obs5,obs6,obs7,obs8,obs9,obs10,obs11,obs12]) cut = 'yes' dist_eye = None foc = None if dic['Observer_0:ProjectionMode']==1: persp = 'on' else: persp = 'off' center = (0.0, 0.0, 0.0) frsp = 0.0 space = 'pos' mode = 'm' params.set("obs",obs) params.set("x0",None) params.set("xp",None) params.set("alpha",None) params.set("view",None) params.set("r_obs",dist) params.set("clip",(gwinClip1,gwinClip2)) params.set("cut",cut) params.set("eye",None) params.set("dist_eye",None) params.set("foc",foc) params.set("persp",persp) params.set("shape",(gwinShapeX,gwinShapeY)) params.set("center",center) params.set("size",(gwinPerspectiveRight,gwinPerspectiveTop)) params.set("frsp",frsp) params.set("space",space) params.set("mode",mode) params.save(file) ################################# def write_cinit_stats(file,stats): ################################# """ write stats from cinit """ R = stats['R'] sr = stats['sr'] sp = stats['sp'] sz = stats['sz'] vc = stats['vc'] vm = stats['vm'] kappa = stats['kappa'] omega = stats['omega'] nu = stats['nu'] Sdens = stats['Sdens'] Q = stats['Q'] f = open(file,'w') line = "# R sr sp sz vc vm kappa omega nu Sdens Q\n" f.write(line) for i in range(len(R)): line = "%g %g %g %g %g %g %g %g %g %g %g\n"%(R[i],sr[i],sp[i],sz[i],vc[i],vm[i],kappa[i],omega[i],nu[i],Sdens[i],Q[i]) f.write(line) f.close() ################################# def write_cinit_stats_new(file,stats): ################################# """ write stats from cinit """ R = stats['R'] sr = stats['sr'] sp = stats['sp'] sz = stats['sz'] vc = stats['vc'] vm = stats['vm'] kappa = stats['kappa'] omega = stats['omega'] nu = stats['nu'] Sdens = stats['Sdens'] Sdensd= stats['Sdensd'] Q = stats['Q'] Ar = stats['Ar'] f = open(file,'w') line = "# R sr sp sz vc vm kappa omega nu Sdens Sdensd Q Ar\n" f.write(line) for i in range(len(R)): line = "%g %g %g %g %g %g %g %g %g %g %g %g %g\n"%(R[i],sr[i],sp[i],sz[i],vc[i],vm[i],kappa[i],omega[i],nu[i],Sdens[i],Sdensd[i],Q[i],Ar[i]) f.write(line) f.close() ################################# def read_cinit_stats(file): ################################# """ read stats from cinit """ def toFloatList(l): return map(float,l) # read the file f = open(file,'r') header = f.readline() if header[0] != '#': raise "NotCinitStatsFileError","file %s not a cinit stats file."%(file) # create dict from header header = string.strip(header[2:]) elts = string.split(header) iobs = {} i = 0 for elt in elts: iobs[elt]=i i = i + 1 # read end of file lines = f.readlines() f.close() # remove trailing lines = map(string.strip, lines) # split lines = map(string.split, lines) # convert into float lines = map(toFloatList, lines) # convert into array lines = array(map(array, lines)) # convert into array lines = transpose(lines) vals = {} for key in iobs.keys(): vals[key] = lines[iobs[key]] return vals ################################# def write_crv(file,r,v): ################################# """ write crv file """ f = open(file,'w') line = "# R Vc\n" f.write(line) for i in range(len(r)): line = "%g %g\n"%(r[i],v[i]) f.write(line) f.close() ################################# def read_crv(file): ################################# """ read crv from cinit """ def toFloatList(l): return map(float,l) # read the file f = open(file,'r') header = f.readline() if header[0] != '#': raise "NotCrvFileError","file %s not a crv stats file."%(file) # create dict from header header = string.strip(header[2:]) elts = string.split(header) iobs = {} i = 0 for elt in elts: iobs[elt]=i i = i + 1 # read end of file lines = f.readlines() f.close() # remove trailing lines = map(string.strip, lines) # split lines = map(string.split, lines) # convert into float lines = map(toFloatList, lines) # convert into array lines = array(map(array, lines)) # convert into array lines = transpose(lines) vals = {} for key in iobs.keys(): vals[key] = lines[iobs[key]] return vals['R'],vals['Vc'] ################################# def write_crv2(file,r,v): ################################# """ write crv2 file """ f = open(file,'w') line = "# R Vc2\n" f.write(line) for i in range(len(r)): line = "%g %g\n"%(r[i],v[i]) f.write(line) f.close() ################################# def read_crv2(file): ################################# """ read crv2 from cinit """ def toFloatList(l): return map(float,l) # read the file f = open(file,'r') header = f.readline() if header[0] != '#': raise "NotCrv2FileError","file %s not a crv2 stats file."%(file) # create dict from header header = string.strip(header[2:]) elts = string.split(header) iobs = {} i = 0 for elt in elts: iobs[elt]=i i = i + 1 # read end of file lines = f.readlines() f.close() # remove trailing lines = map(string.strip, lines) # split lines = map(string.split, lines) # convert into float lines = map(toFloatList, lines) # convert into array lines = array(map(array, lines)) # convert into array lines = transpose(lines) vals = {} for key in iobs.keys(): vals[key] = lines[iobs[key]] return vals['R'],vals['Vc2'] ################################# def write_sden(file,r,v): ################################# """ write crv2 file """ f = open(file,'w') line = "# R sden\n" f.write(line) for i in range(len(r)): line = "%g %g\n"%(r[i],v[i]) f.write(line) f.close() ################################# def read_sden(file): ################################# """ read crv2 from cinit """ def toFloatList(l): return map(float,l) # read the file f = open(file,'r') header = f.readline() if header[0] != '#': raise "NotSden2FileError","file %s not a sden stats file."%(file) # create dict from header header = string.strip(header[2:]) elts = string.split(header) iobs = {} i = 0 for elt in elts: iobs[elt]=i i = i + 1 # read end of file lines = f.readlines() f.close() # remove trailing lines = map(string.strip, lines) # split lines = map(string.split, lines) # convert into float lines = map(toFloatList, lines) # convert into array lines = array(map(array, lines)) # convert into array lines = transpose(lines) vals = {} for key in iobs.keys(): vals[key] = lines[iobs[key]] return vals['R'],vals['sden'] ################################# def read_olddiagnostic(file): ################################# """ read DIAGNOSTIC file created from treeasph """ def read_header(f): for i in range(20): f.readline() def read_step(f): line = f.readline() # Time line = f.readline() Time = float(line[7:20]) f.readline() # nstot, min, max, avg line = f.readline() nstot = int(line[30:38]) min = int(line[43:48]) max = int(line[53:58]) avg = int(line[63:68]) f.readline() f.readline() # mtot, n_stars line = f.readline() mtot = float(line[14:26]) n_stars = int(line[38:45]) # Energy_tot, Energy_kin line = f.readline() Energy_tot = float(line[20:37]) Energy_kin = float(line[52:69]) # Energy_pot, Energy_the line = f.readline() Energy_pot = float(line[20:37]) Energy_the = float(line[52:69]) # Entropy_T, Virial_rat line = f.readline() Entropy_T = float(line[20:37]) Virial_rat = float(line[52:69]) # Disipation, Disipat_sf line = f.readline() Disipation = float(line[20:37]) Disipat_sf = float(line[52:69]) # Delta_instant_Energy line = f.readline() Delta_instant_Energy = float(line[30:47]) # Delta_total_Energy line = f.readline() Delta_total_Energy = float(line[30:47]) f.readline() # amx, amy, amz line = f.readline() amx = float(line[23:40]) amy = float(line[41:57]) amz = float(line[58:74]) # cmpos line = f.readline() cmposx = float(line[15:32]) cmposy = float(line[33:49]) cmposz = float(line[50:66]) # cmvel line = f.readline() cmvelx = float(line[15:32]) cmvely = float(line[33:49]) cmvelz = float(line[50:66]) f.readline() # cpu time per step line = f.readline() cpu_time_per_step = float(line[30:42]) f.readline() f.readline() return Time,nstot,min,max,avg,mtot,n_stars,Energy_tot,Energy_kin,Energy_pot,Energy_the,Entropy_T,Virial_rat,Disipation,Disipat_sf,Delta_instant_Energy,Delta_total_Energy,amx,amy,amz,cmposx,cmposy,cmposz,cmvelx,cmvely,cmvelz,cpu_time_per_step vTime = array([],float) vnstot = array([],float) vmin = array([],float) vmax = array([],float) vavg = array([],float) vmtot = array([],float) vn_stars = array([],float) vEnergy_tot = array([],float) vEnergy_kin = array([],float) vEnergy_pot = array([],float) vEnergy_the = array([],float) vEntropy_T = array([],float) vVirial_rat = array([],float) vDisipation = array([],float) vDisipat_sf = array([],float) vDelta_instant_Energy = array([],float) vDelta_total_Energy = array([],float) vamx = array([],float) vamy = array([],float) vamz = array([],float) vcmposx = array([],float) vcmposy = array([],float) vcmposz = array([],float) vcmvelx = array([],float) vcmvely = array([],float) vcmvelz = array([],float) vcpu_time_per_step = array([],float) f = open(file) read_header(f) while 1: try: output = read_step(f) except ValueError: break Time,nstot,min,max,avg,mtot,n_stars,Energy_tot,Energy_kin,Energy_pot,Energy_the,Entropy_T,Virial_rat,Disipation,Disipat_sf,Delta_instant_Energy,Delta_total_Energy,amx,amy,amz,cmposx,cmposy,cmposz,cmvelx,cmvely,cmvelz,cpu_time_per_step=output vTime = concatenate((vTime,[Time])) vnstot = concatenate((vnstot,[nstot])) vmin = concatenate((vmin,[min])) vmax = concatenate((vmax,[max])) vavg = concatenate((vavg,[avg])) vmtot = concatenate((vmtot,[mtot])) vn_stars = concatenate((vn_stars,[n_stars])) vEnergy_tot = concatenate((vEnergy_tot,[Energy_tot])) vEnergy_kin = concatenate((vEnergy_kin,[Energy_kin])) vEnergy_pot = concatenate((vEnergy_pot,[Energy_pot])) vEnergy_the = concatenate((vEnergy_the,[Energy_the])) vEntropy_T = concatenate((vEntropy_T,[Entropy_T])) vVirial_rat = concatenate((vVirial_rat,[Virial_rat])) vDisipation = concatenate((vDisipation,[Disipation])) vDisipat_sf = concatenate((vDisipat_sf,[Disipat_sf])) vDelta_instant_Energy = concatenate((vDelta_instant_Energy,[Delta_instant_Energy])) vDelta_total_Energy = concatenate((vDelta_total_Energy,[Delta_total_Energy])) vamx = concatenate((vamx,[amx])) vamy = concatenate((vamy,[amy])) vamz = concatenate((vamz,[amz])) vcmposx = concatenate((vcmposx,[cmposx])) vcmposy = concatenate((vcmposy,[cmposy])) vcmposz = concatenate((vcmposz,[cmposz])) vcmvelx = concatenate((vcmvelx,[cmvelx])) vcmvely = concatenate((vcmvely,[cmvely])) vcmvelz = concatenate((vcmvelz,[cmvelz])) vcpu_time_per_step = concatenate((vcpu_time_per_step,[cpu_time_per_step])) vals = {} vals['Time'] = vTime vals['EnergyTot'] = vEnergy_tot vals['EnergyKin'] = vEnergy_kin vals['EnergyPot'] = vEnergy_pot vals['EnergyInt'] = vEnergy_the - vDisipation - vDisipat_sf vals['EnergyRadSph'] = vDisipation vals['EnergySfr'] = vDisipat_sf f.close() return vals ################################# def read_diagnostic(file): ################################# """ read DIAGNOSTIC file created from treeasph """ def read_header(f): for i in range(20): f.readline() def read_step(f): line = f.readline() # Time line = f.readline() Time = float(line[7:20]) f.readline() # nstot, min, max, avg line = f.readline() nstot = int(line[30:38]) min = int(line[43:48]) max = int(line[53:58]) avg = int(line[63:68]) f.readline() f.readline() # mtot, n_stars line = f.readline() mtot = float(line[14:26]) n_stars = int(line[38:45]) # Energy_tot, Energy_kin line = f.readline() Energy_tot = float(line[20:37]) Energy_kin = float(line[52:69]) # Energy_pot, Energy_the line = f.readline() Energy_pot = float(line[20:37]) Energy_the = float(line[52:69]) # Entropy_T, Virial_rat line = f.readline() Entropy_T = float(line[20:37]) Virial_rat = float(line[52:69]) # Disipation, Disipat_sf line = f.readline() Disipation = float(line[20:37]) Disipat_sf = float(line[52:69]) # Disipat_fb line = f.readline() Disipat_fb = float(line[20:37]) # Delta_instant_Energy line = f.readline() Delta_instant_Energy = float(line[30:47]) # Delta_total_Energy line = f.readline() Delta_total_Energy = float(line[30:47]) f.readline() # amx, amy, amz line = f.readline() amx = float(line[23:40]) amy = float(line[41:57]) amz = float(line[58:74]) # cmpos line = f.readline() cmposx = float(line[15:32]) cmposy = float(line[33:49]) cmposz = float(line[50:66]) # cmvel line = f.readline() cmvelx = float(line[15:32]) cmvely = float(line[33:49]) cmvelz = float(line[50:66]) f.readline() # cpu time per step line = f.readline() cpu_time_per_step = float(line[30:42]) f.readline() f.readline() return Time,nstot,min,max,avg,mtot,n_stars,Energy_tot,Energy_kin,Energy_pot,Energy_the,Entropy_T,Virial_rat,Disipation,Disipat_sf,Disipat_fb,Delta_instant_Energy,Delta_total_Energy,amx,amy,amz,cmposx,cmposy,cmposz,cmvelx,cmvely,cmvelz,cpu_time_per_step vTime = array([],float) vnstot = array([],float) vmin = array([],float) vmax = array([],float) vavg = array([],float) vmtot = array([],float) vn_stars = array([],float) vEnergy_tot = array([],float) vEnergy_kin = array([],float) vEnergy_pot = array([],float) vEnergy_the = array([],float) vEntropy_T = array([],float) vVirial_rat = array([],float) vDisipation = array([],float) vDisipat_sf = array([],float) vDisipat_fb = array([],float) vDelta_instant_Energy = array([],float) vDelta_total_Energy = array([],float) vamx = array([],float) vamy = array([],float) vamz = array([],float) vcmposx = array([],float) vcmposy = array([],float) vcmposz = array([],float) vcmvelx = array([],float) vcmvely = array([],float) vcmvelz = array([],float) vcpu_time_per_step = array([],float) f = open(file) read_header(f) while 1: try: output = read_step(f) except ValueError: break Time,nstot,min,max,avg,mtot,n_stars,Energy_tot,Energy_kin,Energy_pot,Energy_the,Entropy_T,Virial_rat,Disipation,Disipat_sf,Disipat_fb,Delta_instant_Energy,Delta_total_Energy,amx,amy,amz,cmposx,cmposy,cmposz,cmvelx,cmvely,cmvelz,cpu_time_per_step=output vTime = concatenate((vTime,[Time])) vnstot = concatenate((vnstot,[nstot])) vmin = concatenate((vmin,[min])) vmax = concatenate((vmax,[max])) vavg = concatenate((vavg,[avg])) vmtot = concatenate((vmtot,[mtot])) vn_stars = concatenate((vn_stars,[n_stars])) vEnergy_tot = concatenate((vEnergy_tot,[Energy_tot])) vEnergy_kin = concatenate((vEnergy_kin,[Energy_kin])) vEnergy_pot = concatenate((vEnergy_pot,[Energy_pot])) vEnergy_the = concatenate((vEnergy_the,[Energy_the])) vEntropy_T = concatenate((vEntropy_T,[Entropy_T])) vVirial_rat = concatenate((vVirial_rat,[Virial_rat])) vDisipation = concatenate((vDisipation,[Disipation])) vDisipat_sf = concatenate((vDisipat_sf,[Disipat_sf])) vDisipat_fb = concatenate((vDisipat_fb,[Disipat_fb])) vDelta_instant_Energy = concatenate((vDelta_instant_Energy,[Delta_instant_Energy])) vDelta_total_Energy = concatenate((vDelta_total_Energy,[Delta_total_Energy])) vamx = concatenate((vamx,[amx])) vamy = concatenate((vamy,[amy])) vamz = concatenate((vamz,[amz])) vcmposx = concatenate((vcmposx,[cmposx])) vcmposy = concatenate((vcmposy,[cmposy])) vcmposz = concatenate((vcmposz,[cmposz])) vcmvelx = concatenate((vcmvelx,[cmvelx])) vcmvely = concatenate((vcmvely,[cmvely])) vcmvelz = concatenate((vcmvelz,[cmvelz])) vcpu_time_per_step = concatenate((vcpu_time_per_step,[cpu_time_per_step])) vals = {} vals['Time'] = vTime vals['EnergyTot'] = vEnergy_tot vals['EnergyKin'] = vEnergy_kin vals['EnergyPot'] = vEnergy_pot vals['EnergyInt'] = vEnergy_the - vDisipation - vDisipat_sf - vDisipat_fb vals['EnergyRadSph'] = vDisipation vals['EnergySfr'] = vDisipat_sf vals['EnergyFeedback'] = vDisipat_fb f.close() return vals +##################################################### +def write_dmp(file,data): +##################################################### + ''' + Write a dmp (pickle) file. In other word, + dump the data object. + + Parameters + ---------- + file : the path to a file + data : a pickable python object + + Examples + -------- + >>> x = {'a':1,'b':2} + >>> io.write_dmp('/tmp/afile.dmp',x) + ''' + f = open(file,'w') + pickle.dump(data, f) + f.close() + + +##################################################### +def read_dmp(file): +##################################################### + ''' + Read a dmp (pickle) file. + + Parameters + ---------- + file : the path to a file + + Returns + ------- + data : a python object + + Examples + -------- + >>> x = {'a':1,'b':2} + >>> io.write_dmp('/tmp/afile.dmp',x) + >>> y = io.read_dmp('/tmp/afile.dmp') + >>> y + {'a': 1, 'b': 2} + ''' + + f = open(file,'r') + data = pickle.load(f) + f.close() + return data + +