Mercurial > public > atmospheric_lidar_ghk / file comparison

comparison: lidar_correction_ghk_pollyxt.py

lidar_correction_ghk_pollyxt.py

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

EARLINET Repository: Lidar polarisation GHK-parameter and error calculation

mercurial