Polarimetry Processing#

astropop polarimetry computation currently is desgined to work only with dual-beam polarimeters with half-wave plate (HWP) and quarter-wave plate (QWP) retarders, performing both circular and linear polarization computations, according the instrument capabilities and design. Savart or Wollaston analyzers are supported. For example of instruments that can be used with astropop, we already tested with success the following instruments:

IAGPOL [1] - Gaveta Polarimétrica do IAG link at Observatório Pico dos Dias, Brazil
SPARC4 [2] - Simultaneous Polarimeter and Rapid Camera in Four bands link at Observatório Pico dos Dias, Brazil
LE2POL [3] - Leicester dual-beam imaging polarimeter at University of Leicester, UK
FORS2 - FORS2, FOcal Reducer/low dispersion Spectrograph 2 at ESO, Chile, in IPOL mode. link
POLIMA2 - Optical Dual-Beam Polarizer at OAN-SPM, Mexico. link

In all this polarimeters, the light beam passed by a retarder and is decomposed in two, ortogonal polarized, beams. So, the ordinary and extraordinary beams produce separated images and all the stars can appear duplicated in the detector. For example, see the image bellow taken from [3]. astropop is designed to handle this kind of data.

LE2POL image, where all stars are duplicated. So, each star produces one image for each ordinary and extraordinary beams. The image is from [3].#

Description of the Polarimetry Reduction#

The reduction process overview is derived from the work of Magalhães et al. [1] and Rodrigues et al. [4] , using a fitting process better described by Campagnolo (2019) [5]. In a few words, the code computes the normalized flux difference between the ordinary and extraordinary beams ($F_{e}$ and $F_{o}$) for each retarder plate position:

\[Z_{i} = \frac{F_{e} - F_{o}}{F_{e} + F_{o}}\]

And fits a equation for a given set o positions. There are two equations that can be used:

For half-wave plates, only linear polarization is computed, and the equation is:

\[Z_i = Q \cos{4 \psi_i} + U \sin{4 \psi_i}\]
For quarter-wave plates, linear and circular polarization is computed, and the equation is:

\[Z_i = Q \cos^2{2 \psi_i} + U \sin{2 \psi_i} \cos{2 \psi_i} - V \sin{2 \psi_i}\]

Where $Q$, $U$ and $V$ are the linear and circular polarization fractions and $psi_i$ is the retarder plate position.

Normalization#

If the instrument have different sensibility for each beam or if the PSF of the two beams do not match (a common thing in Savart analyzers), a normalization $k$ factor must be included in the equation of the relative flux ($Z_i$). This normalization factor scales one of the beams to compensate the difference in sensibility. So, the flux difference equation becomes:

\[Z_{i} = \frac{F_{e} - F_{o} \cdot k}{F_{e} + F_{o} \cdot k}\]

This $k$ factor is computed in different ways for half-wave and quarter-wave plates.

For half-wave plates, $k$ is computed as the ratio between the sum of the extraordinary and ordinary fluxes along a complete set of 4 consecutive retarder positions (0, 22.5, 45 and 67.5 degrees). The explanation is better described in [1]. This is possible because the equation for half-wave plates is symmetrical in 0. So:

\[k = \frac{\sum_{i=1}^{4} F_{e, i}}{\sum_{i=1}^{4} F_{o, i}}\]
For quarter-wave plates, the term $Q\cos^2{2 \psi_i}$ in the equation is always positive, so the results are not symmetrical in 0. Due to this, the equation above cannot be used to compute $k$. Rodrigues (1998) [4] and Lima (2021) [6] describe a iterative process, where $k$ depends on $Q$. According this process, the Stokes parameters are computed with some $k$, and then a new $k$ is computed using the new Stokes parameters. This process is repeated until the computed Stokes parameters converge within a given tolerance.

Both methods are implemented in astropop.

Angle Zero of the Retarder#

