Installing diffuse or meld - failure, can't install pygtk

Tried to install diffuse today.  Had to first install pygtk.  pygtk requires glib.  glib doesn't install defaultly.

To install glib, downloaded it (which wasn't easy, required xz too) and ran this command:

PKG_CONFIG_PATH=/usr/lib/pkgconfig:/usr/local/lib/pkgconfig/ ./configure

glib install worked, now pygobject fails:

checking for GLIB - version >= 2.8.0... yes (version 2.35.2)
checking for PYGOBJECT... no
configure: error: Package requirements (pygobject-2.0 >= 2.21.3) were not met:

Requested 'pygobject-2.0 >= 2.21.3' but version of PyGObject is 2.18.0

Consider adjusting the PKG_CONFIG_PATH environment variable if you
installed software in a non-standard prefix.

pygobject requires gobject-introspection:

http://ftp.gnome.org/pub/GNOME/sources/gobject-introspection/1.31/

do this again:

PKG_CONFIG_PATH=/usr/lib/pkgconfig:/usr/local/lib/pkgconfig/ ./configure

make -j 4'

failure:

  CC     glib_print-glib-print.o

make[2]: *** No rule to make target `/glib-2.0/include/glibconfig.h', needed by `GLib-2.0.gir'.  Stop.

make[2]: *** Waiting for unfinished jobs....

examples/glib-print.c: In function 'main':

examples/glib-print.c:11: warning: 'g_type_init' is deprecated (declared at /usr/local/include/glib-2.0/gobject/gtype.h:669)

girepository/gi-dump-types.c: In function 'main':

girepository/gi-dump-types.c:12: warning: 'g_type_init' is deprecated (declared at /usr/local/include/glib-2.0/gobject/gtype.h:669)

make[1]: *** [all-recursive] Error 1

make: *** [all] Error 2

make -> failure again:

  GEN    g-ir-doc-tool

make[2]: *** No rule to make target `/glib-2.0/include/glibconfig.h', needed by `GLib-2.0.gir'.  Stop.

make[1]: *** [all-recursive] Error 1

make: *** [all] Error 2

well, I give up.

Installing qt on 10.7

Worked following these instructions: http://www.noktec.be/python/how-to-install-pyqt4-on-osx with important exceptions.

1. I installed sip and pyqt against my virtualenv, not against the system python
2. I had to get g++-4.2 to precede g++-4.6 on my path, which was hard.

This command:
PATH=/usr/bin:$PATH CC=g++-4.2 CXX=g++-4.2 ~/virtual-python/bin/python configure.py -d ~/virtual-python/lib/python2.7/site-packages/ --verbose

worked, whereas this one:

CC=g++-4.2 CXX=g++-4.2 ~/virtual-python/bin/python configure.py -d ~/virtual-python/lib/python2.7/site-packages/ --verbose

did not.

Similarly, needed same prefixes when doing make.

test

Using the numbers we measure for $L(H\alpha)$, I'll re-derive key quantities about SN 2009ip.

This analysis implicitly assumes that H$\alpha$ is optically thick at $10^4$ K. I believe we should make that assumption explicit in the text.

In [1]:

from astropy import units, constants
pc = units.pc.in_units('cm')
Iha = 2.9e-14 # erg s^-1 cm^-2 A^-1
D = 24e6 * pc
Lha = Iha * D**2
print Iha," erg/s/cm^2/A"
print D," cm"
print Lha," erg/s/A"

2.9e-14  erg/s/cm^2/A
7.40562619552e+25  cm
1.59045568109e+38  erg/s/A

In [2]:

# Use the wavelength form of the Planck function
h = constants.cgs.h
k = constants.cgs.k_B
c = constants.cgs.c
def planck(lam, temp):
    """ erg s^-1 cm^-2 A^-1 sr^-1 """
    return 2*h*c**2 / lam**5 * (exp(h*c/(k*temperature*lam)) - 1)**-1

In [3]:

angstrom = units.angstrom.in_units('cm') # angstroms
wave = 6562.8*angstrom
temperature = 1e4 # K
print "sanity check: wave * 1e8 = 6562.8 = ",wave*1e8
print "planck(6562.8,1e4) = ",planck(wave,temperature)," erg/s/cm^2/cm"
print "planck(6562.8,1e4) = ",planck(wave,temperature)/1e8," erg/s/cm^2/A"

sanity check: wave * 1e8 = 6562.8 =  6562.8
planck(6562.8,1e4) =  1.22966749446e+15  erg/s/cm^2/cm
planck(6562.8,1e4) =  12296674.9446  erg/s/cm^2/A

We're trying to solve an equation of the form:
$$ L(H\alpha) = \pi R^2 B_\lambda(T) d\lambda $$
Since our luminosity is in erg s$^{-1}$ $\dot{A}^{-1}$, we don't have to worry about the bandwidth.

$$ R = (L(H\alpha) / \pi / B_\lambda(T))^{1/2}$$

In [4]:

Aemitting = Lha / (planck(wave,temperature)*pi/1e8)
print "Emitting area = ",Aemitting," cm^2"
R = (Aemitting/pi)**0.5
print "R = %e" % R, " cm"
au = units.AU.in_units('cm')
print "R = ",R/au," AU"
print "R = ",R/pc, " pc"

Emitting area =  4.11702975893e+30  cm^2
R = 1.144767e+15  cm
R =  76.5229413217  AU
R =  0.000370993688718  pc

In [5]:

# Light crossing time
tc = R / constants.cgs.c
day = units.day.in_units('s')
print "Light Crossing time at %i AU = %f days" % (R/au, tc/day)

Light Crossing time at 76 AU = 0.441960 days

In [6]:

# What if we assume it's a sphere instead of a disk?
# area -> 4 * pi * r^2
R = (Aemitting/(4*pi))**0.5
print "R = %e" % R, " cm"
au = units.AU.in_units('cm')
print "R = ",R/au," AU"
print "R = ",R/pc, " pc"
tc = R / constants.cgs.c
print "Light Crossing time at %i AU = %f days" % (R/au, tc/day)

R = 5.723835e+14  cm
R =  38.2614706609  AU
R =  0.000185496844359  pc
Light Crossing time at 38 AU = 0.220980 days

Because the brightening was sudden, and because a uniform disk is somewhat implausible, we should come up with a different measurement for a disk with an inner hole. I'll use fiducial numbers of 13000 $km s^{-1}$ and $\Delta t \sim 40$ days.

In [14]:

kms = (units.km/units.s).in_units('cm/s')
vejecta = 13e3 * kms
deltat = 40 * day
r_collision = vejecta*deltat
print "Inner disk radius = ",r_collision," or ",r_collision/au," AU"
print "Light crossing time at inner disk radius = ",r_collision / constants.cgs.c / day, "days"
H = Aemitting/r_collision/2/pi
print "Height needed if this is ALL coming from the inner part of a disk: ",H," cm or ",H/au," AU"
print "The speed required to achieve that height in 2 days is ", H / (2*day) / 1e5, "km/s"

Inner disk radius =  4.4928e+15  or  300.325130229  AU
Light crossing time at inner disk radius =  1.73453329503 days
Height needed if this is ALL coming from the inner part of a disk:  1.45843491139e+14  cm or  9.74903522735  AU
The speed required to achieve that height in 2 days is  8440.01684832 km/s

We then want $$\int_{r_{inner}}^{r_{outer}} 2 \pi r dr = A_{emitting}$$ or $$ r_{outer}^2 - r_{inner}^2 = A_{emitting} / \pi$$

In [8]:

router = (Aemitting/pi+r_collision**2)**0.5
print "Outer radius = ",router," or ",router/au," AU"
print "Delta-R = ",router-r_collision," or ",(router-r_collision)/au," AU"
tcross = (router-r_collision) / vejecta
print "The shock crossing time from inner->outer at %g km/s = %g days" % (vejecta/kms, tcross/day)

Outer radius =  4.6363501932e+15  or  309.920867958  AU
Delta-R =  1.43550193199e+14  or  9.59573772859  AU
The shock crossing time from inner->outer at 13000 km/s = 1.27805 days

I think this effectively leaves us with the conclusion that there is a 30 AU radius torus at 300 AU from the central source. That strikes me as somewhat implausible; 30 AU is enormous, particularly at densities $\sim10^{13}$. I'll try computing the mass in the disk assuming it's quite thin, $\sim0.1AU$.

In [9]:

height = 0.1 * au
volume = Aemitting*height
density = 1e13 
nparticles = volume * density
msun = units.msolMass.in_units('g')
mass = nparticles * constants.cgs.m_p * 1.4
print "mass = ",mass,"g = ",mass/msun,"msun"

mass =  1.44223224375e+32 g =  72.5067741062 msun

That's ridiculous. I'll invert the assumption: say the disk/torus is 1 msun; how thin must it be?

In [10]:

mass = 1 * msun
nparticles = mass / (constants.cgs.m_p*1.4)
volume = nparticles/density
height = volume/Aemitting
print "height of disk = ",height,"cm = ",height/au," AU"

height of disk =  20632261267.2 cm =  0.0013791814797  AU

These numbers don't feel right to me. Why such an ultrathin disk? It would HAVE to puff up; the pressures are enormous.

What about the hypothesis that the brightening results from optically thin emission? I'll work under the assumption that H$\alpha$ is "just barely" optically thick (i.e. I'll cheat and use the optically thin approach);
in order to get the high H$\beta$ / H$\alpha$ ratio we actually need it to be pretty optically thick and under the Case B regime (Case C just makes it harder to get the observed $\alpha$/$\beta$ ratio).

