Tuesday 30 October 2012

Carlson's Quadrature Set A

When the radiative transfer equation is solved in three dimensions, it is necessary to discretize the angular dependence of the specific intensity. In principle, the continuous integral over solid angle is replaced by a discrete sum over selected rays. This is the often used technique known as the discrete ordinates method (DOM). However, there is no unique solution for the discretization of the solid angle.

The state-of-art radiative magnetohydrodynamical codes (eg. MURaM and BIFROST) use a quadrature known as Set A by Carlson (1963). The advantage of this particular quadrature is that it is invariant with respect to 90o axis rotation.


The quadrature is defined by directional cosines $\mu_l$ and their weights $w_l$. They are combined into sets that define rays $\Omega_m = (\mu_x, \mu_y, \mu_z)$ and their weights $w_m$. The order of quadrature $n$ specifies the number of $\mu$ levels in the [-1, 1] interval. Since the quadrature is symmetric in respect to axis, it is sufficient to find values of the angles in half of that range. Thus the number of $\mu$ levels is $l_{max} = n/2$. The number of rays per octant in this scheme is $n(n+2)/8$. Each ray is labeled by three numbers indexing $\mu$ level along one of the axes. The sum of indices of a ray is the same of all rays, $l_{max}+2$. For example, when $n=8$ and $l_{max}=4$, possible indices are 1, 2, 3 and 4. Their possible variations are those that have sum of 6. The number of these variations (and the rays) is $n(n+2)/8 = 10$. The number of different classes (variation without repetition with the specified sum) is $l_{max}-1 = 3$. In our example: (4, 1, 1), (3, 2, 1) and (2, 2, 2).

Wednesday 24 October 2012

Talk at IAC solar seminar

At the IAC solar seminar I gave a talk on fast horizontal flows in solar granulation and their spectroscopic signatures. The pdf file with the presentation is here. Most of that work was done in collaboration with my Utrecht colleagues Catherine Fischer, Alexander Vogler and Christoph Keller and it has already been published in A&A (Vitas et al, 2011). Thanks those who attended the talk for a nice and lively discussion!

Figure shows profiles of Fe I 6301.5 A computed at the location of a large supersonic event that produced a shock front. They are computed at $\mu = 0.4$ (a position away from the center of the solar disk) and for 8 azimuthal directions. It's like we've been moving around the shock. The profiles are computed at high resolution, but degraded in both spatial and spectral domain to mimic the capabilities of the SP/SOT instrument at Hinode spacecraft. It is interesting to note how the asymmetry of the line changes with azimuth. When we see the flow is moving toward us (0 degrees), the blue wing of the profile is extended and distorted. When it is moving away from us (180 degrees), the red wing is affected. That behavior is due to the well-known Doppler effect. The figure tells us that the high velocity is not sufficient for the supersonic flows to be identified through line deformation. The velocity had to be roughly aligned with the line of sight of the observer.

Saturday 20 October 2012

LaTeX in Blogger

Recently I found about MathJax, a tool that enables use of LaTeX in blogs and other web platforms. It's an excellent tool! At the Irreducible representations blog there is a nice and simple explanation how to add this feature to your own blog. I copy the explanation here:

To get MathJax to work in Blogger, just go to your Blogger account, click "Design" (top right of the page), and then "Edit HTML". After the first <head> you see, paste
<script src='http://cdn.mathjax.org/mathjax/latest/MathJax.js' type='text/javascript'>    
        HTML: ["input/TeX","output/HTML-CSS"],
        TeX: { extensions: ["AMSmath.js","AMSsymbols.js"], 
               equationNumbers: { autoNumber: "AMS" } },
        extensions: ["tex2jax.js"],
        jax: ["input/TeX","output/HTML-CSS"],
        tex2jax: { inlineMath: [ ['$','$'], ["\\(","\\)"] ],
                   displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
                   processEscapes: true },
        "HTML-CSS": { availableFonts: ["TeX"],
                      linebreaks: { automatic: true } }

And that's all! Now you can use standard LaTeX commands to format mathematical text. For example, the radiative transfer equation along a ray (e.g. see the lecture notes of Rob Rutten, Eq.2.24) can be written as:

\frac{\mathrm{d}I_{\nu}}{\alpha_{\nu}\mathrm{d} s} = S_{\nu} - I_{\nu}

and appears as: $$ \frac{\mathrm{d}I_{\nu}}{\alpha_{\nu}\mathrm{d} s} = S_{\nu} - I_{\nu} $$

Thursday 11 October 2012

Cubic splines: Interpolation and smoothing

Cubic splines are widely used technique for piecewise smooth interpolation. Between any two points of a discretely given function, the function is approximated by a third order polynomial. However, splines can also be used for smoothing or fitting of measurements. The theory behind the splines is nicely explained by Pollock (1999). I implemented that algorithm in IDL procedure splinecoeff.pro.

To try it yourself, here is a small example. Download the text file with the values of the sunspot area data, daily_area.txt from solarscience.msfc.nasa.gov, and run the following:

; Load the data
n = FILE_LINES('daily_area.txt')-1
OPENR, 1, 'daily_area.txt'
a = ''
READF, 1, a
area = FLTARR(6, n)
READF, 1, area

; Extract data between 1977-01-01 and 2011-12-31
initial = WHERE(area[0, *] EQ 1977 AND area[1, *] EQ  1 AND area[2, *] EQ  1)
final   = WHERE(area[0, *] EQ 2011 AND area[1, *] EQ 12 AND area[2, *] EQ 31)
time = JULDAY(area[1,*], area[2,*], area[0,*])
x = time[initial:final]
y = REFORM(area[3, initial:final])
nx = N_ELEMENTS(x)

; Compute coefficients of the cubic splines
coeff = SPLINECOEFF(x, y, lambda = 1.d3)

; Plot the original data and the spline function for given lambda
y1 = FLTARR(nx-1)
x1 = x[0:nx-2]
FOR i = 0, N_ELEMENTS(y)-2 DO y1[i] = coeff.d[i] + $
                                      coeff.c[i] * (x[i+1]-x[i]) + $
                                      coeff.b[i] * (x[i+1]-x[i])^2 + $
                                      coeff.a[i] * (x[i+1]-x[i])^3
PLOT, x, y, psym = 3
OPLOT, x1, y1
The only difference between the interpolation, smoothing and fitting of the measurement is in the free parameter $\lambda$. If $\lambda$ is very small, than the result is the cubic spline interpolation. If $\lambda$ is very large, the measurement is smoothed by cubic splines. In the extreme case, the resulting splines become a linear function. The following animation show how the result change with $\lambda$:

Note: See paper of Djurovic and Paquet (1999) for an example where the cubic splines are used for filtering and smoothing of data.

Wednesday 10 October 2012

Vacuum Tower Telescope (VTT), Teide Observatory (Tenerife, Spain)

This Wednesday I had a chance to assist an observation at the famous VTT (with E.Khomenko and M.Collados). The telescope is one of the most productive instruments in solar physics. The diameter of the primary mirror is 700 mm. The left picture shows the heliostat - it follows the Sun during the day and directs the light into the telescope. The right one shows the tower from the ground. It's an 11-storeys building!

It was a lot of fun to be there. Can't wait the next opportunity!