It is expected that the angle zero of the retarder plate do not match perfectly with the analyzer and the equatorial system of the sky. Due to this, the observation of standard stars may be needed to compute the correction to be applied to fix the instrumental computed position angle and the real position angle of the sky. As this is a standard, post-reduction, data processing step, we will not go deep into this topic here. An offset zero do not affect the results for linear polarization in half-wave plates.

However, for circular polarization, the zero offset is important due to the assymmetry of the equation. So, with this configuration, there is an additional offset to be computed within the Stokes computation process. This instrumental offset, if wrong, leads to wrong results and with high errors. Additionally, the best wy to compute it is by reducing data from a star with strong circular polarization level. Once this offset is determined, it can be used in other reductions, since it is commonly a instrument-dependent value.

Polarimetry Reduction in `astropop`#

There are 2 main steps in the polarimetry reduction process: the pairs matching and the Stokes parameters computation.

To compute the Stokes parameters of a given source, astropop uses arrays of fluxes and retarder positions. So, it is not directly dependent on the photometry routines to process the data and you can use your own already extracted photometry for this process. Also, our helper functions to match the pairs of ordinary and extraordinary beams also use just arrays of x and y positions, and are not dependent on the photometry routines too. So, no part of the polarimetry reduction depend on the use of |astropop| for source extraction and photometry and you can use your own data or code in this process.

For example, in this documentation we will use totally fake photometry and random pairs positions to illustrate the polarimetry reduction process.

Stokes-Least-Squares Polarimetry Processing#

The Stokes-Least-Squares (SLS) method is a method to compute the Stokes parameters by fitting the modulation equations described above to the relative flux difference of the data using a least-squares method with Trust Region Reflective algorithm. For a better description and validation of this method, see [5].

First of all, lets create the fake data to be used here for example.

In [1]: import numpy as np

In [2]: from astropop.polarimetry import halfwave_model, quarterwave_model

In [3]: flux = 1e7  # total flux of the source, in electrons

# retarder positions, in degrees
In [4]: psi_i = np.arange(0, 360, 22.5)

# half-wave plate with only linear polarization
# q=0.05, u=-0.07
In [5]: z_i_half = halfwave_model(psi_i, 0.05, -0.07)

In [6]: z_i_half = np.random.normal(loc=z_i_half, scale=0.005)  # add some noise

In [7]: fo_half = flux*(1+z_i_half)/2  # ordinary beam fluxes

In [8]: fe_half = flux*(1-z_i_half)/2  # extraordinary beam fluxes

# quarter-wave plate with linear and circular polarization
# q=0.05, u=-0.07, v=0.08
In [9]: z_i_quarter = quarterwave_model(psi_i, 0.05, -0.07, 0.08)

In [10]: z_i_quarter = np.random.normal(loc=z_i_quarter, scale=0.005)  # add some noise

In [11]: fo_quarter = flux*(1+z_i_quarter)/2  # ordinary beam fluxes

In [12]: fe_quarter = flux*(1-z_i_quarter)/2  # extraordinary beam fluxes

../_images/polarimetry-3.png — Simulated data for the Stokes-Least-Squares polarimetry processing.#

Create an `SLSDualBeamPolarimetry`#

The Stokes parameters computation is entirely done by the SLSDualBeamPolarimetry class. This class stores some configuration parameters and process a given set of data using the compute method. During the instance creation you must choose the type of the retarder used in the observations, between halfwave and quarterwave options. Aditional arguments can be passed, like normaliation computation, zero position of retarder, and others. Check the class documentation for more information.

First, you must create a SLSDualBeamPolarimetry instance with your options. Here are some examples:

In [13]: from astropop.polarimetry import SLSDualBeamPolarimetry

# half-wave plate
In [14]: s = SLSDualBeamPolarimetry('halfwave')

# quarter-wave plate
In [15]: s = SLSDualBeamPolarimetry('quarterwave')

# half-wave plate and tell the class to compute the k parameter
In [16]: s = SLSDualBeamPolarimetry('halfwave', compute_k=True)