In [18]:

density = 1e6 # assume pretty dense, but still plausibly case B
alphaB = 2.6e-13 # cm^3 s^-1
Eha = h * c / wave
nphot_perang = Lha / Eha # phot/s/A
nphot = nphot_perang * 1 # phot/s
# solve for path length
# nphot = density**2 * alphaB * area * path_length
path_length = nphot / (density**2*alphaB*Aemitting)
print "Emitting path length for density = %0.1e" % density,"cm^-3 is ",path_length," = ",path_length/pc," pc"
# use the rise time of ~2 days as an upper limit on the emitting size scale
maxsize = 2*day*c
print "This is %i times larger than the largest allowed scale" % (path_length/maxsize)
print "Densities must be at least %e for optically thin emission to be plausible" % (density*sqrt(path_length/maxsize))

Emitting path length for density = 1.0e+06 cm^-3 is  4.90881531249e+19  =  15.9083870006  pc
This is 9475 times larger than the largest allowed scale
Densities must be at least 9.734331e+07 for optically thin emission to be plausible

So at densities $>5\times10^7$, if H$\alpha$ starts becoming optically thick but H$\beta$ remains optically thin, we could explain the brightening by a filled sphere of just barely optically thick emission.

Despite all this, we still don't have a great explanation for why the P Cygni profile disappears. We should revisit that if possible.

