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 imffigsize(12,8)

reg_gal2cel requires coords & pyregionRegion Photometry requires pyregioncubes.py requires pyregion for getspec_regcubes.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 massmrange = logspace(1,5,1000) # from 10 msun to million msun clustersschec = imf.modified_schechter(mrange,m0=1e4,m1=100) # cutoff at low massfigure()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>

(figure removed because it was breaking blog renderer)

In [3]:

nclusters = 50000clusters = []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>
:align: center :alt:

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 = 50000Sclusters = []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 = 50000S4clusters = []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 coursenclusters = 100000allclusters = []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'