estimates for dinflo


We begin by summarizing some symbol definitions, and typical values for \zeta Pup.

  • \alpha = 0.65 is the CAK power law exponent
  • \overline{Q} = 2000 is the quality factor of the wind
  • L_* = 8 \times 10^5 L_{\odot} is the luminosity of our massive star
  • M_* = 50 M_{\odot} is the mass of our star
  • R_* = 20 R_{\odot} is the radius of our star
  • G is the gravitational constant
  • c is the speed of light
  • \kappa_e = 0.34 \text{ cm}^2 / \text{g} is the free electron opacity
  • k = 1.38 \times 10^{-16} \text{ eV}^2 / \text{K} is the Boltzmann constant
  • T = 4 \times 10^4 \text{ K} is the star surface temperature
  • m_H = 1.67 \times 10^{-24} \text{ g} is the Boltzmann constant
  • \mu = 1.3 is the mean particle mass divided by the hydrogen mass

We now also introduce the Eddington parameter \Gamma_e, the ratio of the local acceleration from radiation to the gravitational acceleration,

    \[\Gamma_e = \frac{ \kappa_e L_* }{ 4 \pi G M_* c }.\]

For example, we can write the escape velocity (accounting for radiation) as

    \[v_{esc} = \sqrt{ \frac{2 G M_* ( 1 - \Gamma_e )}{R_*} }.\]

Then the mass loss rate predicted by vanilla CAK is

    \[\dot{M}_{cak} = \frac{L_*}{c^2} \frac{\alpha}{1-\alpha} \left( \frac{\overline{Q} \Gamma_e}{1 - \Gamma_e} \right)^{(1-\alpha)/\alpha}\]

The finite disk correction factor is (1+\alpha)^{-1/\alpha}, so

    \[\dot{M}_{fd} = \frac{1}{ (1 + \alpha)^{1/\alpha} } \dot{M}_{cak}.\]

Finally, the finite sound speed introduces yet another correction factor,

    \[\dot{M}_{fd,fs} = \left( 1 + \frac{4 \sqrt{1-\alpha}}{\alpha} \frac{\alpha}{v_{esc}} \right) \dot{M}_{fd}.\]

Now, we calculate the sonic point density \rho_{s}. The sonic point is the location in the wind where the velocity is equal to the sound speed a,

    \[v(R_s) = a,\]

where the isothermal sound speed is

    \[a = \sqrt{ \frac{k T}{\mu m_H} }.\]

Here comes the spooky hydro stuff.


We choose the inflow density to be the density at the sonic point, multiplied by some factor around 5-10. We call this factor d_{fac}, so that

    \[\rho_{inflow} = d_{fac} \rho_s.\]

We define the scale height H = a^2 / g. Then we estimate the sonic point,

    \[R_s = \log (d_{fac}) H + R_{lb}.\]

Here, R_{lb} is the location of the first grid point in the simulation. We can then finally write

    \[\rho_{s} = \frac{ \dot{M} }{4 \pi R_s^2 a}.\]

We can then calculate \rho_{inflow}, and use it in our simulations as the lower boundary condition.

movie making in Python

A picture might be worth a thousand words, but nothing beats a movie to convey time-domain information. In this post, I’ll walk through how you can turn your astro data into a cinematic masterpiece, using Python.

The final product of this post.
The final product of this post.

The first step is to have some video library installed like libav or ffmpeg. If you’re using Ubuntu/Debian, you can run the following command:

sudo apt-get install libav-tools

If you’re using OSX, you can either compile it yourself, install it via homebrew, or get mencoder/ffmpeg via some external package. The easiest way is the last one — installing MPlayer should give you both mencoder and ffmpeg, both of which will work.

If you’re interested in additional pain or you already have homebrew, you can follow these instructions instead. It will do the same thing. 

Once you have this, we can begin with the Python code! Let’s start by importing the various libraries we’ll need.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as manimation

Let’s generate some synthetic data to plot. We’re going to make a Python list, where each object in the list is really a pair of two other lists — namely, the lists of X-values and Y-values of our plot.

# we generate some synthetic data.
#we will make a list of X, Y lists.
frames = []
frameNumber = 240 # we generate 240 different frames