In [22]:

# limb-brightened shell model
half_sph_area = 2*pi*r_collision**2
print "Fraction of sphere that must be emitting: ",Aemitting/half_sph_area

Fraction of sphere that must be emitting:  0.0324616032628

In [1]:

from IPython.core.display import HTML

def blackbg():
	"""Black Background code mirror theme."""
	html = """
<style type="text/css">
.cm-s-ipython { background-color: black; color: lightblue; }
.cm-s-ipython span.cm-keyword {color: #00ff00; font-weight: bold;}
.cm-s-ipython span.cm-number {color: #ee88ee;}
.cm-s-ipython span.cm-operator {color: lime; font-weight: bold;}
.cm-s-ipython span.cm-meta {color: white;}
.cm-s-ipython span.cm-comment {color: cyan; font-style: italic;}
.cm-s-ipython span.cm-string {color: red;}
.cm-s-ipython span.cm-error {color: darkred;}
.cm-s-ipython span.cm-builtin {color: pink; font-weight: bold;}
.CodeMirror pre.CodeMirror-cursor {color: white; border-left: 1px solid white;}
.cm-s-ipython span.cm-variable {color: white;}
div.text_cell_input  {   background-color: black;    color: white;}
div.text_cell {   background-color: black;    color: white;}
div#notebook {    background-color: black;     background-image: none;     color: white;}
div.code_cell { background-color: black;    color: white;}
div.metaedit .maintoolbar{    background-color: black;    color: white;}
body{    background-color: black;    color: white;}
#notebook_name {  background-color: #333333;    color: white; background-image:none;}
span.ui-widget-content{  background-color: #333333;    color: white; background-image:none;}
div.ui-widget-content{  background-color: #333333;    color: white; background-image:none;}
div#ui-button ui-widget ui-state-default ui-corner-all ui-button-text-only output_collapsed vbox{  
  background-color: black;    color: white;}
div.input_area{    background-color: black;    color: white;}
div.ui-widget-content { background-color: #333333; background-image: none;}
div.input_prompt {color:#6666CC; font-weight: bold;}
div.output_prompt {color:#CC6666; font-weight: bold;}
div.output_text {color:white;}
div.text_cell_render {color:white;}
span.item_name {color:white;}
.ui-widget-content a {color:white;}
div.tooltip .ui-corner-all .tooltiptext .smalltooltip {background-color: black; color:white;}
.pln {color:white} 
.typ {color:#f9f}
.lit {color:#BBB}
.kwd {color:#3e3}
.metaedit {background: transparent; border-color: transparent; color:white;}
.ui-button-text {background: #112; color:white; border-color:666;}
.ui-menubar  {background-color: transparent; background: transparent; background-image:none;}
.ui-widget-header  {background-color: transparent; background: transparent; background-image:none;}
.ui-helper-clearfix {background-color: transparent; background: transparent; background-image:none;}
.ui-corner-all {background: #112; color: white;}
.ui-widget  {background: #112; color: white;}
.ui-widget-content {background: #112; color: white;}
div.highlight {background:black;} /* needed for nbviewer */
</style>"""
	return HTML(html)

