Carlson's Quadrature Set B

Beside the quadrature Set A that is described and computed in the previous blog, Carlson (1963) described an alternative set of quadratures, Set B. This set corresponds to odd approximations. It "relates to Set A in the way Double-P relates to Single-P." (ibid). The principles on which this set is constructed are same as in the case of Set A. However, the equations slightly differ. First we determine the $W_l$ coefficients from the system:
$$ \sum_{l=1}^{n/2-1} W_l = \frac{n-3}{3}, $$ $$ W_l^2 = W_1^2 + (l-1)\Delta^\prime, $$ where $\Delta^\prime = 2/n$. This system is again solved by Newton-Raphson formula. The $W_l$ values are then used to get $w_l$ (as with Set A).

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 90^o axis rotation.

Terminology

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).

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.

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'>    
    MathJax.Hub.Config({
        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 } }
    });
</script>

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} $$

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
CLOSE, 1

; 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.

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!

Solar telescope at the Terskol observatory (Caucasus, Russia)

These pictures were taken in March 2006 during a solar-eclipse campaign. The top one shows the heliostat of the large horizontal solar telescope of the Peak Terskol Observatory (Kabardino-Balkaria, Caucasus, Russia) at 3100 m of altitude. The entire building in the middle picture (behind the heliostat) is the actual telescope. The bottom photography shows the interior of the telescope with the primary mirror (D=650 mm). More about this amazing instrument you can find on the website of the observatory and in Burlov-Vasilev et al (1996). However, behind this telescope is a great human power of my dear friend Olexa Andriyenko. Working with him was one of the best experiences of my life.

Dutch Open Telescope (DOT), La Palma (Spain)

The Dutch Open Telescope is very easy to recognize. With its look of an extraterrestrial space ship, it is so different than all the other solar telescopes (well, not only solar) in the world. The minimalistic construction is very important concept as it enables an efficient natural cooling. The diameter of the primary mirror is 450 mm and the focal length only 2 m . These pictures were taken in October, 2007 during an USP-SP workshop at the DOT. The dome of the telescope is retractable (closed in the bottom picture). The DOT website.

New coordinates, new challenges

This spring I changed my job. After a nice and inspiring year that I spent doing the Earth's observation at SRON Netherlands Institute for Space Research, I am back to the solar physics. My new affiliation is Instituto de Astrofísica de Canarias (IAC) in La Laguna (Tenerife, Spain). Here I will be working on an exciting and challenging project: "Solar Partially Ionized Atmosphere (ERC-2011-StG 277829-SPIA)", financed by the European Research Council in the frame of the FP7 Specific Program IDEAS. The PI of the project is Dr. Elena Khomenko, while the other team members are Dr. Angel de Vicente, Dr. Manuel Luna Bennasar, Dr. Antonio Díaz Medina and Prof. Dr. Manuel Collados Vera. More details about the project are available here.

Address: Instituto de Astrofísica de Canarias (IAC) C/ Vía Láctea, s/n E38205 - La Laguna (Tenerife), España Phone: (+34) 922 60 5746 Email: n.vitas@iac.es

The final days of the Utrecht astronomy

Many streets in Utrecht are named after celestial bodies or astronomical terms - Meteoorhof, Jupiterstraat, Eclipshof, Komeetstraat, Venuslaan, etc. There is also Zonstraat - a cozy street dedicated to the Sun that goes all the way through one of the most beautiful neighborhoods to the old observatory Sonnenborgh. Some people in the street keep "solar portraits" on their facades and in the windows. It's a real shame that only memories are left of astronomy in this town.

 

More on the fallout of Utrecht astronomy: 

C. U. Keller: Sterrekundig Instituut Utrecht: The Last Years

R. J. Rutter: Sterrekundig Instituut Utrecht

Nature 

Good Night and Good Luck, Utrecht Astronomers

Austerity bites in Utrecht