import numpy as np
class Galaxy(object):
[docs] def __init__(self,
cmx, cmy, cmz,
cmvx, cmvy, cmvz,
gas_mass, stellar_mass,
HMR, FMR):
self.cm = np.array([cmx,cmy,cmz])
self.vel = np.array([cmvx,cmvy,cmvz])
self.gas_mass = gas_mass
self.stellar_mass = stellar_mass
self.HMR = HMR
self.FMR = FMR
self.central = 1
def deleteAttribute(self,attrib):
[docs] if hasattr(self, attrib):
delattr(self,attrib)
def getProgenGalaxy(self,makeCopy = 1):
[docs] if makeCopy:
import copy
progenGal = copy.deepcopy(self)
progenGal.deleteAttribute('glist')
progenGal.deleteAttribute('slist')
if hasattr(progenGal, 'halo'):
progenGal.halo.deleteAttribute('dmlist')
progenGal.halo.deleteAttribute('glist')
progenGal.halo.deleteAttribute('slist')
progenGal.halo.galaxies = []
progenGal.halo.sub_halos = []
return progenGal
else:
self.deleteAttribute('glist')
self.deleteAttribute('slist')
if hasattr(self, 'halo'):
self.halo.deleteAttribute('dmlist')
self.halo.deleteAttribute('glist')
self.halo.deleteAttribute('slist')
self.halo.galaxies = []
self.halo.sub_halos = []
return self
def rotator(self,obj, DM = 0):
[docs] import rotator as rotator
from scipy.spatial import cKDTree
if not hasattr(obj.particleData, 'gpos'):
obj.particleData.getParticleData()
pd = obj.particleData
gpos = pd.gpos[self.glist]
spos = pd.spos[self.slist]
pos = np.concatenate((gpos,spos))
pos[:,0] -= self.cm[0]
pos[:,1] -= self.cm[1]
pos[:,2] -= self.cm[2]
gvel = pd.gvel[self.glist]
svel = pd.svel[self.slist]
vel = np.concatenate((gvel,svel))
gmass = pd.gmass[self.glist]
smass = pd.smass[self.slist]
mass = np.concatenate((gmass,smass))
pos, d1, d2, ex, ey, ez, newvel = rotator.rotategal(pos,vel,mass)
## project into the xy plane
r = np.sqrt(pos[:,0]**2 + pos[:,1]**2)
vphi = (newvel[:,0] * -1. * pos[:,1] + newvel[:,1] * pos[:,0]) / r
vr = (newvel[:,0] * pos[:,0] + newvel[:,1] * pos[:,1]) / r
self.r_gas_xy = r[0:len(self.glist)]
self.vphi_gas = vphi[0:len(self.glist)]
self.vr_gas = vr[0:len(self.glist)]
self.r_star_xy = r[len(self.glist)::]
self.vphi_star = vphi[len(self.glist)::]
self.vr_star = vr[len(self.glist)::]
self.rotated_gpos = pos[0:len(self.glist)] + self.cm
self.rotated_spos = pos[len(self.glist)::] + self.cm
if DM:
if not hasattr(pd, 'dmTREE'):
print 'constructing DM particle KDTree'
pd.dmTREE = cKDTree(pd.dmpos)
dmindexes = pd.dmTREE.query_ball_point(self.cm, np.max((self.FMR,20)))
pos = pd.dmpos[dmindexes]
pos[:,0] -= self.cm[0]
pos[:,1] -= self.cm[1]
pos[:,2] -= self.cm[2]
vel = pd.dmvel[dmindexes]
newvel = np.zeros((len(vel),3))
for i in range(0,len(pos)):
pos[i] = rotator.rotate_stuff(pos[i,0],pos[i,1],pos[i,2],
self.ALPHA,self.BETA,0,returnval=1)
vx,vy,vz = rotator.rotate_stuff(vel[i,0],vel[i,1],vel[i,2],
self.ALPHA,self.BETA,0)
newvel[i,0] = vx
newvel[i,1] = vy
newvel[i,2] = vz
## project into the xy plane
r = np.sqrt(pos[:,0]**2 + pos[:,1]**2)
vphi = (newvel[:,0] * -1. * pos[:,1] + newvel[:,1] * pos[:,0]) / r
vr = (newvel[:,0] * pos[:,0] + newvel[:,1] * pos[:,1]) / r
self.r_dm_xy = r
self.vphi_dm = vphi
self.vr_dm = vr
self.rotated_dmpos = pos
def calculate_HMR(self,obj,mode):
[docs] pd = obj.particleData
cm = self.cm
if mode == 'gas' or mode == 'H2':
if len(self.glist) < 4:
self.gas_HMR = 0.
self.gas_FMR = 0.
self.H2_HMR = 0.
return
pos = pd.gpos[self.glist]
mass = pd.gmass[self.glist]
if mode == 'H2':
mass *= 0.76 * pd.gfH2[self.glist]
elif mode == 'stellar':
pos = pd.spos[self.slist]
mass = pd.smass[self.slist]
else:
obj.logOutput('wrong mode for calculate_HMR!!')
return
total_mass = np.sum(mass)
half_mass = 0.5 * total_mass
## no H2
if total_mass == 0:
self.H2_HMR = 0.
return
particle_radii = np.sqrt( (pos[:,0] - cm[0])**2 +
(pos[:,1] - cm[1])**2 +
(pos[:,2] - cm[2])**2 )
s_indexes = np.argsort(particle_radii)
s_radii = particle_radii[s_indexes]
s_mass = mass[s_indexes]
m = 0.
r = 0.
for i in range(0,len(s_radii)):
m += s_mass[i]
r = s_radii[i]
if m >= half_mass:
break
if mode == 'gas':
self.gas_HMR = r
self.gas_FMR = s_radii[-1]
elif mode == 'stellar':
self.stellar_HMR = r
elif mode == 'H2':
if total_mass == 0:
self.H2_HMR = 0.
else:
self.H2_HMR = r
return r
class Color(object):
[docs] def __init__(self,name,wavelength_range,index):
self.name = name
self.index = index
self.wavelength_range = wavelength_range
#self.mag = np.zeros(3)
#self.mag.fill(np.nan)
#self.emag = np.zeros(3)
#self.emag.fill(np.nan)
self.mag = {}
self.emag = {}
@property
def flux(self):
tmp = np.zeros(3)
for i in range(0,3):
if self.mag[i] == np.nan:
tmp[i] = 0
else:
tmp[i] = self.mag2flux(self.mag[i])
return tmp
@property
def eflux(self):
tmp = np.zeros(3)
for i in range(0,3):
if self.emag[i] == np.nan:
tmp[i] = 0
else:
tmp[i] = self.mag2flux(self.emag[i])
return tmp
def mag2flux(self,MAG):
[docs] return 10**(MAG / -2.5)
def flux2mag(self,FLUX):
[docs] return -2.5 * np.log10(FLUX)