blackbg()

Out[1]:

ipython notebooks - preserve figure status between cells

To preserve figures between cells (rather than initializing a new one each time), do:
%config InlineBackend.close_figures=False

stackoverflow post answering this that needs more upvotes

netpbm install

I needed to install netpbm in order to use astrometry.net, but the install sucked.

First, had to build with a bunch of flags:
CC=/usr/bin/gcc-4.2 CFLAGS="-arch i386 -arch x86_64" FFLAGS="-m32 -m64" LDFLAGS="-Wall -undefined dynamic_lookup -bundle -arch i386 -arch x86_64"  MACOSX_DEPLOYMENT_TARGET=10.6 CFLAGS_SHLIB="-fno-common" NETPBMLIBTYPE=dylib NETPBMLIBSUFFIX=dylib LDSHLIB="--shared" TIFFLIB=-ltiff JPEGLIB=-ljpeg PNGLIB=-lpng ZLIB=-lz make  (some of those probably aren't necessary, and oops my deployment target should've been 10.7 but I don't think it matters.

Then, I got error messages that looked like this:
strip: symbols referenced by indirect symbol table entries that can't be stripped in: /private/tmp/netpbm/bin/#inst.36429# (for architecture i386)
It was caused by this line:
 /Users/adam/repos/netpbm-10.35.86/buildtools/install.sh -c -s -m 755          atktopbm /tmp/netpbm/bin
I thought that looked a little suspicious, so I tried removing the "-s" (for "stripprg").  I had to remove it in the buildtools/install.sh directory.  Felt like a hack, but apparently worked.  Except that now I get dyld errors:

dyld: Library not loaded: /libnetpbm.10.dylib
  Referenced from: /usr/local/netpbm/bin/pnmfile
  Reason: image not found
Trace/BPT trap: 5

Solution:
export DYLD_LIBRARY_PATH=/usr/local/netpbm/lib/
in ~/.bashrc

UPDATE: Installing on 10.6.8 was trivial, ironically.

Guider Pointing & Blog Motion

There's a new post on a new idea for guiding with, e.g., APO at adamginsburg.blogspot.com, which is where I intend to move all my blogging.

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'