for i in range(0,frameNumber):
xData = np.linspace(0,2*np.pi)
yData = np.sin(xData - i * 2.0*np.pi / frameNumber )
frames.append( [xData, yData] )

Now let’s use the movie writer functionality we imported. We need to initialize the writer with various movie parameters. In the following code, make sure to replace “avconv” with your encoding library!

# now we put them together to make a movie! let's set up the movie writer
framerate = 60 # 24 frames per second
FFMpegWriter = manimation.writers['avconv']
metadata = dict(title='Wave Data', artist='Isaac Newton', comment='')
writer = FFMpegWriter(fps=framerate, metadata=metadata)

Now, let’s set up the matplotlib figure. Here you can customize the title, axis labels, the kind of plot, etc. Here we use a line plot.

# figure setup
fig = plt.figure()
firstFrameX, firstFrameY = frames[0]
l, = plt.plot(firstFrameX, firstFrameY, '-')

plt.ylabel( r'$x$ axis' )
plt.xlabel( r'$y$ axis' )
plt.title( r'$\lambda = 2 \pi$' )

The final step is just to start saving images to the writer, while changing the data of the plot to match that of our frames.

# let's write to the file!
with writer.saving(fig, "anim.mp4", 100):
for i in range(frameNumber):
x, y = frames[i]
l.set_data( x, y)


We’re done! There should be a video file in your directory, once you run this script.

Note that it’s not necessary to do exactly what we did here, and change the data of an existing plot. Try changing the l.set_data call to a plt.plot. You’ll see a very different video! Now add a line after your plotting that clears your plot after every frame. You’ll see that you have restored the video to the previous animation.

LDI update – I think I broke something?

I had run a standard run of Jon’s VH-1 code, and I was expecting to plot the output and then assemble a short video. Unfortunately, the output is filled with NaN values!

This is really disappointing, but it seems like restarting research will not be without hiccup. I’ll have to consult Jon’s old emails and reset my VH-1 parameters. It appears that I must have tried something strange with my old code, the last time I edited it, and then did not revert the changes.

Lesson: I should use git.

reference — ifort

I’m actually enjoying the use of these reference posts. A lot of stuff I figure out how to do, and then a year later, I have to reread the documentation and relearn how to do it.

When writing Fortran code, there seems to be two choices,

  1. gfortran – GNU Fortran project, free as in freedom
  2. ifort – Intel Fortran compiler, free as in beer (for non-commercial use)

On a Linux machine, it’s usually very easy to get gfortran installed and configured. However, there are a few minor advantages to ifort. From some naive tests I ran as a freshman, it appears ifort optimizes more aggressively, leading to ~30% speed improvement. My laptop has an Intel CPU, so some processor-based optimization techniques may be at play here. More importantly, people (I think Jon?) have suggested I change the aging radiative cooling functions in our group’s LDI code. About a year ago Dylan Kee sent me some code by Rich Townsend that implemented the cooling method described here. It appears that his code makes use of some helper functions unique to ifort. I’m super lazy I’ve been very busy with schoolwork, but when I get to this, I’m sure it’ll be easier to use ifort.

The Intel Fortran compiler now comes packaged with Intel® Parallel Studio XE, which for students/academics can be obtained on the Intel website. They’ve improved the installation process dramatically over the last four years — I remember having to do some pretty arcane things to get it installed correctly on Ubuntu, but now there’s a GUI installer which makes it all very easy. I add the following line to my .bashrc,

alias intel="source /opt/intel/bin/ intel64"

This adds the intel compiler environment to my PATH whenever I run the command “intel” in bash.

reference — loading atomDB FITS files

Analyzing data in FITS files is a pretty common astronomical task, and it’s easy with Python. In this reference post, I’ll put together a script that loads the AtomDB emissivity data into Python. I like to use a text editor and console, but this can certainly also be done in IPython. If you don’t have Python installed yet, I recommend Anaconda Python, a batteries-included Python distribution for science.

A hot glue gun.
A hot glue gun.

