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]:

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.

GINSBURG

Sunday, September 22, 2013

Blog moved!

Wednesday, December 12, 2012

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

Tuesday, December 04, 2012

Installing qt on 10.7

Wednesday, November 21, 2012

test

Thursday, November 08, 2012

ipython notebooks - preserve figure status between cells

Sunday, November 04, 2012

netpbm install

Friday, October 05, 2012

Guider Pointing & Blog Motion

Old Posts

Categories