# quarter-wave plate with normalization and zero position at 45 degrees
In [17]: s = SLSDualBeamPolarimetry('quarterwave', compute_k=True, zero=45)

For normalization, there is two (incompatible) arguments: compute_k and k. When using k=value, you are setting a fixed value for the normalization constant, so, it will not be computed by the code. For example, if you set k=1.05, it will use this factor in the nromalization. for compute_k=True, the normalization constant will be computed by the code. In general, all instruments we tested have this $k$ value very close to 1.0, but using normalization may improve some cases.
For the zero computation, you can choose between a fixed value of zero, using zero=value in degrees. So, setting zero=30.2 will force the code to use 30.2° as the zero point of the retarder. When the zero is not specified, there are two behaviors you can expect: - for halfwave retarders, the code will use 0.0 as the zero position, sice it is not possible to compute it from the data. - for `quarterwave` retarders, the code will try to fit the zero together with the Stokes parameters. The code will seek the value the zero inside the limits defined by zero_range argument. So, if you want the code to find a zero between 0 and 45 degrees, you can set zero_range=(0, 45). By default, the code will use zero_range=(0, 90).

Note

It is important to mention that quarterwave retarders find the zero position better when circular polarization is present. So, it is recommended that you use a circular polarized star to compute the zero of the instrument and use it for the other stars in the observation.
For stars that only have linear polarization, the code can compute the zero only if compute_k is set to False. This happens because the zero and k parameters are degenerated in quarterwave retarders due to the $Q\cos^2{2 \psi_i}$. So, using compute_k=True and zero=None will not converge to a solution in this case. For stars with circuar polarimetry, these arguments can be used together and the code will be able to find the zero position and the normalization constant at the same time.

The iteration process is controlled by 2 arguments: iter_tolerance and max_iters. iter_tolerance sets the limit between runs where the code assumes the solution is converged. Like, if the difference between the Stokes parameters of the last two runs is less than iter_tolerance, the code will stop. The max_iters sets the maximum number of iterations the code will run. By default, iter_tolerance=1e-5 and max_iters=100, which give very good results and are enough for most cases.

Stokes Parameters Computation#

The compute method is the method to be used to compute the Stokes parameters from the flux data. It recieves the list of retarder positions (in degrees) and the lists of ordinary and extraordinary fluxes for a source in each retarder position. Optionally, you can pass the lists of errors in the fluxes, that will be used to weight the least squares fit.

For example, computing the Stokes parameters for the half-wave plate data without errors:

In [18]: s = SLSDualBeamPolarimetry('halfwave')

In [19]: result = s.compute(psi_i, fo_half, fe_half)

In [20]: print(result)
q=0.052+-0.002, u=-0.069+-0.002

And with errors:

In [21]: s = SLSDualBeamPolarimetry('quarterwave', zero=0)

# let's assume the errors are the square root of the fluxes
In [22]: fo_quarter_err = np.sqrt(fo_quarter)

In [23]: fe_quarter_err = np.sqrt(fe_quarter)

In [24]: result = s.compute(psi_i, fo_quarter, fe_quarter,
   ....:                    fo_quarter_err, fe_quarter_err)
   ....: 

In [25]: print(result)
q=0.051+-0.002, u=-0.070+-0.003, v=0.077+-0.002

`StokesParameters` class#

The result of compute is a StokesParameters instance. This class is just a container that stores relevant parameters of the computation. Both final Stokes parameters ($Q$, U and V) and intermediate data, like $\psi_i$ retarder positions, $Z_i$ relative diference of fluxes, zero position and total flux in each position. It also compute some important properties, like position angle of linear polarization $\theta$, linear polarization degree $P$ and theoretical error in Stokes parameters. Most of the properties are stored as QFloat containers.

The stored parameters are:

In [26]: print(f'q={result.q}')
q=0.051+-0.002

In [27]: print(f'u={result.u}')
u=-0.070+-0.003

In [28]: print(f'v={result.v}')
v=0.077+-0.002

In [29]: print(f'k={result.k}')
k=1.0

