Adam Ginsburg's blog

• 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'
  

---

• Lognormal Density / IPython Notebook

  An experiment with log normal density distributions (simple, pretty obvious) and a test of including ipython notebooks in blog posts. To do this, I created a gist (https://gist.github.com/3750781) and used the ipython notebook viewer (http://nbviewer.ipython.org/3750781/) to render the notebook in HTML. It worked pretty well; all I had to do was c&p the source and then change 2 lines of the style (background color -> black, foreground color -> white for the markup... maybe I should do it for the rest too).

  Comparison of the derived power spectra and density PDFs for lognormal power distributions. Obviously (if you think about it), the power spectra are always flat for randomly distributed signal. The PDFs change significantly.

---

Page 1 / 1