tag:blogger.com,1999:blog-11690511877546145392024-03-23T06:24:14.401-07:00Nikola VitasNikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.comBlogger56125tag:blogger.com,1999:blog-1169051187754614539.post-59631434334332839082023-06-07T07:58:00.002-07:002023-06-07T07:58:14.085-07:00Opacity for realistic 3D MHD simulations of cool stellar atmospheres<p>The first paper of <b>Andrea Perdomo Garcia</b> is just submitted for publication in <i>Astronomy & Astrophysics</i>, and out on <a href="https://arxiv.org/abs/2306.03744">arxiv.org/abs/2306.03744</a>. The paper is all about computing the opacities for realistic modelling of cool stellar atmospheres. It is divided in three unities. First (Section 3) it describes the computation of detailed <i>monochromatic opacity</i> including millions of atomic and molecular spectral lines and millions of wavelength points. For this the code <a href="http://tlusty.oca.eu/Synspec49/synspec.html">SYNSPEC</a> (Hubeny and Lanz, 2011, 2017a, b) is used. Then (Section 4) the monochromatic opacities are used to construct <i>opacity distribution function</i> which reduces the number of wavelength points from millions to thousands. The results are compared in detail with ones produced by Kurucz. Some striking similarities and some warning differences are found. Finally (Section 5), the opacity distribution function to construct <i>opacity bins</i>. This method, originally proposed by Nordlund (1982) is the key ingredient for realistically simulating stellar atmospheres in 3D as it reduced the problem further, from thousands of wavelength points to only a few. However, the method depends on a choice of some free parameters. In our paper the possible choices are carefully analyzed and some interesting conclusions are offered. </p><p>In Sect.3 there are two figures (Figs.2 and 3) that I find very useful and illustrative. The monochromatic opacity (Fig.2) and the radiative heating rate (Fig.3) are shown as 2D functions of wavelength (X-axis) and height in the atmosphere (Y-axis) for four different cool stars (all with solar metalicity). Optical depths in the continuum and continuum+lines are overplotted.
</p><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgENuzI3McBm6k4RtHWURpiDc9SxgmFKRCZKpkMvCYkdpXcG9X9U0JBiR9WKPyFeTQgr9Oejqhfs1hfjRN2DqKltyv073zBekiNINN50mZPDe6FxHsEW-ii-SzR8vU6gJBS5pQQYMVtwEBvh-YsCJWjPKllf4FVrtW-e4Z_Atprv9CiXdC0u9IYFI_i/s918/fig2_apg2023.png" style="display: block; padding: 1em 0px; text-align: center;"><img border="0" data-original-height="733" data-original-width="918" height="511" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgENuzI3McBm6k4RtHWURpiDc9SxgmFKRCZKpkMvCYkdpXcG9X9U0JBiR9WKPyFeTQgr9Oejqhfs1hfjRN2DqKltyv073zBekiNINN50mZPDe6FxHsEW-ii-SzR8vU6gJBS5pQQYMVtwEBvh-YsCJWjPKllf4FVrtW-e4Z_Atprv9CiXdC0u9IYFI_i/w640-h511/fig2_apg2023.png" width="640" /></a></div><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXBXgBDtEauosL9XeOgmLdeeIKUYBnBp8wapNvjMvaHTw133gusKBZTft2e-MEvPrmCYdmv927uO8xQHfiVeNhUvWn-TanOKBFenZYou5-Z2EkgvmZ1Qa-DBicVxSeWd8JUoitIEf5YkDql1iD2rFN6Vrn4jZWrsYzGxlLfpdN10MFpjB0dHraY11q/s918/fig3_apg2023.png" style="display: block; padding: 1em 0px; text-align: center;"><img border="0" data-original-height="685" data-original-width="918" height="477" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiXBXgBDtEauosL9XeOgmLdeeIKUYBnBp8wapNvjMvaHTw133gusKBZTft2e-MEvPrmCYdmv927uO8xQHfiVeNhUvWn-TanOKBFenZYou5-Z2EkgvmZ1Qa-DBicVxSeWd8JUoitIEf5YkDql1iD2rFN6Vrn4jZWrsYzGxlLfpdN10MFpjB0dHraY11q/w640-h477/fig3_apg2023.png" width="640" /></a></div>
(Andrea is the final year PhD student at Instituto de Astrofisica de Canarias and Univeridad de La Laguna, supervised by Manolo Collados Vera and myself. Stay tuned, more cool stuff is coming out from her research this year.)Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184670.15339576382115538 -51.4080967 56.773863436178843 18.9044033tag:blogger.com,1999:blog-1169051187754614539.post-18425514912677816402023-03-25T02:46:00.000-07:002023-11-24T02:06:00.574-08:00Charles Hermite (1822 - 1901)<p>Sunday morning in Paris was opportunity to walk to the Montparnasse cemetery and pay respect to some of my personal heroes buried there. While the graves of Beckett, Cortázar and Poincaré attract quite many attention and visitors, it is less known that the great French mathematician <a href="https://en.wikipedia.org/wiki/Charles_Hermite">Charles Hermite</a> is buried there as well. Not only that he was a mentor to Henri Poincaré, Henri Padé, Thomas Stieltjes and Mihajlo Petrovic Alas, but his work on interpolation and function approximation is at the very core of the modern numerical methods used in computational fluid dynamics and radiative transfer (even when this is not so obvious or properly acknowledged). His work from 1878 ("Sur la formule d'interpolation de Lagrange") should be read by anyone interested in function approximation. The name on his grave stone are barely readable nowadays. <a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiwvJp9Ef7KOH2F5QfIYADBFhfy8eXTzszDZqqeuB83t6SQlxod1bv4lO7LaPV56KulxzJmaiX7tcAou5qA0OyimK-gjo1z6D1as5vediDAMmoLAev8lOfUqwZlIXcMCIIsOBmO1g6VAxCnoOkZOqwP1kJF4IraHzcOM2YJeIdH_Wv60InsW8JKtzkyoF8/s4160/hermite_grave.jpg" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="4160" data-original-width="3120" height="600" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiwvJp9Ef7KOH2F5QfIYADBFhfy8eXTzszDZqqeuB83t6SQlxod1bv4lO7LaPV56KulxzJmaiX7tcAou5qA0OyimK-gjo1z6D1as5vediDAMmoLAev8lOfUqwZlIXcMCIIsOBmO1g6VAxCnoOkZOqwP1kJF4IraHzcOM2YJeIdH_Wv60InsW8JKtzkyoF8/s600/hermite_grave.jpg" /></a><img alt="" border="0" data-original-height="3120" data-original-width="4160" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhbD27YiMXUvoeG5n2cbZZhPlh5H8MmVOPcHtROxngVMDZXGerFQOcf7pLa3RWqlrva1WhEWlFXCemKyxVqF1aT4Bm-p2PGmfAkdWbhCdn7RB5eRn_12uNxLM8KncJZuApS14l-M2Fee8ucSnjGhKhLA9KQAbwzqfFL2FiekoF2qRmzC_Nz7f8RHz93zR4/s600/hermite_grave2.jpg" width="600" /></p>Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Montparnasse Cemetery, 11 Av. Transversale, 75014 Paris, France48.83846 2.325650248.837753838825407 2.3245773163940431 48.839166161174589 2.3267230836059571tag:blogger.com,1999:blog-1169051187754614539.post-35483625427669826142022-10-03T05:01:00.010-07:002023-09-25T05:04:48.584-07:00Robert Jelle Rutten (1942 - 2022)So long, my teacher and friend. You will be hugely missed.
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhK1sJJyL_RIBHaD0-8-uJcIkr5UkdOy4EsGx7Z8Xdvkc1RKzPpTxQy2zz9nmaLkYKQ8sGrY4msAq_bNhjGDr0SGpksToBxUI1CS6tebzIO6VbJw3_aZDG3tOAP7ch1fJ8zkdhTb0x0MsKrU7pZv1V-ju-oN556-ufwf7WdOxdeJ0DU4pHBThrmZH_VYPc/s1280/robrutten.jpg" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="600" data-original-height="853" data-original-width="1280" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhK1sJJyL_RIBHaD0-8-uJcIkr5UkdOy4EsGx7Z8Xdvkc1RKzPpTxQy2zz9nmaLkYKQ8sGrY4msAq_bNhjGDr0SGpksToBxUI1CS6tebzIO6VbJw3_aZDG3tOAP7ch1fJ8zkdhTb0x0MsKrU7pZv1V-ju-oN556-ufwf7WdOxdeJ0DU4pHBThrmZH_VYPc/s600/robrutten.jpg"/></a></div>
Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.comSanta Cruz de Tenerife, Spain28.4636296 -16.251846713.025238857272351 -33.8299717 43.902020342727653 1.3262782999999985tag:blogger.com,1999:blog-1169051187754614539.post-58334334433489370752022-09-07T08:07:00.004-07:002022-09-07T08:07:35.009-07:001D Solar Atmosphere Models in IDL: Standard Solar Model (from "The Sun", by Stix)Stix, in his book ("<a href="https://link.springer.com/book/10.1007/978-3-642-56042-2" target="_blank">The Sun: An Introduction</a>", Springer, 2002) in Sect.2.4 introduces a standard solar model from the core to the surface defined as $\tau = 2/3$. His Table 2.4 (on p.56) lists the values of various quantities of this model versus the column mass or the height. Check the book for the details of the model. The table in ascii format is here: <a href="https://drive.google.com/file/d/1F_To15ExxaoJPiRUm-qrjkS7lkozopcN/view?usp=sharing">stix_tab.2.4.txt</a>
The columns are: $m/m_\odot$, $r/r_\odot$, $p \mathrm[Pa]$, $T \mathrm[K]$, $\rho \mathrm[kg/m^3]$, $L/L_\odot$, $X$, $\mu$, $\Gamma_1$.
Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0tag:blogger.com,1999:blog-1169051187754614539.post-46443843554470131982022-07-13T09:18:00.003-07:002022-07-13T09:20:29.092-07:00NGC 7319: GTC vs JWST<div class="flourish-embed flourish-photo-slider" data-src="visualisation/10624332"><script src="https://public.flourish.studio/resources/embed.js"></script></div>
The first pictures from James Webb Space Telescope are simply mind blowing! Here is a comparison of NGC 7319 (from Stephan's quintet of galaxies) from <a href="https://en.wikipedia.org/wiki/Gran_Telescopio_Canarias" target="_blank">Gran Telescopio Canarias</a> (still the largest single aperture optical telescope) and from JWST. One should keep in mind that these are different instruments observing and different wavelengths. The alignment is also not perfect. Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0tag:blogger.com,1999:blog-1169051187754614539.post-29586394678791892492020-04-01T02:38:00.009-07:002023-04-12T15:56:27.406-07:00Quadratic splines: Part I<p>Some basic properties of parabola in the context of the second-order Hermite interpolation equation. <br /></p><p>There is a one-dimensional dataset with $n$ points: $x = (x_0, \dots, x_{n-1})$ and $y = (y_0, \dots, y_{n-1})$, such that $x_{i} > x_{i-1}$ for any $i$. We introduce the following labels:
$$h_{i-1} = x_i - x_{i-1}$$
$$\delta_{i-1} = \frac{y_i - y_{i-1}}{x_i - x_{i-1}}$$
$h$ is always positive, while $\delta$ can have any value.
On each segment, we approximate the unknown function with a parabola. On the $[i-1, i]$ segment:
\begin{equation}f(x) = a_0 + a_1 (x - x_{i-1}) + a_2 (x - x_{i-1})^2\end{equation}
The first and the second derivative of the parabola are:
\begin{equation}f'(x) = a_1 + 2 a_2 (x - x_{i-1})\end{equation}
\begin{equation}f''(x) = 2 a_2 \end{equation}
The Hermite interpolation theorem provides us with a tool to determine the $a$ coefficients on the $[i-1, i]$ segment if we know $y_{i-1}$, $y_i$ and one of the first derivative in the end points, either $y'_{i-1}$ or $y'_i$. Let's assume that we know ("know" in the sense "we can compute") $y'_{i-1}$. Then there are three equations that the parabola must obey:
\begin{equation}f(x_{i-1}) = y_{i-1} = a_0\end{equation}
\begin{equation}f'(x_{i-1}) = y'_{i-1} = a_1\end{equation}
\begin{equation}f(x_{i}) = y_i = a_0 + a_1 h_{i-1} + a_2 h_{i-1}^2\end{equation}
The solution is then trivial:
\begin{equation}a_0 = y_{i-1}\end{equation}
\begin{equation}a_1 = y'_{i-1}\end{equation}
\begin{equation}a_2 = \frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2}\end{equation}
Once the $a$ coefficients are known, the parabola is fully determined. That also means that the derivative in the point $i$ is fixed:
$$f'(x_{i}) = y'_i = a_1 + 2 a_2 (x_i - x_{i-1}) = y'_{i-1} + 2 \frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2} h_{i-1}$$
$$y'_i = \frac{2(y_{i} - y_{i-1}) - y'_{i-1} h_{i-1}}{h_{i-1}} $$
\begin{equation}y'_i = 2\delta_{i-1} - y'_{i-1}\end{equation}
<b>Local extremum</b><br /><br />
The extremum of this parabola is located at $x_m$:
$$f'(x_m) = 0 = a_1 + 2 a_2 (x_m - x_{i-1})$$
\begin{equation}x_m = x_{i-1} - \frac{a_1}{2a_2}\end{equation}
<b>Monotonicity</b><br /><br />
The function $f(x)$ is monotone on the $[i-1, i]$ segment if it does not have an extremum on it, i.e. if
$$x_m < x_{i-1}\;\;\;\;\;\;\;\;\;\;\mathrm{or}\;\;\;\;\;\;\;\;\;\; x_m > x_i$$
The first condition readily translates into:
$$x_{i-1} - \frac{a_1}{2a_2} < x_{i-1}$$
$$ \frac{a_1}{2a_2} > 0$$
or, when we substitute $a_1$ and $a_2$:
$$ \frac{y'_{i-1} h_{i-1}^2}{2(y_{i} - y_{i-1} - y'_{i-1} h_{i-1})} > 0$$
$$ \frac{y'_{i-1} h_{i-1}}{2(\delta_{i-1} - y'_{i-1})} > 0$$
As by definition $h_{i-1} > 0$
$$ \frac{y'_{i-1}}{2(\delta_{i-1} - y'_{i-1})} > 0$$
If $y'_{i-1} > 0$ then $\delta_{i-1} - y'_{i-1} > 0$ and $\delta_{i-1} > y'_{i-1}$. If $y'_{i-1} < 0$, then $\delta_{i-1} < y'_{i-1}$. Therefore, the first condition ($x_m < x_{i-1}$) reduces to:
\begin{equation}\boxed{|\delta_{i-1}| > |y'_{i-1}|}\end{equation}
The second condition leads to:
$$x_{i-1} - \frac{a_1}{2a_2} > x_{i}$$
$$-x_{i} + x_{i-1} - \frac{a_1}{2a_2} > 0$$
$$h_{i-1} + \frac{a_1}{2a_2} < 0$$
$$h_{i-1} + \frac{y'_{i-1} h_{i-1}}{2(\delta_{i-1} - y'_{i-1})} < 0$$
$$\frac{2\delta_{i-1} - y'_{i-1}}{2(\delta_{i-1} - y'_{i-1})} < 0$$
There are two cases. First, if $2\delta_{i-1} - y'_{i-1} > 0$, i.e. $2\delta_{i-1} > y'_{i-1}$. then it must be $\delta_{i-1} - y'_{i-1}< 0$, i.e.$\delta_{i-1} < y'_{i-1}$. Therefore, the first case leads to $2\delta_{i-1} > y'_{i-1} > \delta_{i-1}$. Secondly, if $2\delta_{i-1} - y'_{i-1} < 0$, i.e. $2\delta_{i-1} < y'_{i-1}$, it must be $\delta_{i-1} - y'_{i-1}> 0$, i.e. $\delta_{i-1} > y'_{i-1}$. So the second case leades to $\delta_{i-1} > y'_{i-1}> 2\delta_{i-1}$. The first case is possible only when both $\delta_{i-1}$ and $y'_{i-1}$ are positive and the second case is possible only when both are negative. Therefore, together, the two cases reduce the second condition ($x_m > x_{i}$) to:
\begin{equation}\boxed{|2\delta_{i-1}| > |y'_{i-1}| > |\delta_{i-1}|}\end{equation}
Together, the two conditions combine into
\begin{equation}\boxed{|2\delta_{i-1}| > |y'_{i-1}|}\end{equation}
A much faster way to reach the same result is to realize that the the derivatives $y'_{i-1}$ and $y'_{i}$ must be of the same sign if the function is monotone on the $[i-1, i]$ segment, i.e.
$$y'_{i-1} \cdot (2\delta_{i-1} - y'_{i-1}) > 0$$
which obviously leads to $|2\delta_{i-1}| > |y'_{i-1}|$.<br /><br />
<i>Therefore, for a Hermite parabola defined by $y_{i-1}$, $y_i$ and $y'_{i-1}$ to be monotone on the interval $(i-1, i)$, the derivative $y'_{i-1}$ must be of the same sign as the linear slope $\delta_{i-1}$ and smaller in magnitude than $2\delta_{i-1}$.</i><br /><br />
<b>Convexity</b><br /><br />
For the function $f$ to be convex on $(i-1, i)$ (if parabola is convex on a segment, it's convex everywhere) the condition is $f''(x) > 0$. In the case of our parabola:
$$f''(x) = 2a_2 > 0$$.
$$\frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2} > 0$$
or
$$y_{i} - y_{i-1} - y'_{i-1} h_{i-1} > 0$$
\begin{equation}\delta_{i-1} > y'_{i-1}\end{equation}
This result is obvious $\delta_{i-1} = y'_{i-1}$ is the straight line going through the end-points of the segment. For the parabola to be concave, it must hold $f''(x) < 0$, i.e.
\begin{equation}\delta_{i-1} < y'_{i-1}\end{equation} </p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhN68k9wU7O6Vdtwru4FPu_v8vuVLgeN0n2tOMFdN5jgNPDyxmcZTnBHVuiMLqhxVf5u2rtJaspMmfKDitf-gzvkmiQg_vuqT6ixsQC_n4EU6gzVe8in-gXKmGBTaDdmAlg6HRkxJgxyBVkM9srxouDHxVx3mT5diw5tYoZX7Oo8SpVlDJqxSxMBl8C/s640/fig_convex.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="480" data-original-width="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhN68k9wU7O6Vdtwru4FPu_v8vuVLgeN0n2tOMFdN5jgNPDyxmcZTnBHVuiMLqhxVf5u2rtJaspMmfKDitf-gzvkmiQg_vuqT6ixsQC_n4EU6gzVe8in-gXKmGBTaDdmAlg6HRkxJgxyBVkM9srxouDHxVx3mT5diw5tYoZX7Oo8SpVlDJqxSxMBl8C/s16000/fig_convex.png" /></a></div><br /><p><br /> </p>Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184670.15339576382115538 -51.4080967 56.773863436178843 18.9044033tag:blogger.com,1999:blog-1169051187754614539.post-50161411999843122692019-12-07T12:04:00.006-08:002023-09-12T12:17:52.913-07:00XXXI Canary Islands Winter School: Computational Fluid Dynamics in Astrophysics (Photos)<p>Two weeks of interesting lectures, intensive hands-on exercises, discussions and meeting new people... A couple of photos by Claudio Dalla Vecchia that perfectly capture the spirit of the school. </p><p><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgkrHmlgcZBR1xqI-WVfFRTy1bklRYdtOwS36RW2RF-QpBPnizjWvPt5HC5k53C5Zlbc-iEXCNa6UorNJ3wmHW4fwny5zyI-KBG88mE852ge4DtDx-cQHt1kU518PwBiwkY6t05-_xU-olsckwxQHfMLic2nZQHMZbKlYdArenju_AT6tJxdxdhkyk5Puw/s4096/ws2019_pic3.jpg" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="2867" data-original-width="4096" height="448" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgkrHmlgcZBR1xqI-WVfFRTy1bklRYdtOwS36RW2RF-QpBPnizjWvPt5HC5k53C5Zlbc-iEXCNa6UorNJ3wmHW4fwny5zyI-KBG88mE852ge4DtDx-cQHt1kU518PwBiwkY6t05-_xU-olsckwxQHfMLic2nZQHMZbKlYdArenju_AT6tJxdxdhkyk5Puw/w640-h448/ws2019_pic3.jpg" width="640" /></a> <br /></p><span><a name='more'></a></span><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjEmK1L3ae4xRq-sb7WnkWclEzY4pSuemOO6J6anK6eTkgdbrFrSnU56as6iCdG3WpUQD32DnEN8z-z-1g4QmwUBTeq61Vz1XD0zMZRVg7Pib7KVlrhXRrwjZiG753ZAckxsZ-d7cH77hjHfbSBa_sLkPedyBWbb4qOcZeamFOhAvCpnqrGwBuQrWHpBY/s4096/ws2019_pic1.jpg" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="2874" data-original-width="4096" height="450" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjEmK1L3ae4xRq-sb7WnkWclEzY4pSuemOO6J6anK6eTkgdbrFrSnU56as6iCdG3WpUQD32DnEN8z-z-1g4QmwUBTeq61Vz1XD0zMZRVg7Pib7KVlrhXRrwjZiG753ZAckxsZ-d7cH77hjHfbSBa_sLkPedyBWbb4qOcZeamFOhAvCpnqrGwBuQrWHpBY/w640-h450/ws2019_pic1.jpg" width="640" /></a></div><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgs1M8jGSa-M0_uM1eYJZh5IfbQ9GeS0GYpGWHIDJleJJaEVsVOShVKvcWLWJsqPKJkl4513BGMGBLF4kZiaoIiqtSL7uznU-rW_UAGOR_JSGZDh0vEJXOJqtcZb4wuUx9xtGmKs84MaLyv5M3ZyrVHvko6tRalowlr0DnBjXYuDjWMiiaak-SMpv3vfgc/s4096/ws2019_pic4.jpg" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="2867" data-original-width="4096" height="448" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgs1M8jGSa-M0_uM1eYJZh5IfbQ9GeS0GYpGWHIDJleJJaEVsVOShVKvcWLWJsqPKJkl4513BGMGBLF4kZiaoIiqtSL7uznU-rW_UAGOR_JSGZDh0vEJXOJqtcZb4wuUx9xtGmKs84MaLyv5M3ZyrVHvko6tRalowlr0DnBjXYuDjWMiiaak-SMpv3vfgc/w640-h448/ws2019_pic4.jpg" width="640" /></a></div><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhenEjByzrFaZR0nPtV0vkZRFa2sCFp4lTuQrEeUyPlE7ZcREeNwy9BCPQWQwaQOPzqstNjwD0xwDXwq0wH5EwyV4oCzUElMnrBc5rkoyL3bl_PL_dWeapnCMMMw5zrWxF9hWksDSjfd7Z0w6mBwzQMmPRirz9kcxHcQp3noA7nxwaLsMnUxkPitMa9sLQ/s4096/ws2019_pic5.jpg" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="2867" data-original-width="4096" height="448" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhenEjByzrFaZR0nPtV0vkZRFa2sCFp4lTuQrEeUyPlE7ZcREeNwy9BCPQWQwaQOPzqstNjwD0xwDXwq0wH5EwyV4oCzUElMnrBc5rkoyL3bl_PL_dWeapnCMMMw5zrWxF9hWksDSjfd7Z0w6mBwzQMmPRirz9kcxHcQp3noA7nxwaLsMnUxkPitMa9sLQ/w640-h448/ws2019_pic5.jpg" width="640" /></a></div>Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.comSanta Cruz de Tenerife, Spain28.4636296 -16.25184670.15339576382115538 -51.4080967 56.773863436178843 18.9044033tag:blogger.com,1999:blog-1169051187754614539.post-72618503482746486952019-04-26T04:40:00.002-07:002019-04-26T07:04:13.421-07:00XXXI Canary Islands Winter School: Computational Fluid Dynamics in Astrophysics<div class="separator" style="clear: both; text-align: left;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj8SvuASVi9_-qO5VbAp3xifRmaQ8jRYQo-mgbDjNoFFof6Mf_HVXxLV3K5JLjSIZx1OahDPH0DUUwGSUtcqztOYmlfReBzTxwIgYzOdyeNBbBPTcg2X0f2Xl20yZCrmCw9qsEsJ9rebJg/s1600/ws.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj8SvuASVi9_-qO5VbAp3xifRmaQ8jRYQo-mgbDjNoFFof6Mf_HVXxLV3K5JLjSIZx1OahDPH0DUUwGSUtcqztOYmlfReBzTxwIgYzOdyeNBbBPTcg2X0f2Xl20yZCrmCw9qsEsJ9rebJg/s1600/ws.png" width="100%" /></a></div>
<br />
Registration is now open for the XXXI Canary Islands Winter School of
Astrophysics, which will be held this year in La Laguna, Tenerife, Spain
from 19-28 November, on the topic of "<b>Computational Fluid Dynamics in
Astrophysics</b>". Confirmed invited lecturers are<br />
<br />
Fernando <b>Moreno-Insertis</b> (Fluid dynamics, conservation laws and waves)<br />
Ake <b>Nordlund</b> (Stellar and planetary formation)<br />
Maarit <b>Käpylä</b> (Dynamics of solar and stellar convection zones)<br />
Matthias <b>Rempel</b> (Solar Atmosphere)<br />
Sascha <b>Husa</b> (General relativity and compact objects)<br />
Stefanie
<b>Walch-Gassner</b> (Interstellar matter)<br />
Tom <b>Theus</b> (Galaxy evolution)<br />
<br />
The school will consist of lectures, seminars tutorials, and social events
including visits to the observatories in Tenerife and La Palma. It is
aimed at advanced graduate students in astrophysics, as well as
postdoctoral researchers, with a strong interest in the computational
fluid dynamics and it applications astrophysics.<br />
<br />
For further information and to register, please visit<br />
<a href="https://www.blogger.com/goog_125125456"><br /></a>
<a href="http://www.iac.es/winterschool/2019/">http://www.iac.es/winterschool/2019/</a><br />
<br />
Claudio Dalla Vecchia<br />
Nikola Vitas<br />
Elena KhomenkoNikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0tag:blogger.com,1999:blog-1169051187754614539.post-87513726341645552482018-10-17T14:29:00.000-07:002018-10-17T14:32:30.935-07:00Three-dimensional simulations of solar magneto-convection including effects of partial ionization<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgH_F3s-58gWUuvLP_JqIOOnMUg75OzIyPZWtyRpQXAMdFTGDdMJJ6nn5Iea8pT5qtwTbrgZeH8lf2YiJLAVmGjPxIxtqO4VIhDskJWJbIV5toqaO5eytSb0TRkxzZZdoIjzHD8IcaHD2I/s1600/big_photocouv618-2.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1600" data-original-width="1484" height="600" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgH_F3s-58gWUuvLP_JqIOOnMUg75OzIyPZWtyRpQXAMdFTGDdMJJ6nn5Iea8pT5qtwTbrgZeH8lf2YiJLAVmGjPxIxtqO4VIhDskJWJbIV5toqaO5eytSb0TRkxzZZdoIjzHD8IcaHD2I/s400/big_photocouv618-2.jpg" width="555" /></a></div>
<br />
Astronomy & Astrophysics Volume 618 (October 2018) is just out with a figure from our paper on the cover!<br />
<br />
The paper (Khomenko, Vitas, Collados & de Vicente 2018: <a href="https://www.aanda.org/articles/aa/full_html/2018/10/aa33048-18/aa33048-18.html">A&A</a>, <a href="http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2018arXiv180701061K&db_key=PRE&link_type=ABSTRACT">ADS</a>) describes the effects of the partial ionization on the structure, dynamics and energy balance of the low chromosphere. Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-42988124760199213692017-07-11T08:40:00.000-07:002024-01-31T13:32:14.123-08:00Ideal EOS for partially ionized hydrogen plasmaPlasma in the solar photosphere can be safely approximated by the ideal gas (non-interacting, randomly moving point-like particles) in the local thermodynamical equilibrium (TE). That means that<br />
<br />
(1) the equation of state can be approximated by the ideal gas law:<br />
\begin{equation}p = n\,kT,\end{equation}<br />
where $p$ is the gas pressure, $n$ is the total number density, $T$ is the temperature and $k$ is the Boltzmann constant.<br />
<br />
For partially ionized gas: <br />
<br />
(2) the ionization fraction can be approximated by its equilibrium value that is set by Saha's ionization formula:<br />
\begin{equation}\frac{n_{i+1}n_{\mathrm{e}}}{n_i} = \frac{1}{\Phi(T)},\end{equation}<br />
$$\frac{1}{\Phi(T)} = \frac{2 U_{i+1}}{U_i}
\left(\frac{2 \pi m_e k T}{h^2}\right)^{3/2}
\mathrm{e}^{-\chi/kT},$$<br />
where $U_i$ and $U_{i+1}$ are the partition functions of the ionization stages $i$ and $i+1$, $m_\mathrm{e}$ is the electron mass, $h$ is Planck's constant and $\chi$ is the ionization energy. <br />
<br />
For a mixture of gases, it means that each component can be approximated by ideal gas, so that the partial pressure of the $i$ component is:<br />
\begin{equation}p_j = n_j\,kT, \end{equation} <br />
and the total pressure is equal to the sum of the partial pressures (Dalton's law):<br />
\begin{equation}p = \sum p_\mathrm{j},\end{equation}<br />
where the summation goes over all components. The components in this case are the free electrons and all the species of the neutral and charged atoms and molecules that are present in the plasma. In the general case the system of the equations is non-linear and contains large number of equations (one for each atomic or molecular specie in every relevant ionization stage). In the theory of stellar atmosphere it is often necessary to solve this system when two parameters are know. <br />
<br />
Nevertheless, there are some trivial solutions. The most notorious is the example of the atmosphere composed of pure hydrogen when the H2 molecules and the negative hydrogen ions are neglected. In that case there are only three types of particles present: protons, neutral H-atoms and free electrons. The rest of this blog refers to that case only.<br />
<br />
<br />
<b>Solution for given temperature and pressure</b><br />
<br />
This simple exercise is described in Mihalas (1970, p.73, <a href="http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1970stat.book.....M&db_key=AST&link_type=ABSTRACT&high=53bed1f8e332292">1970stat.book.....M</a>). The total number of particles in this case is:<br />
\begin{equation}n = n_\mathrm{e} + n_\mathrm{H} + n_\mathrm{H^+} = n_\mathrm{e} + n_\mathrm{H}^\mathrm{tot},\label{eq:pureh1}\end{equation}<br />
where $n_\mathrm{H}^\mathrm{tot} = n_\mathrm{H} + n_\mathrm{H^+}$ is the total number of neutral and ionized H atoms (the total number of H nuclei). The ratio of $n_\mathrm{H^+}$ and $n_\mathrm{H}$ is given by Saha's equation:<br />
\begin{equation}\frac{n_\mathrm{H^+} n_{\mathrm{e}}}{n_\mathrm{H}} = \frac{1}{\Phi(T)}. \label{eq:pureh2}\end{equation}<br />
In addition we know that the number of electrons must be equal to the number of protons: \begin{equation}n_\mathrm{H^+} = n_\mathrm{e}.\label{eq:pureh3}\end{equation}<br />
If we eliminate $n_\mathrm{H}$ and $n_\mathrm{H^+}$ from Eq.$\ref{eq:pureh2}$ and Eq.$\ref{eq:pureh3}$, then Eq.$\ref{eq:pureh1}$ becomes quadratic equation for the electron pressure:<br />
\begin{equation}n = 2 n_\mathrm{e} + n_\mathrm{e}^2 \Phi(T),\label{eq:pureh4}\end{equation}<br />
with the solution:<br />
\begin{equation}n_{\mathrm{e}} = \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)}.\label{eq:pureh5}\end{equation}<br />
To express this solution in terms of the total pressure $p$ and the electron pressure $p_\mathrm{e}$, we assume that the ideal gas law is valid not only for the mixture, but for all of its components as well. Thus: $p_\mathrm{e} = n_\mathrm{e}\,kT$ and the solution for the electron pressure becomes:<br />
\begin{equation}p_{\mathrm{e}} = \left(\sqrt{1 + \frac{p}{kT}\,\Phi(T)}-1 \right)\,\frac{kT}{\Phi(T)}.\label{eq:pureh6}\end{equation}<br />
(In Mihalas' book there is a typo in this equation.) The partial pressure of the hydrogen atoms (neutral and ionized) is obviously:<br />
$$p_\mathrm{H}^\mathrm{tot} = p - p_\mathrm{e}.$$Once we know the electron pressure, it is easy to find other variables. <br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhxP5oPXzx6WxobJDoOATwa0iH-qR_0C0y3vBuhzIOQufTDKxHow_xU_8rXBixH816BUzDgJUBiDblV3Crw3FByqsXbRPSYj1pvHW1ZzLjuuZgG8CT3QZBa9fT8Vm6YJNUftMmqpr7jj6g/s1600/example_eos_pureh_5_lognumberdensities.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhxP5oPXzx6WxobJDoOATwa0iH-qR_0C0y3vBuhzIOQufTDKxHow_xU_8rXBixH816BUzDgJUBiDblV3Crw3FByqsXbRPSYj1pvHW1ZzLjuuZgG8CT3QZBa9fT8Vm6YJNUftMmqpr7jj6g/s1600/example_eos_pureh_5_lognumberdensities.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmnrlSwKIu2QGJnFEp95Dv4RMucf6j_YNuOG3AGISmc3hv9kvmxN4sg0o0kvftQTAKXtKazMfrh1_6mkOaPgcerOEUzZqZm8HRyPVrIvefvzHHXOJ7WDdgbybRfJC8u9FHZWrcO5rmesg/s1600/example_eos_pureh_4_numberdensities.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmnrlSwKIu2QGJnFEp95Dv4RMucf6j_YNuOG3AGISmc3hv9kvmxN4sg0o0kvftQTAKXtKazMfrh1_6mkOaPgcerOEUzZqZm8HRyPVrIvefvzHHXOJ7WDdgbybRfJC8u9FHZWrcO5rmesg/s1600/example_eos_pureh_4_numberdensities.png" width="400" /> </a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<b>Ionization fraction</b><br />
<br />
The ionization fraction is defined as:<br />
\begin{equation}x \equiv \frac{n_\mathrm{H^+}}{n_\mathrm{H}^\mathrm{tot}}=
\frac{n_\mathrm{e}}{n_\mathrm{H}^\mathrm{tot}}.\label{eq:x}\end{equation}<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjZBstWeZjtJzn10_4wjkK6xgDZxmZG7r-eT-jZf_SBXfqivLsDf8sUYHuYlzan_c3yiYAJPQE34R92dreSactQulu9mTaZFwJ_B-2pFl7afCrurlOQzn2lzDUiVIuLMOLm6rWp4a5ktcs/s1600/example_eos_pureh_1_ionizationfraction.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjZBstWeZjtJzn10_4wjkK6xgDZxmZG7r-eT-jZf_SBXfqivLsDf8sUYHuYlzan_c3yiYAJPQE34R92dreSactQulu9mTaZFwJ_B-2pFl7afCrurlOQzn2lzDUiVIuLMOLm6rWp4a5ktcs/s1600/example_eos_pureh_1_ionizationfraction.png" width="400" /></a></div>
<br />
<b>Density</b><br />
<br />
The mass density of a mixture specified by the number densities of its components is:<br />
\begin{equation}\rho \equiv \sum_i n_i\,m_i = n_\mathrm{H}\,m_\mathrm{H}+n_\mathrm{H^+}\,m_\mathrm{H^+}+n_\mathrm{e}\,m_\mathrm{e},\label{eq:rho}\end{equation}<br />
where $m_\mathrm{H}$, $m_\mathrm{H^+}$ and $m_\mathrm{e}$ are masses of one neutral hydrogen atom, one proton and one electron. Since $n_\mathrm{H^+} = n_\mathrm{e}$ and $m_\mathrm{H} \approx m_\mathrm{H^+} \gg m_\mathrm{e}$, the density is<br />
$$\rho \approx m_\mathrm{H}\,n_\mathrm{H}^\mathrm{tot}.$$<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnTeC4x3YGGN2WB22MQe-nGZviKFjXjN1pdG3crs2kosgf57azEuf9VMQEgioevOEqkfa_nEKJ98fzBPwqtCqEcNGbDF1lsWC9tMIaWGkVZzWQLf3hPvGT7ZFp5aFwQTXFYszLX5YmmJ8/s1600/example_eos_pureh_2_massdensity.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnTeC4x3YGGN2WB22MQe-nGZviKFjXjN1pdG3crs2kosgf57azEuf9VMQEgioevOEqkfa_nEKJ98fzBPwqtCqEcNGbDF1lsWC9tMIaWGkVZzWQLf3hPvGT7ZFp5aFwQTXFYszLX5YmmJ8/s1600/example_eos_pureh_2_massdensity.png" width="400" /></a></div>
<br />
<br />
<b>Alternatives</b><br />
<br />
If the density is given instead of the total pressure, then the system becomes:<br />
\begin{eqnarray}\rho &=& n_\mathrm{H}\,m_\mathrm{H} + n_\mathrm{e}(m_\mathrm{H} + m_\mathrm{e}),\nonumber\\ \Phi\,n_\mathrm{e}^2 &=& n_\mathrm{H},\nonumber\end{eqnarray}<br />
yielding again to a quadratic equation for the electron density:<br />
$$\Phi\,m_\mathrm{H}\,n_\mathrm{e}^2 +(m_\mathrm{H} + m_\mathrm{e}) n_\mathrm{e} - \rho = 0.$$<br />
<br />
If the electron pressure is given instead of the density of the pressure, then the solution for the total pressure is $$p = p_\mathrm{e}(\Phi\,p_\mathrm{e} + 2kT),$$ what is equivalent to Eq.$\ref{eq:pureh4}$.<br />
<br />
<br />
<b>Mean molecular weight</b><br />
<br />
The mean molecular weight is defined as<br />
$$\mu = \frac{\rho}{n\,m_\mathscr{A}},$$<br />
where $m_\mathscr{A}$ is the atomic mass unit in grams. If all gas in neutral,<br />
$$\mu = \frac{m_\mathrm{H}\,n_\mathrm{H}}{n_\mathrm{H}\,m_\mathscr{A}} = A_\mathrm{H},$$<br />
with $A_\mathrm{H}$ being the atomic mass of hydrogen ($\approx 1.008$). If all atoms are ionized, the density remains unchanged, but the number of particles is doubled and, therefore, $\mu \approx 0.5$. (The mean molecular weight is always roughly equal to the total number of nuclei per particle.)<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgz1TR7cxWV_u9bTJuOKmih-ROB55N8rr1oWB9u842kd4X_Hi_Y-gNuVWJ-10EI30ldpd11GpKkiZbCH5tZLUF79r2AaMB0p2Rr4KnC0-0YsJukbZ_zKG-hpklB4ear8yk0-IZZJgDdOuo/s1600/example_eos_pureh_3_meanmolecularweight.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgz1TR7cxWV_u9bTJuOKmih-ROB55N8rr1oWB9u842kd4X_Hi_Y-gNuVWJ-10EI30ldpd11GpKkiZbCH5tZLUF79r2AaMB0p2Rr4KnC0-0YsJukbZ_zKG-hpklB4ear8yk0-IZZJgDdOuo/s1600/example_eos_pureh_3_meanmolecularweight.png" width="400" /></a></div>
<br />
<br />
<b>Specific internal energy</b><br />
<br />
The internal energy is equal to the sum of the kinetic part and the part due to the ionization (in principle, there should be another term due to the excitation here, but we neglect it in this example):<br />
$$E_\mathrm{int} = \frac{3}{2}N\,kT + N_\mathrm{e}\,\chi.$$<br />
The specific internal energy pre volume is then $$\varepsilon_V = \frac{E_\mathrm{int}}{V} = \frac{3}{2}n\,kT + n_\mathrm{e}\,\chi,$$ and the specific internal energy per mass is simply<br />
$$\varepsilon_m = \frac{\varepsilon_V}{\rho}.$$<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhEhNItTGzX_RhKKNpiaZl0vFkOT3XxC3y28g-ntdWlGYKghdPZsmcfmmCqQ_T_eWTjkI6k_cg9lrWEulC5Cc3RUQXKZ7R7Er0SmUh1c4YBgzbns2nqSTTebirTpN01k9VJVUGF-lWbd1w/s1600/example_eos_pureh_11_specinternalenergypervolume.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhEhNItTGzX_RhKKNpiaZl0vFkOT3XxC3y28g-ntdWlGYKghdPZsmcfmmCqQ_T_eWTjkI6k_cg9lrWEulC5Cc3RUQXKZ7R7Er0SmUh1c4YBgzbns2nqSTTebirTpN01k9VJVUGF-lWbd1w/s1600/example_eos_pureh_11_specinternalenergypervolume.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgUAXaC-Jp918eiVxlHwO38EtBaWTVeFQNQSoPEhfAxktFEOgYsWBmrCweyF1eIgsJMG2t9wfbE5N9ag-_fJZwjNzYWP-0euyr1_a4YmN8mlCyDbB7CVQMSTu3XQ9xn6eY1ZN4ZX12sVNk/s1600/example_eos_pureh_10_specinternalenergypermass.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgUAXaC-Jp918eiVxlHwO38EtBaWTVeFQNQSoPEhfAxktFEOgYsWBmrCweyF1eIgsJMG2t9wfbE5N9ag-_fJZwjNzYWP-0euyr1_a4YmN8mlCyDbB7CVQMSTu3XQ9xn6eY1ZN4ZX12sVNk/s1600/example_eos_pureh_10_specinternalenergypermass.png" width="400" /></a></div>
<br />
<br />
<b>Specific heat capacity at constant volume (density)</b><br />
<br />
The heat capacity at constant volume $C_V$ is defined as:
$$C_V = \left(\frac{\partial E_\mathrm{int}}{\partial T}\right)_V.$$
The heat capacity per mass at constant volume $c_V$ is defined as:
\begin{equation}c_V = \left(\frac{\partial \varepsilon_m}{\partial T}\right)_\rho = \left[\frac{\partial }{\partial T}\left(\frac{3}{2}\frac{p}{\rho} + \frac{n_\mathrm{e}}{\rho}\,\chi\right) \right]_\rho = \frac{3}{2}\frac{1}{\rho}\left(\frac{\partial p }{\partial T}\right)_\rho + \frac{\chi}{\rho} \left(\frac{\partial n_\mathrm{e}}{\partial T}\right)_\rho. \label{eq:cv}\end{equation}
From the defition of the ionization fraction (Eq.$\ref{eq:x}$), it directly follows:
$$\left(\frac{\partial n_\mathrm{e}}{\partial T}\right)_\rho = \left(\frac{\partial x n_\mathrm{H}^\mathrm{tot}}{\partial T}\right)_\rho = n_\mathrm{H}^\mathrm{tot}\left(\frac{\partial x }{\partial T}\right)_\rho.$$ <br />
The gas pressure is now \begin{equation}p = n\,kT = (n_\mathrm{H}^\mathrm{tot} + n_\mathrm{e})\,kT = n_\mathrm{H}^\mathrm{tot}(1 + x)\,kT = \frac{\rho}{m_\mathrm{H}}(1 + x)\,kT,\label{eq:p}\end{equation} and thus \begin{equation}\left(\frac{\partial p }{\partial T}\right)_\rho = n_\mathrm{H}^\mathrm{tot} k \left[(1 + x) + T \left(\frac{\partial x}{\partial T}\right)_\rho \right].\label{eq:dpdt}\end{equation}The heat capacity for the constant volume becomes:\begin{eqnarray}c_V &=&
\frac{n_\mathrm{H}^\mathrm{tot}}{\rho}\left[\frac{3}{2} k \left((1 + x) + T \left(\frac{\partial x}{\partial T}\right)_\rho \right) + \chi \left(\frac{\partial x }{\partial T}\right)_\rho\right].\nonumber\\
&=&
\frac{n_\mathrm{H}^\mathrm{tot}}{\rho}\left[\frac{3}{2} k (1 + x) + \left(\frac{3}{2} kT + \chi\right) \left( \frac{\partial x }{\partial T}\right)_\rho\right].\label{eq:cv0}\end{eqnarray} <br />
<br />
To find $\partial x/\partial T$, let's rewrite Saha's equation in terms of $x$:<br />
\begin{equation}\frac{n_\mathrm{H^+} n_{\mathrm{e}}}{n_\mathrm{H}} = \frac{x^2}{1-x}\,n_\mathrm{H}^\mathrm{tot}\end{equation} \begin{equation}\frac{ x^2}{1-x} = \frac{m_\mathrm{H}}{\rho} \left( \frac{2\pi m_\mathrm{e} \,kT}{h^2}\right)^{3/2} \mathrm{e}^{-\frac{\chi}{kT}} = q_\rho\,T^{3/2} \mathrm{e}^{-\frac{\chi}{kT}},\label{eq:sahax}\end{equation} where $q_\rho$ contains all the constants. Now the derivative of Eq.$\ref{eq:sahax}$ with respect to $T$ and constant $\rho$ is\begin{eqnarray}\frac{2x(1-x)+x^2}{(1-x)^2}\,\left(\frac{\partial x}{\partial T}\right)_\rho &=& q_\rho \left[\frac{3}{2}T^{1/2} + T^{3/2}\left(\frac{\chi}{kT^2}\right)\right]\,\mathrm{e}^{-\frac{\chi}{kT}},\nonumber\\\frac{x(2-x)}{(1-x)^2}\,\left(\frac{\partial x}{\partial T}\right)_\rho &=& \frac{1}{T} \left(q_\rho
T^{3/2} \,\mathrm{e}^{-\frac{\chi}{kT}}\right) \left(\frac{3}{2} + \frac{\chi}{kT}\right),\nonumber\\<br />
\left(\frac{\partial x}{\partial T}\right)_\rho &=& \frac{(1-x)^2}{x(2-x)}\,\frac{x^2}{1-x}\,\frac{1}{T} \left(\frac{3}{2} + \frac{\chi}{kT}\right),\nonumber\\ \left(\frac{\partial x}{\partial T}\right)_\rho &=& \frac{x(1-x)}{(2-x)}\, \frac{1}{T}\, \left(\frac{3}{2} + \frac{\chi}{kT}\right).\label{eq:dxdt}\end{eqnarray} At this point it is worth to introduce to two new variables to simplify the notation:\begin{equation}D = D(x) \equiv \frac{x\,(1-x)}{(1+x)\,(2-x)},\end{equation}and\begin{equation}H = H(T) \equiv \frac{3}{2} + \frac{\chi}{kT}.\end{equation}<br />
Note that $D\rightarrow 0$ when $x\rightarrow 0$ or $x\rightarrow 1$. Note as well that $1 - D$ reduces to:<br />
$$1-D = \frac{2}{(1+x)\,(2-x)}.$$<br />
Equation $\ref{eq:dxdt}$ we now write simply as: \begin{equation}\left(\frac{\partial x}{\partial T}\right)_\rho = (1 + x)\,\frac{1}{T}\,D\,H.\label{eq:dxdt1}\end{equation}<br />
Substituting Eq.$\ref{eq:dxdt1}$ into Eq.$\ref{eq:cv0}$ yields to:<br />
\begin{equation}c_V =
\frac{n_\mathrm{H}^\mathrm{tot}}{\rho}\left[\frac{3}{2} k (1 + x) +
\left(\frac{3}{2} kT + \chi\right)\, (1 + x)\,\frac{1}{T}\,D\,H\right]\end{equation} <br />
\begin{equation} c_V=
\frac{k}{m_\mathrm{H}}\,(1 + x)\, \left[\frac{3}{2} + D\,H^2\right],\label{eq:cv7}\end{equation}where we used $n_\mathrm{H}^\mathrm{tot}/\rho = 1/m_\mathrm{H}$. Explicitly in terms of $x$ the specific heat at the constant density is:<br />
\begin{equation} c_V=
\frac{k}{m_\mathrm{H}}\,(1 + x)\, \left[\frac{3}{2} + \frac{x\,(1-x)}{(1+x)\,(2-x)}\,\left(\frac{3}{2} + \frac{\chi}{kT}\right)^2\right].\label{eq:cv8}\end{equation}<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4cYhA7osnu6tN_biV09aLKkggwcWZDAw1B9QIoS7TD30dRP-np9IhZ7bYZlAfdhRBBjrwAOwLFEMx60yMBsJ5ECnZU7DUffPAkx1T-ixIQmDR0qRC6Z1m_wJYk_QwCLRlAavTpZ5LUyM/s1600/example_eos_pureh_6_specificheat.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh4cYhA7osnu6tN_biV09aLKkggwcWZDAw1B9QIoS7TD30dRP-np9IhZ7bYZlAfdhRBBjrwAOwLFEMx60yMBsJ5ECnZU7DUffPAkx1T-ixIQmDR0qRC6Z1m_wJYk_QwCLRlAavTpZ5LUyM/s1600/example_eos_pureh_6_specificheat.png" width="400" /></a></div>
<br />
<br />
<b>Specific heat capacity at constant pressure</b><br />
<br />
The next thing to derive is an equivalent expression for the heat capacity per mass at constant pressure. We start from a relation that can easily be derived from the first principle of thermodynamics:<br />
\begin{equation}c_p = c_V + \frac{T}{\rho^2}\,\left(\frac{\partial
p}{\partial T}\right)_\rho^2\left(\frac{\partial p}{\partial
\rho}\right)_T^{-1}.\end{equation}<br />
<br />
First we rewrite it as:<br />
\begin{eqnarray}c_p &=& c_V + \frac{T}{\rho^2}\,\left[\frac{p}{T}
\left(\frac{\partial \ln p}{\partial \ln
T}\right)_\rho\right]^2\,\left[\frac{p}{\rho}\,\left(\frac{\partial \ln
p}{\partial \ln \rho}\right)_T\right]^{-1}.\nonumber\\ &=& c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\left(\frac{\partial
\ln p}{\partial \ln T}\right)_\rho^2\,\left(\frac{\partial \ln p}{\partial
\ln \rho}\right)_T^{-1}.\end{eqnarray}<br />
<br />
The two derivatives we find by differentiating the logarithm of Eq.$\ref{eq:p}$:<br />
$$\ln p = \ln (k/m_\mathrm{H}) + \ln\rho + \ln T + \ln (1 + x),$$ <br />
$$\mathrm{d} \ln p = \mathrm{d} \ln\rho + \mathrm{d} \ln T + \mathrm{d} \ln (1 + x),$$<br />
yielding to:<br />
\begin{equation}\left(\frac{\partial \ln p}{\partial \ln T}\right)_\rho = 1 + \frac{1}{1+x}\left(\frac{\partial x}{\partial \ln T}\right)_\rho = 1 + D\,H,\label{eq:dlnp_rho}\end{equation}<br />
\begin{equation}\left(\frac{\partial \ln p}{\partial \ln \rho}\right)_T = 1
+ \frac{1}{1+x}\left(\frac{\partial x}{\partial \ln \rho}\right)_T = 1 -
D.\label{eq:dlnp_t}\end{equation}<br />
Therefore, for $c_p$ we can write:<br />
\begin{equation}<br />
c_p = c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D}.<br />
\end{equation}<br />
This expression can be written explicitly in terms of $x$ as:<br />
\begin{eqnarray}<br />
c_p &=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{3}{2} + \frac{1 + 2 D\,H + D\,H^2}{1 - D} \right]\nonumber\\<br />
&=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{3}{2} + \frac{(2-x)(1+x) + 2(1-x)x \,H + (1-x)x\,H^2}{2} \right]\nonumber\\<br />
&=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{5}{2} +
\frac{(1-x)x}{2}( 1 + 2 \,H + x\,H^2) \right]\nonumber\\ &=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{5}{2} +
\frac{(1-x)x}{2}(1 +H)^2 \right]\nonumber\\<br />
&=& \frac{k}{m_\mathrm{H}}\,(1+x) \left[\frac{5}{2} +
\frac{(1-x)x}{2}\left(\frac{5}{2} +\frac{\chi}{kT}\right)^2 \right]<br />
\end{eqnarray} <br />
<br />
The ratio of the heat capacities (sometimes labeled with $\gamma$) for the partially ionized pure hydrogen is:<br />
\begin{equation}\gamma \equiv \frac{c_p}{c_V} = 1 + \frac{(1+ D\,H)^2}{(1 - D)\,(3/2 + D\,H^2)}<br />
.\label{eq:cpcv}\end{equation}<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhWczMq-5rsFbAUvDwwjmAeG2TfbBTkTB3GsmF9s4ADXCgS95r0TPVJYoqFKj462bcuTL5uwkTAo0m72P0my5voVLzhb0crB-ut0v1qCegovPtmbfBt4urdI-FCJ5pOjRT1Rgy4Jtjf98Y/s1600/example_eos_pureh_7_heatcapacityratio.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhWczMq-5rsFbAUvDwwjmAeG2TfbBTkTB3GsmF9s4ADXCgS95r0TPVJYoqFKj462bcuTL5uwkTAo0m72P0my5voVLzhb0crB-ut0v1qCegovPtmbfBt4urdI-FCJ5pOjRT1Rgy4Jtjf98Y/s1600/example_eos_pureh_7_heatcapacityratio.png" width="400" /></a></div>
<br />
<br />
<b>Mayer's relation </b><br />
<br />
Mayer's relation for ideal gas states that $c_p = c_v + nk$. In a gas mixture, this should apply to each component of free particles that contribute to the internal energy. Before we had:
\begin{equation}
c_p = c_V + \frac{k}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D}.<br />
\end{equation}
Therefore, for Mayer's equation to be true, we have to show that:
$$n = \frac{1}{m_\mathrm{H}}\,(1+x)\,\frac{(1 + D\,H)^2}{1-D} $$
<br />
<br />
<b>Adiabatic exponents </b><br />
<br />
In addition we derive the explicit expression of the adiabatic exponents $\Gamma$ of the partially ionized pure hydrogen gas. The exponents are defined as:<br />
\begin{equation}\gamma \equiv \Gamma_1 \equiv \left(\frac{\partial \ln p}{\partial \ln \rho}\right)_S,\end{equation}<br />
\begin{equation}\nabla_\mathrm{ad}\equiv\frac{\Gamma_2 -1}{\Gamma_2} \equiv \left(\frac{\partial \ln T}{\partial \ln p}\right)_S,\label{eq:G1}\end{equation}<br />
\begin{equation}\gamma\,\nabla_\mathrm{ad} \equiv \Gamma_3 -1 \equiv \left(\frac{\partial \ln T}{\partial \ln \rho}\right)_S.\end{equation}<br />
Comment on the meaining of $\Gamma$'s. Note that both $\gamma$ in Eq.$\ref{eq:cpcv}$ and in Eq.$\ref{eq:G1}$ is not the same. It is easy to express the three adiabatic coefficients as function of the gradients in Eqs.$\ref{eq:dlnp_rho}$ and $\ref{eq:dlnp_t}$:<br />
\begin{equation}\Gamma_1 =\left( \frac{\partial \ln p}{\partial \ln \rho}\right)_T,\label{eq:G1_1}\end{equation}<br />
\begin{equation}\frac{\Gamma_2 -1}{\Gamma_2} = \frac{k}{m_\mathrm{H}}\,\frac{1}{c_p} \left( \frac{\partial \ln p}{\partial \ln T}\right)_\rho \left( \frac{\partial \ln p}{\partial \ln \rho}\right)_T^{-1},\label{eq:G2_1}\end{equation}<br />
\begin{equation}\Gamma_3 - 1 = \frac{k}{m_\mathrm{H}}\,\frac{1}{c_v} \left( \frac{\partial \ln p}{\partial \ln T}\right)_\rho.\label{eq:G3_1}\end{equation}<br />
After substituting the results for the gradients Eqs.$\ref{eq:dlnp_rho}$ and $\ref{eq:dlnp_t}$, the final expression of the adiabatic exponents for the partially ionized pure hydrogen plasma are:<br />
\begin{equation}\Gamma_1= \frac{c_p}{c_v}\,(1 - D),\end{equation}<br />
\begin{equation}\frac{\Gamma_2 -1}{\Gamma_2} = \frac{c_p - c_V}{c_p}\,\frac{1}{1 + D\,H},\end{equation}<br />
\begin{equation}\Gamma_3 - 1= \frac{c_p - c_V}{c_V}\,\frac{1-D}{1 + D\,H}.\end{equation}<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi7SCGwdZBM0Cff0PJwK-irjuColQIyLpx54y94T3zhhUzPfF3pqcqN3QGd5C7LjW8RwhIKeFJBjJKr3SWLv3m46MQJc13_zHKyJ9dEAl-tAS82BMK6cE6tvH0C5b4kXQY2OSOjtZuKOBU/s1600/example_eos_pureh_8_gammas.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi7SCGwdZBM0Cff0PJwK-irjuColQIyLpx54y94T3zhhUzPfF3pqcqN3QGd5C7LjW8RwhIKeFJBjJKr3SWLv3m46MQJc13_zHKyJ9dEAl-tAS82BMK6cE6tvH0C5b4kXQY2OSOjtZuKOBU/s1600/example_eos_pureh_8_gammas.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEieayCgijPhks1Zryfa0ungROSwPeTim6yXieNUZZzCJJkaddwk1LjrWzos9NixK17G1iZIQIpdvlhiqEie8F47c0u5pKOQ9p8UafyOXOA2Zc2e27N6hXSjp2lBF2C23UDSsEwsolGLEqg/s1600/example_eos_pureh_9_nablaad.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEieayCgijPhks1Zryfa0ungROSwPeTim6yXieNUZZzCJJkaddwk1LjrWzos9NixK17G1iZIQIpdvlhiqEie8F47c0u5pKOQ9p8UafyOXOA2Zc2e27N6hXSjp2lBF2C23UDSsEwsolGLEqg/s1600/example_eos_pureh_9_nablaad.png" width="400" /></a></div>
<br />
<b>Entropy</b><br />
<br />
Let's now compute the specific entropy per mass $s$. It can be shown that (see my notes on computing entropy for $T$ and $\rho$ as independent variables:
$$ds = \frac{1}{T} \left[\left(\frac{\partial \varepsilon_m}{\partial \rho}\right)_T - \frac{p}{\rho^2}\right] d\rho + \frac{1}{T}\left(\frac{\partial \varepsilon_m}{\partial T}\right)_\rho dT $$
To solve for these derivatives analytically, we need to find $\varepsilon_m = \varepsilon_m(T, \rho)$:
$$\varepsilon_m = \frac{1}{\rho} \left[\frac{3}{2}n\,kT + n_\mathrm{e}\,\chi\right]$$
and its derivatives:
$$\left(\frac{\partial \varepsilon_m}{\partial \rho}\right)_T = -\frac{1}{\rho^2} \left[\frac{3}{2}n\,kT + n_\mathrm{e}\,\chi\right] + \frac{3}{2} \frac{kT}{\rho} \left(\frac{\partial n}{\partial \rho}\right)_T + \frac{\chi}{\rho}\left(\frac{\partial n_\mathrm{e}}{\partial \rho}\right)_T$$
and
$$\left(\frac{\partial \varepsilon_m}{\partial T}\right)_\rho =$$
Where we still need to express $n$ and $n_e$ as functions of $T$, $\rho$. By definition:
$$\rho = n_\mathrm{H}\,m_\mathrm{H} + n_\mathrm{e}(m_\mathrm{H} + m_\mathrm{e})$$
and
$$n = n_\mathrm{H} + 2 n_\mathrm{e}$$
while Saha gives us:
$$\Phi\,n_\mathrm{e}^2 = n_\mathrm{H}$$<br />
As we saw before, the last two equations combine into a quadratic equation with a solution:
$$n_{\mathrm{e}} = \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)}.$$
And for $\rho(n)$ we get explicitly:
$$\rho = (n - n_\mathrm{e}) \,m_\mathrm{H} + n_\mathrm{e} m_\mathrm{e}$$
$$\rho = \left(n - \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)}\right) \,m_\mathrm{H} + \left(\sqrt{1 + n\,\Phi(T)}-1 \right)\,\frac{1}{\Phi(T)} m_\mathrm{e}$$
$$\rho = n m_\mathrm{H} - \frac{\sqrt{1 + n\,\Phi(T)}-1 }{\Phi(T)} (m_\mathrm{H} -m_\mathrm{e})$$
Therefore:
$$\left(\frac{\partial \rho }{\partial n}\right)_T = m_\mathrm{H} - (m_\mathrm{H} -m_\mathrm{e}) \frac{1}{2 \sqrt{1 + n\,\Phi(T)}}$$
$$\left(\frac{\partial n }{\partial \rho}\right)_T = \frac{2 \sqrt{1 + n\,\Phi(T)}}{2 m_\mathrm{H}\, \sqrt{1 + n\,\Phi(T)} - (m_\mathrm{H} -m_\mathrm{e})}$$
<br />
<b>References</b><br />
<br />
The general thermodynamical expressions and their derivation can be found in any textbook on the thermodynamics or on the stellar structure, e.g. in Chandrasekhar, S. (1967, <a href="http://adsabs.harvard.edu/abs/1967aits.book.....C">1967aits.book.....C</a>). There is a brief mentioning of the equation of state for the pure hydrogen in Mihalas, D. (1970, p.73, <a href="http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=1970stat.book.....M&db_key=AST&link_type=ABSTRACT&high=53bed1f8e332292">1970stat.book.....M</a>). Biermann, L. (1942, <a href="http://adsabs.harvard.edu/abs/1942ZA.....21..320B">1942ZA.....21..320B</a>) was the first to derive the explicit expressions for the $c_p$ and $c_V$ (see Lobel, A. (2001, <a href="http://adsabs.harvard.edu/abs/2001ApJ...558..780L">2001ApJ...558..780L</a>)). Knapp, G (<a href="http://www.princeton.edu/astro/people/faculty/gk/" target="_self">lecture notes</a>) gives most of the derivations above neatly and tidily. Hansen, C.J. et al (2004, <a href="http://adsabs.harvard.edu/abs/2004sipp.book.....H">2004sipp.book.....H</a>) derive several cases including the pure hydrogen and discuss the solution.<br />
<br />
<br />
<br />
<br />
<b>IDL implementation</b><br />
<br />
IDL function that evaluates the equation of state for pure hydrogen for given temperature and gas pressure can be downloaded from here: eos_pureh.pro.Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0C. Vía Láctea, 1, 38108 La Laguna, Santa Cruz de Tenerife, Spain28.4744413 -16.308832128.470668869171917 -16.313123634423828 28.478213730828081 -16.304540565576172tag:blogger.com,1999:blog-1169051187754614539.post-29562199983919035302017-07-04T04:50:00.001-07:002017-07-04T04:54:13.003-07:00How to display magnetic field lines in IDL?It is a common problem in visualization of magnetic fields. If we assume that the potential $A$ of the magnetic field is knows, the field lines are, by definition, iso-$A$ curves. The quickest way to display them is by using the CONTOUR procedure. If the potential field is given as a 2D variable potential, then:<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); font-family: 'Andale Mono', fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>IDL> CONTOUR, potential, levels = levels, /xs, /ys</code></pre>
gives:
<div class="separator" style="clear: both; text-align: center;">
<code><br /></code>
<code></code><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUc7HZCBAVH9nuyxQKTaifLn6VtUO6HELxP5VwWBE0KHxqLnZ5Pn0KtKjfAfDHIcP2HQy0ZJ8hrlLaDsFBFliPo2HcEx-kyDep0Y8tYaRdiH-b-q6gR6Cjo4VjpaSpQixTPkgkNikHgh4/s1600/magnetic_1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="400" data-original-width="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiUc7HZCBAVH9nuyxQKTaifLn6VtUO6HELxP5VwWBE0KHxqLnZ5Pn0KtKjfAfDHIcP2HQy0ZJ8hrlLaDsFBFliPo2HcEx-kyDep0Y8tYaRdiH-b-q6gR6Cjo4VjpaSpQixTPkgkNikHgh4/s1600/magnetic_1.png" /></a></div>
However, the information here is not complete without showing the actual direction of the field along the lines. It is easy to do it in IDL:<br />
<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); font-family: 'Andale Mono', fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>IDL> CONTOUR, potential, levels = levels, /xs, /ys
IDL> CONTOUR, potential, levels = levels, /xs, /ys, path_xy = c
IDL> FOR i = 1, N_ELEMENTS(c)/2-1, 50 DO $</code><code><code> ARROW, c[0, i-1], c[1, i-1], c[0, i], c[1, i], /norm, /solid, hsize = 5 </code>
</code></pre>
<div>
<code><br /></code>
<code></code>produces:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjejXBEtJkXZolVs_usqxJB6cCCtYryd9nF0sIlPBDuJAXyLdV4mZK7JKmwrWEWSPQNjeMPk4CrBsqHnL57DUSVD4lyL7DIkVlOBVplYYUzJzUb0lG7hQfG7p5jhXuONCUTJpa_B5AZxh4/s1600/magnetic_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="400" data-original-width="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjejXBEtJkXZolVs_usqxJB6cCCtYryd9nF0sIlPBDuJAXyLdV4mZK7JKmwrWEWSPQNjeMPk4CrBsqHnL57DUSVD4lyL7DIkVlOBVplYYUzJzUb0lG7hQfG7p5jhXuONCUTJpa_B5AZxh4/s1600/magnetic_2.png" /></a></div>
It is much better now, but the arrow heads from different field lines make it crowded. If we are interested in individual field lines, we can add color to this:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg8_m83A3iARpuBjIDreijZZsOtSKvF72ndUUJ9PtNy_iGDCkOij5QmwilwszvPXYKa8jCNm5-9fWsnmfvgxyGfs3c4qEi48dWkW9DCqU0erFrCqLbYfRVT7ojPUrcAlIVt_L-ABDfdX4Y/s1600/magnetic_3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="400" data-original-width="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg8_m83A3iARpuBjIDreijZZsOtSKvF72ndUUJ9PtNy_iGDCkOij5QmwilwszvPXYKa8jCNm5-9fWsnmfvgxyGfs3c4qEi48dWkW9DCqU0erFrCqLbYfRVT7ojPUrcAlIVt_L-ABDfdX4Y/s1600/magnetic_3.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
</div>
Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-34271598481889714632017-04-11T06:48:00.000-07:002017-04-20T09:29:36.921-07:00Deep-learning about horizontal velocities at the solar surfaceThe velocity fields are of great importance for understanding dynamics and structure of the solar atmosphere. The line of sight velocities are coded in the wavelength shifts of the spectral lines, thanks to the Doppler effect, and relatively easy to measure. On the other hand, the orthogonal ("horizontal") components of the velocity vector are impossible to measure directly. <br />
<br />
The most popular method for estimating the horizontal velocities is so-called local correlation tracking (LCT, November & Simon, 1988). It is based on comparing successive images of the solar surface in the continuum light and transforming their differences into information about the horizontal fields. However, the LCT algorithm suffers from several limitations.<br />
<br />
In a paper by Andres Asensio Ramos and Iker S. Requerey (with a small contribution from my side) accepted by A&A and published on Arxiv some weeks ago (<a href="http://adsabs.harvard.edu/abs/2017arXiv170305128A">2017arXiv170305128A</a>) this problem is tackled by the deep-learning approach. A deep fully convolutional neural network is trained on synthetic observations from 3D MHD simulations of the solar photosphere and then applied to the real observation with the IMaX instrument on board the SUNRISE balloon (Martinez Pillet et al, 2011; Solanki, 2010). The method is validated using simulation snapshots of the quiet sun produced with the MANCHA code that I have been developing in the last couple of years.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhLpo0laPEDR2jw8cSsWqgoKzvBqtQyYNswgjbpJNaVvwOSYAeU_-H5EtegYiPGcSw71zpI_Bc3Z78XxpgT5e5wdYJkBGgH4YyPOsBhC_VvTjxUMtdc337mC5EyS-EmP40APGwdnWhQCIQ/s1600/deepvel_fig2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhLpo0laPEDR2jw8cSsWqgoKzvBqtQyYNswgjbpJNaVvwOSYAeU_-H5EtegYiPGcSw71zpI_Bc3Z78XxpgT5e5wdYJkBGgH4YyPOsBhC_VvTjxUMtdc337mC5EyS-EmP40APGwdnWhQCIQ/s1600/deepvel_fig2.png" /></a></div>
<br />Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-60442234384683166442017-03-20T16:03:00.000-07:002017-03-20T16:03:49.756-07:00K-means clusteringThe 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)?<br />
<br />
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).<br />
<br />
<b>k-means</b><br />
<br />
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]$. <br />
<br />
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:<br />
<br />
1. Choose randomly $k$ measurements as initial cluster centers: $c_0, ..., c_{k-1}$. Obviously, each of the clusters is also $n$-dimensional vector.<br />
<br />
2. Compute Euclidean distance $D_{i, l}$ between every measurement $x_i$ and every cluster center $c_l$:<br />
$$D_{i, l} = \sqrt{\sum_{j=0}^{n-1} (x_{i, j} - c_{l, j})^2}.$$<br />
3. Assign every measurement $x_i$ to the cluster represented by the closest cluster center $c_l$.<br />
<br />
4. Now compute new cluster centers by simply averaging all the measurements in each cluster.<br />
<br />
5. Go back to 2. and keep iterating until none of the measurements changes its cluster in two successive iterations.<br />
<br />
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.<br />
<br />
<a name='more'></a><br /><br />
<b>Furthest point</b> <br />
<br />
In the furthest cluster variant, only the very fist cluster center $c_0$ is obtained randomly. Then the distances of all other measurements from $c_0$ are computed and the most distant measurement is selected as $c_1$. The next cluster is selected as the most distant from both $c_0$ and $c_1$, and so on until all initial cluster centers are set. An issue with this method of initialization is that it is very sensitive to outliers especially those on the edges of the measurement space.<br />
<br />
<b>k-means++</b><br />
<b> </b><br />
<a href="http://ilpubs.stanford.edu:8090/778/1/2006-13.pdf">Arthur and Vassilvitskii (2006)</a> proposed an alternative method for the initial randomization. The first cluster center $c_0$ is purely random. Then the distances of all other measurements from $c_0$ are computed as before. Now this distances are used to define probability distribution,<br />
$$p = \frac{D_{i, 0}^2}{\sum_{i=0}^{m-1} D_{i, 0}^2}$$,<br />
and the next cluster is chosen again randomly, but now we used the probability distribution $p$ (instead of uniform). The idea is again to draw new clusters so that they are far from the previously chosen. The difference is that now the new clusters will not necessarily be at the very edges of the measurement space, so it is less likely that isolated outliers would be picked.<br />
<br />
<b>Example of initialization</b><br />
<br />
I constructed a simple test data with 4 clearly separated sets with $m=2$ and then run the three variants for the k-means realization assuming $k=4$. The result is shown in the figure below. The first selected cluster-center is yellow. <br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6W7_5N-d5qDx9-W7mnG7gTdjh1CMYPkbFiFunR-ox9z7Y4lhx_v7oKY4YcavZB2Y1EMckRwqZ30I6Pe1bMCYZSijjPGvBbej2iKfxq4CF32qtVZAQMjGeSB8ZpZudEZpoj3JjtKDnXr4/s1600/fig_initialization_kmeans.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6W7_5N-d5qDx9-W7mnG7gTdjh1CMYPkbFiFunR-ox9z7Y4lhx_v7oKY4YcavZB2Y1EMckRwqZ30I6Pe1bMCYZSijjPGvBbej2iKfxq4CF32qtVZAQMjGeSB8ZpZudEZpoj3JjtKDnXr4/s1600/fig_initialization_kmeans.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
In this trivial case, the 4 clusters are correctly identified already in the first iteration for the furthest-point and k-means++ initializations. In the random initialization, it takes 4 iteration steps. The following figure shows how does the solution propagates step-by-step.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiQ0N9urLhWrJ9Fbl83MjXtPf4WNf51x8yocdqXRjr7jlFCt7iHw4gyKpMONflLH2Clit7iN2umc5qm4xKSUqbxf6zCyKrdE6jWB0Grjmdqlvykg48wu6pCAWNKn_foHLIisi6C7Pg3EiE/s1600/fig_kmeans_animation.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiQ0N9urLhWrJ9Fbl83MjXtPf4WNf51x8yocdqXRjr7jlFCt7iHw4gyKpMONflLH2Clit7iN2umc5qm4xKSUqbxf6zCyKrdE6jWB0Grjmdqlvykg48wu6pCAWNKn_foHLIisi6C7Pg3EiE/s1600/fig_kmeans_animation.gif" /></a></div>
<br />
<b>Silhouette</b><br />
<br />
Rousseew (1987) proposed the concept of silhouette as an estimate how consistent are the clusters and how appropriate is clustering for a given dataset. He introduced a measure of dissimilarity as the mean distance between a measurement and all measurements in one cluster. Let's label with $a_j$ the dissimilarity between a measurement and the cluster to which this measurement is assigned, and with $b_j$ the mean dissimilarity between the same measurement and all other clusters. Then the silhouette is defined for every measurement as:
$$s_j =\frac{b_j - a_j}{\mathrm{max}(a_j, b_j)}.$$
The value of silhouette averaged over the entire dataset tells us how well the data have been clustered. If the value is closed to 1, the clustering is well done. If it is close to 0, the clusters found by the method are not appropriate for the given data.<br />
<br />
<br />
<b>IDL implementation</b><br />
<br />
IDL has its own implementation of k-means (CLUSTER, before it was called KMEANS). To be able to experiment and play with the method, I wrote my KMEANS function that includes the three initializations mentioned above. The result of the function is structure with three tags: clusters, frequencies (percentage of the measurements assigned to each cluster) and silhouette (for every measurement).<br />
<br />
Example:<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); font-family: 'Andale Mono', fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>IDL> n = 100 ; Number of pointsIDL> m = 2 ; Number of dimensions, ie. length of vectors
IDL> k = 4 ; Number of clusters
IDL> x = REFORM(RANDOMU(10, m*n)/4., [m, n])
IDL> x[0, 0:24] += 0.2
IDL> x[1, 0:24] += 0.2
IDL> x[0, 25:49] += 0.8
IDL> x[1, 25:49] += 0.8
IDL> x[0, 50:74] += 0.2
IDL> x[1, 50:74] += 0.8
IDL> x[0, 75:99] += 0.8
IDL> x[1, 75:99] += 0.2
IDL> y = KMEANS(x, k = k, initial = 'kmeanspp')
; Show the initial data points
IDL> PLOT, x[0, *], x[1, *], psym = 4
; Show data points which are assigned to the cluster 2
IDL> id = WHERE(y.clusters EQ 2)
IDL> OPLOT, x[0, id], x[1, id], psym = -1
</code></pre>
<div>
<code><br /></code>
<code><br /></code></div>
<b>Download:</b><a href="https://www.dropbox.com/s/2t76ku43gqts1uz/kmeans.pro?dl=0"> kmeans.pro </a><br />
<br />
<br />
References:<br />
<a href="http://ac.els-cdn.com/0377042787901257/1-s2.0-0377042787901257-main.pdf?_tid=aec52370-0cfd-11e7-998c-00000aacb361&acdnat=1489967115_e241b403c976118fd125bb264e505c09">Rousseeuw (1987), Silhouettes: a graphical aid to the interpretation and validation of cluster analysis</a> <br />
<a href="https://pdfs.semanticscholar.org/d853/79ef05fd9b02f71f947f1b58df474773767f.pdf">Su & Dy (2004), A Deterministic Method for Initializing K-means Clustering</a><br />
<a href="http://ilpubs.stanford.edu:8090/778/1/2006-13.pdf">Arthur & Vassilvitskii (2007), k-means++ : The Advantages of Careful Seeding</a><br />
<a href="http://vldb.org/pvldb/vol5/p622_bahmanbahmani_vldb2012.pdf">Bahmani et al (2013), Scalable K-Means++</a><br />
<br />
Application in solar physics:<br />
<a href="http://adsabs.harvard.edu/abs/2011A%26A...530A..14V">Viticchié & Sánchez Almeida (2011), Asymmetries of the Stokes V profiles observed by HINODE SOT/SP in the quiet Sun</a>Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com1Santa Cruz de Tenerife, Spain28.461749499520266 -16.26628875732421928.447789999520268 -16.286458757324219 28.475708999520265 -16.246118757324219tag:blogger.com,1999:blog-1169051187754614539.post-43722838571565058862017-02-01T16:40:00.000-08:002017-03-20T16:42:50.691-07:00First observation of linear polarization in the forbidden [OI] 630.03 nm lineIn a new paper (<a href="http://adsabs.harvard.edu/abs/2017ApJ...836...29D">de Wijn, Socas-Navarro & Vitas, 2017, ApJ, 836, 29D</a>) 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.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjaxgNcSYFEwU45cCUbigJzSuLFxmGGJYoO0_asPOCX8nK3s7bD2ON-5M0BWAuPwkzkOTYd85beE9n7VxN2jbk_f4h35Kt3AfLb-46is8E1ObVJUAbKdZ0VsY2LSszjvD9YqVkkwg3znzM/s1600/fig2_apj_linpol.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjaxgNcSYFEwU45cCUbigJzSuLFxmGGJYoO0_asPOCX8nK3s7bD2ON-5M0BWAuPwkzkOTYd85beE9n7VxN2jbk_f4h35Kt3AfLb-46is8E1ObVJUAbKdZ0VsY2LSszjvD9YqVkkwg3znzM/s1600/fig2_apj_linpol.jpg" /> </a> </div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
More details of this unique observation will appear soon in a follow-up publication.</div>
Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-19546767211148977012016-11-30T04:16:00.000-08:002017-04-04T07:26:51.674-07:001D 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", <a href="http://adsabs.harvard.edu/abs/1989A%26A...225..204D">1989A&A...225..204D</a>). 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.<br />
<br />
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 <a href="http://adsabs.harvard.edu/abs/1981ApJS...45..635V">Yun et al. (1981)</a>. However, between the optical depths -3 and -4 it becomes close to the VALC model (<a href="http://adsabs.harvard.edu/abs/1981ApJS...45..635V">Vernazza et al, 1981</a>).<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiT53L3XQas2Fv38BHOQ2d4E02Q7Rn_9xNI_tjxR9YJ7RogFaCqoUwve738N2G_iMbJ-5wDX4roc92DXovPfmPcLHr97E_NI466DuhaJl2SM4cw1T8YfDWyjgkFHSAZQMBywXXwkkqKl8c/s1600/fig2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiT53L3XQas2Fv38BHOQ2d4E02Q7Rn_9xNI_tjxR9YJ7RogFaCqoUwve738N2G_iMbJ-5wDX4roc92DXovPfmPcLHr97E_NI466DuhaJl2SM4cw1T8YfDWyjgkFHSAZQMBywXXwkkqKl8c/s1600/fig2.png" /></a></div>
<br />
<a name='more'></a><br />
<br />
This model had been used to construct a simple sunspot model for helioseismology (<a href="http://adsabs.harvard.edu/abs/2011SoPh..268..293C">Cameron et al, 2010</a>) and in many other studies.<br />
<br />
Files: <a href="https://drive.google.com/file/d/0B_BtkJGGslsVZUxwX0NvcGhZUDA/view?usp=sharing">IDL save version</a>, <a href="https://drive.google.com/file/d/0B_BtkJGGslsVYVhMdGIyd0NwLXc/view?usp=sharing">ascii version</a><br />
<br />
<b>Disclaimer</b>: This data is publicly available in Ding & Fang (1989, <a href="http://adsabs.harvard.edu/abs/1989A%26A...225..204D">1989A&A...225..204D</a>).
If you use the data in a publication please cite it properly. I
double-checked that the digitized version here matches the original,
still typos are possible no matter how unlikely. If you have any doubts,
compare the numbers to those in the paper. Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-30723566192740671272016-10-21T04:45:00.002-07:002016-10-21T04:48:19.999-07:00Installing IDL v7.1 on Kubuntu 16.04 (on MacBook pro)<br />
The IDL installation went smooth as usual, but when I tried to run IDL, there was an error:<br />
<br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;"> error while loading shared libraries: libXp.so.6: cannot open shared object file: No such file or directory </span></span><br />
<br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;">In the<a href="http://www.harrisgeospatial.com/Company/PressRoom/Blogs/TabId/836/ArtMID/2928/userid/-1/ArticleID/13759/IDL-fails-to-install-on-Linux-What-to-do.aspx"> IDL help pages ( IDL fails to install on Linux: What to do, scroll down to the Ubuntu section)</a> there is a comment on that that turned only half useful. Two libraries I got installed with no problem:</span></span><br />
<br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo apt-get install libxmu-dev<br />
sudo apt-get install libxmu6</span></span><br />
<br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;">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.</span></span><br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span></span>
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo apt-get install autoconf autogen intltool</span></span><br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo apt-get install git</span></span><br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo apt-get install xutils-dev libtool libx11-dev </span></span><br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;"> x11proto-xext-dev x11proto-print-dev</span></span><br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">git clone https://cgit.freedesktop.org/xorg/lib/libXp/</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;">Then cd to libXp directory and execute</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo ./autogen.sh</span></span></span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;">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:</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo ./configure</span></span><br />
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">sudo make install</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;">And finally I add the path to the library to .bashrc (all in one line):</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-size: x-small;"><span style="font-family: "courier new" , "courier" , monospace;">echo 'export LD_LIBRARY_PATH="/usr/local/lib/:$LD_LIBRARY_PATH"' >> ~/.bashrc</span></span><br />
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;"><br /></span></span>
<span style="font-family: "courier new" , "courier" , monospace;"><span style="font-family: "arial" , "helvetica" , sans-serif;">After that IDL worked normally. </span></span>
<!-- https://afni.nimh.nih.gov/afni/community/board/read.php?1,151673,151673 -->Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com1San Cristóbal de La Laguna, Santa Cruz de Tenerife, Spain28.4874009 -16.31590610000000728.4315669 -16.396587100000005 28.5432349 -16.235225100000008tag:blogger.com,1999:blog-1169051187754614539.post-33390581614387719422016-09-23T05:31:00.000-07:002017-04-05T05:47:22.073-07:00How 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 <a href="http://www.idlcoyote.com/tips/point_in_polygon.html">here</a> and get the source of the <a href="http://www.idlcoyote.com/tip_examples/inside.pro">inside.pro</a>).<br />
<br />
Example:<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); color: black; font-family: Andale Mono,fixed,monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>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
</code></pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhdMgUW24s0O-bimmmIHDs__yksjtgPB9ASe_pLWrSiMIT82EhCRDhAVbKZ1M9O4FQ8pBm-8jVmqE0pPwM2956DvAr8xYZ6DrmpEiNRO-WTYJCSwJ9M6ytiiPbksVqUMQ-uhcj23YIq31s/s1600/blog_inside.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhdMgUW24s0O-bimmmIHDs__yksjtgPB9ASe_pLWrSiMIT82EhCRDhAVbKZ1M9O4FQ8pBm-8jVmqE0pPwM2956DvAr8xYZ6DrmpEiNRO-WTYJCSwJ9M6ytiiPbksVqUMQ-uhcj23YIq31s/s1600/blog_inside.png" /></a></div>
<br />Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-70663010599065755312016-07-16T10:06:00.004-07:002022-07-22T17:16:24.779-07:00Equation of state: Vardya - Mihalas - Wittmann<p><i>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.</i><br /></p><p>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.<br />
<br />
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.</p><p><br />
<br />
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.<br />
<br />
<b>Definitions</b><br />
<br />
Let's first define the pressures:<br />
$p_{\mathrm{H}}$ - partial pressure of the neutral H atoms;<br />
$p_{\mathrm{H^+}}$ - partial pressure of the positive H ions (protons);<br />
</p><span><a name='more'></a></span><span><!--more--></span><div>
$p_{\mathrm{H^-}}$ - partial pressure of the negative H ions;</div>
<div>
$p_{\mathrm{H_2}}$ - partial pressure of the H2 molecules;</div>
<div>
$p_{\mathrm{H_2^+}}$ - partial pressure of the H2+ molecules;</div>
<div>
$p_{\mathrm{e}}$ - electron pressure.</div>
<div>
<br /></div>
<div>
These pressure are proportional to the corresponding number densities:</div>
<div>
$$p_{\mathrm{H}} = n_{\mathrm{H}} k T;$$</div>
<div>
$$p_{\mathrm{H^+}} = n_{\mathrm{H^+}} k T;$$</div>
<div>
$$p_{\mathrm{H^-}} = n_{\mathrm{H^-}} k T;$$</div>
<div>
$$p_{\mathrm{H_2}} = n_{\mathrm{H_2}} k T;$$</div>
<div>
$$p_{\mathrm{H_2^+}} = n_{\mathrm{H_2^+}} k T;$$</div>
<div>
$$p_{\mathrm{e}} = n_{\mathrm{e}} k T.$$<br />
<br />
<b>Particle conservation </b></div>
<div>
The total pressure due to the presence of H in the plasma is then:</div>
\begin{equation}
p_{\mathrm{H^{tot}}} = p_{\mathrm{H}} + p_{\mathrm{H^+}} + p_{\mathrm{H^-}}+p_{\mathrm{H_2}}+p_{\mathrm{H_2^+}}.
\end{equation}
<br />
<div>
However, if we write this equation in terms of the number density, it becomes: </div>
\begin{equation}
[n_{\mathrm{H^{tot}}}]_\mathrm{particles} = n_{\mathrm{H}} + n_{\mathrm{H^+}} + n_{\mathrm{H^-}}+n_{\mathrm{H_2}}+n_{\mathrm{H_2^+}}.
\end{equation}
<br />
<div>
On the other hand, the particle (nuclei) conservation requires that </div>
<div>
\begin{equation}[n_{\mathrm{H^{tot}}}]_\mathrm{nuclei} = n_{\mathrm{H}} + n_{\mathrm{H^+}} + n_{\mathrm{H^-}}+2 n_{\mathrm{H_2}}+2 n_{\mathrm{H_2^+}},
\end{equation}
</div>
<div>
and</div>
<div>
\begin{equation}
\pi_{\mathrm{H^{tot}}} = p_{\mathrm{H}} + p_{\mathrm{H^+}} + p_{\mathrm{H^-}}+2 p_{\mathrm{H_2}}+ 2 p_{\mathrm{H_2^+}},
\end{equation} where $\pi_{\mathrm{H^{tot}}}$ is the total pressure due to the hydrogen particles as if they were all free (atoms or ions), i.e. as if there were no molecules formed. This pressure is clearly fictitious when that conditions is not fulfilled. </div>
<div>
<br /></div>
The <b>particle conservation</b> in the pressure form can be rewritten as (Note that this is different from the first of Eq.(34) in Mihalas' article.):<br />
\begin{equation}f_1 + f_2 + f_3 + 2 f_4 + 2 f_5 = 1, \end{equation}<br />
or as<br />
\begin{equation}f_1\left(1 + \frac{f_2}{f_1} + \frac{f_3}{f_1} + 2\frac{f_4}{f_1} + 2\frac{f_5}{f_1}\right) = 1, \end{equation}<br />
<div>
where</div>
\begin{eqnarray}<br />
f_1&=&\frac{p_{\mathrm{H}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\<br />
f_2&=&\frac{p_{\mathrm{H^+}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\<br />
f_3&=&\frac{p_{\mathrm{H^-}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\<br />
f_4&=&\frac{p_{\mathrm{H_2^+}}}{\pi_{\mathrm{H^{tot}}}},\nonumber\\<br />
f_5&=&\frac{p_{\mathrm{H_2}}}{\pi_{\mathrm{H^{tot}}}}.<br />
\end{eqnarray}<br />
In addition we define:<br />
\begin{equation}<br />
f_{\mathrm{e}}=\frac{p_{\mathrm{e}}}{\pi_{\mathrm{H^{tot}}}}.<br />
\end{equation}<br />
<br />
<b>Ionization and dissociation constants</b><br />
<br />
The relations between the partial pressures are set by the Saha
ionization formula and by the instantaneous chemical equilibrium. These
relations are specified by the corresponding coefficients:<br />
\begin{eqnarray}<br />
K(H)&=&\frac{p_{\mathrm{H^+}}p_{\mathrm{e}}}{p_{\mathrm{H}}},<br />
\nonumber\\<br />
K(H^-)&=&\frac{p_{\mathrm{H}}p_{\mathrm{e}}}{p_{\mathrm{H^-}}},\nonumber\\<br />
K(H_2)&=&\frac{p_{\mathrm{H}^2}}{p_{\mathrm{H_2}}},\nonumber<br />
\nonumber\\<br />
K(H_2^+)&=&\frac{p_{\mathrm{H}} p_{\mathrm{H^+}}}{p_{\mathrm{H_2^+}}}.<br />
\end{eqnarray}<br />
For the sake of simplicity we introduce new variables:<br />
\begin{eqnarray}<br />
g_2 &=& \frac{K(\mathrm{H})}{p_{\mathrm{e}}} ,\nonumber\\<br />
g_3 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H})} ,\nonumber\\<br />
g_4 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2^+})} ,\nonumber\\<br />
g_5 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2}).} <br />
\end{eqnarray}<br />
The new variables can be easily related to the previously introduced $f_i$ ratios:<br />
\begin{eqnarray}<br />
g_2&=& \frac{K(H)}{p_{\mathrm{e}}}=\frac{p_{\mathrm{H^+}}}{p_{\mathrm{H}}} =\frac{f_2}{f_1}\nonumber\\<br />
g_3&=&\frac{p_{\mathrm{e}}}{K(H^-)} = \frac{p_{\mathrm{H^-}}}{p_{\mathrm{H}}}=\frac{f_3}{f_1},\nonumber\\<br />
g_4 &=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2^+})} = \frac{p_\mathrm{e}}{p_\mathrm{H}}\;\frac{p_\mathrm{H_2^+}}{p_\mathrm{H^+}} = \frac{f_\mathrm{e}}{f_1}\;\frac{f_4}{f_2}= \frac{f_\mathrm{e}}{f_1}\;\frac{f_4}{f_1}\;\frac{f_1}{f_2} = \frac{f_\mathrm{e}}{f_1}\;\frac{f_4}{f_1}\;\frac{1}{g_2},\nonumber\\<br />
&\Rightarrow& \frac{f_4}{f_1} = \frac{f_1}{f_\mathrm{e}}\,g_2\,g_4,\nonumber\\<br />
g_5&=& \frac{p_{\mathrm{e}}}{K(\mathrm{H_2})} = \frac{p_\mathrm{e}}{p_\mathrm{H}}\;\frac{p_\mathrm{H_2^+}}{p_\mathrm{H}} = \frac{f_\mathrm{e}}{f_1}\;\frac{f_5}{f_1},\nonumber\\<br />
&\Rightarrow& \frac{f_5}{f_1} = \frac{f_1}{f_\mathrm{e}}\,g_5.\nonumber\\<br />
\end{eqnarray}<br />
<br />
<br />
Now we can write the particle conservation as a function of only two variables $f_1$ and $f_{\mathrm{e}}$:<br />
\begin{equation}f_1(1 + g_2 + g_3) + 2 \frac{f_1^2}{f_{\mathrm{e}}} (g_2 g_4 + g_5) = 1, \end{equation}<br />
<br />
<b>Charge conservation</b><br />
<br />
Another equation required to close the system we find in the charge conservation law:<br />
\begin{equation}<br />
n_{\mathrm{e}} = n_{\mathrm{H^+}}-n_{\mathrm{H^-}}+n_{\mathrm{H_2^+}} + Q,<br />
\end{equation}<br />
where $Q$ stands for all the electrons freed by all single- and double-ionised elements other than hydrogen (index $i$ runs over the atomic elements other than hydrogen; index $j$ runs over the ionization stages):<br />
\begin{equation}<br />
Q = \sum_{z = 2}^{n_i} \alpha_i \sum_{i = j}^{2}n_{ji}.<br />
\end{equation}<br />
<div>
If we rewrite that in terms of pressures:</div>
\begin{equation}<br />
p_{\mathrm{e}} = p_{\mathrm{H^+}}-p_{\mathrm{H^-}}+p_{\mathrm{H_2^+}} + Q,<br />
\end{equation}<br />
or after dividing both sides with $\pi_{\mathrm{H^{tot}}}$<br />
\begin{equation}<br />
f_{\mathrm{e}} = f_2 - f_3 + f_4 + Q.<br />
\end{equation}<br />
<div>
Similar to the particle conservation law, the charge conservation can be rewritten as a function solely of $f_1$ and $f_{\mathrm{e}}$:</div>
<div>
\begin{equation} f_{\mathrm{e}} = f_1 (g_2 - g_3) + \frac{f_1^2}{f_{\mathrm{e}}} g_2 g_4 + Q.\end{equation}</div>
<div>
<br />
<b>Solution</b> <br />
<br />
To solve the system we introduced the following substitutes:<br />
\begin{eqnarray}
a&=&(1 + g_2 + g_3),\nonumber\\
b&=&2\left(1 + \frac{g_2 g_4}{g_5}\right),\nonumber\\
c&=&g_5,\nonumber\\
d&=&g_2 - g_3,\nonumber\\
e&=&\frac{g_2 g_4}{g_5}.
\end{eqnarray}
</div>
<div>
To simplify further derivation, the particle conservation we rewrite as </div>
<div>
\begin{equation}a\, f_1 + c\,b\, \frac{f_1^2}{f_{\mathrm{e}}} = 1, \end{equation}<br />
and the charge conservation as<br />
\begin{equation}f_{\mathrm{e}} = d\,f_1 + e\,c\, \frac{f_1^2}{f_{\mathrm{e}}} + Q.\end{equation}</div>
<div>
Solve Eqs.(19...) and (20) for $f_{\mathrm{e}}$:<br />
\begin{equation} f_{\mathrm{e}} = \frac{c\,b\,f_1^2}{1-a\, f_1},
\end{equation}
\begin{equation} f_{\mathrm{e}}^2 = d\,f_1 f_{\mathrm{e}} + e\,c\,f_1^2 + Q f_{\mathrm{e}},
\end{equation}
and substitute $F_{\mathrm{e}}$ from Eq.(22..) into Eq.(23):<br />
\begin{equation}
\frac{c^2\,b^2\,f_1^4}{(1-a\, f_1)^2} = \frac{cbd\,f_1^3}{(1-a\, f_1)} + e\,c\,f_1^2 + \frac{cb\,Qf_1^2}{(1-a\, f_1^2)}
\end{equation}
We multiply the last expression with $(1-a\, f_1)^2/(c\,f_1^2)$:
\begin{eqnarray}
&\,& cb^2\,f_1^2 = bd\,f_1(1-a\, f_1) + e(1-a\, f_1)^2 + b\,Q(1-a\, f_1),\nonumber\\
&\Rightarrow& cb^2\,f_1^2 = bd\,f_1 -abd\, f_1^2 + e - 2aef_1 + ea^2 f_1^2 + bQ -abQ\,f_1,\nonumber\\
&\Rightarrow& (cb^2+abd-ea^2)\,f_1^2 + (2ae-bd+abQ)\,f_1 - (e + bQ) = 0,\nonumber\\
&\Rightarrow& c_1\,f_1^2 + c_2\,f_1 + c_3 = 0,\end{eqnarray}
</div>
<div>
where<br />
\begin{eqnarray}<br />
c_1 &=& cb^2+abd-ea^2,\nonumber\\<br />
c_2 &=& 2ae-bd+abQ,\nonumber\\<br />
c_3 &=&- (e + bQ).<br />
\end{eqnarray}<br />
The solution of the quadratic equation gives the value of $f_1$:<br />
\begin{equation}<br />
f_1 = \frac{-c_2 \pm \sqrt{c_2^2 - 4c_1c_3}}{2c_1}.<br />
\end{equation}<br />
<br />
<i><u><b>Which solution do we take?</b></u></i><br />
<br />
<b>Backward substitution</b><br />
<br />
Once we know $F_1$, we can evaluate the other $F$ ratios:<br />
\begin{equation}<br />
f_2 = g_2 f_1,<br />
\end{equation}<br />
\begin{equation}<br />
f_3 = g_3 f_1,<br />
\end{equation}<br />
Equation (6) can be written as:<br />
\begin{equation}<br />
a f_1 + (1 + e) f_5 = 1,<br />
\end{equation}<br />
<div>
thus for $f_5$ we have:</div>
<div>
\begin{equation}<br />
f_5 = \frac{1 - a f_1}{1 + e},<br />
\end{equation}</div>
<div>
and for $f_4$:</div>
\begin{equation}<br />
f_4 = \frac{g_2 g_4}{g_5} f_5 = e\,f_5.<br />
\end{equation}<br />
Finally, it is possible to compute $f_{\mathrm{e}}$:<br />
\begin{equation}<br />
f_{\mathrm{e}} = f_2 - f_3 + \frac{1}{2}f_4 + Q.<br />
\end{equation}<br />
<div>
<br /></div>
<div>
The total pressure of hydrogen-nuclei is then (Eq.(9)):</div>
<div>
\begin{equation}
\pi_{\mathrm{H^{tot}}} = \frac{p_{\mathrm{e}}}{f_{\mathrm{e}}},
\end{equation}</div>
<div>
and the partial pressures are</div>
<div>
\begin{eqnarray}
p_{\mathrm{H}} &=& f_1\pi_{\mathrm{H^{tot}}},\nonumber\\
p_{\mathrm{H^+}} &=& f_2\pi_{\mathrm{H^{tot}}},\nonumber\\
p_{\mathrm{H^-}} &=& f_3\pi_{\mathrm{H^{tot}}},\nonumber\\
p_{\mathrm{H_2^+}} &=& \frac{1}{2}\,\pi_{\mathrm{H_2^+}}= \frac{1}{2}f_4\pi_{\mathrm{H^{tot}}},\nonumber\\
p_{\mathrm{H_2}} &=& \frac{1}{2}\,\pi_{\mathrm{H_2}}= \frac{1}{2}f_5\pi_{\mathrm{H^{tot}}}.
\end{eqnarray}</div>
<div>
<br />
<b>Total pressure</b><br />
<br /></div>
<div>
The total pressure of hydrogen particles is therefore</div>
<div>
\begin{eqnarray} p_{\mathrm{H^{tot}}} &=& p_{\mathrm{H}} + p_{\mathrm{H^+}} + p_{\mathrm{H^-}}+p_{\mathrm{H_2}}+p_{\mathrm{H_2^+}} =<br />
\nonumber\\ &=&\left[f_1 + f_2 + f_3 + \frac{1}{2}(f_4 + f_5)\right]\pi_{\mathrm{H^{tot}}}. \end{eqnarray}
</div>
<div>
<br /></div>
Once we know $p_{\mathrm{H^{tot}}}$, the total gas pressure is<br />
\begin{eqnarray}<br />
p_{\mathrm{g}} &=& p_{\mathrm{e}} + p_{\mathrm{a}} = p_{\mathrm{e}} + p_{\mathrm{H^{tot}}}+ p_{\mathrm{nonH}}=p_{\mathrm{e}}\left[1 + \frac{p_{\mathrm{H^{tot}}} + p_{\mathrm{nonH}}}{p_{\mathrm{e}}}\right] .<br />
\end{eqnarray}</div>
Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0tag:blogger.com,1999:blog-1169051187754614539.post-8664995181368441732016-05-03T02:10:00.000-07:002017-04-19T02:11:39.399-07:00Effective 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). <br />
<br />
For each of the levels (u for upper, l for lower) in LS coupling holds:<br />
$$g_\mathrm{LS} =\frac{3}{2}+\frac{S (S+1) - L (L+1)}{2J(J+1)}, $$ <br />
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:<br />
$$\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)).$$ <br />
<br />
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.<br />
<br />
Example:<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); font-family: 'Andale Mono', fixed, monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>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
</code></pre>
<br />
<br />
Download: <a href="https://www.dropbox.com/s/lmp1nm30wli7y4v/lande_factor.pro?dl=0">lande_factor.pro</a>Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-37024363098682935082016-01-15T06:45:00.000-08:002017-04-19T06:55:00.024-07:00Kurucz' Atlas in GNU LinuxI have just found about great effort made by Sbordone, Bonifacio & Castelli to make the Atlas code of Bob Kurucz available at GNU Linux. <br />
<br />
There is a webpage: <a href="http://atmos.obspm.fr/index.php">http://atmos.obspm.fr/</a> providing <a href="http://atmos.obspm.fr/index.php/download">the code</a> and <a href="http://atmos.obspm.fr/index.php/documentation">the documentation</a>.<br />
<br />
<br />
<br />
<br />Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0tag:blogger.com,1999:blog-1169051187754614539.post-31611140922466501852015-07-12T21:08:00.000-07:002015-07-25T14:49:42.461-07:00Notes on some basic quantities, units and constants: Part II<b>Mass fractions</b><br />
<br />
The mass fraction of H, He and
the metals is a quick way to specify the chemical composition of the
stellar plasma. The mass fraction of hydrogen is defined as ratio of
the mass of hydrogen particles and the total mass of all particles (in a
given volume):
$$
X = \frac{M_\mathrm{H}}{M}= \frac{\rho_\mathrm{H}}{\rho},
$$
where $\rho_\mathrm{H}$ and $\rho$ are the respective (mass) densities.
The hydrogen mass density is equal to the hydrogen number density
($n_\mathrm{H}$) times the mass of one hydrogen particle $m_{\mathrm{H}}
= A_{\mathrm{H}}\,m_{\mathscr{A}}$, thus:
$$
X = \frac{M_\mathrm{H}}{M}= \frac{n_\mathrm{H}\,m_{\mathrm{H}}}{\rho}=
\frac{n_\mathrm{H}\,A_{\mathrm{H}}\,m_{\mathscr{A}}}{\rho}.
$$
Similar to that, for helium and the metals we define:
$$
Y = \frac{M_\mathrm{He}}{M}=
\frac{n_\mathrm{He}\,m_{\mathrm{He}}}{\rho}=
\frac{n_\mathrm{He}\,A_{\mathrm{He}}\,m_{\mathscr{A}}}{\rho},
$$
$$
Z = \frac{M_\mathrm{metals}}{M}= \frac{n_\mathrm{metals}\,m_{\mathrm{metals}}}{\rho} =
\frac{\sum_{i=3} n_\mathrm{i}\,A_{\mathrm{i}}\,m_{\mathscr{A}}}{\rho},.
$$<br />
<a name='more'></a>
On the other hand, in the same volume there is the total number of $N$
particles (free and bound atoms and free electrons):
$$
N = N_\mathrm{e} + N_{\mathrm{a}} = N_\mathrm{e} + N_\mathrm{H} +
N_\mathrm{He} + N_\mathrm{metals},
$$
where the indices $\mathrm{e}$ and $\mathrm{a}$ stand for the electrons
and for the atoms.
The same equation for the number densities reads:
$$
n = n_\mathrm{e} + N_{\mathrm{a}} = n_\mathrm{e} + n_\mathrm{H} +
n_\mathrm{He} + n_\mathrm{metals}.
$$
For the mass of one particle of the species $i$ ($i \ge 1$) we have:
$$
m_{i} = A_{i}\,m_{\mathscr{A}},
$$
or expressed relative to hydrogen:
$$
m_{i} =\frac{A_{i}\,m_{\mathrm{H}}}{A_{\mathrm{H}}}.
$$
Since $m_{\mathrm{H}} \approx m_{\mathscr{A}}$, sometimes in the
literature the last equation appears as:
$$
m_{i} =\frac{A_{i}\,m_{\mathscr{A}}}{A_{\mathrm{H}}}.
$$
<br />
<br />
<b>Abundances</b><br />
<br />
Abundances specify
detailed chemical composition of a stellar atmosphere. There are three
common ways to express them in astrophysics: <br />
<br />
1)
$\alpha_i$ - relative number of the atoms (nuclei) of the element Z
relative to the number of hydrogen (nuclei) in the linear scale (for hydrogen, $\alpha_\mathrm{H} \equiv 1$): \begin{equation}\alpha_i = \frac{N_i}{N_\mathrm{H}} = \frac{n_i}{n_\mathrm{H}}\end{equation}; <br />
<br />
2)
$\varepsilon_{i}$ - abundances relative to the amount of hydrogen in the
logarithmic scale normalized so that the abundance of hydrogen
$\varepsilon_{\mathrm{H}} = 12$:\begin{equation}\varepsilon_i = 12 + \log\alpha_i;\end{equation}<br />
<br />
3) $x_i$ - relative number of the nuclei of the element Z relative to the total number of nuclei in the linear scale: \begin{equation}
x_i = \frac{N_i}{N_\mathrm{a}} = \frac{N_i}{\sum N_j} = \frac{10^{\varepsilon_i}}{\sum 10^{\varepsilon_j}}. \end{equation}
For hydrogen for the solar composition, $x_\mathrm{H} \approx 0.9$.<br />
<br />
The
relations between $x_i$ and $\alpha_i$ are the following:
$$
x_i = \frac{N_i}{\sum N_j} = \frac{\alpha_i N_{\mathrm{H}}}{\sum
\alpha_j N_{\mathrm{H}}}= \frac{\alpha_i}{\sum \alpha_j},
$$
$$
\alpha_i = \frac{N_i}{N_\mathrm{H}} = \frac{x_i
N_\mathrm{n}}{N_\mathrm{H}}.
$$
Finally, sometimes there is a need to express the individual abundances
of non-hydrogen species relative to their total number:
$$
(x_i)_{\mathrm{nonH}} = \frac{N_i}{\sum_{j\gt 1} N_j} =
\frac{N_i}{-N_\mathrm{H} + \sum_{j\gt 1} N_j} = \frac{\alpha_i
N_\mathrm{H}}{-N_\mathrm{H} + \sum_{j\gt 1} \alpha_j N_\mathrm{H}} =
\frac{\alpha_i }{-1 + \sum_{j\gt 1} \alpha_j }.
$$
<br />
<br />
<b>Mass fractions in terms of abundances</b><br />
<br />
For
any $X_i$ ($X_i \in [X, Y, Z]$) we can write:
$$
X_i = \frac{n_i m_i}{\rho} = \frac{\alpha_i m_i n_\mathrm{H}}{\rho}=
\frac{\alpha_i m_i n_\mathrm{H}}{\sum \alpha_j m_j n_\mathrm{H}}=
\frac{\alpha_i m_i}{\sum \alpha_j m_j}.
$$
From $m_i = A_i m_\mathrm{H}/ A_\mathrm{H}$ follows:
\begin{equation}
X_i = \frac{\alpha_i A_i}{\sum \alpha_j A_j}.
\label{eq:xyzabund}
\end{equation}
Therefore:
$$
X = \frac{\alpha_\mathrm{H} A_\mathrm{H}}{\sum \alpha_j
A_j}\;\;\;\;\;\;\;Y = \frac{\alpha_\mathrm{He} A_\mathrm{He}}{\sum
\alpha_j A_j}\;\;\;\;\;\;\;Z = \frac{\sum_{j>2} \alpha_j A_j}{\sum
\alpha_j A_j}.
$$
<br />
<br />
<b>Pressure and density in terms of number densities</b><br />
<br />
The
total gas pressure may be represented as a sum of the partial gas
pressures of the particles of different species (possible species are
now atoms (neutrals and ions), electrons and molecules). It is Dalton's
law:
$$
p \equiv p_{\mathrm{tot}} = p_{\mathrm{e}} +
p_{\mathrm{a}^{\mathrm{free}}} + p_{\mathrm{mol}}.
$$
If the ideal gas approximation is applicable, for each component of the
mixture it holds:
$$
p = n\,kT = (n_{\mathrm{e}} + n_{\mathrm{a}}^{\mathrm{free}} +
n_{\mathrm{mol}})\,kT,
$$
where $n_{\mathrm{a}}^{\mathrm{free}}$ represents the sum over all
chemical elements and
$n_{\mathrm{mol}}$ the sum over all molecules.
<br />
<br />
If there is no molecules formed ($n_{\mathrm{a}}^{\mathrm{free}}=n_{\mathrm{a}}$):
$$
p = n\,kT = (n_{\mathrm{e}} + n_{\mathrm{a}})\,kT.
$$
Therefore the total number density of the atoms may be found when the total gas pressure and the electron pressure are known:
$$
n_{\mathrm{a}} = (p - p_{\mathrm{e}})/kT.
$$
<br />
The total mass density is
$$
\rho = \rho_{\mathrm{e}} + \rho_{\mathrm{a}}.
$$
This definition hold no matter if the molecules are formed or not. Since
$\rho_{\mathrm{e}} \ll \rho_{\mathrm{a}}$ it is commonly taken that
$$
\rho \approx \rho_{\mathrm{a}} = \rho_{\mathrm{H}} + \rho_{\mathrm{He}} +
\rho_{\mathrm{m}}.
$$
<br />
Note that pressure is proportional to the number of the particles
while the density is proportional to their mass. As a clear
consequence, the increasing ionization rate (the increasing number of
free electrons) increases the total pressure but not the density.<br />
<br />
<br />
<b>Mean molecular weight</b><br />
<br />
The mean molecular weight is primarily used in chemistry where it is known as <i>molar mass</i>. In chemistry it is defined as <i>the mass of a given substance divided by its amount of substance</i>:<br />
$$
\mu = \frac{M [\mathrm{g}]}{\widetilde{n} [\mathrm{mol}]}.
$$ <br />
If there is, for an example, one mole ($\widetilde{n} = 1 $)
of neutral hydrogen gas, there
is $N_\mathscr{A}$ hydrogen atoms in it. Total mass of these particles
is $
M = N_\mathscr{A} m_{\mathrm{H}}$, where $m_{\mathrm{H}} =
A_{\mathrm{H}} m_\mathscr{A}$ is the mass of one H atom in grams. Numerically,
$N_\mathscr{A}\, m_\mathscr{A}= 1$, and therefore the total mass $M$ is equal
to the value of $A_{\mathrm{H}}$ (but the unit is still gram, not a.m.u.!). The molar mass is simply
that mass divided by the amount of substance, $\widetilde{n} =
1\,\mathrm{mol}$, thus $\mu = 1\, \mathrm{g/mol}$. If the gas is
completely ionized, its mass remains the same, while the amount of
substance (i.e. the number of particles) doubles, so that $\mu = 1/ 2
\,A_{\mathrm{H}}$. <br />
<br />
Equivalently, the $\mu$ can be defined as the mean mass of a
particle in atomic mass units. Let's again consider neutral hydrogen gas. There is only one type of particles in that case,
they are all equal and, therefore, the mean molecular weight is equal to
the mass of one particle in a.m.u., $\mu = A_{\mathrm{H}} \approx 1\,
\mathrm{a.m.u.}$. If the gas gets fully ionized, the total mass in unchanged, the
number of particles is doubled and, clearly, $\mu = 1/2
A_{\mathrm{H}}$. Formally this can be written as: $$
\mu = \frac{\sum N_i\,A_i [\mathrm{a.m.u.}]}{N},
$$ where $N_i$ is the number of particles of type $i$, $N$ is the
total number of particles and $A_i$ is the atomic weight of one particle of type $i$
(in a.m.u.).<br />
<br />
There is often need to relate $\mu$ to the density. The following relations hold:<br />
<br />
\begin{eqnarray}
\mu &=& \frac{M [\mathrm{g}]}{\widetilde{n} [\mathrm{mol}]} =\\
&=& \frac{M [\mathrm{g}]}{N}\,N_\mathscr{A}
[\mathrm{mol^{-1}}] = \\ &=& \frac{M [\mathrm{g}] / V
[\mathrm{cm^3}]}{N / V [\mathrm{cm^3}]}\,N_\mathscr{A}
[\mathrm{mol^{-1}}]= \\&=& \frac{\rho [\mathrm{g\,cm^{-3}}]}{n
[\mathrm{cm^{-3}}]}\,N_\mathscr{A} [\mathrm{mol^{-1}}].\end{eqnarray}<span style="font-size: large;">☞</span>a.m.u. is (in terms of writing equations and crunching numbers) not much
more special then any imperial unit, pound or ounce for example.
Avogadro's number can be taken as the conversion factor between grams
(or kilograms) and a.m.u. <br />
<br />
<b>Mean molecular weight of a mixture</b><br />
<br />
Now
that we understand all about $\mu$ and its units, we take the last of
the the relations from the previous section and apply it to the gas that consists of different
atoms, neutral and ionized, free electrons and molecules. That relation
says that $\mu$ depends upon the density and the number of particles.
The density remains unchanged in the basic micro-processes in the
stellar atmospheres (ionization, recombination, molecular association
and disassociation). On the other hand, the number of the free particles
(atoms, electrons, molecules) changes in these processes and $\mu$
follows that change as well. The actual chemical composition enters
$\mu$ through the density. Therefore, $\mu$ is a single number that
tells us about the chemical composition if the ionization fractions are
known or about the ionization fraction if the chemical composition is
fixed. An obvious question is if we can split the two parts of $\mu$ -
one dependent only on the number of the free electrons and other
dependent only on the chemical composition.<br />
<br />
Let's
assume that there is no molecules in the plasma and write the density
as a sum per specie:
$$
\rho = \sum_{i = 1}^{i_\mathrm{max}} \rho_i = \rho_\mathrm{H} +
\rho_\mathrm{He} + \rho_\mathrm{m},
$$
where we neglected the density of the free electrons ($\rho_\mathrm{e}
\ll \rho_i$ for ant $i$). The total number density explicitly written
is:
$$
n= n_\mathrm{e} + \sum_{i = 1}^{i_\mathrm{max}} n_i = n_\mathrm{e} +
n_\mathrm{H} + n_\mathrm{He} + n_\mathrm{metals}.
$$
The mean molecular weight is then:
$$
\mu = \frac{\rho_\mathrm{H} + \rho_\mathrm{He} +
\rho_\mathrm{metals}}{(n_\mathrm{e} + n_\mathrm{H} + n_\mathrm{He} +
n_\mathrm{metals})\,m_\mathscr{A}},
$$
and its inverse:
$$
\frac{1}{\mu} = \frac{n_\mathrm{e}\,m_\mathscr{A}}{\rho} +
\frac{n_\mathrm{H}\,m_\mathscr{A}}{\rho} +
\frac{n_\mathrm{He}\,m_\mathscr{A}}{\rho} +
\frac{n_\mathrm{metals}\,m_\mathscr{A}}{\rho}.
$$
In the three atomic terms we recognize expressions proportional to the
mass fractions defined above:
$$
\frac{1}{\mu} = \frac{n_\mathrm{e}\,m_\mathscr{A}}{\rho} + X
\frac{m_\mathscr{A}}{m_\mathrm{H}} +
Y\frac{m_\mathscr{A}}{m_\mathrm{He}} +
Z\frac{m_\mathscr{A}}{m_\mathrm{m}}.
$$
The first term:
$$
\frac{1}{\mu_\mathrm{e}} = \frac{n_\mathrm{e}\,m_\mathscr{A}}{\rho}
$$
depends only on the number of the free electrons. The second term
depends only on the chemical composition:
$$
\frac{1}{\mu_\mathrm{a}} = X \frac{m_\mathscr{A}}{m_\mathrm{H}} +
Y\frac{m_\mathscr{A}}{m_\mathrm{He}} +
Z\frac{m_\mathscr{A}}{m_\mathrm{m}}.
$$
If we replace $m_\mathrm{H}$ in the first term with
$A_\mathrm{H}\,m_\mathscr{A}$ and approximate $A_\mathrm{H}\approx 1$,
that term reduces to $X$ only. In the other two terms we use $m_i =
A_i\,m_\mathrm{H}/A_\mathrm{H}$ to get:
$$
\frac{1}{\mu_\mathrm{a}} = X + Y\frac{A_\mathrm{H}}{A_\mathrm{He}} +
Z\frac{A_\mathrm{H}}{A_\mathrm{m}}.
$$
A similar expression employing the abundances instead of the mass
fractions can easily be derived using the relation Eq.\ref{eq:xyzabund}
between the mass fractions and the abundances that we can write for
separately for any element under the consideration:
$$
X_i = \frac{\alpha_i A_i}{\sum \alpha_j A_j}, \;\;\;\;\;\;\;\mathrm{1
\le i \le i_\mathrm{max}}.
$$
For $1/\mu_\mathrm{a}$ now we have:
$$
\frac{1}{\mu_\mathrm{a}} = \sum_i X_i \frac{A_\mathrm{H}}{A_i} = \sum_i
\frac{\alpha_i A_i}{\sum \alpha_j A_j} \frac{A_\mathrm{H}}{A_i} =
\frac{A_\mathrm{H} \sum_i \alpha_i }{\sum \alpha_j A_j},
$$
and finally:
$$
\mu_\mathrm{a} = \frac{\sum \alpha_j A_j}{A_\mathrm{H} \sum_j \alpha_j
}.
$$
This expression, for example, appears in the important article by
Mihalas (1967, in Methods in Computational Physics, eds. Alder, Fernbach
& Rotenberg, p.16, eq.54).
<br />
<br />
<br />
<b>Different forms of the ideal gas law</b><br />
<br />
The
ideal gas law is one of the most useful and widely applicable
approximations in astrophysics. I already used it to introduce the total
and the partial pressure:
$$
p = n\,kT = (n_{\mathrm{e}} + n_{\mathrm{a}})\,kT.
$$
This form comes from statistical mechanics. The constant $k$ is the
Boltzmann constant that is equal to the universal gas constant ($R$)
divided by Avogardo's number ($N_\mathscr{A}$). Other useful forms can
be easily derived:
$$
p = n\,kT = \frac{N}{V}\,\frac{R}{N_\mathscr{A}}\, T =
\frac{N}{N_\mathscr{A}}\,\frac{RT}{V} = \frac{\widetilde{n}RT}{V},
$$
or
$$
p = \frac{\widetilde{n}RT}{V} = \frac{M}{\mu}\,\frac{RT}{V} = \frac{\rho
RT}{\mu}.
$$
Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-64017155134895838682015-07-11T14:40:00.000-07:002015-07-25T14:50:19.278-07:00Notes on some basic quantities, units and constants: Part IThere is a definition of mole as a unit that every student learns at very elementary level. Although mole is one of the
seven base units of the International System, it is a bit specific and sometimes creates confusion. Here is my attempt
to clarify the concept of mole and to derive in one place some useful
relations used in the radiative transfer and atmospheric modeling. The first version of these notes I wrote for a group of students at the University of Belgrade many years ago. I still use them as a personal reminder. <br />
<br />
<b>Counting "Elementary particles"</b><br />
<br />
For the solar/stellar plasma, the "elementary" particles are atoms, ions (positive or negative), free electrons and molecules. In the very cool atmospheres there are dust particles as well. In the solar atmosphere dust can be completely neglected. Regarding the chemical composition, the atmospheric plasma is made out of hydrogen, helium and the metals (all other elements). There is no nuclear reactions and thus the total number density of nuclei per atomic specie is constant with time.<br />
<br />
It is important to distinguish between the number of free atoms and the total number of atoms including those bound in the molecules. The former I denote as $N_{\mathrm{a}}^{\mathrm{free}}$, the latter as $N_{\mathrm{a}}^{\mathrm{tot}}$. The total number of atoms is identical to the number of atomic nuclei. The total number of molecules is $N_{\mathrm{m}}$.<br />
<br />
The total number of particles $N$ is therefore:<br />
$$N = N_\mathrm{e} + N_{\mathrm{a}}^{\mathrm{free}} + N_{\mathrm{m}},$$
or when there is no molecules $$N = N_\mathrm{e} + N_{\mathrm{a}}^{\mathrm{free}} = N_\mathrm{e} + N_\mathrm{H} + N_\mathrm{He} + N_\mathrm{\mathrm{metals}},$$ where $\mathrm{e}$, $\mathrm{m}$, $\mathrm{H}$ and $\mathrm{He}$ stand for the electrons, the molecules, hydrogen, helium and $\mathrm{metals}$ refer to all other elements together. The contribution of the metals can be further divided into the contributions of the individual elements.<br />
<br />
<a name='more'></a><br />
When there is no molecules formed in the plasma, the number of the free atoms (neutral and ions) is equal to the total number of atoms $N_{\mathrm{a}} \equiv N_{\mathrm{a}}^{\mathrm{free}} = N_{\mathrm{a}}^{\mathrm{tot}}$. Whenever this is the case I will simply use $N_{\mathrm{a}}$. In addition I will use $N_{\mathrm{a}}$ instead of $N_{\mathrm{a}}^{\mathrm{tot}}$ when the presence of molecules is not relevant (e.g. the total mass is the same if all the atoms are bound in the molecules as if they were all free).<br />
<br />
With small $n$ I will denote all respective number densities ($V$ is the unit density):
$$n = \frac{N}{V}.$$<b> </b><br />
<b>A brief but important reminder</b><br />
<br />
Any quantity in physics is characterized by its value, its dimension and its unit. For example, dimension of the force in $m \times l \times t^{-2}$, where $m$, $l$ and $t$ are mass, length and time. The standard unit for the force in the SI system is therefore $\mathrm{kg} \times \mathrm{m} \times \mathrm{s}^{-2}$ and it has its own name - Newton. There can be several different units for the same quantity. Different units correspond to different basic systems of units (e.g. CGS and SI units for the force, dyne and Newton) or they just represent smaller or larger versions of the basic units (e.g. in the SI system, mega-meter, pico-meter and so on). However, there are also dimensionless quantities (or quantities with dimension 1), either because they are defined as a ratio of two quantities with the same dimension (e.g. Reynolds number) or because they count discrete features (e.g. number of particles, number of dogs). Dimensionless quantities are, in many cases, also unitless. In principle, we can assign units to them (e.g. 'particles', 'dogs'), but these are not physical units and they should be ignored in the dimensional analysis.<br />
<br />
If a quantity $q$ can be expressed in units $[U_1]$ and $[U_2]$ and if, in these units, it takes values $v_1$ and $v_2$, respectively, then: $$q = v_1\,[U_1] = v_2\,[U_2].$$Since both sides of the equation have the same dimension, there must be a dimensionless number $c_{12}$ such that $[U_1] = c_{12}\,[U_2]$. It is the conversion factor between the two units. Formally, its unit is $[U_1]/[U_2]$. For example, length can be expressed in miles and in kilometers:<br />
$$l = 3\, \mathrm{miles} = 4.82803\,\mathrm{km}$$.<br />
The conversion factor is<br />
$$c = 1.60934 \,\frac{\mathrm{km}}{\mathrm{mile}}.$$ <br />
<b><b>Avogadro's number</b></b><br />
<br />
The simplest definition of a mole says that mole stands for a number. Not for any number, but for a very specific one - <b>Avogadro's number</b>: $$N_{\mathscr{A}} = 6.022140857(74) \times 10^{23}.$$ One mole of hydrogen atoms contains $N_{\mathscr{A}}$ particles of hydrogen atoms. One mole of electrons contains $N_{\mathscr{A}}$ electrons. One mole of dogs contains $N_{\mathscr{A}}$ dogs, etc. <i>The number of moles</i> (also known as<i> </i><b>the amount of substance</b>) $\widetilde{n}$ is
thus defined as the number of free particles divided by Avogadro's number:<br />
$$
\widetilde{n} = \frac{N}{N_{\mathscr{A}}} = \frac{N_{\mathrm{e}} +
N_{\mathrm{a}}^{\mathrm{free}} +
N_{\mathrm{m}}}{N_{\mathscr{A}}}.
$$
The unit of $\widetilde{n}$ is <b>mole</b>. To make it consistent, the unit of Avogadro's number is formally mole<sup>-1</sup>. Since $\widetilde{n}$ is the ratio of two dimensionless numbers, the mole is dimensionless as well. However, note that although both quantities are dimensionless, the amount of substance has a unit, while number of particles is unitless. In principle, we could introduce "particle" as a formal unit of the number of particles. If we do so, the unit of $N_{\mathscr{A}}$ would become particle/mole. In the context of the previous section, $N_{\mathscr{A}}$ is the conversion factor between the number of particles expressed in moles and the number of particles expressed in particles.<br />
<br />
The main purpose of $N_{\mathscr{A}}$ is to help us work with the
large numbers that are typical when we solve equations with quantities
of both micro and macro world.<br />
<br />
<b>Atomic mass unit</b> <br />
<br />
The reciprocal value of $N_{\mathscr{A}}$ is $$m_{\mathscr{A}} = \frac{1}{N_{\mathscr{A}}} = 1.660538921(73)\times
10^{-24}.$$ It gives us the fraction of one mole that corresponds to one particle. For example, if the mass of one mole of particles is 3 g, than the mass of one particle is obviously $m_{\mathscr{A}} \times 3\,\mathrm{g}$. This leads us to the concept of the <b>atomic mass unit</b> (a.m.u.). The a.m.u. is a unit for the mass used in atomic physics. It is not essentially different from the pound used in the imperial system of units. One pound is 453.59237 grams. One a.m.u. is $m_{\mathscr{A}}$ grams:$$1\,\mathrm{a.m.u.} = 1.660538921(73)\times
10^{-24}\,\mathrm{g}.$$ The quantity measured in a.m.u.'s is known as <b>atomic weight</b>, while the quantity measured in grams is the mass itself. The two are simply related as $$m = A\, m_{\mathscr{A}}.$$Therefore, $m_{\mathscr{A}}$ is the conversion factor between mass and atomic weight. It is a dimensionless number; its formal unit is a.m.u./gram. <br />
<br />
Both $m_{\mathscr{A}}$ and $N_{\mathscr{A}}$ are conversion factors but between two quantities of different dimension (mass and number of particles, respectively), and thus the formal dimensions of $m_{\mathscr{A}}$ and $N_{\mathscr{A}}$ differ as well.<br />
<br />
Some additional notes:<br />
<br />
1. In chemistry, atomic weight is often used as it is unitless.<br />
<br />
2. Atomic mass unit (with label a.m.u.) had been replaced in 1961 by unified atomic mass unit (with label u). However, the difference between the two is only in the value of $m_{\mathscr{A}}$ with the rest of the concept unchanged.<br />
<br />
3. Another name for the unified atomic mass unit is Dalton (Da). <br />
<br />
4. When the numbers are dimensionless but still have units, one has to be very careful with dimensional analysis. For example, all the following equations are numerically correct: $$1\,\mathrm{a.m.u.} = m_{\mathscr{A}} \times 1\, \mathrm{g} = \frac{1\,\mathrm{g}}{N_{\mathscr{A}}} = \frac{1\,\mathrm{g}}{1\,\mathrm{mol}}.$$For example for hydrogen all the following forms are common and correct:$$A_\mathrm{H} = 1.008 = 1.008\,\mathrm{a.m.u.} = 1.008\,\frac{\mathrm{g}}{\mathrm{mol}}.$$<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjM3YLYG3rWdYa6veuFZZllEu3SSIEoqgcf7EEqie72IgqzeXWsns3bHe31E4xl_z541JY0zCLGBHUYYdwRZ3X15QrRlZq_55oyrTvguu8c9iUTXGD5LaJzxw-9mW7xZDu96QFwTLlgfy4/s1600/fontaine1.jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjM3YLYG3rWdYa6veuFZZllEu3SSIEoqgcf7EEqie72IgqzeXWsns3bHe31E4xl_z541JY0zCLGBHUYYdwRZ3X15QrRlZq_55oyrTvguu8c9iUTXGD5LaJzxw-9mW7xZDu96QFwTLlgfy4/s1600/fontaine1.jpeg" /></a></div>
<br />
<b>The formal definition</b><br />
<br />
The lecture-book <b>definition of the mole </b>says that: 1 mole is the amount of substance in 12 grams of <sup>12</sup>C (the isotope of carbon with the atomic mass $A_{\mathrm{C}} = 12\;\mathrm{a.m.u.}$).<br />
<br />
This definition specifies the standard particle (one <sup>12</sup>C atom) which mass is the unified atomic mass unit (pre-1961 atomic mass unit was defined based on oxygen) and which number in 12 grams is exactly Avogadro's number. In other words, in 12 grams of <sup>12</sup>C there is exactly
$N_{\mathscr{A}}$ particles with atomic mass $12 \times \mathrm{a.m.u.}
=$ $12 \times (1 / N_{\mathscr{A}})$ grams.Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-74667246383493534622015-05-15T09:30:00.000-07:002020-02-14T05:52:50.257-08:001D Solar Models in IDL: Spruit's Convection ModelThe semi-empirical models of the solar atmosphere rely on the observed intensities and, therefore, they cannot say much about the invisible convection zone below the surface. Spruit (<a href="http://adsabs.harvard.edu/abs/1974SoPh...34..277S">1974SoPh...34..277S</a>) constructed a 1D model of the solar convection zone using the mixing-length theory with 4 free parameters. The model is constructed so that it matches the HSRA model (Gingerich et al, <a href="http://adsabs.harvard.edu/abs/1971SoPh...18..347G">1971SoPh...18..347G</a>) of the solar atmosphere. The model of Spruit is sometimes used to initiate the convection simulations with 3D numerical codes.<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj07QM1Sln_8YVZiXplzfeCzzLQvDuFKReoLdykgR5hd9nRvqXDDZeX2C1MUnF2RX0AMxEmUDY5JczE2CUcZHqQSbefk-Tt-yonclmLg1YKV38gCSuxEiKBSi4q3XwJoFkbv2m13HbA4O8/s1600/spruit_hsra.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj07QM1Sln_8YVZiXplzfeCzzLQvDuFKReoLdykgR5hd9nRvqXDDZeX2C1MUnF2RX0AMxEmUDY5JczE2CUcZHqQSbefk-Tt-yonclmLg1YKV38gCSuxEiKBSi4q3XwJoFkbv2m13HbA4O8/s400/spruit_hsra.png" width="400"></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Fig 1. The HSRA model atmosphere (dashed) and Spruit's model of convective zone (solid). This is reproduced after Fig.3 of Spruit's paper.</td></tr>
</tbody></table>
<br />
<a name='more'></a><br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhM-cuE6B2A1oNQOdURCyRjIre5dPD7Jz-ddgK2DTgk67grc36gZdkajjKHm4K93NGq-mniS-uEfmgHU33SQZWTfC7tXbjQ-s1cs8-7q2afbwPck1oFOiEC1UmWyt64-AJyYwGwlijfr8c/s1600/spruit_hsra_all.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhM-cuE6B2A1oNQOdURCyRjIre5dPD7Jz-ddgK2DTgk67grc36gZdkajjKHm4K93NGq-mniS-uEfmgHU33SQZWTfC7tXbjQ-s1cs8-7q2afbwPck1oFOiEC1UmWyt64-AJyYwGwlijfr8c/s400/spruit_hsra_all.png" width="400"></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Fig 2. Same as Fig1, just the entire models are shown. </td></tr>
</tbody></table>
<br />
<br />
The
data in the files below is Table II of Spruit's paper scanned, OCR-ed and manually double checked. I did some small changes to match the format of other models. The geometrical height scale is in km with zero at $tau_{500}=1$. All values are negative because this model covers only the layers below the surface. The variables in the original table are (the names in parentheses refer to the variable names in my files):<br />
<ul>
<li>depth relative to $tau_{500}=1$ (h),</li>
<li>column mass (m),</li>
<li>gas pressure (pgas),</li>
<li>temperature (t),</li>
<li>electron pressure (pel),</li>
<li>density (rho),</li>
<li>specific heat per mass (cp), </li>
<li>adiabatic gradient (agrad),</li>
<li>mean molecular weight (mu),</li>
<li>Rosseland mean opacity (kappa),</li>
<li>radiative gradient (rgrad),</li>
<li>mean gradient (mgrad),</li>
<li>convective velocity (v),</li>
<li>electrical conductivity (s33).</li>
</ul>
There are two
versions: .dat (ascii) and .sav (IDL). In the IDL files the
model is stored in a structure variable called spruit with the following tags:<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); color: black; font-family: Andale Mono,fixed,monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>spruit = {h:h, m:m, pgas:pgas, pel:pel, t:t, rho:rho, cp:cp, agrad:agrad,
mgrad:mgrad, rgrad:rgrad, mu_mu, kappa:kappa, v:v, s33:s33,
variables:variables, reference:reference}
</code></pre>
where the variables tag has two self-explanatory tags, varnames and varunits. All units are in CGS, except depth that is in km. Both files are all in a tarball:<br />
<br />
<a href="https://db.tt/wQmRehMD">spruit.tar</a><br />
<br />
<b>Disclaimer</b>: This data is publicly available in Spruit (1974, <a href="http://adsabs.harvard.edu/abs/1974SoPh...34..277S">1974SoPh...34..277S</a>).
If you use the data in a publication please cite it appropriately. I
double-checked that the digitized version here matches the original,
still typos are possible no matter how unlikely. If you have any doubts,
compare the numbers to those in the paper. Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com1Santa Cruz de Tenerife, Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-80747783813635451912015-04-12T14:05:00.000-07:002015-04-16T04:11:24.318-07:00Atlantis over Tenerife Japanese astronaut Sochi Noguchi took this <a href="http://spaceflight.nasa.gov/gallery/images/station/crew-23/hires/iss023e044581.jpg">iconic picture</a> of the <a href="http://en.wikipedia.org/wiki/Space_Shuttle_Atlantis">Space Shuttle Atlantis</a> over the Atlantic Ocean and the island of Tenerife. The picure was taken from the <a href="http://en.wikipedia.org/wiki/International_Space_Station">International Space Station</a> on Sunday, May 16 2010 at 10:28 am EDT (1428 GMT) while the Atlantis was getting ready for docking. At that time the ISS was about 350 km above the ground. This was the flight before the last one for the shuttle. In the highest resolution, the white towers of the solar telescopes at the <a href="http://en.wikipedia.org/wiki/Teide_Observatory">Observatorio del Teide</a> (<a href="http://www.iac.es/eno.php?op1=3&lang=en">IAC</a>) are visible between the clouds. To help your eyes, I took a snapshot from the Google Maps: look for the distinctive dark patches in the NASA picture, the towers are tiny white dots just above them.<br />
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiRrlel6xxW6oSTfflBFOrAmuWmd5dEIdzWfl0u7Dz3e83iRzBmClNN4UtBVhdIQ4etoi9QI8XMNY7XOowBMCQjbo6VXBe8bPjNw6VeA8vLDRtUzTn6wxl8B828nmBrk1QQ6cSSzxMl-ew/s1600/iss023e044581.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiRrlel6xxW6oSTfflBFOrAmuWmd5dEIdzWfl0u7Dz3e83iRzBmClNN4UtBVhdIQ4etoi9QI8XMNY7XOowBMCQjbo6VXBe8bPjNw6VeA8vLDRtUzTn6wxl8B828nmBrk1QQ6cSSzxMl-ew/s1600/iss023e044581.jpg" width="550" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><a href="http://spaceflight.nasa.gov/gallery/images/station/crew-23/hires/iss023e044581.jpg">Credit: NASA/Sochi Noguchi (click for hi-res)</a></td></tr>
</tbody></table>
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh9OBvgmEDorecSoGOSyr1RT_jXvs3rHtvM8o9K-RPwPgLrjXOoIa7sijCTBbib5RVLUaObvTTpaewDKPSGJp8g-e_FEUMwjDlEO3rxXHwu6513G5MsCVUPLPAiEKVHIa7qAb-uqRFl6Kk/s1600/izana_google_maps.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh9OBvgmEDorecSoGOSyr1RT_jXvs3rHtvM8o9K-RPwPgLrjXOoIa7sijCTBbib5RVLUaObvTTpaewDKPSGJp8g-e_FEUMwjDlEO3rxXHwu6513G5MsCVUPLPAiEKVHIa7qAb-uqRFl6Kk/s1600/izana_google_maps.jpg" width="550" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><a href="https://www.google.es/maps/@28.2896725,-16.5161694,5525m/data=!3m1!1e3?hl=en">Credit: Google Maps (click to go to the interactive map)</a></td></tr>
</tbody></table>
<br />Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com0Santa Cruz de Tenerife, Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989tag:blogger.com,1999:blog-1169051187754614539.post-22892196635757592232015-02-28T12:25:00.000-08:002016-09-16T08:15:14.122-07:001D Solar Atmosphere Models in IDL: Gingerich et al (HSRA)The Harvard-Smithsonian Reference Atmosphere (HSRA) by Gingerich et al (1971, <a href="http://adsabs.harvard.edu/abs/1971SoPh...18..347G">1971SoPh...18..347G</a>) is another widely used semi-empirical model of the idealized plane-parallel solar atmosphere in hydrostatic equilibrium. This model also includes the chromosphere where the hydrogen ionization is solved using the statistical equilibrium equations. The model is still often used as the reference solar atmosphere in 1D and its number of citations steadily grows (> 820 so far). The paper is very clearly written and it's a very recommendable reading especially if you use this model. Beside the model itself, there is several tables with computed intensities, brightness temperature, optical depths at different wavelengths, etc.<br />
<br />
<a name='more'></a><br />
The depth-dependent variables are listed in two tables. The first one (Table IV, p.359) contains the variables:<br />
<ul>
<li>optical depth,</li>
<li>geometrical depth, </li>
<li>temperature,</li>
<li>pressure, </li>
<li>electron pressure,</li>
<li>opacity,</li>
<li>density, </li>
<li>hydrogen ionization ratio,</li>
<li>distribution of the main electron donors (Mg, Si, Fe, C),</li>
<li>radiative and adiabatic gradients.</li>
</ul>
All these variables are also in my files below. The other table (Table V, p.360) contains the depth distribution of the following variables:<br />
<ul>
<li>Rayleigh scattering contribution in the total absorption,</li>
<li>electron scattering in the total absorption,</li>
<li>opacity per negative hydrogen ion and per neutral hydrogen,</li>
<li>number density of CO,</li>
<li>departure coefficient for hydrogen ground state, $b_1$,</li>
<li>ionization fraction of the metals that are the main electron donors. </li>
</ul>
The data from this table I did not digitized.<br />
<br />
The
data in my files is OCR-ed and manually double checked. There are two
versions of the model: .dat (ascii) and .sav (IDL). In the IDL files the
model is stored in the structure variables hsra that has the following tags:<br />
<pre style="background-color: #eeeeee; border: 1px dashed rgb(153, 153, 153); color: black; font-family: Andale Mono,fixed,monospace; font-size: 12px; line-height: 14px; overflow: auto; padding: 5px; width: 100%;"><code>hsra = {h:h, tau5:tau5, t:t, pgas:pgas, pel:pel, kappa:kappa, rho:rho, hion:hion,
el_metals:el_metals, el_si:el_si, el_mg:el_mg, el_fe:el_fe, el_c:el_c,
rgrad:rgrad, agrad:agrad, variables:variables, reference:reference}</code></pre>
where variables has two self-explanatory tags, varnames and varunits.<br />
<br />
<br />
The files are all in a tarball:<br />
<br />
<a href="https://db.tt/5adimNso">hsra.tar</a><br />
<br />
<br />
<b>Disclaimer</b>: This data is publicly available in Gingerich et al (1971, <a href="http://adsabs.harvard.edu/abs/1971SoPh...18..347G">1971SoPh...18..347G</a>).
If you use the data in a publication please cite it properly. I
double-checked that the digitized version here matches the original,
still typos are possible no matter how unlikely. If you have any doubts,
compare the numbers to those in the paper. Nikola Vitashttp://www.blogger.com/profile/11652414242984631263noreply@blogger.com2Santa Cruz de Tenerife, Santa Cruz de Tenerife, Spain28.4636296 -16.25184669999998728.4077896 -16.332527699999986 28.5194696 -16.171165699999989