In [30]: print(f'psi={result.psi}')
psi=[0.0+-0.0, 22.5+-0.0, 45.0+-0.0, 67.5+-0.0, 90.0+-0.0, 112.5+-0.0,
 135.0+-0.0, 157.5+-0.0, 180.0+-0.0, 202.5+-0.0, 225.0+-0.0, 247.5+-0.0,
 270.0+-0.0, 292.5+-0.0, 315.0+-0.0, 337.5+-0.0] unit=deg

In [31]: print(f'z_i={result.zi}')
z_i=[0.0560+-0.0003, -0.0624+-0.0003, -0.0883+-0.0003, 0.0066+-0.0003,
 0.0565+-0.0003, 0.0392+-0.0003, 0.0806+-0.0003, 0.1152+-0.0003,
 0.0477+-0.0003, -0.0609+-0.0003, -0.0752+-0.0003, 0.0077+-0.0003,
 0.0439+-0.0003, 0.0461+-0.0003, 0.0742+-0.0003, 0.1141+-0.0003]

In [32]: print(f'zero={result.zero}')
zero=0.0+-0.0 deg

In [33]: print(f'flux={result.flux}')
flux=[10000000+-3000, 10000000+-3000, 10000000+-3000, 10000000+-3000,
 10000000+-3000, 10000000+-3000, 10000000+-3000, 10000000+-3000,
 10000000+-3000, 10000000+-3000, 10000000+-3000, 10000000+-3000,
 10000000+-3000, 10000000+-3000, 10000000+-3000, 10000000+-3000]

The linear polarization level is obtained from $P=\sqrt{Q^2+U^2}$ and is accessed by p property:

In [34]: result.p
Out[34]: 
<QFloat at 0x7f62656be0d0>
0.087+-0.003

The position angle of linear polarization is obtained from $\theta=\frac{1}{2}\tan^{-1}\left(\frac{U}{Q}\right)$ and is accessed by theta property:

In [35]: result.theta
Out[35]: 
<QFloat at 0x7f626688d710>
153+-1 deg

The RMS between the model using computed parameters and the $Z_i$ values stored can be accessed directly by rms property. It is only available if zi is stored in the class.

In [36]: result.rms
Out[36]: 0.004408583300385898

A theoretical expected error based on the fluxes SNR can be also computed by StokesParameters using theor_sigma method. It is only available if flux is stored in the class.

In [37]: result.theor_sigma
Out[37]: 
{'q': 0.00012562972690740153,
 'u': 0.00020661862653352698,
 'v': 0.00011693227466553227,
 'p': 0.00018273659819396636}

The computation of this theoretical error is not trivial. For halfwave retarders it is:

\[\sigma_P = \sigma_Q = \sigma_U = \frac{\sqrt{2}}{SNR}\]

and for quarterwave retarders it is:

\[\sigma_Q = \frac{1}{\sqrt{0.396} \cdot SNR}\]

\[\sigma_U = \frac{1}{\sqrt{0.1464} \cdot SNR}\]

\[\sigma_V = \frac{1}{\sqrt{0.4571} \cdot SNR}\]

\[\sigma_P = \frac{\sqrt{Q^2 \cdot \sigma_Q^2 + U^2\sigma_U^2}}{P}\]

Also, a callable model with fixed Stokes parameters can be generated by model method. Useful for plotting or data/fitting comparison. The only argument of the callable model is the psi retarder position.

In [38]: plt.figure(figsize=(6, 4));

In [39]: plt.errorbar(result.psi.nominal, result.zi.nominal,
   ....:              yerr=result.zi.std_dev, fmt='o', label='data');
   ....: 

In [40]: line = np.linspace(0, 360, 300);

In [41]: model = result.model(line*result.psi.unit).nominal

In [42]: plt.plot(line, model, '--', label='model');

In [43]: plt.xlabel(r'$\psi$');

In [44]: plt.ylabel(r'$Z_i$');

In [45]: plt.title('Example of data and the callabel model comparison');

In [46]: plt.legend();

