K-means clustering

The problem of clustering is a rather general one: If one has $m$ observations or measurements in $n$ dimensional space, how to identify $k$ clusters (classes, groups, types) of measurements and their centroids (representatives)?

The k-means method is extremely simple, rather robust and widely used in it numerous variants. It is essentially very similar (but not identical) to Lloyd's algorithm (aka Voronoi relaxation or interpolation used in computer sciences).

k-means

Let's use the following indices: $i$ counts measurements, $i \in [0, m-1]$; $j$ counts dimensions, $j \in [0, n-1]$; $l$ counts clusters, $l \in [0, k-1]$.

Each measurement in $n$-dimensional space is represented by a vector $x_i = \{x_{i, 0}, \dots x_{i, n-1}\}$, where index $i$ is counting different measurements ($i = 0, \dots, m-1$). The algorithm can be summarized as:

1. Choose randomly $k$ measurements as initial cluster centers: $c_0, ..., c_{k-1}$. Obviously, each of the clusters is also $n$-dimensional vector.

2. Compute Euclidean distance $D_{i, l}$ between every measurement $x_i$ and every cluster center $c_l$:
$$D_{i, l} = \sqrt{\sum_{j=0}^{n-1} (x_{i, j} - c_{l, j})^2}.$$
3. Assign every measurement $x_i$ to the cluster represented by the closest cluster center $c_l$.

4. Now compute new cluster centers by simply averaging all the measurements in each cluster.

5. Go back to 2. and keep iterating until none of the measurements changes its cluster in two successive iterations.

This procedure is initiated randomly and the result will be slightly different in every run. The result of clustering (and the actual number of necessary iteration) significantly depends on the initial choice of cluster centers. The easiest way to improve the algorithm is to improve the initial choice, i.e. to alter only the step 1. and then to iterate as before. There are to simple alternatives for the initialization.

First observation of linear polarization in the forbidden [OI] 630.03 nm line

