Many books on radiative transfer and spectroscopy introduce the partition function ($U$, PF) as a statistical weight or probability that an atom is in the ionization stage $r$ at the temperature $T$. It is defined as the sum over all excitation states ($i$) of their statistical weights ($g_i$) multiplied by an

exponential factor (cf. Rutten, 2003, p.30, Eq.289):

$$U_r \equiv \sum_i g_i \mathrm{e}^{-\frac{\chi_{r, i}}{kT}},$$
where $\chi_{r, i}$ is the excitation potential of the state $i$ in the ionization stage $r$, $k$ is the Boltzmann constant) and $T$ is the temperature.

The trouble with PF is that there is an infinite number of the energy levels. No matter how small their contribution becomes (the exponential factor is decreasing with the excitation potential), PF computed after this definition goes to infinity. Farther, that means that the probability of an atom to occupy any finite state would be zero. Since this is clearly unphysical, thus there is a need to define a reliable procedure to cut off the partition function. There are many approached described in literature. However, they are usually rather non-trivial and, therefore, computing partition function even for a limited range of temperature and for a specific atomic element can be a very tedious job. The results of these computations are often published as the coefficients of the best polynomial fits for a range of temperature values. The IDL procedure that you can find below is written to load these tables from several different sources (see below) and to interpolate them for given a temperature, an atomic specie and an ionization stage.

## Wednesday 25 December 2013

## Saturday 12 October 2013

### Use of "some useful atomic data"

To illustrate how the routines in the previous post can be used, I made a couple of plots.

Fig.1 The logarithmic solar abundances relative to hydrogen A (A(H) = 12) after Asplund et al (2009, 2009ARA&A..47..481A) versus the atomic number Z. The data are loaded using load_abundances.pro function.

Fig.1 The logarithmic solar abundances relative to hydrogen A (A(H) = 12) after Asplund et al (2009, 2009ARA&A..47..481A) versus the atomic number Z. The data are loaded using load_abundances.pro function.

Labels:
abundances,
data,
IDL,
sun

Location:
Santa Cruz de Tenerife, Spain

## Saturday 5 October 2013

### Some useful atomic data in IDL

Here I list a couple of my IDL routines that load some useful atomic data
(atomic numbers up to 92 are supported).

This routine returns symbols and names of atomic elements for a given atomic number Z.

load_list_of_elements.pro

Example:

This routine loads the atomic weights for a given atomic number Z. The data comes from

http://www.nist.gov/pml/data/comp.cfm

The original source is Wieser & Berglund (2009, Pure Appl. Chem., Vol.81, No.11, p.2131-2156), published as an IUPAC Tecnical Report available at:

http://pac.iupac.org/publications/pac/pdf/2009/pdf/8111x2131.pdf

load_atomic_weights.pro

Example:

**List of elements**This routine returns symbols and names of atomic elements for a given atomic number Z.

load_list_of_elements.pro

Example:

```
IDL> sym = LOAD_LIST_OF_ELEMENTS([6, 7, 8], element = names)
IDL> PRINT, sym
C N O
IDL> PRINT, names
Carbon Nitrogen Oxygen
```

**Atomic weights (aka relative atomic masses)**This routine loads the atomic weights for a given atomic number Z. The data comes from

http://www.nist.gov/pml/data/comp.cfm

The original source is Wieser & Berglund (2009, Pure Appl. Chem., Vol.81, No.11, p.2131-2156), published as an IUPAC Tecnical Report available at:

http://pac.iupac.org/publications/pac/pdf/2009/pdf/8111x2131.pdf

load_atomic_weights.pro

Example:

```
IDL> a = LOAD_ATOMIC_WEIGHTS([1, 2, 6, 7, 8])
IDL> PRINT, a
1.00790 4.00260 12.0107 14.0067 15.9994
```

Labels:
abundances,
data,
IDL,
tools

Location:
Santa Cruz de Tenerife, Spain

## Sunday 29 September 2013

### Cubic spline interpolation: Periodic splines

If the unknown function is periodic, then the periodicity may be used to specify the boundary conditions. Before (in the case of the natural splines) we had $n+1$ value of x and y and as a result we obtained values of the cubic polynomial at $n$ intervals.