In [47]: plt.tight_layout();

In [48]: plt.show()

Helper Tools#

As seen above, the dual-beam polarimeters produce images where all sources appear duplicated, due to ordinary and extraordinary beams presence. So, how can we separate the two beams? How to match the pairs to correctly process the polarimetry?

The polarimetry module has some helper tools to help with this task. The estimate_dxdy that computes the displacement, in pixels, between the ordinary and extraordinary beams and the match_pairs function that match the pairs of sources in the image, based on a given displacement.

# random positions for the pairs
# 80 pairs of sources in a image of 1024x1024 pixels
# dx=25, dy=35
In [49]: x = np.random.uniform(0, 1024, 80)

In [50]: x = np.concatenate((x, x+25))  # x positions of the rouces

In [51]: y = np.random.uniform(0, 1024, 80)

In [52]: y = np.concatenate((y, y+35))  # y positions of the rouces

# plot the sources
In [53]: plt.figure(figsize=(6, 6));

In [54]: plt.plot(x, y, 'ko', ms=5);

In [55]: plt.title('Example of the pairs of sources');

Displacement Computation#

To find the displacement in $x$ and $y$ directions ($dx$ and $dy$), the estimate_dxdy uses a iterative clipped histogram around the mode of the distribution.

The code compute the distance of eache combination of two sources in the image, considering only the stars where $y_i > y_j$, where $y_i$ and $y_j$ are the $y$ positions of the first and the second sources in the pair, respectively. This filtering eliminates reduncaties in the histogram of positive and negative values. Also, only pairs of stars closer than a dist_limit in each axes are considered.

After, the histogram is computed and the mode is estimated using the maximum value in this histogram. Than the histogram is clipped around this mode, and a new mode is estimated from a new, rebined, histogram. This iterative process repeats for all steps set for the function. For example, setting steps=[50, 5, 3, 1], the histogram will clipped in 50 pixels around the mode in the first iteration, 5 pixels in the second, 3 pixels in the third and 1 pixel in the last iteration. The number of bins of the histogram also interferes in the final result.

For example, let’s estimate the displacement between the ordinary and extraordinary beams of our fake sources.

In [56]: from astropop.polarimetry import estimate_dxdy

In [57]: dx, dy = estimate_dxdy(x, y, dist_limit=50, steps=[50, 5, 3, 1], bins=100)

In [58]: print(f'dx = {dx:.2f} pixels')
dx = -25.00 pixels

In [59]: print(f'dy = {dy:.2f} pixels')
dy = -35.01 pixels

The result of $dx \approx 25$ pixels and $dy \approx 35$ pixels is the expected for the fake sources, since they where used to generate the coordinates. The final precision, in general, depends on the error of the positions of the sources, the number of beams and the size of the clipping in the final step.

Pairs Matching#

Once you have a good value for the displacement, that is instrument-deppendent and may not need to be computed everytime, you can use the match_pairs function to match the pairs of ordinary and extraordinary sources in the image.

This function uses cKDTree to find the nearest neighbors of the ordinary sources in the extraordinary sources. The tolerance defines the maximum distance limit, in pixels, to consider a matching.

The result of the code is a structured ndarray with the indexes of ordinary and extraordinary sorces in the original $x$ and $y$ arrays. The ordinary sources will appear in the o column and the extraordinary sources in the e column.

Note

The use of ordinary and extraordinary names for the beams is just for a convention. This do not represent what beam, in reality, is the ordinary or the extraordinary.

In [60]: from astropop.polarimetry import match_pairs

# let's sort the sources by y position just to mix them
In [61]: sorting = np.argsort(y)

In [62]: x = x[sorting]

In [63]: y = y[sorting]

In [64]: pairs = match_pairs(x, y, 25, 35, tolerance=2.0)

