fast windtabs in Mathematica

I’ve been working on better error estimates in the Np(T) calculations, and I needed a way to quickly calculate the T_w transmission fraction produced by the windtabs model. The point of windtabs is to have a big table full of data (hence “tabs”) to use for model fitting, but I didn’t have access to that table (it’s in a weird XSPEC file format).

The answer is pretty straightforward–make my own table, and interpolate over it. The \beta=1 windtabs transmissions are pretty smooth, so I didn’t need that many points. Using a 30 by 30 grid evenly spaced across [1.01,2] for R_0 and [0,10] for \tau_*, we get a pretty good fit. The plot below is for R_0 = 1.5 R_*, with an exponential absorption model plotted for comparison in red, the full double integral in blue, and the interpolation function in gold. Note the deviation because I only used \tau_* < 10.

If we explore this parameter space a bit, we find that it was indeed pretty smooth (below is \ln T_w).

The end result is that I can really quickly generate T_w values for a particular \tau_* and R_0 combination, for use in a Monte Carlo error estimate. Below is a T_w histogram for normally distributed \tau_* and R_0.

Unfortunately, I’m now stuck! I need an input probability distribution to model the asymmetric error on \tau_*, but I have no idea what it is! Some talking with David suggests that the base distribution (and upper/lower errors) comes from \chi^2 fitting, so some exploration of what XSPEC is doing might be fruitful.

You can get the (rather messy) Mathematica notebook here.

 

Posted in exploration, NpT | Leave a comment

possible trends in spectral type vs. Np(T)

All five of our program stars in the Np(T) paper share a particular Mg XI line at 8.421 Å. If we compare the shock probability on this line for each star, we get an interesting trend, where earlier spectral subtype seems to imply a higher shock probability.

Posted in exploration, NpT | Leave a comment

Finishing the Np(T) Paper

Various tasks remain so that we can publish the \bar{N}p(T) results.

  • Fig 1 – change the color and label from the weak N VII line to the stronger and lower-temperature N VI line complex that you just added. Done.
  • In both Fig 1 and Fig 2, the y-axis label should have “ergs” not “erg.” Done.
  • Update gray box plots with no high temperature cutoff. Done.
  • Update rainbow plot with no high temperature cutoff. Done.
  • Calculate error estimates on N VI. Awaiting Monte Carlo Work.
  • Double check abundance adjustments (zeta Pup and xi Per) Done.
  • Recompute using \alpha = 5/2. Done.

Figure 1. Changed highlighting and label from N VII to N VI.

Figure 2. Changed “erg” to “ergs”.

Posted in NpT | Leave a comment

effects of the 2 * 10^7 K cutoff

In my Np(T) work, I’ve used an upper integration limit of 2 \times 10^7 \text{ K}. This produces the gray box plot of zeta Puppis below.

However, taking this integration limit to 2 \times 10^8 \text{ K} (effectively infinity), has minimal effects but causes the gray box plots to move closer to power law.

Posted in figure, NpT | Leave a comment

working on multi-D: step 1, existing code in 1D

Jon sent me his LDI implementation in the F90 VH-1 version. I had some initial issues with getting it functioning, as I wasn’t familiar with the configuration settings. I thus had some issues with the mass loss rate:

It’s possible that this the fact that I had not used the isothermal option. Jon sent me an updated INDAT file, and we get the expected results:

Posted in multiD, simulation | Leave a comment

research progress for the prior week

  • Wrote some Fortran code to add to VH-1 for parcel tracking, but it doesn’t work. I think at very small time-steps (2.5 s, which is 40 times smaller than the post-simulation tracking) parcels do not have continuous decreases–probably small fluctuations.
  • I made some animations showing density as color on a velocity/distance plot of the simulation output. They’re basically what you expect–reverse shocks are low density but the post-shock gas is very dense.
Posted in update | Leave a comment

parcel tracking animated

I plotted the parcels I track via the Lagrangian mass coordinate method on top of the advecting flow.

NOTE: The X-axis is mislabeled here! This is a plot of velocity vs. height, not mass coordinate.

Parcels seem to get stuck in the shock fronts. This makes sense if we consider where the wind is most dense (in the clumps, which are after the shock fronts!)

Here, line darkness represents density on this velocity vs. height plot.

Posted in exploration, simulation | Leave a comment

varying the kappa cutoff – shock tracking

This calculation was performed with the shock tracking method described in detail here.

One of the annoying open problems in this LDI business is that every simulation has been running with an artificially capped \kappa_0 line driving constant. If we let \kappa_0 be its full value, our simulations crash (so many NaN outputs…)

