Density profiles (constant vs powerlaw)

This example illustrates the usage of the radexp paramter that allows for a powerlaw-like density profile. We compare radexp=-0.5 and radexp=0.0 in two otherwise identical runs.

First we let the density fall of as \(1/\sqrt{R}\).

xstar cfrac=1 temperature=1000 density=1.E+8  spectrum='pow ' trad=-1. rlrad38=1.E-10 column=1.E+17 rlogxi=5. lcpres=0 abundtbl='xdef' modelname='densexp-0.5' niter=99 npass=1 critf=1.E-07 nsteps=6 xeemin=0.04 emult=0.5 taumax=5.  ncn2=9999 radexp=-0.5 vturb=1.

Here the density is held constant.

xstar cfrac=1 temperature=1000 density=1.E+8  spectrum='pow ' trad=-1. rlrad38=1.E-10 column=1.E+17 rlogxi=5. lcpres=0 abundtbl='xdef' modelname='filled sphere' niter=99 npass=1 critf=1.E-07 nsteps=6 xeemin=0.04 emult=0.5 taumax=5.  ncn2=9999 radexp=0.0 vturb=1.

First we look at the density as a function of radius.

>>> import matplotlib.pyplot as plt
>>> from astropy.io import fits
>>> import numpy as np
>>>
>>> a1=fits.open('files/0.5/xout_abund1.fits')[1].data
>>> a2=fits.open('files/0.0/xout_abund1.fits')[1].data
>>>
>>> plt.plot(a1['radius'], a1['n_p'])
>>> plt.plot(a2['radius'], a2['n_p'])
>>> plt.xscale("log")
>>> plt.xlabel("Radius")
>>> plt.ylabel("Density")
>>>
>>> plt.show();

The inner radius is determined by the ionization parameter, density, and luminosity, in this example \(R_\text{in}=10^{17.5}\,\text{cm}\). The outer radius is determined by the column density, so the two clouds have different outer radii. As the ionization structure inside the cloud is also different.

>>> plt.plot(a1['radius'], a1['ion_parameter'])
>>> plt.plot(a2['radius'], a2['ion_parameter'])
>>> plt.xscale("log")
>>> plt.xlabel("Radius")
>>> plt.ylabel("logxi")
>>>
>>> plt.show();

The geometrical dilution of the radiation field is proportional to \(1/R^{2}\) while the density is \(1/\sqrt{R}\). The ionization parameter falls faster for the constant density case, but does not reach values below 2. Consequently, this leads to different charge state distributions. As an example, we compare the total ion columns for silicon.

Note, that the fully ionized ions are not listed in the abundance file. Their column/fraction can be calculated from the total colum density, the elemental abundance and the sum of all other ions.

>>> c1= fits.open('files/0.5/xout_abund1.fits')[2].data
>>> c2= fits.open('files/0.0/xout_abund1.fits')[2].data
>>> si_bare=1e17*3.5e-5-sum([c1.si_xiv[0], c1.si_xiii[0], c1.si_xii[0], c1.si_xi[0], c1.si_x[0], c1.si_ix[0], c1.si_viii[0], c1.si_vii[0], c1.si_vi[0], c1.si_v[0], c1.si_iv[0], c1.si_iii[0], c1.si_ii[0], c1.si_i[0]])
>>> labels = ['Si XIV', 'Si XIII', 'Si XII', 'Si XI', 'Si X', 'Si IX', 'Si VIII',  'Si VII', 'Si VI', 'Si V', 'Si IV', 'Si III', 'Si II', 'Si I', 'stripped']
>>> sizes = [c1.si_xiv[0], c1.si_xiii[0], c1.si_xii[0], c1.si_xi[0], c1.si_x[0], c1.si_ix[0], c1.si_viii[0], c1.si_vii[0], c1.si_vi[0], c1.si_v[0], c1.si_iv[0], c1.si_iii[0], c1.si_ii[0], c1.si_i[0], si_bare]
>>>
>>> fig, ax = plt.subplots()
>>> ax.pie(sizes, labels=labels, autopct='%1.1f%%')

Repeating the same exercise for the constant density case gives.