Now our function is periodic with period $p$, thus $S(x) = S(x+p)$. In practice, we have an additional data point $x_{n+1}$, such that

$$x_{n+1} = x_0 + p,$$

and

$$y_{n+1} = y_0.$$

The total number of data points is $n+2$ and the total number of segments (and cubic splines) between them is $n+1$.

Now our function is periodic with period $p$, thus $S(x) = S(x+p)$. In practice, we have an additional data point $x_{n+1}$, such that

$$x_{n+1} = x_0 + p,$$

and

$$y_{n+1} = y_0.$$

The total number of data points is $n+2$ and the total number of segments (and cubic splines) between them is $n+1$.

Labels:
IDL,
interpolation,
numerics,
tools

Location:
Santa Cruz de Tenerife, Spain

## Thursday 26 September 2013

### Cubic spline interpolation: Natural splines

Last year I published a post on the data interpolation and smoothing using the cubic splines using IDL. (Soon I'll do an update of that post as well.) Here I go a step back and show how the classical spline interpolation works. As before, the derivation follows the book of D.S.G.Pollock (highly recommended to any student interested in the analysis of time-series). I'll first go through some theory general for all the cubic splines, than through the derivation of the

**natural**cubic splines and at the end I'll show some examples. In the next post I'll describe the solution of the**clamped**cubic splines and do the same examples using that variation of the method.
Labels:
IDL,
interpolation,
numerics,
tools

Location:
Santa Cruz de Tenerife, Spain

## Wednesday 21 August 2013

### Three-diagonal periodic system

Three-diagonal periodic system of linear equations can be solved in a similar way as the pure three-diagonal system. The difference is that now we have additional two coefficients in the matrix of the system ($a_0\ne0$, $c_n \ne0$). This form of a system is common when the boundary conditions are periodic (e.g. periodic cubic splines). The solution is again based on the Gaussian elimination, but because of the additional matrix elements, it requires some more steps. The system is strictly diagonal dominant. It can be proved that it has a unique solution.

There is n+1 equation. Our goal is to solve it for n+1 value of $x_i$. The matrix of the system (size $n+1\times n+1$) is of the following form (additional coefficients are shown in bold font):

\begin{equation}\left(\begin{matrix} b_0 & c_0 & 0 & \dots & 0 & 0 & \mathbf{a_0} \\

a_1 & b_1 & c_1 & \dots & 0 & 0 & 0\\

\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\

0 & 0 & 0 & \dots & a_{n-1} & b_{n-1} & c_{n-1} \\

\mathbf{c_n} & 0 & 0 & \dots & 0 & a_{n} & b_n \end{matrix}\right)

\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \\x_{n}\end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1} \\d_{n}\end{matrix}\right)

\label{eq:matrix}\nonumber\end{equation}

or in the form of the individual equations ($n+1$ equation with $n+1$ unknown $x_i$):

\begin{eqnarray}

b_0 x_0 + c_0 x_1 + a_0 x_n &=& d_0,\nonumber\\

a_1 x_0 + b_1 x_1 + c_1 x_2 &=& d_1,\nonumber\\

\dots &=& \dots, \nonumber\\

a_i x_{i-1}+ b_i x_i + c_i x_{i+1} &=& d_i,\nonumber\\

\dots &=& \dots ,\nonumber\\

a_{n-1} x_{n-2} +b_{n-1} x_{n-1} + c_{n-1} x_{n} &=& d_{n-1},\nonumber\\

c_n x_{0} + a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber

\end{eqnarray}

There is n+1 equation. Our goal is to solve it for n+1 value of $x_i$. The matrix of the system (size $n+1\times n+1$) is of the following form (additional coefficients are shown in bold font):

\begin{equation}\left(\begin{matrix} b_0 & c_0 & 0 & \dots & 0 & 0 & \mathbf{a_0} \\

a_1 & b_1 & c_1 & \dots & 0 & 0 & 0\\

\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\

0 & 0 & 0 & \dots & a_{n-1} & b_{n-1} & c_{n-1} \\

\mathbf{c_n} & 0 & 0 & \dots & 0 & a_{n} & b_n \end{matrix}\right)