In a new paper (de Wijn, Socas-Navarro & Vitas, 2017, ApJ, 836, 29D) we present the first results of our observations of a sunspot and an active region using the SP/SOT instrument on board the Hinode satellite. The novelty in our observation is a trick that we used to double the standard wavelength range observed by the instrument. Thanks to that, we were able to see the sun not only in the two iron lines at 630.2 nm, but also in four other lines. One of those is particularly interesting: the forbidden ground-based line of neutral oxygen ([OI] 630.03 nm). It is one of only few oxygen lines in the solar spectrum and probably the best diagnostics of the solar oxygen abundance. For the first time ever we observed the linear polarization in this line! As an M2 (magnetic dipole) transition, it is predicted by the theory (Landi degl'Innocenti and Landi, 2004, Section 6.8) that this line produces the linear polarization signal with the opposite sign to the lines produced by E2 transitions. It is also the first time that linear polarization in M2 and nearby E2 lines is measured simultaneously, so that the flip in sign is obvious (see the left-most spectral line in the red circle in the Figure; in linear polarization it has "W" shape, while all other lines in the wavelength range have "M" shapes). This result may bring new light to the ongoing debate on the solar oxygen crisis.

  

More details of this unique observation will appear soon in a follow-up publication.

1D Solar Atmosphere Models in IDL: Penumbra by Ding & Fang (1989)

Plane-parallel atmosphere in hydrostatic equilibrium published by Ding & Fang (1989, "A semi-empirical model of sunspots penumbra", 1989A&A...225..204D). Statistical equilibrium for hydrogen model-atom with 12 levels plus continuum. The model is produced by fitting observations of penumbra in  2 lines of H and 5 lines of Ca. The observations were carried out on McMath telescope at Kitt Peak National Observatory. The observed sunspot was small, rounded and close to the disk center. The field strength in the umbra was around 1.25 kG and 560 G in the penumbra.

It is interesting to note their Fig.2 (see it below). In the deep photosphere, the temperature in this model is similar to the temperature in the model of Yun et al. (1981). However, between the optical depths -3 and -4 it becomes close to the VALC model (Vernazza et al, 1981).

Installing IDL v7.1 on Kubuntu 16.04 (on MacBook pro)

The IDL installation went smooth as usual, but when I tried to run IDL, there was an error:

error while loading shared libraries: libXp.so.6: cannot open shared object file: No such file or directory 

In the IDL help pages ( IDL fails to install on Linux: What to do, scroll down to the Ubuntu section) there is a comment on that that turned only half useful. Two libraries I got installed with no problem:

sudo apt-get install libxmu-dev
sudo apt-get install libxmu6

but the other three were not found in repositories. Instead of manual installation (two of them are available from http://packages.ubuntu.com/xenial/, but the last one is a bit mysterious. Google finds it only in three pages related to IDL, one of them being the previously mentioned page from the IDL docs). Instead, I got if from Git. As I have just installed Kubuntu there was still some important packages missing, so some of the following steps may be redundant for you.

sudo apt-get install autoconf autogen intltool
sudo apt-get install git
sudo apt-get install xutils-dev libtool libx11-dev 
                          x11proto-xext-dev x11proto-print-dev
git clone https://cgit.freedesktop.org/xorg/lib/libXp/

Then cd to libXp directory and execute

sudo ./autogen.sh

In my case if was complaining about the line 18214 related to XPRINT. I commented out that line from the script and executed it again. After that everything is straight forward:

sudo ./configure
sudo make install

And finally I add the path to the library to .bashrc (all in one line):

echo 'export LD_LIBRARY_PATH="/usr/local/lib/:$LD_LIBRARY_PATH"' >> ~/.bashrc

After that IDL worked normally. 

How to find if point is in or out of polygon? (in IDL)

It's a common problem in visualization and many other applications. There are many ways to solve it, but not all the ways are equally efficient. Baard Krane wrote a very elegant solution in IDL that was later vectorized by William Connelly and popularized in the community - of course - by David Fanning (see his article here and get the source of the inside.pro).

Example:

IDL> px = [0.3, 0.3, 0.7, 0.7]
IDL> py = [0.3, 0.7, 0.7, 0.3]
IDL> x = RANDOMU(seed1, 5000) 
IDL> y = RANDOMU(seed2, 5000) 
IDL> id = INSIDE(x, y, px, py, /ind)          
IDL> PLOT, x, y, psym = 3                     
IDL> OPLOT, [px, 0.3], [py, 0.3]              
IDL> OPLOT, x[id], y[id], psym = 3, col = 255

Equation of state: Vardya - Mihalas - Wittmann

These are my notes on equation of state derived initially by Vardya (1965), described by Mihalas (1967) and popularized by Wittmann (1974). It is widely used (or at least present as an option) in spectral synthesis codes like SIR or NICOLE. It is also prepared for the MANCHA code with several modifications and additionally computed quantities. However, the equation is derived in this particular form to be solvable on computing resources half a century ago. From today's perspective, this formulation is rather obsolete. While the results of the VMW EOS are still largely reliable, various tricks introduced to control numerical stability limit its usability to rather restricted range of the pressure and temperature.

Here I derive a simple equation of state for the solar atmosphere following the classical work of Vardya (1965), Mihalas (1967) and Wittmann (1974). The EOS is based on the Saha ionisation equilibrium and the instanteneous chemical equilibrium for the molecules. The main ingredient is hydrogen. It's included as atomic hydrogen (H), negative hydrogen ion (H-), positive hydrogen ion (H+), and as H2 and H2+ molecules. For all other atoms the neutral and the first two ionisation stages are included.

The equations were first published by Vardya (1967). Mihalas (1967) gave a simple numerical algorithm for an efficient solution of the system. Wittmann (1974) copied the equations and add a corrective factor that provides numerical stability at high temperature.

The derivation here follows Mihalas. However, in the original derivation there is a couple of inconsistencies that obscure the procedure. Here I write the equations in a correct and consistent way.

Definitions

Let's first define the pressures:
$p_{\mathrm{H}}$ - partial pressure of the neutral H atoms;
$p_{\mathrm{H^+}}$ - partial pressure of the positive H ions (protons);

Effective Lande g-factor (in IDL)

The effective Lande g-factor (Shenstone and Blaire, 1929) is defined for a spectral line. It shows how much the $\sigma$ components of the normal Zeeman triplet are separated from the line center.  The higher the Lande factor is, the line is more sensitive to the magnetic field. It is computed as a combination of g-factors for each of the involved atomic levels. Derivation of the equation can be found, for example, in Landi Degl'Innocenti (1982SoPh...77..285L, see it for more details and for an alternative expression useful in the polarized radiative transfer).

For each of the levels (u for upper, l for lower) in LS coupling holds:
$$g_\mathrm{LS} =\frac{3}{2}+\frac{S (S+1) - L (L+1)}{2J(J+1)}, $$ 
where $L$, $S$ and $J$ are the orbitatl, spin and total angular momentum quantum numbers of the atomic level (note that it is undefined for $J=0$). The effective Lande g-factor is then defined as:
$$\bar{g} = \frac{1}{2}(g_\mathrm{l} +g_\mathrm{u}) + \frac{1}{4}(g_\mathrm{l} - g_\mathrm{u})(J_\mathrm{l}(J_\mathrm{l}+1) - J_\mathrm{u}(J_\mathrm{u}+1)).$$

Here is my IDL code for computing the effective Lande g-factor of a spectral line in LS coupling. At the input the code takes the quantum number of the lower and the upper level of a transition. The numbers may be specified either as an array [L, S, J] or in spectroscopic notation.

Example:

IDL> PRINT, LANDE_FACTOR(lower  = [3., 1, 2.], upper = '3p1.0')
  Upper term: 3p1.0 => l = 1.00000s = 1.00000j = 1.00000
          0.250000

IDL> PRINT, LANDE_FACTOR(lower  = '5F1.0', upper = '5d0.0')
  Lower term: 5f1.0 => l = 3.00000s = 2.00000j = 1.00000
  Upper term: 5d0.0 => l = 2.00000s = 2.00000j = 0.00000
          0.00000

Download: lande_factor.pro