We’ll use matplotlib, numpy and astropy — the wrench, hammer, and hot-glue of the Python astro toolbox. I’m using matplotlib 1.4.2 and astropy 1.0.4. Create a directory for your project and put a .py file in it with a descriptive name (like ““). For this example, we’ll be extracting emissivities from FITS files, so download the current version of AtomDB. Extract the files (with “tar -xjvf yourfilehere“) and put it in your project directory.

Let’s write our Python script! Suggested use of this reference: copy and paste the code block by block, running the cumulative script at each step. Each step will output some text, until the end. Play around with it!

We’ll start off by loading our libraries,

# get the libraries we want to use to read FITS files
from import fits

# we'll also want to plot our data and manipulate it
import matplotlib.pyplot as plt
import numpy as np

Now let’s load in the FITS file. I’ve added some debug printing throughout, so you can see what’s going on. You should replace the string in the argument of with whatever FITS file in your AtomDB folder ends with “_linelist.fits“.

# load in the file
hdulist ='atomdb_v3.0.2/apec_v3.0.2_linelist.fits')

# check out the data stuff, we're looking for the emissivities...

Now we get the specific emissivity data we’re looking for, in a table.

# it looks like the emissivities are in the hdulist[1]
db = hdulist[1].data

# let's check out what we've got!

We’re basically done now, since the object db is now a huge table filled with all the AtomDB line data. The rows are specific lines (i.e. db[0] will give you a list with a wavelength, an error, some emissivities, the transition data, see here for more) while the columns are the different types of data. Let’s extract information about the first line in the list.

# Let's display the data for very first line in our linelist
e = db['Emissivity'][0]
t = db['Temperature'][0]
wavelength = db['Lambda'][0]
print e

Now we use matplotlib to make some graphical output,

# plot our two vectors
plt.xlabel('Temperature (K)')
plt.ylabel('Emissivity (blah blah blah units)')
plt.title( str(wavelength) + " Angstrom" )

The end result should be this reassuring plot of the emissivity as a function of temperature.




Some final notes:

  • You might have noticed the strange corner at the end of the plot above. That’s because the temperature and emissivity lists have 0’s at the end of them for some reason. One could maybe trim these lists before plotting.
  • We’re using AtomDB’s data grouped by line, instead of temperature. More information on the formats available are here.
  • This script accesses AtomDB’s data for specific lines. A useful exercise would be to retrieve the emissivity from the continuum via coco. The FITS file should already be in your AtomDB distribution.
  • The db object in this script is an astropy table. Play around with the table. It’d be useful to have a function that searches the table for rows (aka lines) which have a wavelength near a target wavelength. I’ll actually update this post with such a function, but it’s all just application of the things you can do with a table.

delta ionization parameter

I did some implementation of the delta ionization parameter. I’m using this expression,

    \[W = 0.5\left(1-\sqrt{1 - \frac{R^2}{r^2}\right)\]

    \[g_{\delta}(r) = \left(\frac{\rho(r)}{m_p W} \right)^{0.1} ( g_{rad}(r) + g_{scat}(r) )\]

Here, I’m using m_p = 1.67266e-24 grams.

I suspect the line force is somehow strengthened by this (I’m not fully sure if this is making the wind stronger or weaker…), but it’s making the terminal velocity much higher. This simulation used Q/Q_{max} = 0.00395, so I’ll next use try Q/Q_{max} = 0.01 and Q/Q_{max} = 0.1. Hopefully both work.

Multi-D LDI: 2HD+1R via my (flawed) SSF implementation

I have my own implementation of the SSF radiation force with the hydro code VH-1, which has some sort of lingering bug associated with the inner wind. I’ve decided to instead work with Jon Sundqvist’s implementation, so one of my recent projects has been to move the LDI simulation up to 2D.

VH-1 is pretty straightforward in setting the dimension, but Jon’s code has some sort of bug that makes it crash whenever I turn on the second dimension. To be frank, I have a hard time reading other people’s code, so I figured I would first try it on my own SSF code. It works pretty well:

The problem is that my own SSF implementation has a small bug which makes it not have the stable CAK solution very close to the star.

The only hurdle: I originally had a log-gridding function, which was causing crashing. I disabled it, which fixed everything. This is maybe a clue about why Jon’s code fails to run in 2D.