\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \\x_{n}\end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1} \\d_{n}\end{matrix}\right)

\label{eq:matrix}\nonumber\end{equation}

or in the form of the individual equations ($n+1$ equation with $n+1$ unknown $x_i$):

\begin{eqnarray}

b_0 x_0 + c_0 x_1 + a_0 x_n &=& d_0,\nonumber\\

a_1 x_0 + b_1 x_1 + c_1 x_2 &=& d_1,\nonumber\\

\dots &=& \dots, \nonumber\\

a_i x_{i-1}+ b_i x_i + c_i x_{i+1} &=& d_i,\nonumber\\

\dots &=& \dots ,\nonumber\\

a_{n-1} x_{n-2} +b_{n-1} x_{n-1} + c_{n-1} x_{n} &=& d_{n-1},\nonumber\\

c_n x_{0} + a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber

\end{eqnarray}

Location:
Santa Cruz de Tenerife, Spain

## Saturday 17 August 2013

### Three-diagonal system of linear equations

In this post I'll show a fast and accurate solution for the three-diagonal system of linear equations based on Gaussian elimination. In the matrix notation the system has the following form:

\begin{equation}\left(\begin{matrix} b_0 & c_0 & 0 & \dots & 0 & 0 & 0 \\

a_1 & b_1 & c_1 & \dots & 0 & 0 & 0\\

\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\

0 & 0 & 0 & \dots & a_{n-1} & b_{n-1} & c_{n-1} \\

0 & 0 & 0 & \dots & 0 & a_{n} & b_n \end{matrix}\right)

\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \\x_{n}\end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1} \\d_{n}\end{matrix}\right)

\label{eq:matrix00}\nonumber\end{equation}

or in the form of the individual equations ($n+1$ equation with $n+1$ unknown $x_i$):

\begin{eqnarray}

b_0 x_0 + c_0 x_1 &=& d_0,\nonumber\\

a_1 x_0 + b_1 x_1 + c_1 x_2 &=& d_1,\nonumber\\

\dots &=& \dots, \nonumber\\

a_i x_{i-1}+ b_i x_i + c_i x_{i+1} &=& d_i,\nonumber\\

\dots &=& \dots ,\nonumber\\

a_{n-1} x_{n-2} +b_{n-1} x_{n-1} + c_{n-1} x_{n} &=& d_{n-1},\nonumber\\

a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber

\end{eqnarray}

\begin{equation}\left(\begin{matrix} b_0 & c_0 & 0 & \dots & 0 & 0 & 0 \\

a_1 & b_1 & c_1 & \dots & 0 & 0 & 0\\

\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\

0 & 0 & 0 & \dots & a_{n-1} & b_{n-1} & c_{n-1} \\

0 & 0 & 0 & \dots & 0 & a_{n} & b_n \end{matrix}\right)

