diff --git a/pNbody/units.py b/pNbody/units.py index fc33107..cd05cba 100644 --- a/pNbody/units.py +++ b/pNbody/units.py @@ -1,672 +1,711 @@ import copy from pNbody import io,param ################################################################################ class Units: ################################################################################ ''' Units ''' def __init__(self,symbol,factor=1.,power=1,ulist=[]): ''' ''' self.symbol = symbol self.factor = float(factor) self.power = int(power) self.ulist = ulist self.GatherBaseUnits() def set_symbol(self,symbol): self.ulist = [copy.deepcopy(self)] self.power = 1 self.factor = 1.0 self.symbol = symbol self.GatherBaseUnits() def GatherBaseUnits(self): ''' create the self.bases list (recursively) ''' self.bases = [] for unit in self.ulist: if len(unit.ulist)==0: # a base unit self.bases.append(unit) # add the unit else: # go deeper in the tree unit.GatherBaseUnits() self.bases = self.bases + unit.bases self.ReduceBaseUnits() def ReduceBaseUnits(self): allunits = {} factorTot = 1.0 for base in self.bases: factorTot = factorTot * base.factor if allunits.has_key(base.symbol): allunits[base.symbol] = allunits[base.symbol] + base.power else: allunits[base.symbol] = base.power self.bases = [] flag = 1 for key in allunits.keys(): factor = 1.0 if flag: factor = factorTot flag = 0 power = allunits[key] symbol = key ulist = [] self.bases.append(Units(symbol,factor,power,ulist)) def get_basefactor(self): ''' return the factor relative to the base units actually, it is stored as the factor of the first base unit ''' return bases[0].factor #def to_units(self,newunits): # ''' # return a new unit object in new units # ''' # # # check that the units are similar def str_base(self): string = '' if self.factor != 1.0: string = string + self.factor.__str__() if len(string) > 0: string = string + ' ' if self.symbol == None: string = string + '-' else: string = string + self.symbol if self.power != 1: string = string + '^%d'%self.power return string def __str__(self): string = '' # (1) display using symbol if self.factor != 0: string = string+self.str_base() str1 = "[%s]"%string string = '' # (2) display using (base units) for base in self.bases: if len(string) > 0: string = string + ' ' string = string + base.str_base() str2 = "[%s]"%string if str2 != '[]': str1 = "%s\t%s"%(str1,str2) return str1 def __mul__(self,y): x = copy.deepcopy(self) y = copy.deepcopy(y) # multiply with a scalar if type(y)== int or type(y) == float : x.factor = x.factor*y x.mul_list(x.ulist,y) x.GatherBaseUnits() return x # multiply with a unit elif type(y)==type(x): factor = x.factor * y.factor if x.symbol == y.symbol: # same units x.symbol = x.symbol x.power = x.power + y.power x.ulist = x.ulist + y.ulist x.GatherBaseUnits() return x else: # different unit symbol = None power = 1 factor = 1.0 if x.ulist == []: # case of a base unit x.ulist = [copy.deepcopy(x)] if y.ulist == []: # case of a base unit y.ulist = [copy.deepcopy(y)] ulist = x.ulist + y.ulist return Units(symbol,factor,power,ulist) def __div__(self,y): x = copy.deepcopy(self) y = copy.deepcopy(y) # divide with a scalar if type(y) == int or type(y) == float: x.factor = x.factor/y x.div_list(x.ulist,y) x.GatherBaseUnits() return x # divide with a unit elif type(y)==type(x): factor = x.factor / y.factor if x.symbol == y.symbol: # same units x.symbol = None x.power = 0 x.ulist = [] x.GatherBaseUnits() return x else: # different unit symbol = None power = 1 factor = 1.0 if x.ulist == []: # case of a base unit x.ulist = [copy.deepcopy(x)] if y.ulist == []: # case of a base unit y.ulist = [copy.deepcopy(y)] y.pow_list(y.ulist,-1) y.GatherBaseUnits() ulist = x.ulist + y.ulist return Units(symbol,factor,power,ulist) def __rdiv__(self,y): x = copy.deepcopy(self) # divide with a scalar if type(y)==int or type(y)==float: return y*x**-1 else: raise "TypeError","unsupported operand type(s) for / or div()" def __pow__(self,y): x = copy.deepcopy(self) y = copy.deepcopy(y) # power with a scalar if type(y)==int or type(y)==float: x.factor = x.factor**y x.power = x.power*y x.pow_list(x.ulist,y) x.GatherBaseUnits() return x else: raise "TypeError","unsupported operand type(s) for ** or pow()" def mul_list(self,ulist,y): for u in ulist: if len(u.ulist)==0: # a base unit u.factor = u.factor*y else: # go deeper in the tree u.mul_list(u.ulist,y) u.GatherBaseUnits() def div_list(self,ulist,y): for u in ulist: if len(u.ulist)==0: # a base unit u.factor = u.factor/y else: # go deeper in the tree u.mul_list(u.ulist,y) u.GatherBaseUnits() def pow_list(self,ulist,y): for u in ulist: if len(u.ulist)==0: # a base unit u.factor = u.factor**y u.power = u.power*y else: # go deeper in the tree u.pow_list(u.ulist,y) u.GatherBaseUnits() __rmul__ = __mul__ ################################################################################ class UnitSystem: ################################################################################ ''' Units system ''' def __init__(self,UnitSysName,UnitLst): ''' UnitDic = {'lenght':UnitLenght,'mass':UnitMass,'time':UnitTime} The units used to define a system of units must be "basic", i.e., not based on several units - ''' + usefull variables + + self.dic_of_factors : value of base unit in current units + self.UnitDic : contains base units and value of the choosen unit in this base units + ''' + UnitDic = {} for u in UnitLst: u.GatherBaseUnits() if len(u.bases) > 1: raise "UnitSystemError","%s is not a pure Unit."%(u.symbol) elif len(u.bases) == 1: UnitDic[u.bases[0].symbol] = u.bases[0] elif len(u.bases) == 0: # base unit UnitDic[u.symbol] = u - dic_of_factors = {} - dic_of_powers = {} + + + self.set_dics(UnitDic) - for UnitType in UnitDic.keys(): - u = UnitDic[UnitType] - - if u.bases ==[]: # base unit - dic_of_factors[UnitType] = 1/u.factor - dic_of_powers[UnitType] = u.power - else: - dic_of_factors[UnitType] = 1/u.bases[0].factor - dic_of_powers[UnitType] = u.bases[0].power + self.UnitSysName = UnitSysName + self.UnitLst = UnitLst self.UnitDic = UnitDic - self.dic_of_factors = dic_of_factors # value of base unit in current units - self.dic_of_powers = dic_of_powers # not usefull - + self.UnitLength = self.UnitDic['m'] self.UnitMass = self.UnitDic['kg'] self.UnitTime = self.UnitDic['s'] self.UnitVelocity = self.UnitDic['m']/self.UnitDic['s'] self.UnitSpecEnergy = self.UnitVelocity**2 self.UnitEnergy = self.UnitSpecEnergy*self.UnitMass self.UnitDensity = self.UnitMass/self.UnitLength**3 self.UnitSurfaceDensity = self.UnitMass/self.UnitLength**2 self.UnitSurface = self.UnitLength**2 + + def set_dics(self,UnitDic): + + dic_of_factors = {} + dic_of_powers = {} + + for UnitType in UnitDic.keys(): + u = UnitDic[UnitType] + + if u.bases ==[]: # base unit + dic_of_factors[UnitType] = 1/u.factor + dic_of_powers[UnitType] = u.power + else: + dic_of_factors[UnitType] = 1/u.bases[0].factor + dic_of_powers[UnitType] = u.bases[0].power + + self.dic_of_factors = dic_of_factors # value of base unit in current units + self.dic_of_powers = dic_of_powers # not usefull + + + + def get_UnitLength_in_cm(self): f = self.UnitDic['m'].factor c = PhysCte(f,Unit_m) return c.into(cgs) def get_UnitMass_in_g(self): f = self.UnitDic['kg'].factor c = PhysCte(f,Unit_kg) return c.into(cgs) def get_UnitTime_in_s(self): f = self.UnitDic['s'].factor c = PhysCte(f,Unit_s) return c.into(cgs) def get_UnitVelocity_in_cm_per_s(self): return self.get_UnitLength_in_cm()/self.get_UnitTime_in_s() def get_UnitEnergy_in_cgs(self): return self.get_UnitMass_in_g()*self.get_UnitVelocity_in_cm_per_s()**2 def get_UnitDensity_in_cgs(self): return self.get_UnitMass_in_g()/self.get_UnitLength_in_cm()**3 def info(self): ''' print some info ''' print "units info" print " UnitLength_in_cm =%g"%(self.get_UnitLength_in_cm()) print " UnitVelocity_in_cm_per_s =%g"%(self.get_UnitVelocity_in_cm_per_s()) print " UnitMass_in_g =%g"%(self.get_UnitMass_in_g()) def getparam(self): param = {} param["UnitLength_in_cm"] = self.get_UnitLength_in_cm() param["UnitVelocity_in_cm_per_s"] = self.get_UnitVelocity_in_cm_per_s() param["UnitMass_in_g"] = self.get_UnitMass_in_g() return param def convertionFactorTo(self,newUnits): ''' return the conversion factor to obtain the new units ''' f = 1.0 # loop over all base units of the new Units if newUnits.bases==[]: baselist = [newUnits] else: baselist = newUnits.bases for base in baselist: factor = base.factor symbol = base.symbol power = base.power # multiply by the right factor f = f /( factor * (self.dic_of_factors[symbol])**power ) return f def into(self,newUnits): ''' return into the new units ''' return PhysCte(self.convertionFactorTo(newUnits),newUnits) + def multiply_UnitLength(self,f): + self.UnitDic['m'].factor *= f + self.set_dics(self.UnitDic) + + def multiply_UnitMass(self,f): + self.UnitDic['kg'].factor *= f + self.set_dics(self.UnitDic) + + def multiply_UnitTime(self,f): + self.UnitDic['s'].factor *= f + self.set_dics(self.UnitDic) + + + + def CorrectFromHubbleParameter(self,HubbleParam): + print "Units Correction : HubbleParam = %g"%HubbleParam + self.multiply_UnitLength(1.0/HubbleParam) + self.multiply_UnitMass(1.0/HubbleParam) + self.multiply_UnitTime(1.0/HubbleParam) + ################################################################################ class PhysCte: ################################################################################ ''' Physical constant ''' def __init__(self,value,Unit): #super(PhysCte,self).__init__(value) self.value = value self.Unit = Unit def factor_to_base(self): if self.Unit.bases == []: factor = self.Unit.factor else: factor = self.Unit.bases[0].factor return factor def __str__(self): return float.__str__(float(self.value)) + ' %s'%(self.Unit) def into(self,SystemName): # here, we could return an object f = self.value if self.Unit.bases==[]: unit = self.Unit symbol = unit.symbol f = f * unit.factor * (SystemName.dic_of_factors[symbol])**unit.power else: for unit in self.Unit.bases: symbol = unit.symbol f = f * unit.factor * (SystemName.dic_of_factors[symbol])**unit.power return float(f) def __float__(self): return float(self.value) def __mul__(self,y): return PhysCte(self.value*y.value,self.Unit*y.Unit) def __add__(self,y): pass def __sub__(self,y): pass def __div__(self,y): return PhysCte(self.value/float(y.value),self.Unit/y.Unit) def __pow__(self,y): pass ################################################################################ # define some units ################################################################################ # base units Unit_m = Units('m') Unit_kg = Units('kg') Unit_s = Units('s') Unit_mol = Units('mol') Unit_C = Units('C') Unit_K = Units('K') # other lenght units Unit_cm = 0.01*Unit_m Unit_cm.set_symbol('cm') Unit_km = 1e3*Unit_m Unit_km.set_symbol('km') Unit_Mm = 1e6*Unit_m Unit_Mm.set_symbol('Mm') Unit_Gm = 1e9*Unit_m Unit_Gm.set_symbol('Gm') Unit_kpc= 1000*3.085e18*Unit_cm Unit_kpc.set_symbol('kpc') Unit_pc= 3.085e18*Unit_cm Unit_pc.set_symbol('pc') Unit_ua= 1.495978e11*Unit_m Unit_ua.set_symbol('ua') # other mass units Unit_g = 0.001*Unit_kg Unit_g.set_symbol('g') Unit_Ms = 1.9891e33*Unit_g Unit_Ms.set_symbol('Ms') Unit_Msol = 1.9891e33*Unit_g Unit_Msol.set_symbol('Ms') Unit_Mg = 2.23e11*Unit_Ms Unit_Mg.set_symbol('Mg') Unit_Mt = 5.9742e24*Unit_kg Unit_Mt.set_symbol('Mt') Unit_Mj = 317.893*Unit_Mt Unit_Mj.set_symbol('Mj') # other time units Unit_h = 3600*Unit_s Unit_h.set_symbol('h') Unit_yr = 31536000*Unit_s Unit_yr.set_symbol('yr') Unit_kyr= 1e3*Unit_yr Unit_kyr.set_symbol('kyr') Unit_Myr= 1e6*Unit_yr Unit_Myr.set_symbol('Myr') Unit_Gyr= 1e9*Unit_yr Unit_Gyr.set_symbol('Myr') Unit_dy= 86400*Unit_s Unit_dy.set_symbol('days') Unit_hr= 3600*Unit_s Unit_hr.set_symbol('hr') Unit_century= 100*Unit_yr Unit_century.set_symbol('century') # other speed units Unit_kmh= Unit_km/Unit_h Unit_kmh.set_symbol('kmh') Unit_kms= Unit_km/Unit_s Unit_kms.set_symbol('kms') # other units Unit_ms = Unit_m/Unit_s Unit_ms2 = Unit_ms**2 Unit_J = Unit_kg * Unit_ms2 Unit_J.set_symbol('J') Unit_erg = Unit_g * (Unit_cm/Unit_s)**2 Unit_erg.set_symbol('erg') Unit_N = Unit_kg * Unit_m/(Unit_s)**2 Unit_N.set_symbol('N') Unit_Pa = Unit_N/Unit_m**2 Unit_Pa.set_symbol('Pa') Unit_G = Unit_N*Unit_m**2/Unit_kg**2 Unit_d = Unit_kg/Unit_m**3 Unit_Lsol = 3.839 * Unit_erg/Unit_s Unit_Lsol.set_symbol('Lsol') ################################################################################ # define some common unit systems ################################################################################ mks = UnitSystem('mks',[Unit_m, Unit_kg, Unit_s , Unit_K, Unit_mol, Unit_C]) cgs = UnitSystem('cgs',[Unit_cm, Unit_g, Unit_s , Unit_K, Unit_mol, Unit_C]) gal = UnitSystem('gal',[Unit_kpc,Unit_Mg, Unit_Myr, Unit_K, Unit_mol, Unit_C]) ''' constantes are now defined in ctes.py ################################################## # define some const. in mks ################################################## slight = PhysCte(2.99792458e8,Unit_m/Unit_s) G = PhysCte(6.6732e-11,Unit_G) q_electron = PhysCte(1.6022e-19,Unit_C) planck = PhysCte(6.6262e-34,Unit_J*Unit_s) boltzmann = PhysCte(1.3807e-23,Unit_J/Unit_K) m_electron = PhysCte(9.1095e-31,Unit_kg) m_proton = PhysCte(1.6726e-27,Unit_kg) m_neutron = PhysCte(1.6750e-27,Unit_kg) Na = PhysCte(6.0220e+23,Unit_mol) R = PhysCte(8.3144e+00,Unit_J/Unit_mol) Av = PhysCte(6.828e-50 ,Unit_Pa*Unit_m**6) Bv = PhysCte(4.419e-29,Unit_m**3) ''' ################################# def Set_SystemUnits_From_Params(params): ################################# ''' return a system of units from given parameters params is a dictionary that must constains at least params['UnitVelocity_in_cm_per_s'] params['UnitMass_in_g'] params['UnitLength_in_cm'] ''' UnitVelocity_in_cm_per_s = params['UnitVelocity_in_cm_per_s'] UnitMass_in_g = params['UnitMass_in_g'] UnitLength_in_cm = params['UnitLength_in_cm'] UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s # now from the params, define a system of units Unit_lenght = Unit_cm * UnitLength_in_cm Unit_mass = Unit_g * UnitMass_in_g Unit_time = Unit_s * UnitTime_in_s localsystem = UnitSystem('local',[Unit_lenght, Unit_mass, Unit_time , Unit_K, Unit_mol, Unit_C]) return localsystem ################################# def Set_SystemUnits_From_File(unitsfile): ################################# ''' return a system of units from given file coding units unitsfile file is either a gadget parameter file or a pNbody units file ''' params = {} # try to read a gadget file try: gparams = io.read_params(unitsfile) params['UnitLength_in_cm'] = gparams['UnitLength_in_cm'] params['UnitMass_in_g'] = gparams['UnitMass_in_g'] params['UnitVelocity_in_cm_per_s'] = gparams['UnitVelocity_in_cm_per_s'] #params['Omega0'] = gparams['Omega0'] #params['OmegaLambda'] = gparams['OmegaLambda'] #params['OmegaBaryon'] = gparams['OmegaBaryon'] #params['BoxSize'] = gparams['BoxSize'] #params['ComovingIntegrationOn'] = gparams['ComovingIntegrationOn'] params = gparams except: # try to read a pNbody units file try: gparams = param.Params(unitsfile,None) params['UnitLength_in_cm'] = gparams.get("UnitLength_in_cm") params['UnitMass_in_g'] = gparams.get("UnitMass_in_g") params['UnitVelocity_in_cm_per_s'] = gparams.get("UnitVelocity_in_cm_per_s") #params['Omega0'] = gparams.get('Omega0') #params['OmegaLambda'] = gparams.get('OmegaLambda') #params['OmegaBaryon'] = gparams.get('OmegaBaryon') #params['BoxSize'] = gparams.get('BoxSize') #params['ComovingIntegrationOn'] = gparams.get('ComovingIntegrationOn') except: raise IOError(015,'format of unitsfile %s unknown ! Pease check.'%(unitsfile)) localsystem = Set_SystemUnits_From_Params(params) return localsystem