Coverage for src/susi/model/susi_mm.py: 14%
98 statements
« prev ^ index » next coverage.py v7.5.0, created at 2025-06-13 14:15 +0000
1# -*- coding: utf-8 -*-
2import sys
3import logging
4from os import path
5import numpy as np
7__author__ = "Francisco Iglesias"
8__copyright__ = "Copyright 2020, Max Planck Institute for Solar System Research"
9__version__ = "0.01"
10__maintainer__ = "Francisco Iglesias"
11__email__ = "franciscoaiglesias@gmail.com"
13# Definition of common Mueller matrices
16def susi_mmret(pang, ret_in, Cd=0, wl0=0, wl=0, tran=0):
17 """
18 Computes the Mueller matrix of a linear retarder
19 acording to eq 4.49 in "del Toro Iniesta 2003"
21 INPUTS
22 pang: Scalar angle of the optical axis [deg]
23 ret : Array [nx,ny] of retardances at wl0 for each pixel [deg]
24 wl0 : Scalar central wavelength [nm]
25 Cd : Scalar coeficient for the retardence wavelength dependece aproximation.
26 See Gisler et. al. 2003. Use 0 to have no dispersion, in this
27 case wl and wl0 won't be used (default)
28 wl : Scalar current wavelength [nm]
30 OPTIONAL INPUTS:
31 tran: Set to retrive the transpose of the mueller matrix
33 OUTPUTS
34 [nx,ny,4,4] Mueller matrix of the linear retarder defined by Inputs
35 """
36 dim = np.shape(ret_in)
37 ret = ret_in / 360.0 # in frac of wl
38 phi = pang * np.pi / 180.0 # [Rad]
40 c2 = np.cos(2.0 * phi) # just some intermediate constants
41 s2 = np.sin(2.0 * phi)
43 if Cd != 0: # retardance [Rad] at the current wavelength
44 ret = 2.0 * np.pi * (ret * wl0 / wl - Cd / (wl0**2.0 * wl) + Cd / wl**3.0) # with dispersion
45 else:
46 ret = 2.0 * np.pi * ret # no dispersion
48 cr = np.cos(ret)
49 sr = np.sin(ret)
51 one = np.ones((dim))
52 zero = np.zeros((dim))
54 if tran == 0:
55 return np.transpose(
56 np.array(
57 [
58 [one, zero, zero, zero],
59 [zero, c2**2 + s2**2 * cr, c2 * s2 * (1 - cr), -s2 * sr],
60 [zero, c2 * s2 * (1 - cr), s2**2 + c2**2 * cr, c2 * sr],
61 [zero, s2 * sr, -c2 * sr, cr],
62 ]
63 ),
64 (2, 3, 0, 1),
65 )
66 else:
67 return np.transpose(
68 np.array(
69 [
70 [one, zero, zero, zero],
71 [zero, c2**2 + s2**2 * cr, c2 * s2 * (1 - cr), s2 * sr],
72 [zero, c2 * s2 * (1 - cr), s2**2 + c2**2 * cr, -c2 * sr],
73 [zero, -s2 * sr, c2 * sr, cr],
74 ]
75 ),
76 (2, 3, 0, 1),
77 )
80def susi_mmpol(pang, k0, k90):
81 """
82 Computes the Mueller matrix of a partial linear polarizer
83 acording to eq 4.48 in "del Toro Iniesta 2003"
85 INPUTS
86 pang:Scalar angle of the optical axis measured counterclockwise from positive x axis [deg]
87 k0: Array [nx,ny] of atenuation factors of the 0 deg E field component for each pixel
88 k90: Array [nx,ny] of atenuation factors of the 90 deg E field component for each pixel
90 OUTPUTS
91 [nx,ny,4,4] Mueller matrix of the partial linear polarizer defined by Inputs
92 """
93 dim = np.shape(k0)
94 phi = pang * np.pi / 180.0
95 gamma = 2.0 * k0 * k90 / (k0**2.0 + k90**2.0)
96 beta = (k0**2.0 - k90**2.0) / (k0**2.0 + k90**2.0)
97 c2 = np.cos(2.0 * phi)
98 s2 = np.sin(2.0 * phi)
100 one = np.ones((dim))
101 zero = np.zeros((dim))
102 mm = np.array(
103 [
104 [one, beta * c2, beta * s2, zero],
105 [beta * c2, c2**2 + gamma * s2**2, (1 - gamma) * c2 * s2, zero],
106 [beta * s2, (1 - gamma) * c2 * s2, s2**2 + gamma * c2**2, zero],
107 [zero, zero, zero, gamma],
108 ]
109 )
110 return np.transpose(
111 mm * 0.5 * (k0**2 + k90**2),
112 (2, 3, 0, 1),
113 )
116def susi_psg(nx, ny, nmeas, ret, phi, phip, alpha, ba):
117 """
118 Computes the output Stokes vector of a PSG (polarizer followed by retarder) assuming an unpolarized input
119 Stokes equal to [1,0,0,0]. Many outpus can be created if many PSG configurations (nmeas) are passed
120 (position angles).See Mathematica notebook susi_modulation_matrix.nb for details.
122 INPUTS
123 nx, ny, nmeas: size of the ouput
124 ret: PSG retardance [rad] per pixel, and per retarder position (nx,ny,nmeas)
125 phi: Twice the ret pos ang [Rad] (array[nmeas])
126 phip: Twice the pol pos ang [Rad] (array[nmeas])
127 alpha: 0.5*alpha in susi_mmpol (scalar)
128 ba: beta*alpha in susi_mmpol (scalar)
130 OUTPUTS
131 [nx,ny,nmeas,4] Output Stokes images for each nmeas
132 """
133 c2 = np.cos(phi)
134 s2 = np.sin(phi)
135 s4 = np.sin(2.0 * phi)
136 cr = np.cos(ret)
137 sr = np.sin(ret)
138 srr = np.sin(ret / 2.0) ** 2.0
140 c2p = np.cos(phip)
141 s2p = np.sin(phip)
142 c2d = np.cos(phip - phi)
143 s2d = np.sin(phip - phi)
145 # modeled calibration stokes vector for all polcal measurements and pixels [Nx,Ny,nmeas,4]
146 a = cr * np.expand_dims(ba * s2p * c2**2.0, (0, 1))
147 b = np.expand_dims(ba * s2p * s2**2.0, (0, 1))
148 c = srr * np.expand_dims(ba * c2p * s4, (0, 1))
149 return np.stack(
150 [
151 np.full((nx, ny, nmeas), alpha),
152 np.expand_dims(ba * c2d * c2, (0, 1)) - cr * np.expand_dims(ba * s2d * s2, (0, 1)),
153 a + b + c,
154 sr * np.expand_dims(-ba * s2d, (0, 1)),
155 ],
156 axis=3,
157 )
160def susi_mmrot(ang):
161 """
162 Ref.: Stenflo, J. 1994, Solar Magnetic Fields, Cuvillier, eq. (13.3)
163 Counter-clockwise rotation of the Stokes reference system
165 INPUTS:
166 ang: rotation angle (deg). It can be an array [nx.ny]
168 OUTPUTS
169 [nx,ny,4,4] Rotation Mueller matrix
170 """
171 dim = np.shape(ang)
172 i = ang * np.pi / 90.0 # 2*[Rad]
173 c2 = np.cos(i)
174 s2 = np.sin(i)
176 one = np.ones((dim))
177 zero = np.zeros((dim))
179 return np.transpose(
180 np.array([[one, zero, zero, zero], [zero, c2, s2, zero], [zero, -s2, c2, zero], [zero, zero, zero, one]]),
181 (2, 3, 0, 1),
182 )
185def susi_mmirr(n, k, aoi, aor, ref=1, ellip=1):
186 """
187 ; Returns the Mueller matrix of a reflection by a mirror, normalized to its 0,0 element.
188 ; Ref.: Stenflo, J. 1994, Solar Magnetic Fields, Cuvillier, section 13.3.6
189 ; Positive Stokes Q is perpendicular to the plane of incidence
190 ;
191 ; The Mueller matrix is normalized to 1 for incident unpolarized light. In
192 ; ellipsometric mode, and with reflectivities given (see arguments and keywords below)
193 ; the Mueller matrix is normalized to the absolute reflectivity for unpolarized light.
194 ;
195 ; The sign convention for Stokes Q, U, V is consistent with rotations on a Poincare
196 ; sphere with angle delta, and in a right handed (Q, U, V) coordinate system.
197 ;
198 INPUTS (they can be arrays [nx, ny], except ellip that must be a scalar)
200 n: (see ellip) real part of the complex refractive index of the mirror
201 k: (see ellip) imaginary part of the complex refractive index of the mirror
202 aoi: angle of incidence in the mirror [deg]
203 aor: angle of rotation applied to the output matrix [deg]
204 ref: reflectivity to scale the normalized Mueller matrix of the mirror
205 ellip: ret to interpret n and k as ellipsometric parameters psi [deg] and delta [deg], respectively
207 OUTPUTS
208 [nx,ny,4,4] Mueller matrix of the mirror
209 """
210 # if input are scalars, convert to 1x1 array filled with the same value
211 scalar_output = False
212 if not isinstance(n, np.ndarray):
213 n = np.array([[n]])
214 scalar_output = True
215 if not isinstance(k, np.ndarray):
216 k = np.array([[k]])
217 if not isinstance(aoi, np.ndarray):
218 aoi = np.array([[aoi]])
219 if not isinstance(aor, np.ndarray):
220 aor = np.array([[aor]])
221 if not isinstance(ref, np.ndarray):
222 ref = np.array([[ref]])
224 if ellip == 1:
225 rho = np.tan(n * np.pi / 180.0)
226 delta = -(np.pi + k * np.pi / 180.0)
227 else:
228 i = aoi * np.pi / 180.0
229 p = n**2.0 - k**2.0 - np.sin(i) ** 2.0
230 q = 2.0 * n * k
231 a = np.sqrt(p**2.0 + q**2.0)
232 rp = (1.0 / np.sqrt(2.0)) * (p + a) ** 0.5
233 rm = (1.0 / np.sqrt(2.0)) * (-p + a) ** 0.5
234 s = np.sin(i) * np.tan(i)
235 rho = ((a + s**2.0 - 2.0 * s * rp) / (a + s**2.0 + 2.0 * s * rp)) ** 0.5
236 delta = np.arctan(2.0 * s * rm / (a - s**2.0))
238 m11 = 0.5 * (1.0 + rho**2.0)
239 m12 = 0.5 * (1.0 - rho**2.0)
240 m33 = -rho * np.cos(delta)
241 m34 = -rho * np.sin(delta)
242 zero = np.zeros((np.shape(n)))
244 m = np.transpose(
245 np.array([[m11, m12, zero, zero], [m12, m11, zero, zero], [zero, zero, m33, m34], [zero, zero, -m34, m33]]),
246 (2, 3, 0, 1),
247 )
249 if np.any(aor != 0):
250 m = susi_mmrot(-aor) @ m @ susi_mmrot(aor)
252 m /= np.expand_dims(m[:, :, 0, 0], (2, 3)) # norm
254 if ellip == 1 and np.any(ref != 1): # norm to absolute reflectivity
255 m *= np.expand_dims(ref, (2, 3))
257 if scalar_output:
258 return m[0, 0]
260 return m