-
Notifications
You must be signed in to change notification settings - Fork 0
/
DRRP_Functions.py
231 lines (178 loc) · 10.8 KB
/
DRRP_Functions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
# These are the functions needed to run the traditional Dual Rotating Retarder Polarimeter measurements and analysis
# This method uses measurements of intensity with a linear polarizer in the polarization state generator
# One channel of a Wollaston prism can serve as a linear polarizer in this case
import numpy as np
from numpy.linalg import inv
from astropy.io import fits
import os
import re
from scipy.optimize import curve_fit
# Mueller matrix for a linear polarizer, with angle a between transmission axis and horizontal (radians)
def linear_polarizer(a):
M01 = np.cos(2*a)
M02 = np.sin(2*a)
M10 = np.cos(2*a)
M11 = np.cos(2*a)**2
M12 = np.cos(2*a)*np.sin(2*a)
M20 = np.sin(2*a)
M21 = np.cos(2*a)*np.sin(2*a)
M22 = np.sin(2*a)**2
return 0.5*np.array([[1, M01, M02, 0],
[M10, M11, M12, 0],
[M20, M21, M22, 0],
[0, 0, 0, 0]])
# Mueller matrix for a linear retarder (waveplate). Angle of fast axis a, retardance r in radians
def linear_retarder(a, r):
M11 = np.cos(2*a)**2 + np.cos(r)*np.sin(2*a)**2
M12 = np.cos(2*a)*np.sin(2*a)*(1-np.cos(r))
M13 = -np.sin(2*a)*np.sin(r)
M21 = M12
M22 = np.sin(2*a)**2 + np.cos(2*a)**2*np.cos(r)
M23 = np.cos(2*a)*np.sin(r)
M31 = -M13
M32 = -M23
M33 = np.cos(r)
return np.array([[1, 0, 0, 0],
[0, M11, M12, M13],
[0, M21, M22, M23],
[0, M31, M32, M33]])
# Sorting function for extracting filenames based on last number in the filename (the angle of rotation)
def extract_number(filename):
match = re.findall(r'\d+(?:\.\d+)?', filename)
if match:
return float(match[-1])
# Function to subtract dark frames from raw frames. The new reduced images are saved to a different folder
def dark_subtraction(image_file, dark_file, old_directory, new_directory):
# Open the dark image and extract pixel values
fits.open(dark_file)
dark = fits.getdata(dark_file)
dark_median = np.median(dark, axis=0)
# Search through the desired raw data folder
for filename in os.listdir(old_directory):
if filename.startswith(image_file): # Call specific files starting with the desired name
with fits.open(os.path.join(old_directory, filename)) as hdul:
img_data = hdul[0].data
img_median = np.median(img_data, axis=0)
reduced_data = img_median - dark_median
# Save the newly reduced image to a reduced data folder
new_filename = f"Reduced_{filename}"
new_filepath = os.path.join(new_directory, new_filename)
fits.writeto(new_filepath, reduced_data, overwrite=True)
# Get intensity values from each spot in the reduced images. reduced_filename should just be the start of the name (leave out the last number, the angle).
def extract_intensities(reduced_filename, reduced_folder, lcenter, rcenter, maxradius, cutoff=5000):
I_left = np.array([])
I_right = np.array([])
bad_indices = np.array([])
for filename in sorted(os.listdir(reduced_folder), key = extract_number):
if filename.startswith(reduced_filename):
with fits.open(os.path.join(reduced_folder, filename)) as hdul:
reduced_img_data = hdul[0].data
ys, xs, = np.indices(reduced_img_data.shape)
lradius = np.sqrt((ys-lcenter[0])**2+(xs-lcenter[1])**2)
rradius = np.sqrt((ys-rcenter[0])**2+(xs-rcenter[1])**2)
lbackground_mask = (lradius > 20) & (lradius < 26)
rbackground_mask = (rradius > 20) & (rradius < 26) # Index the background around each spot, take the median value
background_lmedian = np.median(reduced_img_data[lbackground_mask])
background_rmedian = np.median(reduced_img_data[rbackground_mask])
lflux = np.sum(reduced_img_data[lradius < maxradius] - background_lmedian) # Now take the flux with the background mask subtracted
rflux = np.sum(reduced_img_data[rradius < maxradius] - background_rmedian)
I_left = np.append(I_left, lflux)
I_right = np.append(I_right, rflux)
if lflux+rflux < cutoff:
print("Warning: low flux detected, check the image " + filename + ", index: " + str(sorted(os.listdir(reduced_folder), key = extract_number).index(filename)))
bad_indices = np.append(bad_indices, sorted(os.listdir(reduced_folder), key = extract_number).index(filename))
else:
continue
return I_left, I_right, bad_indices
# Gives the condition number of eventual Mueller matrix (made by Jaren)
def condition_number(matrix):
minv = np.linalg.pinv(matrix)
# Compute maximum norm
norm = np.linalg.norm(matrix, ord=np.inf)
ninv = np.linalg.norm(minv, ord=np.inf)
return norm*ninv
# Function to compute the Mueller matrix of a sample based on DRRP intensity measurements and calibration parameters
def calibrated_full_mueller_polarimetry(thetas, a1, a2, w1, w2, r1, r2, I_meas=1, LPA_angle=0, return_condition_number=False, M_in=None):
nmeas = len(thetas)
Wmat = np.zeros([nmeas, 16])
Pmat = np.zeros([nmeas])
th = thetas
for i in range(nmeas):
# Mueller Matrix of generator (linear polarizer and a quarter wave plate)
Mg = linear_retarder(th[i]+w1, np.pi/2+r1) @ linear_polarizer(0+a1)
# Mueller Matrix of analyzer (one channel of the Wollaston prism is treated as a linear polarizer. The right spot is horizontal (0) and the left spot is vertical(pi/2))
Ma = linear_polarizer(LPA_angle+a2) @ linear_retarder(th[i]*5+w2, np.pi/2+r2)
# Data reduction matrix. Taking the 0 index ensures that intensity is the output
Wmat[i,:] = np.kron(Ma[0,:], Mg[:,0])
# M_in is some example Mueller matrix. Providing this input will test theoretical Mueller matrix. Otherwise, the raw data is used
if M_in is not None:
Pmat[i] = (Ma[0,:] @ M_in @ Mg[:,0]) * I_meas
else:
Pmat[i] = I_meas[i]
# Compute Mueller matrix using Moore-Penrose pseudo invervse
M = np.linalg.pinv(Wmat) @ Pmat
M = np.reshape(M,[4,4])
if return_condition_number == True:
return M, condition_number(Wmat)
else:
return M
# Define the identity matrix and other matrices which are useful for the Mueller calculus
M_identity = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
A = np.array([1, 0, 0, 0])
B = np.array([[1], [0], [0], [0]])
# This is the full Mueller matrix equation for our setup. The output is a list, useful for curve fitting. Variables with 1 refer to the generator, 2 refers to analyzer.
def calibration_function(t, a1, a2, w1, w2, r1, r2):
prediction = [None]*len(t)
for i in range(len(t)):
prediction[i] = float(A @ linear_polarizer(a2) @ linear_retarder(5*t[i]+w2, np.pi/2+r2) @ M_identity @ linear_retarder(t[i]+w1, np.pi/2+r1) @ linear_polarizer(a1) @ B)
return prediction
# Calibration function designed for data from the left spot, which is the vertial alignment. This changes the angle of the analyzing LP
def vertical_calibration_function(t, a1, a2, w1, w2, r1, r2):
prediction = [None]*len(t)
for i in range(len(t)):
prediction[i] = float(A @ linear_polarizer(a2+np.pi/2) @ linear_retarder(5*t[i]+w2, np.pi/2+r2) @ M_identity @ linear_retarder(t[i]+w1, np.pi/2+r1) @ linear_polarizer(a1) @ B)
return prediction
# Basically the same as above, but with an optional input matrix to simulate data
def output_simulation_function(t, a1, a2, w1, w2, r1, r2, LPA_angle=0, M_in=None):
if M_in is None:
M = M_identity
else:
M = M_in
prediction = [None]*len(t)
for i in range(len(t)):
prediction[i] = float(A @ linear_polarizer(LPA_angle+a2) @ linear_retarder(5*t[i]+w2, np.pi/2+r2) @ M @ linear_retarder(t[i]+w1, np.pi/2+r1) @ linear_polarizer(a1) @ B)
return prediction
# After testing, each of the above functions works individually. Now combine them into one function to rule them all
# Finds the mueller matrix derived from each channel separately, then averages the two retardances found this way
# First three inputs must come from the calibration data, last three inputs correspond to the HWP sample
def ultimate_polarimetry(cal_angles, cal_left_intensity, cal_right_intensity, sample_angles, sample_left_intensity, sample_right_intensity):
initial_guess = [0, 0, 0, 0, 0, 0]
parameter_bounds = ([-np.pi, -np.pi, -np.pi, -np.pi, -np.pi/2, -np.pi/2], [np.pi, np.pi, np.pi, np.pi, np.pi/2, np.pi/2])
# Find parameters from calibration of the left spot
lnormalized_intensity = cal_left_intensity/(2*max(cal_left_intensity))
lpopt, lpcov = curve_fit(vertical_calibration_function, cal_angles, lnormalized_intensity, p0=initial_guess, bounds=parameter_bounds)
print(lpopt, "Left parameters for a1, a2, w1, w2, r1, and r2. 1 for generator, 2 for analyzer")
#print(np.sqrt(np.diag(lpcov)))
# Find parameters from calibration of the right spot
rnormalized_intensity = cal_right_intensity/(2*max(cal_right_intensity))
rpopt, rpcov = curve_fit(calibration_function, cal_angles, rnormalized_intensity, p0=initial_guess, bounds=parameter_bounds)
print(rpopt, "Right parameters for a1, a2, w1, w2, r1, and r2. 1 for generator, 2 for analyzer")
#print(np.sqrt(np.diag(rpcov)))
# Optional print the calibration matrices (should be close to identity) to see how well the parameters compensate
MlCal = calibrated_full_mueller_polarimetry(cal_angles, lpopt[0], lpopt[1], lpopt[2], lpopt[3], lpopt[4], lpopt[5], cal_left_intensity, LPA_angle=np.pi/2)
print(MlCal/MlCal.max(), ' Left calibration')
MrCal = calibrated_full_mueller_polarimetry(cal_angles, rpopt[0], rpopt[1], rpopt[2], rpopt[3], rpopt[4], rpopt[5], cal_right_intensity)
print(MrCal/MrCal.max(), ' Right calibration')
# Use the parameters found above from curve fitting to construct the actual Mueller matrix of the sample
Ml = calibrated_full_mueller_polarimetry(sample_angles, lpopt[0], lpopt[1], lpopt[2], lpopt[3], lpopt[4], lpopt[5], sample_left_intensity, LPA_angle=np.pi/2)
Ml = Ml/Ml.max()
Mr = calibrated_full_mueller_polarimetry(sample_angles, rpopt[0], rpopt[1], rpopt[2], rpopt[3], rpopt[4], rpopt[5], sample_right_intensity)
Mr = Mr/Mr.max()
np.set_printoptions(suppress=True)
# Extract retardance from the last entry of the mueller matrix, which should just be cos(phi)
lretardance = np.arccos(Ml[3,3])/(2*np.pi)
rretardance = np.arccos(Mr[3,3])/(2*np.pi)
print(lretardance, ' This is the retardance found from the left spot')
print(rretardance, ' This is the retardance found from the right spot')
avg_retardance = (lretardance+rretardance)/2
return Ml, Mr, avg_retardance