IMF and number of massive stars per cluster

Some experiments determining what fraction of massive stars form in clusters above a given mass, given a few assumptions about the cluster initial mas function. I'm sure this has been done in papers, but I needed to do the exercise for myself.... and I don't know which papers off the top of my head (but I should add them to this post if/when I find them).

In [1]:

from agpy import imf
figsize(12,8)

reg_gal2cel requires coords & pyregion
Region Photometry requires pyregion
cubes.py requires pyregion for getspec_reg
cubes.py requires pywcs for some subimage_integ,aper_wordl2pix,getspec, and coords_in_image

In [2]:

# From a Schecter mass function with some cutoff, generate a collection of clusters
# Then, determine whether "most" massive stars are in massive clusters or if they're evenly distributed with mass
mrange = logspace(1,5,1000) # from 10 msun to million msun clusters
schec = imf.modified_schechter(mrange,m0=1e4,m1=100) # cutoff at low mass
figure()
loglog(mrange,schec,label="Modified Schecter $m_1=10^2$, $m_0=10^4$")
loglog(mrange,imf.schechter(mrange,m0=1e5),label="Schecter $m_0=10^5$")
loglog(mrange,imf.schechter(mrange,m0=1e4),label="Schecter $m_0=10^4$")
xlabel("Cluster Mass")
ylabel("Fractional Number of Clusters")
legend(loc='best')

Out [2]:

<matplotlib.legend.Legend at 0x10e81b950>

In [3]:

nclusters = 50000
clusters = []
obcount = []
ocount = []
mclusters = []
for rn in np.random.rand(nclusters):
    mcluster = imf.inverse_imf(rn, massfunc=imf.modified_schechter, mmax=1e7, mmin=10, m1=100, m0=1e5)
    cluster = imf.make_cluster(mcluster, silent=True)
    obcount.append((cluster>8).sum())
    ocount.append((cluster>20).sum())
    mclusters.append(mcluster)
    clusters.append(cluster)

In [4]:

close(1)
figure(1,figsize=(18,10)); clf()
plot(sort(mclusters),np.cumsum(np.array(obcount)[argsort(mclusters)])/1./np.sum(obcount),'.',label="$M_0=10^5$ OB Mod Schec")
plot(sort(mclusters),np.cumsum(np.array(ocount)[argsort(mclusters)])/1./np.sum(ocount),'.',label="$M_0=10^5$ O Mod Schec")
semilogx(sort(mclusters),np.array(mclusters)*0+0.5,'r--')
xlabel("Cluster Mass $M_i$")
ylabel("Fraction of Clusters with $M<M_i$")

Out [4]:

<matplotlib.text.Text at 0x10e71b5d0>

In [5]:

imf.inverse_imf(0.99, massfunc=imf.modified_schechter, m1=100, m0=1e5, mmin=10, mmax=1e7)

Out [5]:

7619.788265653388

In [6]:

max(mclusters)

Out [6]:

311567.90346958087

In [7]:

np.cumsum(obcount)

Out [7]:

array([     1,      1,      1, ..., 307506, 307506, 307594])

In [8]:

argsort(mclusters)

Out [8]:

array([40808, 22411, 41334, ..., 23973, 21675, 31040])

In [9]:

nclusters = 50000
Sclusters = []
Sobcount = []
Socount = []
Smclusters = []
for rn in np.random.rand(nclusters):
    mcluster = imf.inverse_imf(rn, massfunc='schechter', mmax=1e7, mmin=10, m0=1e5)
    cluster = imf.make_cluster(mcluster, silent=True)
    Sobcount.append((cluster>8).sum())
    Socount.append((cluster>20).sum())
    Smclusters.append(mcluster)
    Sclusters.append(cluster)

In [10]:

fig = figure(1)
plot(sort(Smclusters),np.cumsum(np.array(Sobcount)[argsort(Smclusters)])/1./np.sum(Sobcount),'.',label="$M_0=10^5$ OB Schechter")
plot(sort(Smclusters),np.cumsum(np.array(Socount)[argsort(Smclusters)])/1./np.sum(Socount),'.',label="$M_0=10^5$ O Schechter")
legend(loc='best')

Out [10]:

<matplotlib.legend.Legend at 0x19f9b3810>

In [11]:

nclusters = 50000
S4clusters = []
S4obcount = []
S4ocount = []
S4mclusters = []
for rn in np.random.rand(nclusters):
    mcluster = imf.inverse_imf(rn, massfunc='schechter', mmax=1e6, mmin=10, m0=1e4)
    cluster = imf.make_cluster(mcluster, silent=True)
    S4obcount.append((cluster>8).sum())
    S4ocount.append((cluster>20).sum())
    S4mclusters.append(mcluster)
    S4clusters.append(cluster)

In [12]:

fig = figure(1)
plot(sort(S4mclusters),np.cumsum(np.array(S4obcount)[argsort(S4mclusters)])/1./np.sum(S4obcount),'.',label="$M_0=10^4$ OB Schechter")
plot(sort(S4mclusters),np.cumsum(np.array(S4ocount)[argsort(S4mclusters)])/1./np.sum(S4ocount),'.',label="$M_0=10^4$ O Schechter")
legend(loc='best')

Out [12]:

<matplotlib.legend.Legend at 0x10e7fead0>

In [13]:

sum(S4mclusters)

Out [13]:

3125486.249456272

In [14]:

# this scenario is where the smallest cluster is a single star, so the exponential turnover is at 1 msun, and the minimum cluster mass is 0.1 msun
# need to go to larger numbers for this one, of course
nclusters = 100000
allclusters = []
allobcount = []
allocount = []
allmclusters = []
for rn in np.random.rand(nclusters):
    mcluster = imf.inverse_imf(rn, massfunc=imf.modified_schechter, mmax=1e8, mmin=0.1, m0=1e4, m1=1)
    cluster = imf.make_cluster(mcluster, silent=True)
    allobcount.append((cluster>8).sum())
    allocount.append((cluster>20).sum())
    allmclusters.append(mcluster)
    allclusters.append(cluster)

In [15]:

fig = figure(1)
plot(sort(allmclusters),np.cumsum(np.array(allobcount)[argsort(allmclusters)])/1./np.sum(allobcount),'.',label="$M_0=10^4 M_1=1$ OB Mod Schechter")
plot(sort(allmclusters),np.cumsum(np.array(allocount)[argsort(allmclusters)])/1./np.sum(allocount),'.',label="$M_0=10^4 M_1=1$ O Mod Schechter")
legend(loc='best')

Out [15]:

<matplotlib.legend.Legend at 0x19f9ddd10>

In [16]:

A = argmin(abs(0.5-np.cumsum(np.array(allobcount)[argsort(allmclusters)])/1./np.sum(allobcount)))
sort(allmclusters)[A]
A = argmin(abs(0.5-np.cumsum(np.array(allocount)[argsort(allmclusters)])/1./np.sum(allocount)))
sort(allmclusters)[A]

Out [16]:

689.71604560136109

In [21]:

%pastebin 0-16 raw=False

Out [21]:

u'https://gist.github.com/3751918'