In [65]: print(pairs)
[(  3,   0) (  4,   0) (  6,   1) (  7,   2) (  8,   3) (  9,   5)
 ( 13,  10) ( 15,  11) ( 16,  12) ( 18,  14) ( 20,  17) ( 24,  19)
 ( 28,  21) ( 29,  22) ( 31,  23) ( 33,  25) ( 34,  26) ( 35,  27)
 ( 36,  30) ( 37,  32) ( 40,  38) ( 41,  39) ( 45,  42) ( 46,  43)
 ( 48,  44) ( 50,  47) ( 53,  49) ( 55,  51) ( 57,  52) ( 59,  54)
 ( 60,  56) ( 61,  58) ( 66,  62) ( 67,  63) ( 68,  64) ( 71,  65)
 ( 75,  69) ( 77,  70) ( 79,  72) ( 82,  73) ( 84,  74) ( 86,  76)
 ( 90,  78) ( 91,  80) ( 92,  81) ( 93,  83) ( 94,  85) ( 96,  87)
 ( 97,  88) ( 98,  89) (103,  95) (106,  99) (107, 100) (108, 101)
 (111, 102) (113, 104) (115, 105) (116, 109) (117, 110) (120, 112)
 (121, 114) (124, 118) (125, 119) (129, 122) (130, 123) (133, 126)
 (134, 127) (135, 128) (137, 131) (139, 132) (142, 136) (144, 138)
 (145, 140) (147, 141) (149, 143) (151, 146) (153, 148) (155, 150)
 (157, 152) (158, 154) (159, 156)]

So, to access the coordinates of the nth pair, you can use the pairs['o'][n] and pairs['e'][n] to get the index.

In [66]: n = 27

# indexes of the stars
In [67]: index_o = pairs['o'][n]

In [68]: index_e = pairs['e'][n]

# coordinates of the ordinary beam
In [69]: x_o = x[index_o]

In [70]: y_o = y[index_o]

# coordinates of the extraordinary beam
In [71]: x_e = x[index_e]

In [72]: y_e = y[index_e]

In [73]: print(f"star: {n}\n"
   ....:       f"ordinary: index={index_o} x={x_o:.2f} y={y_o:.2f}\n"
   ....:       f"extraordinary: index={index_e} x={x_e:.2f} y={y_e:.2f}")
   ....: 
star: 27
ordinary: index=55 x=27.46 y=465.15
extraordinary: index=51 x=2.46 y=430.15

In [74]: plt.figure(figsize=(5, 5));

In [75]: plt.plot(x, y, 'ko', ms=5, alpha=0.5);

In [76]: plt.plot([x_o, x_e], [y_o, y_e], 'r-', lw=2);

In [77]: plt.title(f"Pair n={n} matched");

# plot all matched pairs
In [78]: plt.figure(figsize=(5, 5));

In [79]: plt.plot(x, y, 'ko', ms=5, alpha=0.5);

In [80]: for pair in pairs:
   ....:     index_o = pair['o']
   ....:     index_e = pair['e']
   ....:     x_o = x[index_o]
   ....:     y_o = y[index_o]
   ....:     x_e = x[index_e]
   ....:     y_e = y[index_e]
   ....:     plt.plot([x_o, x_e], [y_o, y_e], 'r-', lw=2)
   ....: 

In [81]: plt.title(f"All pairs matched");

References#

Polarimetry API#

astropop.polarimetry Package#

Compute polarimetry of dual beam polarimeters images.

Functions#

`estimate_dxdy`(x, y[, steps, bins, dist_limit])	Estimate the displacement between the two beams.
`match_pairs`(x, y, dx, dy[, tolerance])	Match the pairs of ordinary/extraordinary points (x, y).
`quarterwave_model`(psi, q, u, v[, zero])	Compute polarimetry z(psi) model for quarter wavelength retarder.
`halfwave_model`(psi, q, u[, zero])	Compute polarimetry z(psi) model for half wavelength retarder.

Classes#

`SLSDualBeamPolarimetry`(retarder, k, zero, ...)	Polarimetry computation for Stokes Least Squares algorithm.
`StokesParameters`(retarder, q, u[, v, k, ...])	Store the Stokes parameters results from dual beam polarimeters.