diff --git a/config/formats/gadget.py b/config/formats/gadget.py index 74f91bb..4d5db54 100644 --- a/config/formats/gadget.py +++ b/config/formats/gadget.py @@ -1,1200 +1,1388 @@ #################################################################################################################################### # # GADGET CLASS # #################################################################################################################################### class Nbody_gadget(NbodyDefault): def __init__(self,p_name=None,pos=None,vel=None,mass=None,num=None,tpe=None,ftype=None,status='old',byteorder=sys.byteorder,pio='no',local=False,log=None,unitsfile=None): NbodyDefault.__init__(self,p_name=p_name,pos=pos,vel=vel,mass=mass,num=num,tpe=tpe,ftype=ftype,status=status,byteorder=byteorder,pio=pio,local=local,log=log,unitsfile=unitsfile) def InitSpec(self): - self.getComovingIntegration() + self.initComovingIntegration() def get_read_fcts(self): return [self.read_particles] def get_write_fcts(self): return [self.write_particles] def get_mxntpe(self): return 6 def get_default_spec_vars(self): ''' return specific variables default values for the class ''' return {'massarr' :array([0,0,self.nbody,0,0,0]), 'atime' :0., 'redshift' :0., 'flag_sfr' :0, 'flag_feedback' :0, 'nall' :array([0,0,self.nbody,0,0,0]), 'flag_cooling' :0, 'num_files' :1, 'boxsize' :0., 'omega0' :0., 'omegalambda' :0., 'hubbleparam' :0., 'flag_age' :0., 'hubbleparam' :0., 'flag_metals' :0., 'nallhw' :array([0,0,0,0,0,0]), 'flag_entr_ic' :0, 'flag_chimie_extraheader':0, 'critical_energy_spec' :0., 'empty' :48*'', 'comovingintegration' :None } def get_massarr_and_nzero(self): """ return massarr and nzero !!! when used in //, if only a proc has a star particle, !!! nzero is set to 1 for all cpu, while massarr has a length of zero !!! """ if self.has_var('massarr') and self.has_var('nzero'): if self.massarr !=None and self.nzero!=None: - print "warning : get_massarr_and_nzero : here we use massarr and nzero",self.massarr,self.nzero + if mpi.mpi_IsMaster(): + print "warning : get_massarr_and_nzero : here we use massarr and nzero",self.massarr,self.nzero return self.massarr,self.nzero massarr = zeros(len(self.npart),float) nzero = 0 # for each particle type, see if masses are equal for i in range(len(self.npart)): first_elt = sum((arange(len(self.npart)) < i) * self.npart) last_elt = first_elt + self.npart[i] if first_elt!=last_elt: c = (self.mass[first_elt]==self.mass[first_elt:last_elt]).astype(int) if sum(c)==len(c): massarr[i] = self.mass[first_elt] else: nzero = nzero+len(c) return massarr.tolist(),nzero - def getComovingIntegration(self): + def initComovingIntegration(self): """ set true if the file has been runned using the comoving integration scheme """ - flag = True + + # do nothing, if the value is already set + if self.comovingintegration!=None: + return + flag = True + # this is not very good, however, there is not other choice... if self.omega0==0 and self.omegalambda==0: flag = False self.comovingintegration = flag - print "comovingintegration=",self.comovingintegration - + #print "init : ComovingIntegration=",self.comovingintegration + def isComovingIntegrationOn(self): """ return true if the file has been runned using the comoving integration scheme """ - - flag = True - - # this is not very good, however, there is not other choice... - if self.omega0==0 and self.omegalambda==0: - flag = False - - self.comovingintegration = flag - - return flag + return self.comovingintegration def setComovingIntegrationOn(self): - self.comovingintegration = True + self.comovingintegration = True def setComovingIntegrationOff(self): - self.comovingintegration = False + self.comovingintegration = False - - + def ComovingIntegrationInfo(self): + if self.isComovingIntegrationOn(): + print "ComovingIntegration" + print " on (a=%5.3f h=%5.3f)"%(self.atime,self.hubbleparam) + else: + print "ComovingIntegration" + print " off" + def read_particles(self,f): ''' read gadget file ''' ########################################## # read the header and send it to each proc ########################################## tpl = (24,48,float,float,int32,int32,24,int32,int32,float,float,float,float,int32,int32,24,int32,int32,float,48) header = io.ReadBlock(f,tpl,byteorder=self.byteorder,pio=self.pio) 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,flag_chimie_extraheader,critical_energy_spec,empty = header if fabs(flag_chimie_extraheader)>1: # if the header is empty, it may be indetermined flag_chimie_extraheader=0 npart = fromstring(npart,int32) massarr = fromstring(massarr,float) nall = fromstring(nall,int32) nallhw = fromstring(nallhw,int32) if sys.byteorder != self.byteorder: npart.byteswap(True) massarr.byteswap(True) nall.byteswap(True) nallhw.byteswap(True) if flag_metals: # gas metal properties NELEMENTS = flag_metals self.NELEMENTS = NELEMENTS else: NELEMENTS = 0 self.NELEMENTS = NELEMENTS ########################################## # computes nzero ########################################## # count number of particles that have non constant masses and then # have masses storded further if self.pio == 'no': npart_tot = npart npart_all = libutil.get_npart_all(npart,mpi.mpi_NTask()) npart = npart_all[mpi.mpi_ThisTask()] # local npart_read = npart_tot nbody_read = sum(npart_read) npart_m_read = npart_tot * (massarr==0) ngas = npart[0] ngas_read = npart_tot[0] nstars = npart[1] nstars_read = npart_tot[1] # compute nzero nzero=0 mass = array([]) for i in range(len(npart_tot)): if massarr[i] == 0: nzero = nzero+npart_tot[i] else: mass = concatenate((mass,ones(npart[i])*massarr[i])) else: npart_tot = mpi.mpi_allreduce(npart) npart_all = None # each proc read for himself npart = npart # local npart_read = None # each proc read for himself nbody_read = sum(npart) npart_m_read = None # each proc read for himself ngas = npart[0] ngas_read = ngas nstars = npart[1] nstars_read = nstars # compute nzero nzero=0 mass = array([]) for i in range(len(npart)): if massarr[i] == 0: nzero = nzero+npart[i] else: mass = concatenate((mass,ones(npart[i])*massarr[i])) nbody = sum(npart) nbody_tot = sum(npart_tot) if npart_m_read !=None: if sum(npart_m_read) != nzero: raise "sum(npart_m) (%d) != nzero (%d)"%(sum(npart_m_read),nzero),npart_m_read ########################################## # optionnally read extra header ########################################## if flag_chimie_extraheader==1: if mpi.mpi_IsMaster(): print "reading chimie extra-header..." tpl = (int32, int(NELEMENTS)*4, 256-4-int(NELEMENTS)*4) - nelts,ChimieSolarAbundances,labels = io.ReadBlock(f,tpl,byteorder=self.byteorder,pio=self.pio) + nelts,ChimieSolarMassAbundances,labels = io.ReadBlock(f,tpl,byteorder=self.byteorder,pio=self.pio) nelts = int(nelts) self.ChimieNelements = nelts ChimieElements = string.split(labels,',')[:nelts] self.ChimieElements = ChimieElements - ChimieSolarAbundances = fromstring(ChimieSolarAbundances,float32) - self.ChimieSolarAbundances = {} + ChimieSolarMassAbundances = fromstring(ChimieSolarMassAbundances,float32) + self.ChimieSolarMassAbundances = {} for i,elt in enumerate(self.ChimieElements): - self.ChimieSolarAbundances[elt] = ChimieSolarAbundances[i] + self.ChimieSolarMassAbundances[elt] = ChimieSolarMassAbundances[i] else: self.ChimieNelements = int(NELEMENTS) self.ChimieElements = ['Fe','Mg','O','Metals'] - self.ChimieSolarAbundances = {} - self.ChimieSolarAbundances['Fe'] = 0.001771 - self.ChimieSolarAbundances['Mg'] = 0.00091245 - self.ChimieSolarAbundances['O'] = 0.0108169 - self.ChimieSolarAbundances['Metals'] = 0.02 + self.ChimieSolarMassAbundances = {} + self.ChimieSolarMassAbundances['Fe'] = 0.001771 + self.ChimieSolarMassAbundances['Mg'] = 0.00091245 + self.ChimieSolarMassAbundances['O'] = 0.0108169 + self.ChimieSolarMassAbundances['Metals'] = 0.02 ########################################## # read and send particles attribute ########################################## if mpi.mpi_IsMaster(): print "reading pos..." pos = io.ReadDataBlock(f,float32,shape=(nbody_read,3),byteorder=self.byteorder,pio=self.pio,npart=npart_read) if mpi.mpi_IsMaster(): print "reading vel..." vel = io.ReadDataBlock(f,float32,shape=(nbody_read,3),byteorder=self.byteorder,pio=self.pio,npart=npart_read) if mpi.mpi_IsMaster(): print "reading num..." num = io.ReadDataBlock(f,int32 ,shape=(nbody_read,) ,byteorder=self.byteorder,pio=self.pio,npart=npart_read) ########################################## # read mass if needed ########################################## if nzero != 0: if mpi.mpi_IsMaster(): print "reading mass..." massnzero = io.ReadDataBlock(f,float32,shape=(nzero,),byteorder=self.byteorder,pio=self.pio,npart=npart_m_read) #mass = concatenate((mass,massnzero)) #if nzero==nbody: # this is maybe needed when pio='yes' Tue Feb 1 21:56:48 CET 2011 if nzero==nbody_tot: mass = massnzero else: mass = array([]) i1 = 0 for i in range(len(npart)): if npart[i]!=0: # particles belong to the class if massarr[i] != 0: mass = concatenate((mass,ones(npart[i])*massarr[i])) else: i2 = i1+npart[i] mass = concatenate((mass,massnzero[i1:i2])) i1 = i2 if i2!=len(massnzero): raise "i2=",i2,"""!=len(massnzero)""" if len(mass)!=nbody_tot: raise "len(mass)=",len(mass),"!=nbody_tot" ''' if massarr[i] == 0: if npart[i]!=0: if len(massnzero)!=npart[i]: raise "this case is not taken into account, sorry !" mass = concatenate((mass,massnzero)) else: mass = concatenate((mass,ones(npart[i])*massarr[i])) ''' # extentions u = None rho = None rsp = None - opt = None + opt1= None opt2= None erd = None dte = None pot = None tstar = None minit = None idp = None metals = None if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading u..." u = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) u = concatenate((u,zeros(nbody-ngas).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading rho..." rho = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) rho = concatenate((rho,zeros(nbody-ngas).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading rsp..." rsp = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) rsp = concatenate((rsp,zeros(nbody-ngas).astype(float32))) # here it is the end of the minimal output if flag_metals: # gas metal properties if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading metals..." metals = io.ReadDataBlock(f,float32,shape=(ngas_read,NELEMENTS),byteorder=self.byteorder,pio=self.pio,npart=None) metals = concatenate((metals,zeros((nbody-ngas,NELEMENTS)).astype(float32))) if flag_age: # stellar properties if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading tstar..." tstar = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None) tstar = concatenate((-1*ones(ngas).astype(float32),tstar,-1*ones(nbody-ngas-nstars).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading minit..." minit = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None) minit = concatenate((0*ones(ngas).astype(float32),minit,0*ones(nbody-ngas-nstars).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading idp..." idp = io.ReadDataBlock(f,int32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None) idp = concatenate((-1*ones(ngas).astype(float32),idp,-1*ones(nbody-ngas-nstars).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading rho_stars..." rho_stars = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None) rho = concatenate((rho[:ngas],rho_stars,zeros(nbody-ngas-nstars).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): if mpi.mpi_IsMaster(): print "reading rsp_stars..." rsp_stars = io.ReadDataBlock(f,float32,shape=(nstars_read,),byteorder=self.byteorder,pio=self.pio,npart=None) rsp = concatenate((rsp[:ngas],rsp_stars,zeros(nbody-ngas-nstars).astype(float32))) if flag_metals: # stars metal properties if not io.end_of_file(f,pio=self.pio,MPI=MPI): + if mpi.mpi_IsMaster(): print "reading metals_stars..." metals_stars = io.ReadDataBlock(f,float32,shape=(nstars_read,NELEMENTS),byteorder=self.byteorder,pio=self.pio,npart=None) metals = concatenate(( metals[:ngas,:], metals_stars, zeros((nbody-ngas-nstars,NELEMENTS)).astype(float32) )) - + # other variables if not io.end_of_file(f,pio=self.pio,MPI=MPI): - opt = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) - if (len(opt)==ngas): - opt = concatenate((opt,zeros(nbody-ngas).astype(float32))) + if mpi.mpi_IsMaster(): print "reading opt1..." + opt1 = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) + if (len(opt1)==ngas): + opt1 = concatenate((opt1,zeros(nbody-ngas).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): + if mpi.mpi_IsMaster(): print "reading opt2..." opt2 = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) if (len(opt2)==ngas): opt2 = concatenate((opt2,zeros(nbody-ngas).astype(float32))) + + if not io.end_of_file(f,pio=self.pio,MPI=MPI): erd = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) if (len(erd)==ngas): erd = concatenate((erd,zeros(nbody-ngas).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): dte = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) dte = concatenate((dte,zeros(nbody-ngas).astype(float32))) if not io.end_of_file(f,pio=self.pio,MPI=MPI): pot = io.ReadDataBlock(f,float32,shape=(ngas_read,),byteorder=self.byteorder,pio=self.pio,npart=None) # make global self.npart = npart self.massarr = massarr self.atime = atime self.redshift = redshift self.flag_sfr = flag_sfr self.flag_feedback = flag_feedback self.nall = nall self.flag_cooling = flag_cooling self.num_files = num_files self.boxsize = boxsize self.omega0 = omega0 self.omegalambda = omegalambda self.hubbleparam = hubbleparam self.flag_age = flag_age self.flag_metals = flag_metals self.nallhw = nallhw self.flag_entr_ic = flag_entr_ic self.flag_chimie_extraheader = flag_chimie_extraheader self.critical_energy_spec = critical_energy_spec self.empty = empty self.nbody = nbody self.pos = pos self.vel = vel self.mass = mass self.num = num self.tpe = array([],int32) for i in range(len(npart)): self.tpe = concatenate( (self.tpe,ones(npart[i])*i) ) self.nzero = nzero self.u = u self.rho = rho self.tstar = tstar self.minit = minit self.idp = idp self.metals = metals - self.opt = opt - self.opt2= opt2 + self.opt1 = opt1 + self.opt2 = opt2 self.erd = erd self.dte = dte self.pot = pot self.rsp = rsp if type(self.massarr) == ndarray: self.massarr = self.massarr.tolist() if type(self.nall) == ndarray: self.nall = self.nall.tolist() if type(self.nallhw) == ndarray: self.nallhw = self.nallhw.tolist() def write_particles(self,f): ''' specific format for particle file ''' # here, we must let the user decide if we creates # the mass block or not, event if all particles have the same mass massarr,nzero = self.get_massarr_and_nzero() if self.pio == 'yes': ''' here, we have to compute also mass for each proc ''' npart = self.npart nall = self.npart_tot num_files = mpi.NTask npart_write = None npart_m_write = None else: npart = self.npart nall = self.npart_tot npart_all = array(mpi.mpi_allgather(npart)) num_files = 1 npart_write = self.npart npart_m_write = array(self.npart) * (array(self.massarr)==0) # compute the global massarr and global nzero nzero_tot = mpi.mpi_sum(nzero) massarr_all = array(mpi.mpi_allgather(massarr)) massarr_tot = zeros(len(npart),float) for i in range(len(npart)): # keep only values where there are particles massarr_all_red = compress(npart_all[:,i]!=0,massarr_all[:,i]) if len(massarr_all_red)>0: if (massarr_all_red == massarr_all_red).all(): massarr_tot[i] = massarr_all[0,i] else: # not equal raise "this case is not implemented" massarr_tot[i] = 0.0 nzero_tot = nzero_tot + sum(npart_write[:,i]) # now, re-compute nzero for the current node massarr = massarr_tot nzero = 0 for i in range(len(npart)): if massarr[i] == 0: nzero = nzero+npart[i] # now that we have the right massarr and nzero, # we can compute massnzero for each node nzero_all = zeros((mpi.NTask,len(self.npart))) if nzero != 0: ni = 0 massnzero=array([],float32) for i in range(len(self.npart)): if npart[i]!=0 and massarr[i]==0.: massnzero = concatenate((massnzero,self.mass[ni:ni+self.npart[i]])) nzero_all[mpi.ThisTask,i] = nzero_all[mpi.ThisTask,i] + npart[i] ni = ni + self.npart[i] nzero_all = mpi.mpi_allreduce(nzero_all) if self.pio == 'yes': if nzero!=0 and len(massnzero)==0: # !!! because zere is a bug see warning in get_massarr_and_nzero nzero = 0 else: npart = self.npart_tot nzero = mpi.mpi_allreduce(nzero) # to ensure that all nodes # will do -> write mass if needed # header if sys.byteorder == self.byteorder: npart = array(npart,int32).tostring() massarr = array(massarr,float).tostring() else: npart = array(npart,int32).byteswap().tostring() massarr = array(massarr,float).byteswap().tostring() atime = self.atime redshift = self.redshift flag_sfr = self.flag_sfr flag_feedback = self.flag_feedback if sys.byteorder == self.byteorder: nall = array(nall,int32).tostring() else: nall = array(nall,int32).byteswap().tostring() flag_cooling = self.flag_cooling num_files = num_files boxsize = self.boxsize omega0 = self.omega0 omegalambda = self.omegalambda hubbleparam = self.hubbleparam flag_age = self.flag_age flag_metals = self.flag_metals if sys.byteorder == self.byteorder: nallhw = array(self.nallhw,float).tostring() else: nallhw = array(self.nallhw,float).byteswap().tostring() flag_entr_ic = self.flag_entr_ic flag_chimie_extraheader = self.flag_chimie_extraheader critical_energy_spec = self.critical_energy_spec empty = self.empty # header tpl = ((npart,24),(massarr,48),(atime,float),(redshift,float),(flag_sfr,int32),(flag_feedback,int32),(nall,24),(flag_cooling,int32),(num_files,int32),(boxsize,float),(omega0,float),(omegalambda,float),(hubbleparam,float),(flag_age,int32),(flag_metals,int32),(nallhw,24),(flag_entr_ic,int32),(flag_chimie_extraheader,int32),(critical_energy_spec,float),(empty,48)) io.WriteBlock(f,tpl,byteorder=self.byteorder) # extra header if self.flag_chimie_extraheader: print "writing chimie extra-header..." - SolarAbundances = zeros(self.ChimieNelements,float32) + SolarMassAbundances = zeros(self.ChimieNelements,float32) labels = "" for i,elt in enumerate(self.ChimieElements): - SolarAbundances[i] = self.ChimieSolarAbundances[elt] + SolarMassAbundances[i] = self.ChimieSolarMassAbundances[elt] labels = labels + "%s,"%elt labels_len = (256-4-self.ChimieNelements*4) labels = labels + (labels_len-len(labels)) *" " - tpl = ((self.ChimieNelements,int32), (SolarAbundances,float32), (labels,len(labels))) + tpl = ((self.ChimieNelements,int32), (SolarMassAbundances,float32), (labels,len(labels))) io.WriteBlock(f,tpl,byteorder=self.byteorder) # positions io.WriteArray(f,self.pos.astype(float32),byteorder=self.byteorder,pio=self.pio,npart=npart_write) # velocities io.WriteArray(f,self.vel.astype(float32),byteorder=self.byteorder,pio=self.pio,npart=npart_write) # id io.WriteArray(f,self.num.astype(int32),byteorder=self.byteorder,pio=self.pio,npart=npart_write) - print npart_write # write mass if needed if nzero != 0: io.WriteArray(f,massnzero.astype(float32),byteorder=self.byteorder,pio=self.pio,npart=npart_m_write) # write extension if self.has_array('u'): if self.u != None: io.WriteArray(f,self.u[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write u" if self.has_array('rho'): if self.rho != None: io.WriteArray(f,self.rho[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write rho" if self.has_array('rsp'): if self.rsp != None: io.WriteArray(f,self.rsp[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write rsp" # this is the end of the minimal output if self.flag_metals: if self.has_array('metals'): io.WriteArray(f,self.metals[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write metals" if self.flag_age: if self.has_array('tstar'): io.WriteArray(f,self.tstar[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write tstar" if self.has_array('minit'): io.WriteArray(f,self.minit[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write minit" if self.has_array('idp'): print io.WriteArray(f,self.idp[self.npart[0]:self.npart[0]+self.npart[1]].astype(int32),byteorder=self.byteorder,pio=self.pio,npart=None) print "write idp" if self.has_array('rho'): data = self.rho[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32) if len(data)>0: io.WriteArray(f,data,byteorder=self.byteorder,pio=self.pio,npart=None) print "write rho (stars)" if self.has_array('rsp'): data = self.rsp[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32) if len(data)>0: io.WriteArray(f,data,byteorder=self.byteorder,pio=self.pio,npart=None) print "write rsp (stars)" if self.flag_metals: if self.has_array('metals'): data = self.metals[self.npart[0]:self.npart[0]+self.npart[1]].astype(float32) if len(data)>0: io.WriteArray(f,data,byteorder=self.byteorder,pio=self.pio,npart=None) print "write metals (stars)" if self.has_array('opt'): if self.opt != None: io.WriteArray(f,self.opt[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) if self.has_array('opt2'): if self.opt2 != None: io.WriteArray(f,self.opt2[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) if self.has_array('erd'): if self.erd != None: io.WriteArray(f,self.erd[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) if self.has_array('dte'): if self.dte != None: io.WriteArray(f,self.dte[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) if self.has_array('pot'): if self.pot != None: io.WriteArray(f,self.pot[:self.npart[0]].astype(float32),byteorder=self.byteorder,pio=self.pio,npart=None) def spec_info(self): """ Write spec info """ infolist = [] infolist.append("") #infolist.append("nzero : %s"%self.nzero) #infolist.append("npart : %s"%self.npart) #infolist.append("massarr : %s"%self.massarr) infolist.append("atime : %s"%self.atime) infolist.append("redshift : %s"%self.redshift) infolist.append("flag_sfr : %s"%self.flag_sfr) infolist.append("flag_feedback : %s"%self.flag_feedback) infolist.append("nall : %s"%self.nall) infolist.append("flag_cooling : %s"%self.flag_cooling) infolist.append("num_files : %s"%self.num_files) infolist.append("boxsize : %s"%self.boxsize) infolist.append("omega0 : %s"%self.omega0) infolist.append("omegalambda : %s"%self.omegalambda) infolist.append("hubbleparam : %s"%self.hubbleparam) infolist.append("flag_age : %s"%self.flag_age) infolist.append("flag_metals : %s"%self.flag_metals) infolist.append("nallhw : %s"%self.nallhw) infolist.append("flag_entr_ic : %s"%self.flag_entr_ic) infolist.append("critical_energy_spec: %s"%self.critical_energy_spec) infolist.append("") if self.has_array('u'): infolist.append("len u : %s"%len(self.u)) infolist.append("u[0] : %s"%self.u[0]) infolist.append("u[-1] : %s"%self.u[-1]) if self.has_array('rho'): infolist.append("len rho : %s"%len(self.rho)) infolist.append("rho[0] : %s"%self.rho[0]) infolist.append("rho[-1] : %s"%self.rho[-1]) if self.has_array('rsp'): infolist.append("len rsp : %s"%len(self.rsp)) infolist.append("rsp[0] : %s"%self.rsp[0]) infolist.append("rsp[-1] : %s"%self.rsp[-1]) if self.has_array('opt'): infolist.append("len opt : %s"%len(self.opt)) infolist.append("opt[0] : %s"%self.opt[0]) infolist.append("opt[-1] : %s"%self.opt[-1]) if self.has_array('opt2'): infolist.append("len opt2 : %s"%len(self.opt2)) infolist.append("opt2[0] : %s"%self.opt2[0]) infolist.append("opt2[-1] : %s"%self.opt2[-1]) if self.has_array('erd'): infolist.append("len erd : %s"%len(self.erd)) infolist.append("erd[0] : %s"%self.erd[0]) infolist.append("erd[-1] : %s"%self.erd[-1]) if self.has_array('dte'): infolist.append("len dte : %s"%len(self.dte)) infolist.append("dte[0] : %s"%self.dte[0]) infolist.append("dte[-1] : %s"%self.dte[-1]) if self.has_array('tstar'): infolist.append("len tstar : %s"%len(self.tstar)) infolist.append("tstar[0] : %s"%self.tstar[0]) infolist.append("tstar[-1] : %s"%self.tstar[-1]) if self.has_array('idp'): infolist.append("len idp : %s"%len(self.idp)) infolist.append("idp[0] : %s"%self.idp[0]) infolist.append("idp[-1] : %s"%self.idp[-1]) return infolist #def select(self,tpe='gas'): def select(self,*arg,**kw): """ Return an N-body object that contain only particles of a certain type, defined by in gadget: gas : gas particles halo : halo particles disk : disk particles bulge : bulge particles stars : stars particles bndry : bndry particles sph : gas with u > u_c sticky : gas with u < u_c """ index = {'gas':0,'halo':1,'disk':2,'bulge':3,'stars':4,'bndry':5,'stars1':1,'halo1':2} # this allows to write nb.select(('gas','disk')) if len(arg)==1: if type(arg[0])==types.TupleType: arg = arg[0] tpes = arg # create the selection vector c = zeros(self.nbody) for tpe in tpes: if type(tpe) == types.StringType: if (tpe=='sph'): c = c+(self.u>self.critical_energy_spec)*(self.tpe==0) elif (tpe=='sticky'): c = c+(self.u0)*(self.tpe==0) elif (tpe=='cg'): c = c+(self.u<0)*(self.tpe==0) elif (tpe=='all'): return self elif not index.has_key(tpe): print "unknown type, do nothing %s"%(tpe) return self else: i = index[tpe] c = c+(self.tpe==i) elif type(tpe) == types.IntType: c = c+(self.tpe==tpe) return self.selectc(c) ''' elif type(tpe) == types.StringType: if (tpe=='sph'): nb = self.select('gas') return nb.selectc((self.u>self.critical_energy_spec)) if (tpe=='sticky'): nb = self.select('gas') return nb.selectc((nb.u<=nb.critical_energy_spec)) if not index.has_key(tpe): print "unknown type %s"%(tpe) return self i = index[tpe] else: i = tpe if self.npart[i]==0: #print "no particle of type %s"%(tpe) return self.selectc(zeros(self.nbody)) n1 = sum(self.npart[:i]) n2 = n1 + self.npart[i]-1 return self.sub(n1,n2) ''' def subdis(self,mode='dd',val=None): """ Equivalent of select """ return self.select(mode) def Z(self): """ total metallicity """ - return log10(self.metals[:,self.NELEMENTS-1] / self.ChimieSolarAbundances['Metals'] + 1.0e-20) + return log10(self.metals[:,self.NELEMENTS-1] / self.ChimieSolarMassAbundances['Metals'] + 1.0e-20) def Fe(self): """ metallicity Fe """ - return log10(self.metals[:,0] / self.ChimieSolarAbundances['Fe'] + 1.0e-20) + return log10(self.metals[:,0] / self.ChimieSolarMassAbundances['Fe'] + 1.0e-20) def Mg(self): """ magnesium """ - return log10(self.metals[:,1] / self.ChimieSolarAbundances['Mg'] + 1.0e-20) + return log10(self.metals[:,1] / self.ChimieSolarMassAbundances['Mg'] + 1.0e-20) def O(self): """ Oxygen """ - return log10(self.metals[:,2] / self.ChimieSolarAbundances['O'] + 1.0e-20) + return log10(self.metals[:,2] / self.ChimieSolarMassAbundances['O'] + 1.0e-20) def Ba(self): """ Barium """ - return log10(self.metals[:,3] / self.ChimieSolarAbundances['Ba'] + 1.0e-20) + return log10(self.metals[:,3] / self.ChimieSolarMassAbundances['Ba'] + 1.0e-20) def MgFe(self): eps = 1e-20 - MgFe = log10((self.metals[:,1]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarAbundances['Mg'] * self.ChimieSolarAbundances['Fe']) + MgFe = log10((self.metals[:,1]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarMassAbundances['Mg'] * self.ChimieSolarMassAbundances['Fe']) return MgFe def BaFe(self): eps = 1e-20 - BaFe = log10((self.metals[:,3]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarAbundances['Ba'] * self.ChimieSolarAbundances['Fe']) + BaFe = log10((self.metals[:,3]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarMassAbundances['Ba'] * self.ChimieSolarMassAbundances['Fe']) return BaFe def SiFe(self): eps = 1e-20 - SiFe = log10((self.metals[:,3]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarAbundances['Si'] * self.ChimieSolarAbundances['Fe']) + SiFe = log10((self.metals[:,3]+eps)/(self.metals[:,0]+eps) / self.ChimieSolarMassAbundances['Si'] * self.ChimieSolarMassAbundances['Fe']) return SiFe def luminosity_spec(self,tnow=None): """ compute specific luminosity, using metalicity u_mass = 1.e10 LuminosityElement*u_mass/1.0d6, ' x 10^6 Lsun' """ # initialize SSP from pNbody.SSP import libvazdekis # vazdekis_kb_mu1.3.txt : krupa 01 revisited self.LObj = libvazdekis.VazdekisLuminosities(os.path.join(OPTDIR,'SSP','vazdekis_kb_mu1.3.txt')) self.LObj.ExtrapolateMatrix(order=1,s=0) self.LObj.CreateInterpolator() self.LObj.Extrapolate2DMatrix() u_mass = 1.e10 u_time = 4.7287e6 SolarAbun_Fe = 0.001771 if self.tstar==None: return array([],float32) if tnow==None: tnow = self.atime Ages = (tnow-self.tstar)*u_time*1.0e-9 # ages in Gyr Zs = self.Z() # compute luminosities using LObj L = self.LObj.Luminosities(Zs,Ages) return L def old_luminosity_spec(self,tnow=None): """ compute specific luminosity, using metalicity u_mass = 1.e10 LuminosityElement*u_mass/1.0d6, ' x 10^6 Lsun' """ u_mass = 1.e10 u_time = 4.7287e6 SolarAbun_Fe = 0.001771 if self.tstar==None: return 0 if tnow==None: tnow = self.atime age = (tnow-self.tstar)*u_time*1.0e-9 Fe = log10(self.metals[:,0] / SolarAbun_Fe + 1.0e-20) met_bins = array([-1.68, -1.28,-0.68]) age_bins = array([0.007,0.02,0.07,0.10,0.11,0.13,0.14,0.16,0.18,0.20,0.22,0.25,0.28, 0.32,0.35,0.40,0.45,0.50,0.56,0.63,0.71,0.79,0.89,1.00,1.12,1.26,1.41,1.58,1.78,2.00, 2.24,2.51,2.82,3.16,3.55,3.98,4.47,5.01,5.62,6.31,7.08,7.94,8.91,10.00,11.22,12.59,14.13,15.85]) L0 = array([1./(0.003*2.0),1./(0.007*2.0),1./(0.054*2.0),1./0.172,1./0.122,1./0.126,1./0.132,1./0.138,1./0.145, 1./0.153,1./0.166,1./0.179,1./0.196,1./0.218,1./0.232,1./0.248,1./0.267,1./0.292,1./0.313,1./0.339, 1./0.371,1./0.397,1./0.429,1./0.464,1./0.471,1./0.461,1./0.597,1./0.640,1./0.681,1./0.746,1./0.813, 1./0.899,1./0.975,1./1.085,1./1.202,1./1.328,1./1.464,1./1.606,1./1.753,1./1.929,1./2.100,1./2.282, 1./2.485,1./2.762,1./2.987,1./3.230,1./3.504,1./3.794,1./4.100]) L1 = array([1./(0.003 * 2.0),1./(0.007 * 2.0),1./(0.054 * 2.0),1./0.154,1./0.126,1./0.134,1./0.141,1./0.149, 1./0.162,1./0.177,1./0.195,1./0.212,1./0.227,1./0.245,1./0.263,1./0.278,1./0.296,1./0.324,1./0.348, 1./0.376,1./0.412,1./0.444,1./0.481,1./0.516,1./0.550,1./0.543,1./0.590,1./0.688,1./0.747,1./0.819, 1./0.896,1./0.985,1./1.077,1./1.197,1./1.317,1./1.444,1./1.573,1./1.735,1./1.888,1./2.066,1./2.304, 1./2.502,1./2.729,1./3.008,1./3.256,1./3.564,1./3.740,1./4.212,1./4.491]) L2 = array([1./(0.003 * 2.0),1./(0.007 * 2.0),1./(0.054 * 2.0),1./0.163,1./0.169,1./0.159,1./0.172,1./0.184, 1./0.198,1./0.214,1./0.230,1./0.253,1./0.270,1./0.291,1./0.318,1./0.345,1./0.371,1./0.401,1./0.432, 1./0.464,1./0.501,1./0.540,1./0.583,1./0.633,1./0.684,1./0.710,1./0.694,1./0.870,1./0.989,1./1.092, 1./1.183,1./1.296,1./1.410,1./1.534,1./1.665,1./1.809,1./2.023,1./2.200,1./2.477,1./2.692,1./2.927, 1./3.214,1./3.490,1./3.783,1./4.087,1./4.411,1./4.742,1./5.115,1./5.515]) L3 = array([1./(0.003 * 2.0),1./(0.007 * 2.0),1./(0.054 * 2.0),1./0.166,1./0.177,1./0.179,1./0.192,1./0.206, 1./0.221,1./0.241,1./0.259,1./0.282,1./0.308,1./0.331,1./0.357,1./0.386,1./0.415,1./0.448,1./0.485, 1./0.524,1./0.562,1./0.614,1./0.671,1./0.725,1./0.785,1./0.770,1./0.937,1./1.068,1./1.189,1./1.318, 1./1.448,1./1.575,1./1.736,1./1.918,1./2.066,1./2.236,1./2.439,1./2.728,1./2.977,1./3.331,1./3.722, 1./4.063,1./4.387,1./4.717,1./5.047,1./5.474,1./5.972,1./6.636,1./7.018]) L = array([L0,L1,L2,L3]) #i = searchsorted(met_bins,Fe) # metallicity index #j = searchsorted(age_bins,age) # age index i = met_bins.searchsorted(Fe) j = age_bins.searchsorted(age) #return take(L,[i,j],[0,1])*self.mass return L[i,j] def luminosity(self,tnow=None): return self.luminosity_spec(tnow)*self.mass def Age(self): '''in Gyrs (for treeasph units)''' u_time = 4.7287e6 age = (self.atime-self.tstar)*u_time*1.0e-9 return age + def age(self): + '''in CU ''' + return (self.atime-self.tstar) + ################################################################# - # physical values + # physical values (with correct unit conversion) ################################################################# + def Rxyz(self,a=None,h=None,units=None): + """ + return the radius of each particles in physical units, i.e. + correct it from the scaling factor and h if necessary (i.e. comoving integration is on) + """ + + print "... compute Rxyz()" + + + # set factor unit + funit=1.0 + if units!=None: + funit = self.localsystem_of_units.convertionFactorTo(units) + print "... factor units = %g"%funit + + + if self.isComovingIntegrationOn(): + print " converting to physical units (a=%5.3f h=%5.3f)"%(self.atime,self.hubbleparam) + return self.rxyz()*self.atime/self.hubbleparam * funit + else: + return self.rxyz() * funit + + + def Rxy(self,a=None,h=None,units=None): + """ + return the radius of each particles in physical units, i.e. + correct it from the scaling factor and h if necessary (i.e. comoving integration is on) + """ + + print "... compute Rxy()" + + + # set factor unit + funit=1.0 + if units!=None: + funit = self.localsystem_of_units.convertionFactorTo(units) + print "... factor units = %g"%funit + + + if self.isComovingIntegrationOn(): + print " converting to physical units (a=%5.3f h=%5.3f)"%(self.atime,self.hubbleparam) + return self.rxy()*self.atime/self.hubbleparam * funit + else: + return self.rxy() * funit + + + + def Rho(self,a=None,h=None,units=None): ''' return the density of particles. a : scaling factor h : hubble parameter units : output units different cases : comoving integration (self.comovingintegration==True) 1) convert into physical coorinates 2) if a=1 -> stay in comoving (in this case, we can also use nb.rho) non comoving integration (self.comovingintegration==False) 1) do not convert 2) if I want to force a behavior : put a=0.1 -> ''' - print "..compute Rho()" - - # !!!!!!!!!!! - # !!!!!!!!!!! - self.hubbleparam=0.704 # this must be set correctely - # !!!!!!!!!!! - # !!!!!!!!!!! + print "... compute Rho()" + # set factor unit funit=1.0 if units!=None: funit = self.localsystem_of_units.convertionFactorTo(units) - - print self.rho.mean() - print self.rho.mean()*funit - print self.rho.mean()/(self.atime**3)*self.hubbleparam**2 *funit,self.atime**3,self.hubbleparam**2 + print "... factor units = %g"%funit - if self.comovingintegration: - print "Rho : converting to physical units (a=%5.3f h=%5.3f)"%(self.atime,self.hubbleparam) - print "mean:",log10((self.rho/self.atime**3*self.hubbleparam**2 *funit).mean()) + if self.isComovingIntegrationOn(): + print " converting to physical units (a=%5.3f h=%5.3f)"%(self.atime,self.hubbleparam) return self.rho/self.atime**3*self.hubbleparam**2 *funit else: - print "Rho : converting to physical units (h=%5.3f)"%(self.hubbleparam) - return self.rho*self.hubbleparam**2 *funit - - - - - - + return self.rho*funit def T(self): ''' u does not depends on a nor h ''' - print "----------------------------------> temperature...." + print "... compute T()" gamma = self.unitsparameters.get('gamma') xi = self.unitsparameters.get('xi') ionisation = self.unitsparameters.get('ionisation') mu = thermodyn.MeanWeight(xi,ionisation) mh = ctes.PROTONMASS.into(self.localsystem_of_units) k = ctes.BOLTZMANN.into(self.localsystem_of_units) thermopars = {"k":k,"mh":mh,"mu":mu,"gamma":gamma} - T = where((self.u>0),thermodyn.Tru(self.Rho(),self.u,thermopars),0) - return T - - + T = where((self.u>0),thermodyn.Tru(None,self.u,thermopars),0) + return T - - def Tcool(self): + + def Tcool(self,units=None): from pNbody import cooling + #print "... compute Tcool()" + if self.metals == None: FeH = zeros(self.nbody).astype(float32) else: FeH = self.metals[:,0] #l = cooling.get_lambda_from_Density_EnergyInt_FeH(self.rho,self.u,FeH) #dudt = l/self.rho #tcool = self.u/dudt + # parameters for the cooling + + from pNbody import cooling + cooling_params = self.localsystem_of_units.getparam() + cooling_params['CoolingFile'] = "/home/epfl/revaz/.pNbody/cooling_with_metals.dat" + cooling.init_cooling(cooling_params) + tcool = cooling.get_cooling_time_from_Density_EnergyInt_FeH(self.rho,self.u,FeH) + + # set factor unit + funit=1.0 + if units!=None: + funit = self.localsystem_of_units.convertionFactorTo(units) + print "... factor units = %g"%funit + tcool = tcool*funit + + if self.isComovingIntegrationOn(): + print "Tcool : isComovingIntegrationOn not implemented" + sys.exit() return tcool.astype(float32) + def StellarAge(self,units=None): + ''' + stellar age + ''' + print "... compute StellarAge()" + + age = (self.atime-self.tstar) + + + # set factor unit + funit=1.0 + if units!=None: + funit = self.localsystem_of_units.convertionFactorTo(units) + print "... factor units = %g"%funit + + if self.isComovingIntegrationOn(): + print "with comoving interation on, you need to convert da in cosmic time" + print "this is not implemented" + sys.exit() + else: + return age*funit + + + + + + + + + + + + + def dLdt(self): + from pNbody import cooling + + if self.metals == None: + FeH = zeros(self.nbody).astype(float32) + else: + FeH = self.metals[:,0] + + l = cooling.get_lambda_from_Density_EnergyInt_FeH(self.rho,self.u,FeH) + dLdt = self.mass * l/self.rho + + return dLdt.astype(float32) + + + def CosmicTime(self): + """ + return cosmic time in Gyrs + """ + from pNbody import cosmo + from pNbody import ctes + + + Hubble = ctes.HUBBLE.into(self.localsystem_of_units) + pars = {"Hubble":Hubble,"HubbleParam":self.hubbleparam,"OmegaLambda":self.omegalambda,"Omega0":self.omega0} + + print "WARNING : THIS FUNCTION MUST BE TESTED !!!" + + return cosmo.CosmicTime_a(self.atime,pars) + + + + def Redshift(self): + """ + return redshift + """ + from pNbody import cosmo + return cosmo.Z_a(self.atime) + + + + def sfr(self,dt): + """ + star formation rate per particle + + all units are in code units + """ + + sfr = where( (self.atime-self.tstar) < dt, self.mass/dt ,0 ) + + return sfr + + def toPhysicalUnits(self,a=None,h=None): + """ + convert from comobile units to physical units + correct from the scaling factor and + from the hubble parameter + """ + + if self.isComovingIntegrationOn(): + + if a==None: + a = self.atime + if h==None: + h = self.hubbleparam + + print " converting to physical units (a=%5.3f h=%5.3f)"%(a,h) + Hubble = ctes.HUBBLE.into(self.localsystem_of_units) + OmegaLambda= self.omegalambda + Omega0 = self.omega0 + print " (HubbleCte =%5.3f)"%Hubble + print " (OmegaLambda=%5.3f)"%OmegaLambda + print " (Omega0 =%5.3f)"%Omega0 + + + pars = {"Hubble":Hubble,"OmegaLambda":OmegaLambda,"Omega0":Omega0} + + Ha = cosmo.Hubble_a(a,pars=pars) + + self.vel = self.pos*Ha*a + self.vel*sqrt(a) + self.pos = self.pos*a/h + self.mass= self.mass/h + + if self.has_array('u'): + self.u = self.u + if self.has_array('rho'): + self.rho = self.rho/a**2 * h**2 + + def TimeStepLevel(self): + """ + return the timestep level in log2 + """ + + return (log10(self.opt1)/log10(2)).astype(int) + diff --git a/examples/ic/addmetals/addmetals.py b/examples/ic/addmetals/addmetals.py index 3094880..6a8831c 100755 --- a/examples/ic/addmetals/addmetals.py +++ b/examples/ic/addmetals/addmetals.py @@ -1,206 +1,206 @@ #!/usr/bin/env python from pNbody import * from optparse import OptionParser from PyChem import chemistry ######################################## # # parser # ######################################## def parse_options(): usage = "usage: %prog [options] file" parser = OptionParser(usage=usage) parser.add_option("-t", action="store", dest="ftype", type="string", default = 'gadget', help="type of the file", metavar=" TYPE") parser.add_option("-o", action="store", dest="outputfile", type="string", default = None, help="outputfile name", metavar=" STING") parser.add_option("--exec", action="store", dest="execline", type="string", default = None, help="give command to execute before", metavar=" STRING") parser.add_option("-p", action="store", dest="chimieparameterfile", type="string", default = "chimie.dat", help="chimie parameter file", metavar=" STING") parser.add_option("--NumberOfTables", action="store", dest="NumberOfTables", type="int", default = 1, help="NumberOfTables", metavar=" INT") parser.add_option("--DefaultTable", action="store", dest="DefaultTable", type="int", default = 0, help="DefaultTable", metavar=" INT") (options, args) = parser.parse_args() if len(args) == 0: print "you must specify a filename" sys.exit(0) file = None else: file = args return file,options ################################################################################ # # MAIN # ################################################################################ file,opt = parse_options() ######################### # open nbody file nb = Nbody(file,ftype=opt.ftype) if opt.execline != None: exec(opt.execline) # set internal energy system_of_units = units.Set_SystemUnits_From_File('params') out_units = units.UnitSystem('local',[units.Unit_kpc,units.Unit_Msol,units.Unit_Gyr,units.Unit_K]) ToGyr = system_of_units.convertionFactorTo(out_units.UnitTime) FromGyr = 1/ToGyr #################################### # convert disk and bulge to stars1 nb0 = nb.select(0) nb1 = nb.select(1) nb2 = nb.select(2) nb3 = nb.select(3) nb4 = nb.select(4) nb5 = nb.select(5) nb3.set_tpe(1) nb4.set_tpe(1) nb5.set_tpe(1) nb = nb0 + nb3 + nb4 + nb5 + nb2 del nb0,nb1,nb2,nb3,nb4,nb5 nb.init() ######################### # open chimie file chemistry.init_chimie(opt.chimieparameterfile,opt.NumberOfTables,opt.DefaultTable) ######################### # set the chimie extra-header nb.flag_chimie_extraheader = 1 nb.ChimieNelements = int(chemistry.get_nelts()) nb.ChimieElements = chemistry.get_elts_labels() -nb.ChimieSolarAbundances = chemistry.get_elts_SolarAbundances() +nb.ChimieSolarMassAbundances = chemistry.get_elts_SolarMassAbundances() NELEMENTS = nb.ChimieNelements nb.flag_metals=NELEMENTS ######################### # global init nb.rho = ones(nb.nbody).astype(float32) nb.rsp = ones(nb.nbody).astype(float32) # must be!=0 with CHIMIE_INPUT_ALL nb.metals = zeros((nb.nbody,NELEMENTS)).astype(float32) ######################### # variables for the gas # Fe Fe = -1 + nb.rxy()*(-1/10.) Fe = Fe.astype(float32) -metals = nb.ChimieSolarAbundances['Fe'] * (10**Fe - 1.0e-20) +metals = nb.ChimieSolarMassAbundances['Fe'] * (10**Fe - 1.0e-20) nb.metals[:,0] = where(nb.tpe==0, metals,nb.metals[:,0] ) ######################### # variables for the stars nb.flag_age=1 # age (formation time) nb.tstar = zeros(nb.nbody).astype(float32) tstar = clip( random.normal(-3,1,nb.nbody) ,-5,0) * FromGyr tstar = tstar.astype(float32) nb.tstar = where(nb.tpe==1, tstar,nb.tstar ) # minit (initial mass) nb.minit = zeros(nb.nbody).astype(float32) minit = nb.mass # !!! minit = minit.astype(float32) nb.minit = where(nb.tpe==1, minit,nb.minit ) # idp (stellar index) nb.idp = zeros(nb.nbody).astype(float32) idp = -1 * ones(nb.nbody) idp = idp.astype(int32) nb.idp = where(nb.tpe==1, idp,nb.idp ) # Fe Fe = -1.5 + nb.rxy()*(-1/10.) Fe = Fe.astype(float32) -metals = nb.ChimieSolarAbundances['Fe'] * (10**Fe - 1.0e-20) +metals = nb.ChimieSolarMassAbundances['Fe'] * (10**Fe - 1.0e-20) nb.metals[:,0] = where(nb.tpe==1, metals,nb.metals[:,0] ) ######################### # write outputs if opt.outputfile!=None: nb.rename(opt.outputfile) nb.write() diff --git a/examples/movies/movie08/Mkgmov.py b/examples/movies/movie08/Mkgmov.py index 020357e..3b48977 100644 --- a/examples/movies/movie08/Mkgmov.py +++ b/examples/movies/movie08/Mkgmov.py @@ -1,530 +1,530 @@ import os from pNbody import * from pNbody.param import Params import copy import types import gzip ####################################################################### # some usefull functions ####################################################################### def ReadNbodyParameters(paramname): ''' read param from a parameter Nbody file ''' gparams = Params(paramname,None) param = {} # create new params for p in gparams.params: param[p[0]]=p[3] return param def gzip_compress(file): f = open(file,'r') content = f.read() f.close() f = gzip.open(file+'.gz', 'wb') f.write(content) f.close() ####################################################################### # # C L A S S D E F I N I T I O N # ####################################################################### class Movie(): def __init__(self,parameterfile='filmparam.py',format=None,imdir=None,timesteps=None,pio=False,compress=True): self.DEFAULT_TIMESTEPS = None self.DEFAULT_FORMAT = "fits" self.DEFAULT_IMDIR = "fits" self.DEFAULT_PIO = False self.DEFAULT_COMPRESS = True self.DEFAULT_SCALE = "log" self.DEFAULT_CD = 0. self.DEFAULT_MN = 0. self.DEFAULT_MX = 0. self.DEFAULT_PALETTE = "light" self.parameterfile=parameterfile # read the parameter file self.read_parameterfile() # use options if format != None: self.film['format'] = format if imdir != None: self.film['imdir'] = imdir if timesteps != None: self.film['timesteps'] = timesteps if pio != None: self.film['pio'] = pio if compress != None: self.film['compress'] = compress self.imdir = self.film['imdir'] self.pio = self.film['pio'] self.compress = self.film['compress'] # deal with timesteps self.set_timesteps() self.ifile=0 if mpi.mpi_IsMaster(): if self.pio: self.pio = "yes" else: self.pio = "no" if self.parameterfile==None: print "you must specify a parameter file" sys.exit() if not os.path.exists(self.parameterfile): print "file %s does not exists"%self.parameterfile sys.exit() if not os.path.exists(self.imdir): os.mkdir(self.imdir) def info(self): print "INFO INFO INFO" #print self.film print self.parameterfile print self.getftype() print "INFO INFO INFO" def read_parameterfile(self): if not os.path.isfile(self.parameterfile): raise IOError(915,'file %s not found ! Pease check the file name.'%(self.parameterfile)) filmparam = __import__(os.path.splitext(self.parameterfile)[0],globals(), locals(), [], -1) self.film = filmparam.film # set some defaults if not self.film.has_key('timesteps'): self.film['timesteps'] = self.DEFAULT_TIMESTEPS if not self.film.has_key('imdir'): self.film['imdir'] = self.DEFAULT_IMDIR if not self.film.has_key('format'): self.film['format'] = self.DEFAULT_FORMAT if not self.film.has_key('pio'): self.film['pio'] = self.DEFAULT_PIO if not self.film.has_key('compress'): self.film['compress'] = self.DEFAULT_COMPRESS self.setftype(self.film['ftype']) # post process for i,frame in enumerate(self.film['frames']): frame['id'] = i # check #for frame in self.film['frames']: # print frame['id'] # for component in frame['components']: # print " ",component['id'] ################################## # time steps stuffs ################################## def set_timesteps(self): """ define self.times (which is a list) based on the value contained in self.film['timesteps'] """ # self.times if self.film['timesteps']=='every': self.times="every" elif type(self.film['timesteps']) == types.StringType: fname = self.film['timesteps'] if not os.path.isfile(fname): raise IOError(916,'file %s not found ! Pease check the file name.'%(fname)) times = io.read_ascii(fname,[0])[0] times = take( times,len(times)-1-arange(len(times))) # invert order times = times.tolist() self.times = times elif type(self.film['timesteps']) == types.ListType: self.times = self.film['timesteps'] elif type(self.film['timesteps']) == types.TupleType: t0 = self.film['timesteps'][0] t1 = self.film['timesteps'][1] dt = self.film['timesteps'][2] times = arange(t0,t1,dt) times = take( times,len(times)-1-arange(len(times))) # invert order times = times.tolist() self.times = times else: self.times=[] def set_next_time(self): if self.times!="every": if len(self.times)>0: self.times.pop() def get_next_time(self): if self.times=="every": return 0.0 if len(self.times)==0: return None else: return self.times[-1] def getftype(self): return self.ftype def setftype(self,ftype): self.ftype = ftype def AppplyFilmParam(self,nb,film): # set time reference for this file exec("nb.tnow = %s"%film['time']) # exec1 if film['exec'] != None: exec(film['exec']) # macro if film['macro'] != None: execfile(film['macro']) return nb def AppplyFrameParam(self,nb,frame): nbf = nb # exec if frame['exec'] != None: exec(frame['exec']) # macro if frame['macro'] != None: execfile(frame['macro']) return nbf def AppplyComponentParam(self,nbf,component): if component['id'][0] == '@': # here, all tasks must have an object containing all particles # ok, but not in the right order !!! nbfc = Nbody(componentid,self.getftype()) nbfc = nbfc.SendAllToAll() nbfc.sort() nbfc.componentid=component['id'][1:] elif component['id'][0] == '#': nbfc = nbf nbfc.componentid = component['id'] else: nbfc = nbf.select(component['id']) nbfc.componentid = component['id'] # exec if component['exec'] != None: exec(component['exec']) # macro if component['macro'] != None: execfile(component['macro']) #print "------------------------" #print min(nbfc.u),max(nbfc.u) #print min(nbfc.rho),max(nbfc.rho) #print min(nbfc.tpe),max(nbfc.tpe) #print "temperature",min(nbfc.T()),max(nbfc.T()) #print nbfc.nbody #print min(nbfc.rsp),max(nbfc.rsp) #print "------------------------" return nbfc def dump(self,dict): # exctract dict atime = dict['atime'] pos = dict['pos'] # create nbody object nb = Nbody(pos=pos,ftype='gadget') nb.atime = atime # add other arrays if dict.has_key('header'): header = dict['header'] for key in header.keys(): setattr(nb, key, header[key] ) if nb.flag_chimie_extraheader: - if not dict.has_key('ChimieSolarAbundances'): - raise KeyError(200,'key ChimieSolarAbundances not found in dict.') + if not dict.has_key('ChimieSolarMassAbundances'): + raise KeyError(200,'key ChimieSolarMassAbundances not found in dict.') if not dict.has_key('nelts'): raise KeyError(201,'key nelts not found in dict.') - nb.ChimieSolarAbundances = dict['ChimieSolarAbundances'] + nb.ChimieSolarMassAbundances = dict['ChimieSolarMassAbundances'] nb.ChimieNelements = dict['nelts'] # add other arrays if dict.has_key('vel'): nb.vel = dict['vel'] if dict.has_key('num'): nb.num = dict['num'] if dict.has_key('mass'): nb.mass = dict['mass'] if dict.has_key('tpe'): nb.tpe = dict['tpe'] if dict.has_key('u'): nb.u = dict['u'] if dict.has_key('rho'): nb.rho = dict['rho'] if dict.has_key('rsp'): nb.rsp = dict['rsp'] if dict.has_key('metals'): nb.metals = dict['metals'] nb.init() #print "################################" #print "writing qq.dat" #print "################################" #nb.rename('qq.dat') #nb.write() self.dumpimage(nb=nb) def dumpimage(self,nb=None,file=None): if file!=None: #print "(task=%04d) reading "%(mpi.ThisTask),file,self.getftype(),self.pio nb = Nbody(file,ftype=self.getftype(),pio=self.pio) else: if nb==None: raise "you must specify at least a file or give an nbody object" film = self.film nb = self.AppplyFilmParam(nb,film) for frame in film['frames']: nbf = self.AppplyFrameParam(nb,frame) for component in frame['components']: nbfc = self.AppplyComponentParam(nbf,component) # find the observer position # 1) from params # 2) from pfile # 3) from tdir # and transform into parameter if frame['tdir']!=None: tfiles = glob.glob(os.path.join(frame['tdir'],"*")) tfiles.sort() bname = os.path.basename(file) tfiles_for_this_file = [] for j in xrange(len(tfiles)): tfile = "%s.%05d"%(os.path.basename(file),j) tmp_tfile = os.path.join(frame['tdir'],tfile) if os.path.exists(tmp_tfile): tfiles_for_this_file.append(tmp_tfile) elif frame['pfile']!=None: if not os.path.isfile(frame['pfile']): print "parameter file %s does not exists(1)..."%(frame['pfile']) # read from pfile defined in frame param = ReadNbodyParameters(frame['pfile']) tfiles_for_this_file = [None] else: # take frame as parameter param = copy.copy(frame) tfiles_for_this_file = [None] # loop over different oberver positions for this file for j,tfile in enumerate(tfiles_for_this_file): if tfile!=None: param = ReadNbodyParameters(tfile) # add parameters defined by user in the parameter file for key in component.keys(): param[key] = component[key] # set image shape using frame param['shape'] = (frame['width'],frame['height']) # compute map mat = nbfc.CombiMap(param) if mpi.mpi_IsMaster(): # save output if self.film["format"]=="fits": output = '%04d_%04d-%s-%06d.fits'%(self.ifile,frame['id'],nbfc.componentid,j) output = os.path.join(self.imdir,output) print nb.atime,output if os.path.exists(output): os.remove(output) header = [('TIME',nb.tnow,'snapshot time')] io.WriteFits(transpose(mat), output, extraHeader = header) # compress if self.compress: gzip_compress(output) os.remove(output) elif self.film["format"]=="png": output = '%04d_%04d-%s-%06d.png'%(self.ifile,frame['id'],nbfc.componentid,j) output = os.path.join(self.imdir,output) print nb.atime,output if not frame.has_key('scale'): frame['scale']=self.DEFAULT_SCALE if not frame.has_key('cd'): frame['cd']=self.DEFAULT_CD if not frame.has_key('mn'): frame['mn']=self.DEFAULT_MN if not frame.has_key('mx'): frame['mx']=self.DEFAULT_MX if not frame.has_key('palette'): frame['palette']=self.DEFAULT_PALETTE matint,mn_opt,mx_opt,cd_opt = set_ranges(mat,scale=frame['scale'],cd=frame['cd'],mn=frame['mn'],mx=frame['mx']) frame['mn'] = mn_opt frame['mx'] = mx_opt frame['cd'] = cd_opt img = get_image(matint,palette_name=frame['palette']) img.save(output) print frame['mn'],frame['mx'],frame['cd'] # increment counter self.ifile+=1 diff --git a/pNbody/Mkgmov.py b/pNbody/Mkgmov.py index 217c26e..797cb57 100644 --- a/pNbody/Mkgmov.py +++ b/pNbody/Mkgmov.py @@ -1,809 +1,809 @@ import Mtools import Mtools as mt import os import glob from pNbody import * from pNbody.param import Params import copy import types import gzip ######################################################################################################## # # Some info on the structure : # # 1) nb : the full nbody object (main loop) # # 2) nbf : normally set the display parameters, position of the observer # but could be used to set the component # # 3) nbfc : normally set the component to display and/or the value to dispplay # # component['id'][0] == '@' # this is used to open simple nbody object (to draw a box, for example) # component['id'][0] == '#' # no selection # component['id'] == 'gas' # selection gas # # # # # a script may be called at level 2) # # frame['ext_cmd'] = [] # frame['ext_cmd'].append("""from plot_P import MakePlot""") # frame['ext_cmd'].append("""MakePlot([nbf],output)""") # # it must save an output output, the latter is opened as an img # and appened to imgs # # # a script may be called at level 3) # # 1) # component['ext_cmd'] = [] # component['ext_cmd'].append("""from plot_P import MakePlot""") # component['ext_cmd'].append("""MakePlot([nbfc],output)""") # # it must save an output output, the latter is opened as an img # and appened to imgs # # # 2) using # component['to_img']='script_name' # # ######################################################################################################## ####################################################################### # some usefull functions ####################################################################### def ReadNbodyParameters(paramname): ''' read param from a parameter Nbody file ''' gparams = Params(paramname,None) param = {} # create new params for p in gparams.params: param[p[0]]=p[3] return param def gzip_compress(file): f = open(file,'r') content = f.read() f.close() f = gzip.open(file+'.gz', 'wb') f.write(content) f.close() ####################################################################### # # C L A S S D E F I N I T I O N # ####################################################################### class Movie(): def __init__(self,parameterfile='filmparam.py',format=None,imdir=None,timesteps=None,pio=False,compress=True): self.DEFAULT_TIMESTEPS = None self.DEFAULT_FORMAT = "fits" self.DEFAULT_IMDIR = "fits" self.DEFAULT_PIO = False self.DEFAULT_COMPRESS = True self.DEFAULT_SCALE = "log" self.DEFAULT_CD = 0. self.DEFAULT_MN = 0. self.DEFAULT_MX = 0. self.DEFAULT_PALETTE = "light" self.parameterfile=parameterfile # read the parameter file self.read_parameterfile() # use options if format != None: self.film['format'] = format if imdir != None: self.film['imdir'] = imdir if timesteps != None: self.film['timesteps'] = timesteps if pio != None: self.film['pio'] = pio if compress != None: self.film['compress'] = compress self.imdir = self.film['imdir'] self.pio = self.film['pio'] self.compress = self.film['compress'] # deal with timesteps self.set_timesteps() self.ifile=-1 self.file_offset=0 if mpi.mpi_IsMaster(): # whith gmkgmov, the init is only runned by the master, this line is not needed... if self.pio: self.pio = "yes" else: self.pio = "no" if self.parameterfile==None: print "you must specify a parameter file" sys.exit() if not os.path.exists(self.parameterfile): print "file %s does not exists"%self.parameterfile sys.exit() if not os.path.exists(self.imdir): os.mkdir(self.imdir) else: print "directory %s exists !!!"%self.imdir files = os.listdir(self.imdir) print "the directory contains %d files"%(len(files)) # png files png_files = glob.glob(os.path.join(self.imdir,"*.png")) n_png_files = len(png_files) print "the directory contains %d png files"%(n_png_files) # fits files fits_files = glob.glob(os.path.join(self.imdir,"*.fits")) n_fits_files = len(fits_files) print "the directory contains %d fits files"%(n_fits_files) if n_png_files > 0: self.file_offset = n_png_files def info(self): print "INFO INFO INFO" #print self.film print self.parameterfile print self.getftype() print "INFO INFO INFO" def read_parameterfile(self): if not os.path.isfile(self.parameterfile): raise IOError(915,'file %s not found ! Pease check the file name.'%(self.parameterfile)) # import the parameter file as a module module_name = os.path.basename(os.path.splitext(self.parameterfile)[0]) module_dir = os.path.dirname(self.parameterfile) if sys.path.count(module_dir) == 0: sys.path.append(module_dir) filmparam = __import__(module_name,globals(), locals(), [], -1) self.film = filmparam.film # set some defaults if not self.film.has_key('timesteps'): self.film['timesteps'] = self.DEFAULT_TIMESTEPS if not self.film.has_key('imdir'): self.film['imdir'] = self.DEFAULT_IMDIR if not self.film.has_key('format'): self.film['format'] = self.DEFAULT_FORMAT if not self.film.has_key('pio'): self.film['pio'] = self.DEFAULT_PIO if not self.film.has_key('compress'): self.film['compress'] = self.DEFAULT_COMPRESS self.setftype(self.film['ftype']) # post process for i,frame in enumerate(self.film['frames']): frame['id'] = i # check #for frame in self.film['frames']: # print frame['id'] # for component in frame['components']: # print " ",component['id'] ################################## # time steps stuffs ################################## def set_timesteps(self): """ define self.times (which is a list) based on the value contained in self.film['timesteps'] """ # self.times if self.film['timesteps']=='every': self.times="every" elif type(self.film['timesteps']) == types.StringType: fname = self.film['timesteps'] if not os.path.isfile(fname): raise IOError(916,'file %s not found ! Pease check the file name.'%(fname)) times = io.read_ascii(fname,[0])[0] times = take( times,len(times)-1-arange(len(times))) # invert order times = times.tolist() self.times = times elif type(self.film['timesteps']) == types.ListType: self.times = self.film['timesteps'] elif type(self.film['timesteps']) == types.TupleType: t0 = self.film['timesteps'][0] t1 = self.film['timesteps'][1] dt = self.film['timesteps'][2] times = arange(t0,t1,dt) times = take( times,len(times)-1-arange(len(times))) # invert order times = times.tolist() self.times = times else: self.times=[] def set_next_time(self): if self.times!="every": if len(self.times)>0: self.times.pop() def get_next_time(self): if self.times=="every": return 0.0 if len(self.times)==0: return None else: return self.times[-1] def getftype(self): return self.ftype def setftype(self,ftype): self.ftype = ftype def ApplyFilmParam(self,nb,film): # set time reference for this file exec("nb.tnow = %s"%film['time']) # exec1 if film.has_key('exec'): if film['exec'] != None: exec(film['exec']) # macro if film.has_key('macro'): if film['macro'] != None: execfile(film['macro']) return nb def ApplyFrameParam(self,nb,frame): nbf = nb # exec if frame.has_key('exec'): if frame['exec'] != None: exec(frame['exec']) # macro if frame.has_key('macro'): if frame['macro'] != None: execfile(frame['macro']) return nbf def ApplyComponentParam(self,nbf,component): if component['id'][0] == '@': # here, all tasks must have an object containing all particles # ok, but not in the right order !!! nbfc = Nbody(componentid,self.getftype()) nbfc = nbfc.SendAllToAll() nbfc.sort() nbfc.componentid=component['id']#[1:] elif component['id'][0] == '#': nbfc = nbf nbfc.componentid = component['id']#[1:] else: nbfc = nbf.select(component['id']) nbfc.componentid = component['id'] # exec if component.has_key('exec'): if component['exec'] != None: exec(component['exec']) # macro if component.has_key('macro'): if component['macro'] != None: execfile(component['macro']) #print "------------------------" #print min(nbfc.u),max(nbfc.u) #print min(nbfc.rho),max(nbfc.rho) #print min(nbfc.tpe),max(nbfc.tpe) #print "temperature",min(nbfc.T()),max(nbfc.T()) #print nbfc.nbody #print min(nbfc.rsp),max(nbfc.rsp) #print "------------------------" return nbfc def dump(self,dict): # exctract dict atime = dict['atime'] pos = dict['pos'] # create nbody object nb = Nbody(pos=pos,ftype='gadget') nb.atime = atime # add other arrays if dict.has_key('vel'): nb.vel = dict['vel'] if dict.has_key('num'): nb.num = dict['num'] if dict.has_key('mass'): nb.mass = dict['mass'] if dict.has_key('tpe'): nb.tpe = dict['tpe'] if dict.has_key('u'): nb.u = dict['u'] if dict.has_key('rho'): nb.rho = dict['rho'] if dict.has_key('rsp'): nb.rsp = dict['rsp'] if dict.has_key('metals'): nb.metals = dict['metals'] #!!! nb.flag_chimie_extraheader=1 nb.ChimieNelements = 5 nb.flag_metals = 5 - nb.ChimieSolarAbundances={} - nb.ChimieSolarAbundances['Fe']=0.00176604 + nb.ChimieSolarMassAbundances={} + nb.ChimieSolarMassAbundances['Fe']=0.00176604 nb.init() #print "################################" #print "writing qq.dat" #print "################################" #nb.rename('qq.dat') #nb.write() self.dumpimage(nb=nb) def dumpimage(self,nb=None,file=None): # increment counter self.ifile+=1 # skip file if needed if self.ifile0: ################################################# # 1) use an outer script to create an img (this is a bit redundant with 2.2, see below ) ################################################# output= "/tmp/%015d.png"%(int(random.random()*1e17)) for cmd in frame['ext_cmd']: exec(cmd) if mpi.mpi_IsMaster(): img = Image.open(output) imgs.append(img) if os.path.exists(output): os.remove(output) continue # composition parameters if frame.has_key('cargs'): if len(frame['cargs'])!=0: frame['compose']=True datas = [] else: frame['compose']=False else: frame['cargs']=[] frame['compose']=False for component in frame['components']: if mpi.mpi_IsMaster(): print "------------------------" print "component",component['id'] print "------------------------" nbfc = self.ApplyComponentParam(nbf,component) # find the observer position # 1) from params # 2) from pfile # 3) from tdir # and transform into parameter if frame['tdir']!=None: tfiles = glob.glob(os.path.join(frame['tdir'],"*")) tfiles.sort() bname = os.path.basename(file) tfiles_for_this_file = [] for j in xrange(len(tfiles)): tfile = "%s.%05d"%(os.path.basename(file),j) tmp_tfile = os.path.join(frame['tdir'],tfile) if os.path.exists(tmp_tfile): tfiles_for_this_file.append(tmp_tfile) elif frame['pfile']!=None: if not os.path.isfile(frame['pfile']): print "parameter file %s does not exists(1)..."%(frame['pfile']) # read from pfile defined in frame param = ReadNbodyParameters(frame['pfile']) tfiles_for_this_file = [None] else: # take frame as parameter param = copy.copy(frame) tfiles_for_this_file = [None] # loop over different oberver positions for this file for iobs,tfile in enumerate(tfiles_for_this_file): if tfile!=None: param = ReadNbodyParameters(tfile) # add parameters defined by user in the parameter file for key in component.keys(): param[key] = component[key] # set image shape using frame param['shape'] = (frame['width'],frame['height']) # compute map mat = nbfc.CombiMap(param) if mpi.mpi_IsMaster(): if frame['compose']: datas.append(mat) if component.has_key('ext_cmd'): ################################################# # 1) use an outer script to create an img ################################################# if len(component['ext_cmd'])>0: output= "/tmp/%015d.png"%(int(random.random()*1e17)) for cmd in component['ext_cmd']: exec(cmd) if mpi.mpi_IsMaster(): img = Image.open(output) imgs.append(img) if os.path.exists(output): os.remove(output) elif self.film["format"]=="fits": ################################# # 1) save fits file ################################# output = '%04d_%04d-%s-%06d.fits'%(self.ifile,frame['id'],component['id'],iobs) output = os.path.join(self.imdir,output) print nb.atime,output if os.path.exists(output): os.remove(output) header = [('TIME',nb.tnow,'snapshot time')] io.WriteFits(transpose(mat), output, extraHeader = header) # compress if self.compress: gzip_compress(output) os.remove(output) elif self.film["format"]=="png": ################################# # 2) output png file or ... ################################# output = '%04d_%04d-%s-%06d.png'%(self.ifile,frame['id'],nbfc.componentid,iobs) output = os.path.join(self.imdir,output) print nb.atime,output # here, we should use component['scale'] ... not frame['scale'], no ? if not frame.has_key('scale'): frame['scale']=self.DEFAULT_SCALE if not frame.has_key('cd'): frame['cd']=self.DEFAULT_CD if not frame.has_key('mn'): frame['mn']=self.DEFAULT_MN if not frame.has_key('mx'): frame['mx']=self.DEFAULT_MX if not frame.has_key('palette'): frame['palette']=self.DEFAULT_PALETTE matint,mn_opt,mx_opt,cd_opt = set_ranges(mat,scale=frame['scale'],cd=frame['cd'],mn=frame['mn'],mx=frame['mx']) frame['mn'] = mn_opt frame['mx'] = mx_opt frame['cd'] = cd_opt img = get_image(matint,palette_name=frame['palette']) img.save(output) print frame['mn'],frame['mx'],frame['cd'] # need to create an img if component.has_key('to_img'): if component['to_img']==True: ########################################## # 2.1) create an img and apply commmands ########################################## # get params if not component.has_key('scale'): component['scale']=self.DEFAULT_SCALE if not component.has_key('cd'): component['cd']=self.DEFAULT_CD if not component.has_key('mn'): component['mn']=self.DEFAULT_MN if not component.has_key('mx'): component['mx']=self.DEFAULT_MX if not component.has_key('palette'): component['palette']=self.DEFAULT_PALETTE matint,mn_opt,mx_opt,cd_opt = set_ranges(mat,scale=component['scale'],cd=component['cd'],mn=component['mn'],mx=component['mx']) img = get_image(matint,palette_name=component['palette']) print mn_opt,mx_opt,cd_opt # here we can add img commands.... if component.has_key('img_cmd'): if len(component['img_cmd'])>0: for cmd in component['img_cmd']: exec(cmd) # append img to list img.atime = nb.atime imgs.append(img) elif type(component['to_img'])==types.StringType: ########################################## # 2.2) use an outer script to create an img from mat ########################################## output= "/tmp/%015d.png"%(int(random.random()*1e17)) # get params if not component.has_key('scale'): component['scale']=self.DEFAULT_SCALE if not component.has_key('cd'): component['cd']=self.DEFAULT_CD if not component.has_key('mn'): component['mn']=self.DEFAULT_MN if not component.has_key('mx'): component['mx']=self.DEFAULT_MX if not component.has_key('palette'): component['palette']=self.DEFAULT_PALETTE component['atime']=nbfc.atime mk = __import__(component['to_img'],globals(), locals(), [], -1) mk.MkImage(mat,output,component) img = Image.open(output) imgs.append(img) os.remove(output) del nbf ####################### # compose components ####################### if frame['compose']: if mpi.mpi_IsMaster(): img,cargs = Mtools.fits_compose_colors_img(datas,frame['cargs']) # save #output = '%04d_%04d.png'%(self.ifile,frame['id']) #output = os.path.join(self.imdir,output) #img.save(output) # append img to list img.atime = nb.atime imgs.append(img) del nb ####################### # compose frames ####################### if mpi.mpi_IsMaster(): if film.has_key('img_cmd'): if len(film['img_cmd'])>0: for cmd in film['img_cmd']: exec(cmd) output = '%04d.png'%(self.ifile) output = os.path.join(self.imdir,output) img.save(output) img.i=0