The current Bartol group’s LDI implementation runs at 0.01 \kappa_0. The true value is 100 times what we’re using, which seems pretty bad!

To investigate this, I ran the shock tracking analysis on six different values of this so-called \kappa_{max} cutoff, linearly spaced from 0.001 to 0.01.

0.0001 0.0028 0.0046 0.0064 0.0082 0.01

Again, I relax the solution for 200 kilseconds and analyze the following 200 kilseconds. For the self-excited LDI, I plot the shock histogram of each \kappa_{max} value below.

As you can guess, increasing \kappa_{max} increases the shock temperatures. How about for the limb-darkened + perturbed LDI? The plot is similar.

Posted in exploration, figure, result, simulation | Leave a comment

shock tracking

Last week, I wrote a shock tracking method which focuses on individual shocks rather than the Lagrangian mass coordinate. It works like this: consider a snapshot of the wind velocity.

Instead of interpolating the velocity of a particle’s Lagrangian mass coordinate, we instead take each of the above shocks (theirs sizes marked in red), and convert the jumps to temperatures.

(1)   \begin{equation*}  T = \frac{3}{16} \frac{\mu m_h}{k} (\Delta v)^2 \end{equation*}

We use mass weighted binning, by adding the mass that flows through this shock in a single time-step to the correct bin. For example, we might start with bins like this:

1 \times 10^6 \text{ K} 2 \times 10^6 \text{ K} 3 \times 10^6 \text{ K} 5 \times 10^6 \text{ K} 6 \times 10^6 \text{ K}
0 0 0 0 0

We suddenly measure a single shock of temperature 3.5 \times 10^6 K, and it has velocity v(r), density \rho(r), and distance from the star r. Then for a simulation time-step \Delta t (in our case, 2.5 seconds), we would have a mass of

    \[\Delta m_{shock} = (\text{mass flux}) \times (\text{time}) = (4 \pi \rho(r) v(r) r^2) \times \Delta t\]

This is the mass flux from right before the shock hits. We add that mass to the appropriate bin.

1 \times 10^6 \text{ K} 2 \times 10^6 \text{ K} 3 \times 10^6 \text{ K} 5 \times 10^6 \text{ K} 6 \times 10^6 \text{ K}
0 0 \Delta m_{shock} 0 0

Now, if we run this at every time-step, we can make a plot that looks like this. I do some summation so that the mass weighted histogram is reverse cumulative like Np(T).

This is the result of analysis on purely self-excited shocks (black) and limb darkening + base perturbations (red), featuring 100 linearly spaced temperature bins. I ran the group’s LDI simulation for 200 kiloseconds to let the solution relax, and then ran my analysis code for another 200 kiloseconds.

Notes

  • This method doesn’t need interpolation or the Lagrangian mass coordinate. In my opinion, this makes it more robust.
  • However, I’m unsure if this is a completely unbiased sample!
  • The total shock rate N_0 is hard to extract from this method, unfortunately.

Finally, these results are not entirely off the mark when compared to the observations.

This plot is simply the Np(T) box plots, with the limb-darkened and perturbed histogram superimposed. The only parameter I adjusted was the normalization N_0, since I can’t calculate one from this method (yet! I have an idea…)

 

Posted in exploration, figure, NpT, result, simulation | Leave a comment

finding 68% intervals (tedious work I’ve been avoiding)

Vero suggested a while ago that instead of using FWHM, we could find the interval at which 68% of the total area of the ratio function \Lambda_{line}(T) / \Lambda_{tot}(T) is covered. However, there are many such intervals — to break the degeneracy, I stipulate \Lambda_l(a) = \Lambda_l(b), for an interval defined as (a, b). For robustness, I also stipulate that the bounds must be on opposite sides of the temperature of peak cooling T_{peak}, so a < T_{peak} < b.

It turns out this is pretty easy–just a lot of root finding. Below is the interval for oxygen VII, with the interval shown. I avoid using a log plot since it distorts areas significantly.

The interesting thing about this method is that it takes into account the long tails. Here’s two silicon and two magnesium lines.

If we look closer at the silicon lines, we note something pretty crazy–the interval extends over an incredibly long temperature range!

This does some bad things to our box plots. Actually, I think we all saw this coming–the box plots now look very similar to my heat plots from a few posts back. I removed the error bars in the boxes, as they were redundant and now visually distracting with the addition of the tails.

If we’re going to use this, I can’t fit just the box centers anymore.

Posted in NpT | Leave a comment