Coverage for src/susi/model/pol_eff.py: 54%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3import os 

4import numpy as np 

5import matplotlib.pyplot as plt 

6from scipy import interpolate 

7from scipy.io import readsav 

8from .. import ROOT_DIR 

9 

10#: File path to the SAV files 

11FILE_PATH = os.path.abspath(os.path.join(ROOT_DIR, "..", "data", "susi", "modeled", "susi_pol_eff_20221105.sav")) 

12 

13#: Stokes names for plotting 

14STKN = ["I", "Q", "U", "V"] 

15 

16 

17def pol_eff(wl): 

18 """ 

19 Returns modeled polarization efficiency ([cam1], [cam2]) for a given wavelength as interpolated from .sav file. 

20 

21 :param wl: Array or scalar with wavelengths at wich you want to get the efficiency [nm]. 

22 Must be in the 300-410 nm range 

23 

24 """ 

25 if np.min(wl) < 300.0 or np.max(wl) > 400.0: 

26 raise Exception("The input wavelength is not int the 300-400 nm range") 

27 

28 eff = readsav(FILE_PATH) 

29 oeff1, oeff2 = [], [] 

30 for i in range(4): 

31 tck = interpolate.splrep(eff["wl"], eff["all_eff"][0, :, i], s=0) 

32 oeff1.append((interpolate.splev(wl, tck))) 

33 tck = interpolate.splrep(eff["wl"], eff["all_eff"][1, :, i], s=0) 

34 oeff2.append((interpolate.splev(wl, tck))) 

35 return np.array([oeff1, oeff2]) 

36 

37 

38def plot_eff(wl): 

39 """ 

40 plots the polarization efficiency for a given wavelength as interpolated from .sav file. 

41 

42 :param wl: Array with wl range to plot [min, max] in nm 

43 """ 

44 fig = plt.figure(figsize=(21, 7)) 

45 for i in range(4): 

46 fig.add_subplot(1, 4, i + 1) 

47 plt.title(STKN[i]) 

48 plt.xlabel(r"$\lambda$ [nm]") 

49 plt.plot(wl, pol_eff(wl)[1, i, :], "--r", label="SUSIM_US550") 

50 plt.plot(wl, pol_eff(wl)[0, i, :], "--b", label="SUSIM_US560") 

51 if i == 0: 

52 plt.ylim([0, 1.1]) 

53 plt.plot([np.min(wl), np.max(wl)], [1, 1], "-", color="black") 

54 plt.legend(loc="best") 

55 else: 

56 plt.ylim([0, 1]) 

57 plt.plot([np.min(wl), np.max(wl)], [1 / np.sqrt(3), 1 / np.sqrt(3)], "-", color="black") 

58 plt.grid(which="both", axis="both", linestyle="--", linewidth=0.5) 

59 plt.show()