\left(\begin{matrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \\x_{n}\end{matrix}\right) = \left(\begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_{n-1} \\d_{n}\end{matrix}\right)

\label{eq:matrix00}\nonumber\end{equation}

or in the form of the individual equations ($n+1$ equation with $n+1$ unknown $x_i$):

\begin{eqnarray}

b_0 x_0 + c_0 x_1 &=& d_0,\nonumber\\

a_1 x_0 + b_1 x_1 + c_1 x_2 &=& d_1,\nonumber\\

\dots &=& \dots, \nonumber\\

a_i x_{i-1}+ b_i x_i + c_i x_{i+1} &=& d_i,\nonumber\\

\dots &=& \dots ,\nonumber\\

a_{n-1} x_{n-2} +b_{n-1} x_{n-1} + c_{n-1} x_{n} &=& d_{n-1},\nonumber\\

a_{n} x_{n-1} +b_{n} x_{n} &=& d_{n}.\nonumber

\end{eqnarray}

Location:
Santa Cruz de Tenerife, Spain

## Tuesday 16 July 2013

### "The Sun is miasma of incandescent plasma"

Labels:
science+art,
sun

Location:
Santa Cruz de Tenerife, Spain

## Sunday 2 June 2013

### Kourganoff graphs

Scharzschild equation provides a tool to easily find and analyze the solution of the radiative transfer cases for some simple source functions. First of all, the lambda operator is clearly linear, so\begin{equation}

\Lambda[a\,f(t) + b\,g(t)] = a\,\Lambda[f(t)] + b\,\Lambda[g(t)].

\end{equation}

\Lambda[a\,f(t) + b\,g(t)] = a\,\Lambda[f(t)] + b\,\Lambda[g(t)].

\end{equation}

Labels:
IDL,
lambda operator,
numerics,
radiative transfer

## Saturday 1 June 2013

### Formal solution and Schwarzschild equation

The radiation transport equation (RTE) for time-independent plane-parallel atmosphere has well known form (e.g. Rutten, 2003):

\begin{equation}

\mu\,\frac{\mathrm{d}I_{\nu}}{\mathrm{d} \tau_{\nu}} = I_{\nu} - S_{\nu},

\end{equation}

where $I_{\nu}$ is specific intensity, $S_{\nu}$ source function, $\tau_{\nu}$ is the radial optical depth and $\mu$ is the cosine of the angle between the line of sight and vertical axis. Mathematically, this is a linear first order differential equation with constant coefficients. Physically, it describes propagation of radiation through an interacting medium. Properties of the medium are hidden in the source function and the optical depth. However, to make the first step toward the solution of RTE, we will assume that the optical depth is independent variable and that the source function is a know function of it.

\begin{equation}

\mu\,\frac{\mathrm{d}I_{\nu}}{\mathrm{d} \tau_{\nu}} = I_{\nu} - S_{\nu},

\end{equation}

where $I_{\nu}$ is specific intensity, $S_{\nu}$ source function, $\tau_{\nu}$ is the radial optical depth and $\mu$ is the cosine of the angle between the line of sight and vertical axis. Mathematically, this is a linear first order differential equation with constant coefficients. Physically, it describes propagation of radiation through an interacting medium. Properties of the medium are hidden in the source function and the optical depth. However, to make the first step toward the solution of RTE, we will assume that the optical depth is independent variable and that the source function is a know function of it.

Labels:
integration,
lambda operator,
radiative transfer,
tools

Location:
Stari Grad, Serbia

## Sunday 7 April 2013

### Gauss quadrature and Legendre polynomials

Numerical integration is among the most common tasks in astrophysics. Simple formulae, like trapezoidal, Newton-Cotes or Simpson, are often not enough accurate. Gauss quadrature provides a more accurate solution, but its implementation is a bit more difficult. Here I will demonstrate the implementation of that algorithm in IDL. The algorithm includes finding the coefficients of Legendre polynomials and their zeros. Both tasks are explicitly solved in the following as well.

But let's go step by step. Our task is to solve numerically the integral:

$$I = \int_a^b f(x)\, \mathrm{d} x.$$

But let's go step by step. Our task is to solve numerically the integral:

$$I = \int_a^b f(x)\, \mathrm{d} x.$$

## Sunday 10 March 2013

### Exponential integrals

The formal solution of the radiative transfer equation can be expressed using exponential integrals. These are mathematical functions defined as:

\begin{equation}

E_n(x) = \int_1^{\infty}\, \exp (-yx)\frac{\mathrm{d}x}{y^n}, \;\;\;\;n = 1, 2, ...

\end{equation}

There is several important properties of these functions:

\begin{equation}

E_n(x) = \int_1^{\infty}\, \exp (-yx)\frac{\mathrm{d}x}{y^n}, \;\;\;\;n = 1, 2, ...

\end{equation}

There is several important properties of these functions:

## Thursday 21 February 2013

### Spinoza's tomb

The human Mind cannot be absolutely destroyed with the Body, but something of it remains which is eternal. (from

*Ethica*, Part V, proposition 23)Baruch (Benedict de) Spinoza

(24 November 1632 - 21 February 1677)

Dutch philosopher, born in Jodenbuurt in Amsterdam, buried in the churchyard of Nieuwe Kerk in Den Haag.

(Photo was taken on a sunny afternoon in May, 2010.)

Labels:
philosophy,
science+history,
spinoza

Location:
Zuidwal, The Hague, The Netherlands

Subscribe to:
Posts (Atom)