Mercurial > public > atmospheric_lidar_ghk / file comparison

comparison: lidar_correction_ghk.py

lidar_correction_ghk.py

changeset 1
5dffdb7caec9
parent 0
d56e90f31b9e
child 5
9d5aa2422c02
equal deleted inserted replaced
0:d56e90f31b9e 1:5dffdb7caec9
1 # -*- coding: utf-8 -*-
2 """
3 Copyright 2016 Volker Freudenthaler
4
5 Licensed under the EUPL, Version 1.1 only (the "Licence").
6
7 You may not use this work except in compliance with the Licence.
8 A copy of the licence is distributed with the code. Alternatively, you may obtain
9 a copy of the Licence at:
10
11 https://joinup.ec.europa.eu/community/eupl/og_page/eupl
12
13 Unless required by applicable law or agreed to in writing, software distributed
14 under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS
15 OF ANY KIND, either express or implied. See the Licence for the specific language governing
16 permissions and limitations under the Licence.
17
18 Equation reference: http://www.atmos-meas-tech-discuss.net/amt-2015-338/amt-2015-338.pdf
19 With equations code from Appendix C
20 Python 3.4.2
21 """
22 #!/usr/bin/env python3
23 from __future__ import print_function
24 #import math
25 import numpy as np
26 import sys
27 import os
28
29 #import seaborn as sns
30 import matplotlib.pyplot as plt
31 from time import clock
32
33 #from matplotlib.backends.backend_pdf import PdfPages
34 #pdffile = '{}.pdf'.format('path')
35 #pp = PdfPages(pdffile)
36 ## pp.savefig can be called multiple times to save to multiple pages
37 #pp.savefig()
38 #pp.close()
39
40 from contextlib import contextmanager
41 @contextmanager
42 def redirect_stdout(new_target):
43 old_target, sys.stdout = sys.stdout, new_target # replace sys.stdout
44 try:
45 yield new_target # run some code with the replaced stdout
46 finally:
47 sys.stdout.flush()
48 sys.stdout = old_target # restore to the previous value
49 '''
50 real_raw_input = vars(__builtins__).get('raw_input',input)
51 '''
52 try:
53 import __builtin__
54 input = getattr(__builtin__, 'raw_input')
55 except (ImportError, AttributeError):
56 pass
57
58 from distutils.util import strtobool
59 def user_yes_no_query(question):
60 sys.stdout.write('%s [y/n]\n' % question)
61 while True:
62 try:
63 return strtobool(input().lower())
64 except ValueError:
65 sys.stdout.write('Please respond with \'y\' or \'n\'.\n')
66
67 #if user_yes_no_query('want to exit?') == 1: sys.exit()
68
69 '''
70 ## {{{ http://code.activestate.com/recipes/577058/ (r2)
71 def query_yes_no(question, default="yes"):
72 valid = {"yes":"yes", "y":"yes", "ye":"yes",
73 "no":"no", "n":"no"}
74 if default == None:
75 prompt = " [y/n] "
76 elif default == "yes":
77 prompt = " [Y/n] "
78 elif default == "no":
79 prompt = " [y/N] "
80 else:
81 raise ValueError("invalid default answer: '%s'" % default)
82
83 while 1:
84 sys.stdout.write(question + prompt)
85 choice = input().lower()
86 if default is not None and choice == '':
87 return default
88 elif choice in valid.keys():
89 return valid[choice]
90 else:
91 sys.stdout.write("Please respond with 'yes' or 'no' "\
92 "(or 'y' or 'n').\n")
93 ## end of http://code.activestate.com/recipes/577058/ }}}
94 '''
95 abspath = os.path.abspath(__file__)
96 dname = os.path.dirname(abspath)
97 fname = os.path.basename(abspath)
98 os.chdir(dname)
99
100 #PrintToOutputFile = True
101
102 sqr05 = 0.5**0.5
103
104 # ---- Initial definition of variables; the actual values will be read in with exec(open('./optic_input.py').read()) below
105 LID = "internal"
106 EID = "internal"
107 # --- IL Laser IL and +-Uncertainty
108 bL = 1. #degree of linear polarization; default 1
109 RotL, dRotL, nRotL = 0.0, 0.0, 1 #alpha; rotation of laser polarization in degrees; default 0
110 # --- ME Emitter and +-Uncertainty
111 DiE, dDiE, nDiE = 0., 0.00, 1 # Diattenuation
112 TiE = 1. # Unpolarized transmittance
113 RetE, dRetE, nRetE = 0., 180.0, 0 # Retardance in degrees
114 RotE, dRotE, nRotE = 0., 0.0, 0 # beta: Rotation of optical element in degrees
115 # --- MO Receiver Optics including telescope
116 DiO, dDiO, nDiO = -0.055, 0.003, 1
117 TiO = 0.9
118 RetO, dRetO, nRetO = 0., 180.0, 2
119 RotO, dRotO, nRotO = 0., 0.1, 1 #gamma
120 # --- PBS MT transmitting path defined with (TS,TP); and +-Uncertainty
121 TP, dTP, nTP = 0.98, 0.02, 1
122 TS, dTS, nTS = 0.001, 0.001, 1
123 TiT = 0.5 * (TP + TS)
124 DiT = (TP-TS)/(TP+TS)
125 # PolFilter
126 RetT, dRetT, nRetT = 0., 180., 0
127 ERaT, dERaT, nERaT = 0.001, 0.001, 1
128 RotaT, dRotaT, nRotaT = 0., 3., 1
129 DaT = (1-ERaT)/(1+ERaT)
130 TaT = 0.5*(1+ERaT)
131 # --- PBS MR reflecting path defined with (RS,RP); and +-Uncertainty
132 RS, dRS, nRS = 1 - TS, 0., 0
133 RP, dRP, nRP = 1 - TP, 0., 0
134 TiR = 0.5 * (RP + RS)
135 DiR = (RP-RS)/(RP+RS)
136 # PolFilter
137 RetR, dRetR, nRetR = 0., 180., 0
138 ERaR, dERaR, nERaR = 0.001, 0.001, 1
139 RotaR,dRotaR,nRotaR = 90., 3., 1
140 DaR = (1-ERaR)/(1+ERaR)
141 TaR = 0.5*(1+ERaR)
142
143 # Parellel signal detected in the transmitted channel => Y = 1, or in the reflected channel => Y = -1
144 Y = -1.
145
146 # Calibrator = type defined by matrix values
147 LocC = 4 # location of calibrator: behind laser = 1; behind emitter = 2; before receiver = 3; before PBS = 4
148
149 TypeC = 3 # linear polarizer calibrator
150 # example with extinction ratio 0.001
151 DiC, dDiC, nDiC = 1.0, 0., 0 # ideal 1.0
152 TiC = 0.5 # ideal 0.5
153 RetC, dRetC, nRetC = 0., 0., 0
154 RotC, dRotC, nRotC = 0.0, 0.1, 0 #constant calibrator offset epsilon
155 RotationErrorEpsilonForNormalMeasurements = False # is in general False for TypeC == 3 calibrator
156
157 # Rotation error without calibrator: if False, then epsilon = 0 for normal measurements
158 RotationErrorEpsilonForNormalMeasurements = True
159
160 # LDRCal assumed atmospheric linear depolarization ratio during the calibration measurements (first guess)
161 LDRCal0,dLDRCal,nLDRCal= 0.25, 0.04, 1
162 LDRCal = LDRCal0
163 # measured LDRm will be corrected with calculated parameters
164 LDRmeas = 0.015
165 # LDRtrue for simulation of measurement => LDRsim
166 LDRtrue = 0.5
167 LDRtrue2 = 0.004
168
169 # Initialize other values to 0
170 ER, nER, dER = 0.001, 0, 0.001
171 K = 0.
172 Km = 0.
173 Kp = 0.
174 LDRcorr = 0.
175 Eta = 0.
176 Ir = 0.
177 It = 0.
178 h = 1.
179
180 Loc = ['', 'behind laser', 'behind emitter', 'before receiver', 'before PBS']
181 Type = ['', 'mechanical rotator', 'hwp rotator', 'linear polarizer', 'qwp rotator', 'circular polarizer', 'real HWP +-22.5°']
182 dY = ['reflected channel', '', 'transmitted channel']
183
184 # end of initial definition of variables
185 # *******************************************************************************************************************************
186
187 # --- Read actual lidar system parameters from ./optic_input.py (must be in the same directory)
188
189 #InputFile = 'optic_input_ver6e_POLIS_355.py'
190 #InputFile = 'optic_input_ver6e_POLIS_355_JA.py'
191 #InputFile = 'optic_input_ver6c_POLIS_532.py'
192 #InputFile = 'optic_input_ver6e_POLIS_532.py'
193 #InputFile = 'optic_input_ver8c_POLIS_532.py'
194 #InputFile = 'optic_input_ver6e_MUSA.py'
195 #InputFile = 'optic_input_ver6e_MUSA_JA.py'
196 #InputFile = 'optic_input_ver6e_PollyXTSea.py'
197 #InputFile = 'optic_input_ver6e_PollyXTSea_JA.py'
198 #InputFile = 'optic_input_ver6e_PollyXT_RALPH.py'
199 #InputFile = 'optic_input_ver8c_PollyXT_RALPH.py'
200 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_2.py'
201 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_3.py'
202 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_4.py'
203 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_5.py'
204 #InputFile = 'optic_input_ver8c_PollyXT_RALPH_6.py'
205 InputFile = 'optic_input_ver8c_PollyXT_RALPH_7.py'
206 #InputFile = 'optic_input_ver8a_MOHP_DPL_355.py'
207 #InputFile = 'optic_input_ver9_MOHP_DPL_355.py'
208 #InputFile = 'optic_input_ver6e_RALI.py'
209 #InputFile = 'optic_input_ver6e_RALI_JA.py'
210 #InputFile = 'optic_input_ver6e_RALI_new.py'
211 #InputFile = 'optic_input_ver6e_RALI_act.py'
212 #InputFile = 'optic_input_ver6e_MULHACEN.py'
213 #InputFile = 'optic_input_ver6e_MULHACEN_JA.py'
214 #InputFile = 'optic_input_ver6e-IPRAL.py'
215 #InputFile = 'optic_input_ver6e-IPRAL_JA.py'
216 #InputFile = 'optic_input_ver6e-LB21.py'
217 #InputFile = 'optic_input_ver6e-LB21_JA.py'
218 #InputFile = 'optic_input_ver6e_Bertha_b_355.py'
219 #InputFile = 'optic_input_ver6e_Bertha_b_532.py'
220 #InputFile = 'optic_input_ver6e_Bertha_b_1064.py'
221
222 '''
223 print("From ", dname)
224 print("Running ", fname)
225 print("Reading input file ", InputFile, " for")
226 '''
227 input_path = os.path.join('.', 'system_settings', InputFile)
228 # this works with Python 2 - and 3?
229 exec(open(input_path).read(), globals())
230 # end of read actual system parameters
231
232 # --- Manual Parameter Change ---
233 # (use for quick parameter changes without changing the input file )
234 #DiO = 0.
235 #LDRtrue = 0.45
236 #LDRtrue2 = 0.004
237 #Y = -1
238 #LocC = 4 #location of calibrator: 1 = behind laser; 2 = behind emitter; 3 = before receiver; 4 = before PBS
239 ##TypeC = 6 Don't change the TypeC here
240 #RotationErrorEpsilonForNormalMeasurements = True
241 #LDRCal = 0.25
242 #bL = 0.8
243 ## --- Errors
244 RotL0, dRotL, nRotL = RotL, dRotL, nRotL
245
246 DiE0, dDiE, nDiE = DiE, dDiE, nDiE
247 RetE0, dRetE, nRetE = RetE, dRetE, nRetE
248 RotE0, dRotE, nRotE = RotE, dRotE, nRotE
249
250 DiO0, dDiO, nDiO = DiO, dDiO, nDiO
251 RetO0, dRetO, nRetO = RetO, dRetO, nRetO
252 RotO0, dRotO, nRotO = RotO, dRotO, nRotO
253
254 DiC0, dDiC, nDiC = DiC, dDiC, nDiC
255 RetC0, dRetC, nRetC = RetC, dRetC, nRetC
256 RotC0, dRotC, nRotC = RotC, dRotC, nRotC
257
258 TP0, dTP, nTP = TP, dTP, nTP
259 TS0, dTS, nTS = TS, dTS, nTS
260 RetT0, dRetT, nRetT = RetT, dRetT, nRetT
261
262 ERaT0, dERaT, nERaT = ERaT, dERaT, nERaT
263 RotaT0,dRotaT,nRotaT= RotaT,dRotaT,nRotaT
264
265 RP0, dRP, nRP = RP, dRP, nRP
266 RS0, dRS, nRS = RS, dRS, nRS
267 RetR0, dRetR, nRetR = RetR, dRetR, nRetR
268
269 ERaR0, dERaR, nERaR = ERaR, dERaR, nERaR
270 RotaR0,dRotaR,nRotaR= RotaR,dRotaR,nRotaR
271
272 LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,nLDRCal
273 #LDRCal0,dLDRCal,nLDRCal=LDRCal,dLDRCal,0
274 # ---------- End of manual parameter change
275
276 RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC = RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0
277 TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR = TP0, TS0, RP0, RS0 , ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0
278 LDRCal = LDRCal0
279 DTa0, TTa0, DRa0, TRa0, LDRsimx, LDRCorr = 0,0,0,0,0,0
280
281 TiT = 0.5 * (TP + TS)
282 DiT = (TP-TS)/(TP+TS)
283 ZiT = (1. - DiT**2)**0.5
284 TiR = 0.5 * (RP + RS)
285 DiR = (RP-RS)/(RP+RS)
286 ZiR = (1. - DiR**2)**0.5
287
288 # --------------------------------------------------------
289 def Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS, RP, RS, ERaT, RotaT, RetT, ERaR, RotaR, RetR, LDRCal):
290 # ---- Do the calculations of bra-ket vectors
291 h = -1. if TypeC == 2 else 1
292 # from input file: assumed LDRCal for calibration measurements
293 aCal = (1.-LDRCal)/(1+LDRCal)
294 # from input file: measured LDRm and true LDRtrue, LDRtrue2 =>
295 #ameas = (1.-LDRmeas)/(1+LDRmeas)
296 atrue = (1.-LDRtrue)/(1+LDRtrue)
297 #atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)
298
299 # angles of emitter and laser and calibrator and receiver optics
300 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon
301 S2a = np.sin(2*np.deg2rad(RotL))
302 C2a = np.cos(2*np.deg2rad(RotL))
303 S2b = np.sin(2*np.deg2rad(RotE))
304 C2b = np.cos(2*np.deg2rad(RotE))
305 S2ab = np.sin(np.deg2rad(2*RotL-2*RotE))
306 C2ab = np.cos(np.deg2rad(2*RotL-2*RotE))
307 S2g = np.sin(np.deg2rad(2*RotO))
308 C2g = np.cos(np.deg2rad(2*RotO))
309
310 # Laser with Degree of linear polarization DOLP = bL
311 IinL = 1.
312 QinL = bL
313 UinL = 0.
314 VinL = (1. - bL**2)**0.5
315
316 # Stokes Input Vector rotation Eq. E.4
317 A = C2a*QinL - S2a*UinL
318 B = S2a*QinL + C2a*UinL
319 # Stokes Input Vector rotation Eq. E.9
320 C = C2ab*QinL - S2ab*UinL
321 D = S2ab*QinL + C2ab*UinL
322
323 # emitter optics
324 CosE = np.cos(np.deg2rad(RetE))
325 SinE = np.sin(np.deg2rad(RetE))
326 ZiE = (1. - DiE**2)**0.5
327 WiE = (1. - ZiE*CosE)
328
329 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4
330 # b = beta
331 IinE = (IinL + DiE*C)
332 QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL))
333 UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL))
334 VinE = (-ZiE*SinE*D + ZiE*CosE*VinL)
335
336 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
337 IinF = IinE
338 QinF = aCal*QinE
339 UinF = -aCal*UinE
340 VinF = (1.-2.*aCal)*VinE
341
342 # receiver optics
343 CosO = np.cos(np.deg2rad(RetO))
344 SinO = np.sin(np.deg2rad(RetO))
345 ZiO = (1. - DiO**2)**0.5
346 WiO = (1. - ZiO*CosO)
347
348 # calibrator
349 CosC = np.cos(np.deg2rad(RetC))
350 SinC = np.sin(np.deg2rad(RetC))
351 ZiC = (1. - DiC**2)**0.5
352 WiC = (1. - ZiC*CosC)
353
354 # Stokes Input Vector before the polarising beam splitter Eq. E.31
355 A = C2g*QinE - S2g*UinE
356 B = S2g*QinE + C2g*UinE
357
358 IinP = (IinE + DiO*aCal*A)
359 QinP = (C2g*DiO*IinE + aCal*QinE - S2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE))
360 UinP = (S2g*DiO*IinE - aCal*UinE + C2g*(WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE))
361 VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
362
363 #-------------------------
364 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue
365 #F11sim = TiO*(IinE + DiO*atrue*A)/IinE
366 #-------------------------
367
368 # For PollyXT
369 # analyser
370 #RS = 1 - TS
371 #RP = 1 - TP
372
373 TiT = 0.5 * (TP + TS)
374 DiT = (TP-TS)/(TP+TS)
375 ZiT = (1. - DiT**2)**0.5
376 TiR = 0.5 * (RP + RS)
377 DiR = (RP-RS)/(RP+RS)
378 ZiR = (1. - DiR**2)**0.5
379 CosT = np.cos(np.deg2rad(RetT))
380 SinT = np.sin(np.deg2rad(RetT))
381 CosR = np.cos(np.deg2rad(RetR))
382 SinR = np.sin(np.deg2rad(RetR))
383
384 DaT = (1-ERaT)/(1+ERaT)
385 DaR = (1-ERaR)/(1+ERaR)
386 TaT = 0.5*(1+ERaT)
387 TaR = 0.5*(1+ERaR)
388
389 S2aT = np.sin(np.deg2rad(h*2*RotaT))
390 C2aT = np.cos(np.deg2rad(2*RotaT))
391 S2aR = np.sin(np.deg2rad(h*2*RotaR))
392 C2aR = np.cos(np.deg2rad(2*RotaR))
393
394 # Aanalyzer As before the PBS Eq. D.5
395 ATP1 = (1+C2aT*DaT*DiT)
396 ATP2 = Y*(DiT+C2aT*DaT)
397 ATP3 = Y*S2aT*DaT*ZiT*CosT
398 ATP4 = S2aT*DaT*ZiT*SinT
399 ATP = np.array([ATP1,ATP2,ATP3,ATP4])
400
401 ARP1 = (1+C2aR*DaR*DiR)
402 ARP2 = Y*(DiR+C2aR*DaR)
403 ARP3 = Y*S2aR*DaR*ZiR*CosR
404 ARP4 = S2aR*DaR*ZiR*SinR
405 ARP = np.array([ARP1,ARP2,ARP3,ARP4])
406
407 DTa = ATP2*Y/ATP1
408 DRa = ARP2*Y/ARP1
409
410 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
411 if LocC == 4: # Calibrator before the PBS
412 #print("Calibrator location not implemented yet")
413
414 #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
415 #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC))
416 S2e = np.sin(np.deg2rad(h*2*RotC))
417 C2e = np.cos(np.deg2rad(2*RotC))
418 # rotated AinP by epsilon Eq. C.3
419 ATP2e = C2e*ATP2 + S2e*ATP3
420 ATP3e = C2e*ATP3 - S2e*ATP2
421 ARP2e = C2e*ARP2 + S2e*ARP3
422 ARP3e = C2e*ARP3 - S2e*ARP2
423 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4])
424 ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4])
425 # Stokes Input Vector before the polarising beam splitter Eq. E.31
426 A = C2g*QinE - S2g*UinE
427 B = S2g*QinE + C2g*UinE
428 #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
429 Co = ZiO*SinO*VinE
430 Ca = (WiO*B - 2*ZiO*SinO*VinE)
431 #C = Co + aCal*Ca
432 #IinP = (IinE + DiO*aCal*A)
433 #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
434 #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
435 #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
436 IinPo = IinE
437 QinPo = (C2g*DiO*IinE - S2g*Co)
438 UinPo = (S2g*DiO*IinE + C2g*Co)
439 VinPo = ZiO*CosO*VinE
440
441 IinPa = DiO*A
442 QinPa = QinE - S2g*Ca
443 UinPa = -UinE + C2g*Ca
444 VinPa = ZiO*(SinO*B - 2*CosO*VinE)
445
446 IinP = IinPo + aCal*IinPa
447 QinP = QinPo + aCal*QinPa
448 UinP = UinPo + aCal*UinPa
449 VinP = VinPo + aCal*VinPa
450 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3
451 #QinPe = C2e*QinP + S2e*UinP
452 #UinPe = C2e*UinP - S2e*QinP
453 QinPoe = C2e*QinPo + S2e*UinPo
454 UinPoe = C2e*UinPo - S2e*QinPo
455 QinPae = C2e*QinPa + S2e*UinPa
456 UinPae = C2e*UinPa - S2e*QinPa
457 QinPe = C2e*QinP + S2e*UinP
458 UinPe = C2e*UinP - S2e*QinP
459
460 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
461 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
462 # parameters for calibration with aCal
463 AT = ATP1*IinP + h*ATP4*VinP
464 BT = ATP3e*QinP - h*ATP2e*UinP
465 AR = ARP1*IinP + h*ARP4*VinP
466 BR = ARP3e*QinP - h*ARP2e*UinP
467 # Correction paremeters for normal measurements; they are independent of LDR
468 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
469 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
470 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
471 GT = np.dot(ATP,IS1)
472 GR = np.dot(ARP,IS1)
473 HT = np.dot(ATP,IS2)
474 HR = np.dot(ARP,IS2)
475 else:
476 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
477 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
478 GT = np.dot(ATPe,IS1)
479 GR = np.dot(ARPe,IS1)
480 HT = np.dot(ATPe,IS2)
481 HR = np.dot(ARPe,IS2)
482 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
483 # parameters for calibration with aCal
484 AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP)
485 BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
486 AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP)
487 BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
488 # Correction paremeters for normal measurements; they are independent of LDR
489 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
490 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
491 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
492 GT = np.dot(ATP,IS1)
493 GR = np.dot(ARP,IS1)
494 HT = np.dot(ATP,IS2)
495 HR = np.dot(ARP,IS2)
496 else:
497 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
498 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
499 GT = np.dot(ATPe,IS1e)
500 GR = np.dot(ARPe,IS1e)
501 HT = np.dot(ATPe,IS2e)
502 HR = np.dot(ARPe,IS2e)
503 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
504 # parameters for calibration with aCal
505 AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP)
506 BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
507 AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP)
508 BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
509 # Correction paremeters for normal measurements; they are independent of LDR
510 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
511 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
512 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
513 GT = np.dot(ATP,IS1)
514 GR = np.dot(ARP,IS1)
515 HT = np.dot(ATP,IS2)
516 HR = np.dot(ARP,IS2)
517 else:
518 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
519 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
520 GT = np.dot(ATPe,IS1e)
521 GR = np.dot(ARPe,IS1e)
522 HT = np.dot(ATPe,IS2e)
523 HR = np.dot(ARPe,IS2e)
524 else:
525 print("Calibrator not implemented yet")
526 sys.exit()
527
528 elif LocC == 3: # C before receiver optics Eq.57
529
530 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC))
531 #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
532 S2e = np.sin(np.deg2rad(2*RotC))
533 C2e = np.cos(np.deg2rad(2*RotC))
534
535 # As with C before the receiver optics (rotated_diattenuator_X22x5deg.odt)
536 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
537 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
538 AF3 = np.array([S2g*DiO,S2g*C2g*WiO,1-C2g**2*WiO,C2g*ZiO*SinO])
539 AF4 = np.array([0,S2g*SinO,-C2g*SinO,CosO])
540
541 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
542 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
543 ATF2 = ATF[1]
544 ATF3 = ATF[2]
545 ARF2 = ARF[1]
546 ARF3 = ARF[2]
547
548 # rotated AinF by epsilon
549 ATF1 = ATF[0]
550 ATF4 = ATF[3]
551 ATF2e = C2e*ATF[1] + S2e*ATF[2]
552 ATF3e = C2e*ATF[2] - S2e*ATF[1]
553 ARF1 = ARF[0]
554 ARF4 = ARF[3]
555 ARF2e = C2e*ARF[1] + S2e*ARF[2]
556 ARF3e = C2e*ARF[2] - S2e*ARF[1]
557
558 ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4])
559 ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4])
560
561 QinEe = C2e*QinE + S2e*UinE
562 UinEe = C2e*UinE - S2e*QinE
563
564 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
565 IinF = IinE
566 QinF = aCal*QinE
567 UinF = -aCal*UinE
568 VinF = (1.-2.*aCal)*VinE
569
570 IinFo = IinE
571 QinFo = 0.
572 UinFo = 0.
573 VinFo = VinE
574
575 IinFa = 0.
576 QinFa = QinE
577 UinFa = -UinE
578 VinFa = -2.*VinE
579
580 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3
581 QinFe = C2e*QinF + S2e*UinF
582 UinFe = C2e*UinF - S2e*QinF
583 QinFoe = C2e*QinFo + S2e*UinFo
584 UinFoe = C2e*UinFo - S2e*QinFo
585 QinFae = C2e*QinFa + S2e*UinFa
586 UinFae = C2e*UinFa - S2e*QinFa
587
588 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
589 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
590 # parameters for calibration with aCal
591 AT = ATF1*IinF + ATF4*h*VinF
592 BT = ATF3e*QinF - ATF2e*h*UinF
593 AR = ARF1*IinF + ARF4*h*VinF
594 BR = ARF3e*QinF - ARF2e*h*UinF
595 # Correction paremeters for normal measurements; they are independent of LDR
596 if (not RotationErrorEpsilonForNormalMeasurements):
597 GT = ATF1*IinE + ATF4*VinE
598 GR = ARF1*IinE + ARF4*VinE
599 HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE
600 HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE
601 else:
602 GT = ATF1*IinE + ATF4*h*VinE
603 GR = ARF1*IinE + ARF4*h*VinE
604 HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE
605 HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE
606 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
607 # p = +45°, m = -45°
608 IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF])
609 IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe])
610 AT = np.dot(ATFe,IF1e)
611 AR = np.dot(ARFe,IF1e)
612 BT = np.dot(ATFe,IF2e)
613 BR = np.dot(ARFe,IF2e)
614
615 # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6
616 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
617 IS1 = np.array([IinE,0,0,VinE])
618 IS2 = np.array([0,QinE,-UinE,-2*VinE])
619 GT = np.dot(ATF,IS1)
620 GR = np.dot(ARF,IS1)
621 HT = np.dot(ATF,IS2)
622 HR = np.dot(ARF,IS2)
623 else:
624 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
625 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
626 GT = np.dot(ATFe,IS1e)
627 GR = np.dot(ARFe,IS1e)
628 HT = np.dot(ATFe,IS2e)
629 HR = np.dot(ARFe,IS2e)
630
631 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
632 # parameters for calibration with aCal
633 IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF])
634 IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe])
635 AT = np.dot(ATFe,IF1e)
636 AR = np.dot(ARFe,IF1e)
637 BT = np.dot(ATFe,IF2e)
638 BR = np.dot(ARFe,IF2e)
639
640 # Correction paremeters for normal measurements; they are independent of LDR
641 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
642 #IS1 = np.array([IinE,0,0,VinE])
643 #IS2 = np.array([0,QinE,-UinE,-2*VinE])
644 IS1 = np.array([IinFo,0,0,VinFo])
645 IS2 = np.array([0,QinFa,UinFa,VinFa])
646 GT = np.dot(ATF,IS1)
647 GR = np.dot(ARF,IS1)
648 HT = np.dot(ATF,IS2)
649 HR = np.dot(ARF,IS2)
650 else:
651 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
652 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
653 #IS1e = np.array([IinFo,0,0,VinFo])
654 #IS2e = np.array([0,QinFae,UinFae,VinFa])
655 GT = np.dot(ATFe,IS1e)
656 GR = np.dot(ARFe,IS1e)
657 HT = np.dot(ATFe,IS2e)
658 HR = np.dot(ARFe,IS2e)
659
660 else:
661 print('Calibrator not implemented yet')
662 sys.exit()
663
664 elif LocC == 2: # C behind emitter optics Eq.57 -------------------------------------------------------
665 #print("Calibrator location not implemented yet")
666 S2e = np.sin(np.deg2rad(2*RotC))
667 C2e = np.cos(np.deg2rad(2*RotC))
668
669 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt)
670 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
671 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
672 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO])
673 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
674
675 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
676 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
677 ATF1 = ATF[0]
678 ATF2 = ATF[1]
679 ATF3 = ATF[2]
680 ATF4 = ATF[3]
681 ARF1 = ARF[0]
682 ARF2 = ARF[1]
683 ARF3 = ARF[2]
684 ARF4 = ARF[3]
685
686 # AS with C behind the emitter
687 # terms without aCal
688 ATE1o, ARE1o = ATF1, ARF1
689 ATE2o, ARE2o = 0., 0.
690 ATE3o, ARE3o = 0., 0.
691 ATE4o, ARE4o = ATF4, ARF4
692 # terms with aCal
693 ATE1a, ARE1a = 0. , 0.
694 ATE2a, ARE2a = ATF2, ARF2
695 ATE3a, ARE3a = -ATF3, -ARF3
696 ATE4a, ARE4a = -2*ATF4, -2*ARF4
697 # rotated AinEa by epsilon
698 ATE2ae = C2e*ATF2 + S2e*ATF3
699 ATE3ae = -S2e*ATF2 - C2e*ATF3
700 ARE2ae = C2e*ARF2 + S2e*ARF3
701 ARE3ae = -S2e*ARF2 - C2e*ARF3
702
703 ATE1 = ATE1o
704 ATE2e = aCal*ATE2ae
705 ATE3e = aCal*ATE3ae
706 ATE4 = (1-2*aCal)*ATF4
707 ARE1 = ARE1o
708 ARE2e = aCal*ARE2ae
709 ARE3e = aCal*ARE3ae
710 ARE4 = (1-2*aCal)*ARF4
711
712 # rotated IinE
713 QinEe = C2e*QinE + S2e*UinE
714 UinEe = C2e*UinE - S2e*QinE
715
716 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
717 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4
718 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE
719 BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe)
720 AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE
721 BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe)
722
723 # Correction paremeters for normal measurements; they are independent of LDR
724 if (not RotationErrorEpsilonForNormalMeasurements):
725 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
726 GT = ATE1o*IinE + ATE4o*h*VinE
727 GR = ARE1o*IinE + ARE4o*h*VinE
728 HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE
729 HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE
730 else:
731 GT = ATE1o*IinE + ATE4o*h*VinE
732 GR = ARE1o*IinE + ARE4o*h*VinE
733 HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE
734 HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE
735
736 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5
737 # p = +45°, m = -45°
738 AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe
739 BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE)
740 AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe
741 BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE)
742
743 # Correction paremeters for normal measurements; they are independent of LDR
744 if (not RotationErrorEpsilonForNormalMeasurements):
745 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
746 GT = ATE1o*IinE + ATE4o*VinE
747 GR = ARE1o*IinE + ARE4o*VinE
748 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
749 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
750 else:
751 D = IinE + DiC*QinEe
752 A = DiC*IinE + QinEe
753 B = ZiC*(CosC*UinEe + SinC*VinE)
754 C = -ZiC*(SinC*UinEe - CosC*VinE)
755 GT = ATE1o*D + ATE4o*C
756 GR = ARE1o*D + ARE4o*C
757 HT = ATE2a*A + ATE3a*B + ATE4a*C
758 HR = ARE2a*A + ARE3a*B + ARE4a*C
759
760 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt
761 # p = +22.5°, m = -22.5°
762 IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE])
763 IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe])
764 ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4])
765 AREe = np.array([ARE1,ARE2e,ARE3e,ARE4])
766 AT = np.dot(ATEe,IE1e)
767 AR = np.dot(AREe,IE1e)
768 BT = np.dot(ATEe,IE2e)
769 BR = np.dot(AREe,IE2e)
770
771 # Correction paremeters for normal measurements; they are independent of LDR
772 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
773 GT = ATE1o*IinE + ATE4o*VinE
774 GR = ARE1o*IinE + ARE4o*VinE
775 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
776 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
777 else:
778 D = IinE + DiC*QinEe
779 A = DiC*IinE + QinEe
780 B = ZiC*(CosC*UinEe + SinC*VinE)
781 C = -ZiC*(SinC*UinEe - CosC*VinE)
782 GT = ATE1o*D + ATE4o*C
783 GR = ARE1o*D + ARE4o*C
784 HT = ATE2a*A + ATE3a*B + ATE4a*C
785 HR = ARE2a*A + ARE3a*B + ARE4a*C
786
787 else:
788 print('Calibrator not implemented yet')
789 sys.exit()
790
791 else:
792 print("Calibrator location not implemented yet")
793 sys.exit()
794
795 # Determination of the correction K of the calibration factor
796 IoutTp = TaT*TiT*TiO*TiE*(AT + BT)
797 IoutTm = TaT*TiT*TiO*TiE*(AT - BT)
798 IoutRp = TaR*TiR*TiO*TiE*(AR + BR)
799 IoutRm = TaR*TiR*TiO*TiE*(AR - BR)
800
801 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
802 Etapx = IoutRp/IoutTp
803 Etamx = IoutRm/IoutTm
804 Etax = (Etapx*Etamx)**0.5
805
806 Eta = (TaR*TiR)/(TaT*TiT) # Eta = Eta*/K Eq. 84
807 K = Etax / Eta
808
809 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet
810 #Eta_test_p = (IoutRp/IoutTp)
811 #Eta_test_m = (IoutRm/IoutTm)
812 #Eta_test = (Eta_test_p*Eta_test_m)**0.5
813
814 # ----- Forward simulated signals and LDRsim with atrue; from input file
815 It = TaT*TiT*TiO*TiE*(GT+atrue*HT)
816 Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR)
817 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file
818 LDRsim = Ir/It # simulated uncorrected LDR with Y from input file
819 # Corrected LDRsimCorr from forward simulated LDRsim (atrue)
820 # LDRsimCorr = (1./Eta*LDRsim*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRsim*(GT-HT))
821 if Y == -1.:
822 LDRsimx = 1./LDRsim
823 else:
824 LDRsimx = LDRsim
825
826 # The following is correct without doubt
827 #LDRCorr = (LDRsim*K/Etax*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
828
829 # The following is a test whether the equations for calibration Etax and normal signal (GHK, LDRsim) are consistent
830 LDRCorr = (LDRsim/Eta*(GT+HT)-(GR+HR))/((GR-HR)-LDRsim*K/Etax*(GT-HT))
831
832 TTa = TiT*TaT #*ATP1
833 TRa = TiR*TaR #*ARP1
834
835 F11sim = 1/(TiO*TiE)*((HR*Etax/K*It/TTa-HT*Ir/TRa)/(HR*GT-HT*GR)) # IL = 1, Etat = Etar = 1
836
837 return (GT, HT, GR, HR, K, Eta, LDRsimx, LDRCorr, DTa, DRa, TTa, TRa, F11sim)
838 # *******************************************************************************************************************************
839
840 # --- CALC truth
841 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
842
843 # --------------------------------------------------------
844 with open('output_' + LID + '.dat', 'w') as f:
845 with redirect_stdout(f):
846 print("From ", dname)
847 print("Running ", fname)
848 print("Reading input file ", InputFile) #, " for Lidar system :", EID, ", ", LID)
849 print("for Lidar system: ", EID, ", ", LID)
850 # --- Print iput information*********************************
851 print(" --- Input parameters: value ±error / ±steps ----------------------")
852 print("{0:8} {1:8} {2:8.5f}; {3:8} {4:7.4f}±{5:7.4f}/{6:2d}".format("Laser: ", "DOLP = ", bL, " rotation alpha = ", RotL0, dRotL, nRotL))
853 print(" Diatt., Tunpol, Retard., Rotation (deg)")
854 print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Emitter ", DiE0, dDiE, TiE, RetE0, dRetE, RotE0, dRotE, nDiE, nRetE, nRotE))
855 print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Receiver ", DiO0, dDiO, TiO, RetO0, dRetO, RotO0, dRotO, nDiO, nRetO, nRotO))
856 print("{0:12} {1:7.4f}±{2:7.4f}/{8:2d}, {3:7.4f}, {4:3.0f}±{5:3.0f}/{9:2d}, {6:7.4f}±{7:7.4f}/{10:2d}".format("Calibrator ", DiC0, dDiC, TiC, RetC0, dRetC, RotC0, dRotC, nDiC, nRetC, nRotC))
857 print("{0:12}".format(" --- Pol.-filter ---"))
858 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("ERT, ERR :", ERaT0, dERaT, nERaT, ERaR0, dERaR, nERaR))
859 print("{0:12}{1:7.4f}±{2:7.4f}/{3:2d}, {4:7.4f}±{5:7.4f}/{6:2d}".format("RotaT , RotaR :", RotaT0, dRotaT, nRotaT, RotaR0,dRotaR,nRotaR))
860 print("{0:12}".format(" --- PBS ---"))
861 print("{0:12}{1:7.4f}±{2:7.4f}/{9:2d}, {3:7.4f}±{4:7.4f}/{10:2d}, {5:7.4f}±{6:7.4f}/{11:2d},{7:7.4f}±{8:7.4f}/{12:2d}".format("TP,TS,RP,RS :", TP0, dTP, TS0, dTS, RP0, dRP, RS0, dRS, nTP, nTS, nRP, nRS))
862 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}, {5:1.0f}".format("DT,TT,DR,TR,Y :", DiT, TiT, DiR, TiR, Y))
863 print("{0:12}".format(" --- Combined PBS + Pol.-filter ---"))
864 print("{0:12}{1:7.4f},{2:7.4f}, {3:7.4f},{4:7.4f}".format("DTa,TTa,DRa,TRa: ", DTa0, TTa0, DRa0, TRa0))
865 print()
866 print("Rotation Error Epsilon For Normal Measurements = ", RotationErrorEpsilonForNormalMeasurements)
867 #print ('LocC = ', LocC, Loc[LocC], '; TypeC = ',TypeC, Type[TypeC])
868 print(Type[TypeC], Loc[LocC], "; Parallel signal detected in", dY[int(Y+1)])
869 # end of print actual system parameters
870 # ******************************************************************************
871
872 #print()
873 #print(" --- LDRCal during calibration | simulated and corrected LDRs -------------")
874 #print("{0:8} |{1:8}->{2:8},{3:9}->{4:9} |{5:8}->{6:8}".format(" LDRCal"," LDRtrue", " LDRsim"," LDRtrue2", " LDRsim2", " LDRmeas", " LDRcorr"))
875 #print("{0:8.5f} |{1:8.5f}->{2:8.5f},{3:9.5f}->{4:9.5f} |{5:8.5f}->{6:8.5f}".format(LDRCal, LDRtrue, LDRsim, LDRtrue2, LDRsim2, LDRmeas, LDRCorr))
876 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
877 #print("{0:6.3f}±{1:5.3f}/{2:2d}|{3:8.5f}->{4:8.5f}->{5:8.5f}".format(LDRCal0, dLDRCal, nLDRCal, LDRtrue, LDRsimx, LDRCorr))
878 #print("{0:8} |{1:8}->{2:8}->{3:8}".format(" LDRCal"," LDRtrue", " LDRsimx", " LDRcorr"))
879 #print(" --- LDRCal during calibration")
880 print("{0:26}: {1:6.3f}±{2:5.3f}/{3:2d}".format("LDRCal during calibration", LDRCal0, dLDRCal, nLDRCal))
881
882 #print("{0:8}={1:8.5f};{2:8}={3:8.5f}".format(" IinP",IinP," F11sim",F11sim))
883 print()
884
885 K0List = np.zeros(3)
886 LDRsimxList = np.zeros(3)
887 LDRCalList = 0.004, 0.2, 0.45
888 for i,LDRCal in enumerate(LDRCalList):
889 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
890 K0List[i] = K0
891 LDRsimxList[i] = LDRsimx
892
893 print("{0:8},{1:8},{2:8},{3:8},{4:9},{5:9},{6:9}".format(" GR", " GT", " HR", " HT", " K(0.004)", " K(0.2)", " K(0.45)"))
894 print("{0:8.5f},{1:8.5f},{2:8.5f},{3:8.5f},{4:9.5f},{5:9.5f},{6:9.5f}".format(GR0, GT0, HR0, HT0, K0List[0], K0List[1], K0List[2]))
895 print('========================================================================')
896
897 print("{0:9},{1:9},{2:9}".format(" LDRtrue", " LDRsimx", " LDRCorr"))
898 LDRtrueList = 0.004, 0.02, 0.2, 0.45
899 for i,LDRtrue in enumerate(LDRtrueList):
900 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
901 print("{0:9.5f},{1:9.5f},{2:9.5f}".format(LDRtrue, LDRsimx, LDRCorr))
902
903
904 file = open('output_' + LID + '.dat', 'r')
905 print (file.read())
906 file.close()
907
908 '''
909 if(PrintToOutputFile):
910 f = open('output_ver7.dat', 'w')
911 old_target = sys.stdout
912 sys.stdout = f
913
914 print("something")
915
916 if(PrintToOutputFile):
917 sys.stdout.flush()
918 f.close
919 sys.stdout = old_target
920 '''
921 # --- CALC again truth with LDRCal0 to reset all 0-values
922 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal0)
923
924 # --- Start Errors calculation
925
926 iN = -1
927 N = ((nRotL*2+1)*
928 (nRotE*2+1)*(nRetE*2+1)*(nDiE*2+1)*
929 (nRotO*2+1)*(nRetO*2+1)*(nDiO*2+1)*
930 (nRotC*2+1)*(nRetC*2+1)*(nDiC*2+1)*
931 (nTP*2+1)*(nTS*2+1)*(nRP*2+1)*(nRS*2+1)*(nERaT*2+1)*(nERaR*2+1)*
932 (nRotaT*2+1)*(nRotaR*2+1)*(nRetT*2+1)*(nRetR*2+1)*(nLDRCal*2+1))
933 print("N = ",N ," ", end="")
934
935 if N > 1e6:
936 if user_yes_no_query('Warning: processing ' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
937 if N > 5e6:
938 if user_yes_no_query('Warning: the memory required for ' + str(N) + ' samples might be ' + '{0:5.1f}'.format(N/4e6) + ' GB. Do you anyway want to proceed?') == 0: sys.exit()
939
940 #if user_yes_no_query('Warning: processing' + str(N) + ' samples will take very long. Do you want to proceed?') == 0: sys.exit()
941
942 # --- Arrays for plotting ------
943 LDRmin = np.zeros(5)
944 LDRmax = np.zeros(5)
945 F11min = np.zeros(5)
946 F11max = np.zeros(5)
947
948 LDRrange = np.zeros(5)
949 LDRrange = 0.004, 0.02, 0.1, 0.3, 0.45
950 #aLDRsimx = np.zeros(N)
951 #aLDRsimx2 = np.zeros(N)
952 #aLDRcorr = np.zeros(N)
953 #aLDRcorr2 = np.zeros(N)
954 aERaT = np.zeros(N)
955 aERaR = np.zeros(N)
956 aRotaT = np.zeros(N)
957 aRotaR = np.zeros(N)
958 aRetT = np.zeros(N)
959 aRetR = np.zeros(N)
960 aTP = np.zeros(N)
961 aTS = np.zeros(N)
962 aRP = np.zeros(N)
963 aRS = np.zeros(N)
964 aDiE = np.zeros(N)
965 aDiO = np.zeros(N)
966 aDiC = np.zeros(N)
967 aRotC = np.zeros(N)
968 aRetC = np.zeros(N)
969 aRotL = np.zeros(N)
970 aRetE = np.zeros(N)
971 aRotE = np.zeros(N)
972 aRetO = np.zeros(N)
973 aRotO = np.zeros(N)
974 aLDRCal = np.zeros(N)
975 aA = np.zeros((5,N))
976 aX = np.zeros((5,N))
977 aF11corr = np.zeros((5,N))
978
979 atime = clock()
980 dtime = clock()
981
982 # --- Calc Error signals
983 #GT, HT, GR, HR, K, Eta, LDRsim = Calc(RotL, RotE, RetE, DiE, RotO, RetO, DiO, RotC, RetC, DiC, TP, TS)
984 # ---- Do the calculations of bra-ket vectors
985 h = -1. if TypeC == 2 else 1
986
987 # from input file: measured LDRm and true LDRtrue, LDRtrue2 =>
988 ameas = (1.-LDRmeas)/(1+LDRmeas)
989 atrue = (1.-LDRtrue)/(1+LDRtrue)
990 atrue2 = (1.-LDRtrue2)/(1+LDRtrue2)
991
992 for iLDRCal in range(-nLDRCal,nLDRCal+1):
993 # from input file: assumed LDRCal for calibration measurements
994 LDRCal = LDRCal0
995 if nLDRCal > 0: LDRCal = LDRCal0 + iLDRCal*dLDRCal/nLDRCal
996
997 GT0, HT0, GR0, HR0, K0, Eta0, LDRsimx, LDRCorr, DTa0, DRa0, TTa0, TRa0, F11sim0 = Calc(RotL0, RotE0, RetE0, DiE0, RotO0, RetO0, DiO0, RotC0, RetC0, DiC0, TP0, TS0, RP0, RS0, ERaT0, RotaT0, RetT0, ERaR0, RotaR0, RetR0, LDRCal)
998 aCal = (1.-LDRCal)/(1+LDRCal)
999 for iRotL, iRotE, iRetE, iDiE \
1000 in [(iRotL,iRotE,iRetE,iDiE)
1001 for iRotL in range(-nRotL,nRotL+1)
1002 for iRotE in range(-nRotE,nRotE+1)
1003 for iRetE in range(-nRetE,nRetE+1)
1004 for iDiE in range(-nDiE,nDiE+1)]:
1005
1006 if nRotL > 0: RotL = RotL0 + iRotL*dRotL/nRotL
1007 if nRotE > 0: RotE = RotE0 + iRotE*dRotE/nRotE
1008 if nRetE > 0: RetE = RetE0 + iRetE*dRetE/nRetE
1009 if nDiE > 0: DiE = DiE0 + iDiE*dDiE/nDiE
1010
1011 # angles of emitter and laser and calibrator and receiver optics
1012 # RotL = alpha, RotE = beta, RotO = gamma, RotC = epsilon
1013 S2a = np.sin(2*np.deg2rad(RotL))
1014 C2a = np.cos(2*np.deg2rad(RotL))
1015 S2b = np.sin(2*np.deg2rad(RotE))
1016 C2b = np.cos(2*np.deg2rad(RotE))
1017 S2ab = np.sin(np.deg2rad(2*RotL-2*RotE))
1018 C2ab = np.cos(np.deg2rad(2*RotL-2*RotE))
1019
1020 # Laser with Degree of linear polarization DOLP = bL
1021 IinL = 1.
1022 QinL = bL
1023 UinL = 0.
1024 VinL = (1. - bL**2)**0.5
1025
1026 # Stokes Input Vector rotation Eq. E.4
1027 A = C2a*QinL - S2a*UinL
1028 B = S2a*QinL + C2a*UinL
1029 # Stokes Input Vector rotation Eq. E.9
1030 C = C2ab*QinL - S2ab*UinL
1031 D = S2ab*QinL + C2ab*UinL
1032
1033 # emitter optics
1034 CosE = np.cos(np.deg2rad(RetE))
1035 SinE = np.sin(np.deg2rad(RetE))
1036 ZiE = (1. - DiE**2)**0.5
1037 WiE = (1. - ZiE*CosE)
1038
1039 # Stokes Input Vector after emitter optics equivalent to Eq. E.9 with already rotated input vector from Eq. E.4
1040 # b = beta
1041 IinE = (IinL + DiE*C)
1042 QinE = (C2b*DiE*IinL + A + S2b*(WiE*D - ZiE*SinE*VinL))
1043 UinE = (S2b*DiE*IinL + B - C2b*(WiE*D - ZiE*SinE*VinL))
1044 VinE = (-ZiE*SinE*D + ZiE*CosE*VinL)
1045
1046 #-------------------------
1047 # F11 assuemd to be = 1 => measured: F11m = IinP / IinE with atrue
1048 #F11sim = (IinE + DiO*atrue*(C2g*QinE - S2g*UinE))/IinE
1049 #-------------------------
1050
1051 for iRotO, iRetO, iDiO, iRotC, iRetC, iDiC, iTP, iTS, iRP, iRS, iERaT, iRotaT, iRetT, iERaR, iRotaR, iRetR \
1052 in [(iRotO,iRetO,iDiO,iRotC,iRetC,iDiC,iTP,iTS,iRP,iRS,iERaT,iRotaT,iRetT,iERaR,iRotaR,iRetR )
1053 for iRotO in range(-nRotO,nRotO+1)
1054 for iRetO in range(-nRetO,nRetO+1)
1055 for iDiO in range(-nDiO,nDiO+1)
1056 for iRotC in range(-nRotC,nRotC+1)
1057 for iRetC in range(-nRetC,nRetC+1)
1058 for iDiC in range(-nDiC,nDiC+1)
1059 for iTP in range(-nTP,nTP+1)
1060 for iTS in range(-nTS,nTS+1)
1061 for iRP in range(-nRP,nRP+1)
1062 for iRS in range(-nRS,nRS+1)
1063 for iERaT in range(-nERaT,nERaT+1)
1064 for iRotaT in range(-nRotaT,nRotaT+1)
1065 for iRetT in range(-nRetT,nRetT+1)
1066 for iERaR in range(-nERaR,nERaR+1)
1067 for iRotaR in range(-nRotaR,nRotaR+1)
1068 for iRetR in range(-nRetR,nRetR+1)]:
1069
1070 iN = iN + 1
1071 if (iN == 10001):
1072 ctime = clock()
1073 print(" estimated time ", "{0:4.2f}".format(N/10000 * (ctime-atime)), "sec ") #, end="")
1074 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r")
1075 ctime = clock()
1076 if ((ctime - dtime) > 10):
1077 print("\r elapsed time ", "{0:5.0f}".format((ctime-atime)), "sec ", end="\r")
1078 dtime = ctime
1079
1080 if nRotO > 0: RotO = RotO0 + iRotO*dRotO/nRotO
1081 if nRetO > 0: RetO = RetO0 + iRetO*dRetO/nRetO
1082 if nDiO > 0: DiO = DiO0 + iDiO*dDiO/nDiO
1083 if nRotC > 0: RotC = RotC0 + iRotC*dRotC/nRotC
1084 if nRetC > 0: RetC = RetC0 + iRetC*dRetC/nRetC
1085 if nDiC > 0: DiC = DiC0 + iDiC*dDiC/nDiC
1086 if nTP > 0: TP = TP0 + iTP*dTP/nTP
1087 if nTS > 0: TS = TS0 + iTS*dTS/nTS
1088 if nRP > 0: RP = RP0 + iRP*dRP/nRP
1089 if nRS > 0: RS = RS0 + iRS*dRS/nRS
1090 if nERaT > 0: ERaT = ERaT0 + iERaT*dERaT/nERaT
1091 if nRotaT > 0:RotaT= RotaT0+ iRotaT*dRotaT/nRotaT
1092 if nRetT > 0: RetT = RetT0 + iRetT*dRetT/nRetT
1093 if nERaR > 0: ERaR = ERaR0 + iERaR*dERaR/nERaR
1094 if nRotaR > 0:RotaR= RotaR0+ iRotaR*dRotaR/nRotaR
1095 if nRetR > 0: RetR = RetR0 + iRetR*dRetR/nRetR
1096
1097 #print("{0:5.2f}, {1:5.2f}, {2:5.2f}, {3:10d}".format(RotL, RotE, RotO, iN))
1098
1099 # receiver optics
1100 CosO = np.cos(np.deg2rad(RetO))
1101 SinO = np.sin(np.deg2rad(RetO))
1102 ZiO = (1. - DiO**2)**0.5
1103 WiO = (1. - ZiO*CosO)
1104 S2g = np.sin(np.deg2rad(2*RotO))
1105 C2g = np.cos(np.deg2rad(2*RotO))
1106 # calibrator
1107 CosC = np.cos(np.deg2rad(RetC))
1108 SinC = np.sin(np.deg2rad(RetC))
1109 ZiC = (1. - DiC**2)**0.5
1110 WiC = (1. - ZiC*CosC)
1111
1112 #For POLLY_XT
1113 # analyser
1114 #RS = 1 - TS
1115 #RP = 1 - TP
1116 TiT = 0.5 * (TP + TS)
1117 DiT = (TP-TS)/(TP+TS)
1118 ZiT = (1. - DiT**2)**0.5
1119 TiR = 0.5 * (RP + RS)
1120 DiR = (RP-RS)/(RP+RS)
1121 ZiR = (1. - DiR**2)**0.5
1122 CosT = np.cos(np.deg2rad(RetT))
1123 SinT = np.sin(np.deg2rad(RetT))
1124 CosR = np.cos(np.deg2rad(RetR))
1125 SinR = np.sin(np.deg2rad(RetR))
1126
1127 DaT = (1-ERaT)/(1+ERaT)
1128 DaR = (1-ERaR)/(1+ERaR)
1129 TaT = 0.5*(1+ERaT)
1130 TaR = 0.5*(1+ERaR)
1131
1132 S2aT = np.sin(np.deg2rad(h*2*RotaT))
1133 C2aT = np.cos(np.deg2rad(2*RotaT))
1134 S2aR = np.sin(np.deg2rad(h*2*RotaR))
1135 C2aR = np.cos(np.deg2rad(2*RotaR))
1136
1137 # Aanalyzer As before the PBS Eq. D.5
1138 ATP1 = (1+C2aT*DaT*DiT)
1139 ATP2 = Y*(DiT+C2aT*DaT)
1140 ATP3 = Y*S2aT*DaT*ZiT*CosT
1141 ATP4 = S2aT*DaT*ZiT*SinT
1142 ATP = np.array([ATP1,ATP2,ATP3,ATP4])
1143
1144 ARP1 = (1+C2aR*DaR*DiR)
1145 ARP2 = Y*(DiR+C2aR*DaR)
1146 ARP3 = Y*S2aR*DaR*ZiR*CosR
1147 ARP4 = S2aR*DaR*ZiR*SinR
1148 ARP = np.array([ARP1,ARP2,ARP3,ARP4])
1149
1150 TTa = TiT*TaT #*ATP1
1151 TRa = TiR*TaR #*ARP1
1152
1153 # ---- Calculate signals and correction parameters for diffeent locations and calibrators
1154 if LocC == 4: # Calibrator before the PBS
1155 #print("Calibrator location not implemented yet")
1156
1157 #S2ge = np.sin(np.deg2rad(2*RotO + h*2*RotC))
1158 #C2ge = np.cos(np.deg2rad(2*RotO + h*2*RotC))
1159 S2e = np.sin(np.deg2rad(h*2*RotC))
1160 C2e = np.cos(np.deg2rad(2*RotC))
1161 # rotated AinP by epsilon Eq. C.3
1162 ATP2e = C2e*ATP2 + S2e*ATP3
1163 ATP3e = C2e*ATP3 - S2e*ATP2
1164 ARP2e = C2e*ARP2 + S2e*ARP3
1165 ARP3e = C2e*ARP3 - S2e*ARP2
1166 ATPe = np.array([ATP1,ATP2e,ATP3e,ATP4])
1167 ARPe = np.array([ARP1,ARP2e,ARP3e,ARP4])
1168 # Stokes Input Vector before the polarising beam splitter Eq. E.31
1169 A = C2g*QinE - S2g*UinE
1170 B = S2g*QinE + C2g*UinE
1171 #C = (WiO*aCal*B + ZiO*SinO*(1-2*aCal)*VinE)
1172 Co = ZiO*SinO*VinE
1173 Ca = (WiO*B - 2*ZiO*SinO*VinE)
1174 #C = Co + aCal*Ca
1175 #IinP = (IinE + DiO*aCal*A)
1176 #QinP = (C2g*DiO*IinE + aCal*QinE - S2g*C)
1177 #UinP = (S2g*DiO*IinE - aCal*UinE + C2g*C)
1178 #VinP = (ZiO*SinO*aCal*B + ZiO*CosO*(1-2*aCal)*VinE)
1179 IinPo = IinE
1180 QinPo = (C2g*DiO*IinE - S2g*Co)
1181 UinPo = (S2g*DiO*IinE + C2g*Co)
1182 VinPo = ZiO*CosO*VinE
1183
1184 IinPa = DiO*A
1185 QinPa = QinE - S2g*Ca
1186 UinPa = -UinE + C2g*Ca
1187 VinPa = ZiO*(SinO*B - 2*CosO*VinE)
1188
1189 IinP = IinPo + aCal*IinPa
1190 QinP = QinPo + aCal*QinPa
1191 UinP = UinPo + aCal*UinPa
1192 VinP = VinPo + aCal*VinPa
1193 # Stokes Input Vector before the polarising beam splitter rotated by epsilon Eq. C.3
1194 #QinPe = C2e*QinP + S2e*UinP
1195 #UinPe = C2e*UinP - S2e*QinP
1196 QinPoe = C2e*QinPo + S2e*UinPo
1197 UinPoe = C2e*UinPo - S2e*QinPo
1198 QinPae = C2e*QinPa + S2e*UinPa
1199 UinPae = C2e*UinPa - S2e*QinPa
1200 QinPe = C2e*QinP + S2e*UinP
1201 UinPe = C2e*UinP - S2e*QinP
1202
1203 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
1204 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
1205 # parameters for calibration with aCal
1206 AT = ATP1*IinP + h*ATP4*VinP
1207 BT = ATP3e*QinP - h*ATP2e*UinP
1208 AR = ARP1*IinP + h*ARP4*VinP
1209 BR = ARP3e*QinP - h*ARP2e*UinP
1210 # Correction paremeters for normal measurements; they are independent of LDR
1211 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
1212 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
1213 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
1214 GT = np.dot(ATP,IS1)
1215 GR = np.dot(ARP,IS1)
1216 HT = np.dot(ATP,IS2)
1217 HR = np.dot(ARP,IS2)
1218 else:
1219 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
1220 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
1221 GT = np.dot(ATPe,IS1)
1222 GR = np.dot(ARPe,IS1)
1223 HT = np.dot(ATPe,IS2)
1224 HR = np.dot(ARPe,IS2)
1225 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
1226 # parameters for calibration with aCal
1227 AT = ATP1*IinP + ATP3e*UinPe + ZiC*CosC*(ATP2e*QinPe + ATP4*VinP)
1228 BT = DiC*(ATP1*UinPe + ATP3e*IinP) - ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
1229 AR = ARP1*IinP + ARP3e*UinPe + ZiC*CosC*(ARP2e*QinPe + ARP4*VinP)
1230 BR = DiC*(ARP1*UinPe + ARP3e*IinP) - ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
1231 # Correction paremeters for normal measurements; they are independent of LDR
1232 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
1233 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
1234 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
1235 GT = np.dot(ATP,IS1)
1236 GR = np.dot(ARP,IS1)
1237 HT = np.dot(ATP,IS2)
1238 HR = np.dot(ARP,IS2)
1239 else:
1240 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
1241 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
1242 GT = np.dot(ATPe,IS1e)
1243 GR = np.dot(ARPe,IS1e)
1244 HT = np.dot(ATPe,IS2e)
1245 HR = np.dot(ARPe,IS2e)
1246 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
1247 # parameters for calibration with aCal
1248 AT = ATP1*IinP + sqr05*DiC*(ATP1*QinPe + ATP2e*IinP) + (1-0.5*WiC)*(ATP2e*QinPe + ATP3e*UinPe) + ZiC*(sqr05*SinC*(ATP3e*VinP-ATP4*UinPe) + ATP4*CosC*VinP)
1249 BT = sqr05*DiC*(ATP1*UinPe + ATP3e*IinP) + 0.5*WiC*(ATP2e*UinPe + ATP3e*QinPe) - sqr05*ZiC*SinC*(ATP2e*VinP - ATP4*QinPe)
1250 AR = ARP1*IinP + sqr05*DiC*(ARP1*QinPe + ARP2e*IinP) + (1-0.5*WiC)*(ARP2e*QinPe + ARP3e*UinPe) + ZiC*(sqr05*SinC*(ARP3e*VinP-ARP4*UinPe) + ARP4*CosC*VinP)
1251 BR = sqr05*DiC*(ARP1*UinPe + ARP3e*IinP) + 0.5*WiC*(ARP2e*UinPe + ARP3e*QinPe) - sqr05*ZiC*SinC*(ARP2e*VinP - ARP4*QinPe)
1252 # Correction paremeters for normal measurements; they are independent of LDR
1253 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
1254 IS1 = np.array([IinPo,QinPo,UinPo,VinPo])
1255 IS2 = np.array([IinPa,QinPa,UinPa,VinPa])
1256 GT = np.dot(ATP,IS1)
1257 GR = np.dot(ARP,IS1)
1258 HT = np.dot(ATP,IS2)
1259 HR = np.dot(ARP,IS2)
1260 else:
1261 IS1e = np.array([IinPo+DiC*QinPoe,DiC*IinPo+QinPoe,ZiC*(CosC*UinPoe+SinC*VinPo),-ZiC*(SinC*UinPoe-CosC*VinPo)])
1262 IS2e = np.array([IinPa+DiC*QinPae,DiC*IinPa+QinPae,ZiC*(CosC*UinPae+SinC*VinPa),-ZiC*(SinC*UinPae-CosC*VinPa)])
1263 GT = np.dot(ATPe,IS1e)
1264 GR = np.dot(ARPe,IS1e)
1265 HT = np.dot(ATPe,IS2e)
1266 HR = np.dot(ARPe,IS2e)
1267 else:
1268 print("Calibrator not implemented yet")
1269 sys.exit()
1270
1271 elif LocC == 3: # C before receiver optics Eq.57
1272
1273 #S2ge = np.sin(np.deg2rad(2*RotO - 2*RotC))
1274 #C2ge = np.cos(np.deg2rad(2*RotO - 2*RotC))
1275 S2e = np.sin(np.deg2rad(2*RotC))
1276 C2e = np.cos(np.deg2rad(2*RotC))
1277
1278 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt)
1279 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
1280 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
1281 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO])
1282 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
1283
1284 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
1285 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
1286 ATF1 = ATF[0]
1287 ATF2 = ATF[1]
1288 ATF3 = ATF[2]
1289 ATF4 = ATF[3]
1290 ARF1 = ARF[0]
1291 ARF2 = ARF[1]
1292 ARF3 = ARF[2]
1293 ARF4 = ARF[3]
1294
1295 # rotated AinF by epsilon
1296 ATF2e = C2e*ATF[1] + S2e*ATF[2]
1297 ATF3e = C2e*ATF[2] - S2e*ATF[1]
1298 ARF2e = C2e*ARF[1] + S2e*ARF[2]
1299 ARF3e = C2e*ARF[2] - S2e*ARF[1]
1300
1301 ATFe = np.array([ATF1,ATF2e,ATF3e,ATF4])
1302 ARFe = np.array([ARF1,ARF2e,ARF3e,ARF4])
1303
1304 QinEe = C2e*QinE + S2e*UinE
1305 UinEe = C2e*UinE - S2e*QinE
1306
1307 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
1308 IinF = IinE
1309 QinF = aCal*QinE
1310 UinF = -aCal*UinE
1311 VinF = (1.-2.*aCal)*VinE
1312
1313 IinFo = IinE
1314 QinFo = 0.
1315 UinFo = 0.
1316 VinFo = VinE
1317
1318 IinFa = 0.
1319 QinFa = QinE
1320 UinFa = -UinE
1321 VinFa = -2.*VinE
1322
1323 # Stokes Input Vector before receiver optics rotated by epsilon Eq. C.3
1324 QinFe = C2e*QinF + S2e*UinF
1325 UinFe = C2e*UinF - S2e*QinF
1326 QinFoe = C2e*QinFo + S2e*UinFo
1327 UinFoe = C2e*UinFo - S2e*QinFo
1328 QinFae = C2e*QinFa + S2e*UinFa
1329 UinFae = C2e*UinFa - S2e*QinFa
1330
1331 # Calibration signals and Calibration correction K from measurements with LDRCal / aCal
1332 if (TypeC == 2) or (TypeC == 1): # rotator calibration Eq. C.4
1333 AT = ATF1*IinF + ATF4*h*VinF
1334 BT = ATF3e*QinF - ATF2e*h*UinF
1335 AR = ARF1*IinF + ARF4*h*VinF
1336 BR = ARF3e*QinF - ARF2e*h*UinF
1337
1338 # Correction paremeters for normal measurements; they are independent of LDR
1339 if (not RotationErrorEpsilonForNormalMeasurements):
1340 GT = ATF1*IinE + ATF4*VinE
1341 GR = ARF1*IinE + ARF4*VinE
1342 HT = ATF2*QinE - ATF3*UinE - ATF4*2*VinE
1343 HR = ARF2*QinE - ARF3*UinE - ARF4*2*VinE
1344 else:
1345 GT = ATF1*IinE + ATF4*h*VinE
1346 GR = ARF1*IinE + ARF4*h*VinE
1347 HT = ATF2e*QinE - ATF3e*h*UinE - ATF4*h*2*VinE
1348 HR = ARF2e*QinE - ARF3e*h*UinE - ARF4*h*2*VinE
1349
1350 elif (TypeC == 3) or (TypeC == 4): # linear polariser calibration Eq. C.5
1351 # p = +45°, m = -45°
1352 IF1e = np.array([IinF, ZiC*CosC*QinFe, UinFe, ZiC*CosC*VinF])
1353 IF2e = np.array([DiC*UinFe, -ZiC*SinC*VinF, DiC*IinF, ZiC*SinC*QinFe])
1354
1355 AT = np.dot(ATFe,IF1e)
1356 AR = np.dot(ARFe,IF1e)
1357 BT = np.dot(ATFe,IF2e)
1358 BR = np.dot(ARFe,IF2e)
1359
1360 # Correction paremeters for normal measurements; they are independent of LDR --- the same as for TypeC = 6
1361 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
1362 IS1 = np.array([IinE,0,0,VinE])
1363 IS2 = np.array([0,QinE,-UinE,-2*VinE])
1364
1365 GT = np.dot(ATF,IS1)
1366 GR = np.dot(ARF,IS1)
1367 HT = np.dot(ATF,IS2)
1368 HR = np.dot(ARF,IS2)
1369 else:
1370 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
1371 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
1372 GT = np.dot(ATFe,IS1e)
1373 GR = np.dot(ARFe,IS1e)
1374 HT = np.dot(ATFe,IS2e)
1375 HR = np.dot(ARFe,IS2e)
1376
1377 elif (TypeC == 6): # diattenuator calibration +-22.5° rotated_diattenuator_X22x5deg.odt
1378 # p = +22.5°, m = -22.5°
1379 IF1e = np.array([IinF+sqr05*DiC*QinFe, sqr05*DiC*IinF+(1-0.5*WiC)*QinFe, (1-0.5*WiC)*UinFe+sqr05*ZiC*SinC*VinF, -sqr05*ZiC*SinC*UinFe+ZiC*CosC*VinF])
1380 IF2e = np.array([sqr05*DiC*UinFe, 0.5*WiC*UinFe-sqr05*ZiC*SinC*VinF, sqr05*DiC*IinF+0.5*WiC*QinFe, sqr05*ZiC*SinC*QinFe])
1381
1382 AT = np.dot(ATFe,IF1e)
1383 AR = np.dot(ARFe,IF1e)
1384 BT = np.dot(ATFe,IF2e)
1385 BR = np.dot(ARFe,IF2e)
1386
1387 # Correction paremeters for normal measurements; they are independent of LDR
1388 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
1389 #IS1 = np.array([IinE,0,0,VinE])
1390 #IS2 = np.array([0,QinE,-UinE,-2*VinE])
1391 IS1 = np.array([IinFo,0,0,VinFo])
1392 IS2 = np.array([0,QinFa,UinFa,VinFa])
1393 GT = np.dot(ATF,IS1)
1394 GR = np.dot(ARF,IS1)
1395 HT = np.dot(ATF,IS2)
1396 HR = np.dot(ARF,IS2)
1397 else:
1398 #IS1e = np.array([IinE,DiC*IinE,ZiC*SinC*VinE,ZiC*CosC*VinE])
1399 #IS2e = np.array([DiC*QinEe,QinEe,-ZiC*(CosC*UinEe+2*SinC*VinE),ZiC*(SinC*UinEe-2*CosC*VinE)])
1400 IS1e = np.array([IinFo+DiC*QinFoe,DiC*IinFo+QinFoe,ZiC*(CosC*UinFoe+SinC*VinFo),-ZiC*(SinC*UinFoe-CosC*VinFo)])
1401 IS2e = np.array([IinFa+DiC*QinFae,DiC*IinFa+QinFae,ZiC*(CosC*UinFae+SinC*VinFa),-ZiC*(SinC*UinFae-CosC*VinFa)])
1402 GT = np.dot(ATFe,IS1e)
1403 GR = np.dot(ARFe,IS1e)
1404 HT = np.dot(ATFe,IS2e)
1405 HR = np.dot(ARFe,IS2e)
1406
1407
1408 else:
1409 print('Calibrator not implemented yet')
1410 sys.exit()
1411
1412 elif LocC == 2: # C behind emitter optics Eq.57
1413 #print("Calibrator location not implemented yet")
1414 S2e = np.sin(np.deg2rad(2*RotC))
1415 C2e = np.cos(np.deg2rad(2*RotC))
1416
1417 # AS with C before the receiver optics (see document rotated_diattenuator_X22x5deg.odt)
1418 AF1 = np.array([1,C2g*DiO,S2g*DiO,0])
1419 AF2 = np.array([C2g*DiO,1-S2g**2*WiO,S2g*C2g*WiO,-S2g*ZiO*SinO])
1420 AF3 = np.array([S2g*DiO, S2g*C2g*WiO, 1-C2g**2*WiO, C2g*ZiO*SinO])
1421 AF4 = np.array([0, S2g*SinO, -C2g*SinO, CosO])
1422
1423 ATF = (ATP1*AF1+ATP2*AF2+ATP3*AF3+ATP4*AF4)
1424 ARF = (ARP1*AF1+ARP2*AF2+ARP3*AF3+ARP4*AF4)
1425 ATF1 = ATF[0]
1426 ATF2 = ATF[1]
1427 ATF3 = ATF[2]
1428 ATF4 = ATF[3]
1429 ARF1 = ARF[0]
1430 ARF2 = ARF[1]
1431 ARF3 = ARF[2]
1432 ARF4 = ARF[3]
1433
1434 # AS with C behind the emitter --------------------------------------------
1435 # terms without aCal
1436 ATE1o, ARE1o = ATF1, ARF1
1437 ATE2o, ARE2o = 0., 0.
1438 ATE3o, ARE3o = 0., 0.
1439 ATE4o, ARE4o = ATF4, ARF4
1440 # terms with aCal
1441 ATE1a, ARE1a = 0. , 0.
1442 ATE2a, ARE2a = ATF2, ARF2
1443 ATE3a, ARE3a = -ATF3, -ARF3
1444 ATE4a, ARE4a = -2*ATF4, -2*ARF4
1445 # rotated AinEa by epsilon
1446 ATE2ae = C2e*ATF2 + S2e*ATF3
1447 ATE3ae = -S2e*ATF2 - C2e*ATF3
1448 ARE2ae = C2e*ARF2 + S2e*ARF3
1449 ARE3ae = -S2e*ARF2 - C2e*ARF3
1450
1451 ATE1 = ATE1o
1452 ATE2e = aCal*ATE2ae
1453 ATE3e = aCal*ATE3ae
1454 ATE4 = (1-2*aCal)*ATF4
1455 ARE1 = ARE1o
1456 ARE2e = aCal*ARE2ae
1457 ARE3e = aCal*ARE3ae
1458 ARE4 = (1-2*aCal)*ARF4
1459
1460 # rotated IinE
1461 QinEe = C2e*QinE + S2e*UinE
1462 UinEe = C2e*UinE - S2e*QinE
1463
1464 # --- Calibration signals and Calibration correction K from measurements with LDRCal / aCal
1465 if (TypeC == 2) or (TypeC == 1): # +++++++++ rotator calibration Eq. C.4
1466 AT = ATE1o*IinE + (ATE4o+aCal*ATE4a)*h*VinE
1467 BT = aCal * (ATE3ae*QinEe - ATE2ae*h*UinEe)
1468 AR = ARE1o*IinE + (ARE4o+aCal*ARE4a)*h*VinE
1469 BR = aCal * (ARE3ae*QinEe - ARE2ae*h*UinEe)
1470
1471 # Correction paremeters for normal measurements; they are independent of LDR
1472 if (not RotationErrorEpsilonForNormalMeasurements):
1473 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
1474 GT = ATE1o*IinE + ATE4o*h*VinE
1475 GR = ARE1o*IinE + ARE4o*h*VinE
1476 HT = ATE2a*QinE + ATE3a*h*UinEe + ATE4a*h*VinE
1477 HR = ARE2a*QinE + ARE3a*h*UinEe + ARE4a*h*VinE
1478 else:
1479 GT = ATE1o*IinE + ATE4o*h*VinE
1480 GR = ARE1o*IinE + ARE4o*h*VinE
1481 HT = ATE2ae*QinE + ATE3ae*h*UinEe + ATE4a*h*VinE
1482 HR = ARE2ae*QinE + ARE3ae*h*UinEe + ARE4a*h*VinE
1483
1484 elif (TypeC == 3) or (TypeC == 4): # +++++++++ linear polariser calibration Eq. C.5
1485 # p = +45°, m = -45°
1486 AT = ATE1*IinE + ZiC*CosC*(ATE2e*QinEe + ATE4*VinE) + ATE3e*UinEe
1487 BT = DiC*(ATE1*UinEe + ATE3e*IinE) + ZiC*SinC*(ATE4*QinEe - ATE2e*VinE)
1488 AR = ARE1*IinE + ZiC*CosC*(ARE2e*QinEe + ARE4*VinE) + ARE3e*UinEe
1489 BR = DiC*(ARE1*UinEe + ARE3e*IinE) + ZiC*SinC*(ARE4*QinEe - ARE2e*VinE)
1490
1491 # Correction paremeters for normal measurements; they are independent of LDR
1492 if (not RotationErrorEpsilonForNormalMeasurements):
1493 # Stokes Input Vector before receiver optics Eq. E.19 (after atmosphere F)
1494 GT = ATE1o*IinE + ATE4o*VinE
1495 GR = ARE1o*IinE + ARE4o*VinE
1496 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
1497 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
1498 else:
1499 D = IinE + DiC*QinEe
1500 A = DiC*IinE + QinEe
1501 B = ZiC*(CosC*UinEe + SinC*VinE)
1502 C = -ZiC*(SinC*UinEe - CosC*VinE)
1503 GT = ATE1o*D + ATE4o*C
1504 GR = ARE1o*D + ARE4o*C
1505 HT = ATE2a*A + ATE3a*B + ATE4a*C
1506 HR = ARE2a*A + ARE3a*B + ARE4a*C
1507
1508 elif (TypeC == 6): # real HWP calibration +-22.5° rotated_diattenuator_X22x5deg.odt
1509 # p = +22.5°, m = -22.5°
1510 IE1e = np.array([IinE+sqr05*DiC*QinEe, sqr05*DiC*IinE+(1-0.5*WiC)*QinEe, (1-0.5*WiC)*UinEe+sqr05*ZiC*SinC*VinE, -sqr05*ZiC*SinC*UinEe+ZiC*CosC*VinE])
1511 IE2e = np.array([sqr05*DiC*UinEe, 0.5*WiC*UinEe-sqr05*ZiC*SinC*VinE, sqr05*DiC*IinE+0.5*WiC*QinEe, sqr05*ZiC*SinC*QinEe])
1512 ATEe = np.array([ATE1,ATE2e,ATE3e,ATE4])
1513 AREe = np.array([ARE1,ARE2e,ARE3e,ARE4])
1514 AT = np.dot(ATEe,IE1e)
1515 AR = np.dot(AREe,IE1e)
1516 BT = np.dot(ATEe,IE2e)
1517 BR = np.dot(AREe,IE2e)
1518
1519 # Correction paremeters for normal measurements; they are independent of LDR
1520 if (not RotationErrorEpsilonForNormalMeasurements): # calibrator taken out
1521 GT = ATE1o*IinE + ATE4o*VinE
1522 GR = ARE1o*IinE + ARE4o*VinE
1523 HT = ATE2a*QinE + ATE3a*UinE + ATE4a*VinE
1524 HR = ARE2a*QinE + ARE3a*UinE + ARE4a*VinE
1525 else:
1526 D = IinE + DiC*QinEe
1527 A = DiC*IinE + QinEe
1528 B = ZiC*(CosC*UinEe + SinC*VinE)
1529 C = -ZiC*(SinC*UinEe - CosC*VinE)
1530 GT = ATE1o*D + ATE4o*C
1531 GR = ARE1o*D + ARE4o*C
1532 HT = ATE2a*A + ATE3a*B + ATE4a*C
1533 HR = ARE2a*A + ARE3a*B + ARE4a*C
1534
1535 else:
1536 print('Calibrator not implemented yet')
1537 sys.exit()
1538
1539 # Calibration signals with aCal => Determination of the correction K of the real calibration factor
1540 IoutTp = TaT*TiT*TiO*TiE*(AT + BT)
1541 IoutTm = TaT*TiT*TiO*TiE*(AT - BT)
1542 IoutRp = TaR*TiR*TiO*TiE*(AR + BR)
1543 IoutRm = TaR*TiR*TiO*TiE*(AR - BR)
1544 # --- Results and Corrections; electronic etaR and etaT are assumed to be 1
1545 #Eta = TiR/TiT # Eta = Eta*/K Eq. 84
1546 Etapx = IoutRp/IoutTp
1547 Etamx = IoutRm/IoutTm
1548 Etax = (Etapx*Etamx)**0.5
1549 K = Etax / Eta0
1550 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f},{4:6.3f},{5:6.3f},{6:6.3f},{7:6.3f},{8:6.3f},{9:6.3f},{10:6.3f}".format(AT, BT, AR, BR, DiC, ZiC, RetO, TP, TS, Kp, Km))
1551 #print("{0:6.3f},{1:6.3f},{2:6.3f},{3:6.3f}".format(DiC, ZiC, Kp, Km))
1552
1553 # For comparison with Volkers Libreoffice Müller Matrix spreadsheet
1554 #Eta_test_p = (IoutRp/IoutTp)
1555 #Eta_test_m = (IoutRm/IoutTm)
1556 #Eta_test = (Eta_test_p*Eta_test_m)**0.5
1557
1558 # *************************************************************************
1559 iLDR = -1
1560 for LDRTrue in LDRrange:
1561 iLDR = iLDR + 1
1562 atrue = (1-LDRTrue)/(1+LDRTrue)
1563 # ----- Forward simulated signals and LDRsim with atrue; from input file
1564 It = TaT*TiT*TiO*TiE*(GT+atrue*HT) # TaT*TiT*TiC*TiO*IinL*(GT+atrue*HT)
1565 Ir = TaR*TiR*TiO*TiE*(GR+atrue*HR) # TaR*TiR*TiC*TiO*IinL*(GR+atrue*HR)
1566
1567 # LDRsim = 1/Eta*Ir/It # simulated LDR* with Y from input file
1568 LDRsim = Ir/It # simulated uncorrected LDR with Y from input file
1569 '''
1570 if Y == 1.:
1571 LDRsimx = LDRsim
1572 LDRsimx2 = LDRsim2
1573 else:
1574 LDRsimx = 1./LDRsim
1575 LDRsimx2 = 1./LDRsim2
1576 '''
1577 # ----- Backward correction
1578 # Corrected LDRCorr from forward simulated LDRsim (atrue) with assumed true G0,H0,K0
1579 LDRCorr = (LDRsim*K0/Etax*(GT0+HT0)-(GR0+HR0))/((GR0-HR0)-LDRsim*K0/Etax*(GT0-HT0))
1580
1581 # -- F11corr from It and Ir and calibration EtaX
1582 Text1 = "F11corr from It and Ir with calibration EtaX: x-axis: F11corr(LDRtrue) / F11corr(LDRtrue = 0.004) - 1"
1583 F11corr = 1/(TiO*TiE)*((HR0*Etax/K0*It/TTa-HT0*Ir/TRa)/(HR0*GT0-HT0*GR0)) # IL = 1 Eq.(64)
1584
1585 #Text1 = "F11corr from It and Ir without corrections but with calibration EtaX: x-axis: F11corr(LDRtrue) devided by F11corr(LDRtrue = 0.004)"
1586 #F11corr = 0.5/(TiO*TiE)*(Etax*It/TTa+Ir/TRa) # IL = 1 Eq.(64)
1587
1588 # -- It from It only with atrue without corrections - for BERTHA (and PollyXTs)
1589 #Text1 = " x-axis: IT(LDRtrue) / IT(LDRtrue = 0.004) - 1"
1590 #F11corr = It/(TaT*TiT*TiO*TiE) #/(TaT*TiT*TiO*TiE*(GT0+atrue*HT0))
1591 # !!! see below line 1673ff
1592
1593 aF11corr[iLDR,iN] = F11corr
1594 aA[iLDR,iN] = LDRCorr
1595
1596 aX[0,iN] = GR
1597 aX[1,iN] = GT
1598 aX[2,iN] = HR
1599 aX[3,iN] = HT
1600 aX[4,iN] = K
1601
1602 aLDRCal[iN] = iLDRCal
1603 aERaT[iN] = iERaT
1604 aERaR[iN] = iERaR
1605 aRotaT[iN] = iRotaT
1606 aRotaR[iN] = iRotaR
1607 aRetT[iN] = iRetT
1608 aRetR[iN] = iRetR
1609
1610 aRotL[iN] = iRotL
1611 aRotE[iN] = iRotE
1612 aRetE[iN] = iRetE
1613 aRotO[iN] = iRotO
1614 aRetO[iN] = iRetO
1615 aRotC[iN] = iRotC
1616 aRetC[iN] = iRetC
1617 aDiO[iN] = iDiO
1618 aDiE[iN] = iDiE
1619 aDiC[iN] = iDiC
1620 aTP[iN] = iTP
1621 aTS[iN] = iTS
1622 aRP[iN] = iRP
1623 aRS[iN] = iRS
1624
1625 # --- END loop
1626 btime = clock()
1627 print("\r done in ", "{0:5.0f}".format(btime-atime), "sec") #, end="\r")
1628
1629 # --- Plot -----------------------------------------------------------------
1630 #sns.set_style("whitegrid")
1631 #sns.set_palette("bright", 6)
1632
1633 '''
1634 fig2 = plt.figure()
1635 plt.plot(aA[2,:],'b.')
1636 plt.plot(aA[3,:],'r.')
1637 plt.plot(aA[4,:],'g.')
1638 #plt.plot(aA[6,:],'c.')
1639 plt.show
1640 '''
1641 # Plot LDR
1642 def PlotSubHist(aVar, aX, X0, daX, iaX, naX):
1643 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
1644 iLDR = -1
1645 for LDRTrue in LDRrange:
1646 iLDR = iLDR + 1
1647
1648 LDRmin[iLDR] = np.min(aA[iLDR,:])
1649 LDRmax[iLDR] = np.max(aA[iLDR,:])
1650 Rmin = LDRmin[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995
1651 Rmax = LDRmax[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005
1652
1653 #plt.subplot(5,2,iLDR+1)
1654 plt.subplot(1,5,iLDR+1)
1655 (n, bins, patches) = plt.hist(aA[iLDR,:],
1656 bins=100, log=False,
1657 range=[Rmin, Rmax],
1658 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled')
1659
1660 for iaX in range(-naX,naX+1):
1661 plt.hist(aA[iLDR,aX == iaX],
1662 range=[Rmin, Rmax],
1663 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5)))
1664
1665 if (iLDR == 2): plt.legend()
1666
1667 plt.tick_params(axis='both', labelsize=9)
1668 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
1669
1670 #plt.title(LID + ' ' + aVar, fontsize=18)
1671 #plt.ylabel('frequency', fontsize=10)
1672 #plt.xlabel('LDRcorr', fontsize=10)
1673 #fig.tight_layout()
1674 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
1675 #plt.show()
1676 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
1677 #plt.close
1678 return
1679
1680 if (nRotL > 0): PlotSubHist("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
1681 if (nRetE > 0): PlotSubHist("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
1682 if (nRotE > 0): PlotSubHist("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
1683 if (nDiE > 0): PlotSubHist("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
1684 if (nRetO > 0): PlotSubHist("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
1685 if (nRotO > 0): PlotSubHist("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
1686 if (nDiO > 0): PlotSubHist("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
1687 if (nDiC > 0): PlotSubHist("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
1688 if (nRotC > 0): PlotSubHist("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
1689 if (nRetC > 0): PlotSubHist("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
1690 if (nTP > 0): PlotSubHist("TP", aTP, TP0, dTP, iTP, nTP)
1691 if (nTS > 0): PlotSubHist("TS", aTS, TS0, dTS, iTS, nTS)
1692 if (nRP > 0): PlotSubHist("RP", aRP, RP0, dRP, iRP, nRP)
1693 if (nRS > 0): PlotSubHist("RS", aRS, RS0, dRS, iRS, nRS)
1694 if (nRetT > 0): PlotSubHist("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
1695 if (nRetR > 0): PlotSubHist("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
1696 if (nERaT > 0): PlotSubHist("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
1697 if (nERaR > 0): PlotSubHist("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
1698 if (nRotaT > 0): PlotSubHist("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
1699 if (nRotaR > 0): PlotSubHist("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
1700 if (nLDRCal > 0): PlotSubHist("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
1701
1702 plt.show()
1703 plt.close
1704
1705 print()
1706 #print("IT(LDRtrue) devided by IT(LDRtrue = 0.004)")
1707 print(Text1)
1708 print()
1709
1710 iLDR = 5
1711 for LDRTrue in LDRrange:
1712 iLDR = iLDR - 1
1713 aF11corr[iLDR,:] = aF11corr[iLDR,:] / aF11corr[0,:] - 1.0
1714
1715 # Plot F11
1716 def PlotSubHistF11(aVar, aX, X0, daX, iaX, naX):
1717 fig, ax = plt.subplots(nrows=1, ncols=5, sharex=True, sharey=True, figsize=(25, 2))
1718 iLDR = -1
1719 for LDRTrue in LDRrange:
1720 iLDR = iLDR + 1
1721
1722 '''
1723 F11min[iLDR] = np.min(aF11corr[iLDR,:])
1724 F11max[iLDR] = np.max(aF11corr[iLDR,:])
1725 Rmin = F11min[iLDR] * 0.995 # np.min(aA[iLDR,:]) * 0.995
1726 Rmax = F11max[iLDR] * 1.005 # np.max(aA[iLDR,:]) * 1.005
1727 '''
1728 #Rmin = 0.8
1729 #Rmax = 1.2
1730
1731 #plt.subplot(5,2,iLDR+1)
1732 plt.subplot(1,5,iLDR+1)
1733 (n, bins, patches) = plt.hist(aF11corr[iLDR,:],
1734 bins=100, log=False,
1735 alpha=0.5, normed=False, color = '0.5', histtype='stepfilled')
1736
1737 for iaX in range(-naX,naX+1):
1738 plt.hist(aF11corr[iLDR,aX == iaX],
1739 bins=100, log=False, alpha=0.3, normed=False, histtype='stepfilled', label = str(round(X0 + iaX*daX/naX,5)))
1740
1741 if (iLDR == 2): plt.legend()
1742
1743 plt.tick_params(axis='both', labelsize=9)
1744 #plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
1745
1746 #plt.title(LID + ' ' + aVar, fontsize=18)
1747 #plt.ylabel('frequency', fontsize=10)
1748 #plt.xlabel('LDRcorr', fontsize=10)
1749 #fig.tight_layout()
1750 fig.suptitle(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]) + ' - ' + aVar, fontsize=14, y=1.05)
1751 #plt.show()
1752 #fig.savefig(LID + '_' + aVar + '.png', dpi=150, bbox_inches='tight', pad_inches=0)
1753 #plt.close
1754 return
1755
1756 if (nRotL > 0): PlotSubHistF11("RotL", aRotL, RotL0, dRotL, iRotL, nRotL)
1757 if (nRetE > 0): PlotSubHistF11("RetE", aRetE, RetE0, dRetE, iRetE, nRetE)
1758 if (nRotE > 0): PlotSubHistF11("RotE", aRotE, RotE0, dRotE, iRotE, nRotE)
1759 if (nDiE > 0): PlotSubHistF11("DiE", aDiE, DiE0, dDiE, iDiE, nDiE)
1760 if (nRetO > 0): PlotSubHistF11("RetO", aRetO, RetO0, dRetO, iRetO, nRetO)
1761 if (nRotO > 0): PlotSubHistF11("RotO", aRotO, RotO0, dRotO, iRotO, nRotO)
1762 if (nDiO > 0): PlotSubHistF11("DiO", aDiO, DiO0, dDiO, iDiO, nDiO)
1763 if (nDiC > 0): PlotSubHistF11("DiC", aDiC, DiC0, dDiC, iDiC, nDiC)
1764 if (nRotC > 0): PlotSubHistF11("RotC", aRotC, RotC0, dRotC, iRotC, nRotC)
1765 if (nRetC > 0): PlotSubHistF11("RetC", aRetC, RetC0, dRetC, iRetC, nRetC)
1766 if (nTP > 0): PlotSubHistF11("TP", aTP, TP0, dTP, iTP, nTP)
1767 if (nTS > 0): PlotSubHistF11("TS", aTS, TS0, dTS, iTS, nTS)
1768 if (nRP > 0): PlotSubHistF11("RP", aRP, RP0, dRP, iRP, nRP)
1769 if (nRS > 0): PlotSubHistF11("RS", aRS, RS0, dRS, iRS, nRS)
1770 if (nRetT > 0): PlotSubHistF11("RetT", aRetT, RetT0, dRetT, iRetT, nRetT)
1771 if (nRetR > 0): PlotSubHistF11("RetR", aRetR, RetR0, dRetR, iRetR, nRetR)
1772 if (nERaT > 0): PlotSubHistF11("ERaT", aERaT, ERaT0, dERaT, iERaT, nERaT)
1773 if (nERaR > 0): PlotSubHistF11("ERaR", aERaR, ERaR0, dERaR, iERaR, nERaR)
1774 if (nRotaT > 0): PlotSubHistF11("RotaT", aRotaT, RotaT0, dRotaT, iRotaT, nRotaT)
1775 if (nRotaR > 0): PlotSubHistF11("RotaR", aRotaR, RotaR0, dRotaR, iRotaR, nRotaR)
1776 if (nLDRCal > 0): PlotSubHistF11("LDRCal", aLDRCal, LDRCal0, dLDRCal, iLDRCal, nLDRCal)
1777
1778 plt.show()
1779 plt.close
1780 '''
1781 # only histogram
1782 #print("******************* " + aVar + " *******************")
1783 fig, ax = plt.subplots(nrows=5, ncols=2, sharex=True, sharey=True, figsize=(10, 10))
1784 iLDR = -1
1785 for LDRTrue in LDRrange:
1786 iLDR = iLDR + 1
1787 LDRmin[iLDR] = np.min(aA[iLDR,:])
1788 LDRmax[iLDR] = np.max(aA[iLDR,:])
1789 Rmin = np.min(aA[iLDR,:]) * 0.999
1790 Rmax = np.max(aA[iLDR,:]) * 1.001
1791 plt.subplot(5,2,iLDR+1)
1792 (n, bins, patches) = plt.hist(aA[iLDR,:],
1793 range=[Rmin, Rmax],
1794 bins=200, log=False, alpha=0.2, normed=False, color = '0.5', histtype='stepfilled')
1795 plt.tick_params(axis='both', labelsize=9)
1796 plt.plot([LDRTrue, LDRTrue], [0, np.max(n)], 'r-', lw=2)
1797 plt.show()
1798 plt.close
1799 '''
1800
1801 # --- Plot LDRmin, LDRmax
1802 fig2 = plt.figure()
1803 plt.plot(LDRrange,LDRmax-LDRrange, linewidth=2.0, color='b')
1804 plt.plot(LDRrange,LDRmin-LDRrange, linewidth=2.0, color='g')
1805
1806 plt.xlabel('LDRtrue', fontsize=18)
1807 plt.ylabel('LDRTrue-LDRmin, LDRTrue-LDRmax', fontsize=14)
1808 plt.title(LID + ' ' + str(Type[TypeC]) + ' ' + str(Loc[LocC]), fontsize=18)
1809 #plt.ylimit(-0.07, 0.07)
1810 plt.show()
1811 plt.close
1812
1813 # --- Save LDRmin, LDRmax to file
1814 # http://stackoverflow.com/questions/4675728/redirect-stdout-to-a-file-in-python
1815 with open('LDR_min_max_ver7_' + LID + '.dat', 'w') as f:
1816 with redirect_stdout(f):
1817 print(LID)
1818 print("LDRtrue, LDRmin, LDRmax")
1819 for i in range(len(LDRrange)):
1820 print("{0:7.4f},{1:7.4f},{2:7.4f}".format(LDRrange[i], LDRmin[i], LDRmax[i]))
1821
1822 '''
1823 # --- Plot K over LDRCal
1824 fig3 = plt.figure()
1825 plt.plot(LDRCal0+aLDRCal*dLDRCal/nLDRCal,aX[4,:], linewidth=2.0, color='b')
1826
1827 plt.xlabel('LDRCal', fontsize=18)
1828 plt.ylabel('K', fontsize=14)
1829 plt.title(LID, fontsize=18)
1830 plt.show()
1831 plt.close
1832 '''
1833
1834 # Additional plot routines ======>
1835 '''
1836 #******************************************************************************
1837 # 1. Plot LDRcorrected - LDR(measured Icross/Iparallel)
1838 LDRa = np.arange(1.,100.)*0.005
1839 LDRCorra = np.arange(1.,100.)
1840 if Y == - 1.: LDRa = 1./LDRa
1841 LDRCorra = (1./Eta*LDRa*(GT+HT)-(GR+HR))/((GR-HR)-1./Eta*LDRa*(GT-HT))
1842 if Y == - 1.: LDRa = 1./LDRa
1843 #
1844 #fig = plt.figure()
1845 plt.plot(LDRa,LDRCorra-LDRa)
1846 plt.plot([0.,0.5],[0.,0.5])
1847 plt.suptitle('LDRcorrected - LDR(measured Icross/Iparallel)', fontsize=16)
1848 plt.xlabel('LDR', fontsize=18)
1849 plt.ylabel('LDRCorr - LDR', fontsize=16)
1850 #plt.savefig('test.png')
1851 #
1852 '''
1853 '''
1854 #******************************************************************************
1855 # 2. Plot LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue
1856 LDRa = np.arange(1.,100.)*0.005
1857 LDRsima = np.arange(1.,100.)
1858
1859 atruea = (1.-LDRa)/(1+LDRa)
1860 Ita = TiT*TiO*IinL*(GT+atruea*HT)
1861 Ira = TiR*TiO*IinL*(GR+atruea*HR)
1862 LDRsima = Ira/Ita # simulated uncorrected LDR with Y from input file
1863 if Y == -1.: LDRsima = 1./LDRsima
1864 #
1865 #fig = plt.figure()
1866 plt.plot(LDRa,LDRsima)
1867 plt.plot([0.,0.5],[0.,0.5])
1868 plt.suptitle('LDRsim (simulated measurements without corrections = Icross/Iparallel) over LDRtrue', fontsize=10)
1869 plt.xlabel('LDRtrue', fontsize=18)
1870 plt.ylabel('LDRsim', fontsize=16)
1871 #plt.savefig('test.png')
1872 #
1873 '''

EARLINET Repository: Lidar polarisation GHK-parameter and error calculation

mercurial