7207821764

Committed 14 Dec 2023 10:49AM UTC coverage: 90.688% (-0.7%) from 91.426%

Build # 7207821764

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #145 from LSSTDESC/144-run_module_suite-deprecated

test remove run_module_suite

Run Details

2 of 7 new or added lines in 2 files covered. (28.57%)

93 existing lines in 21 files now uncovered.

7226 of 7968 relevant lines covered (90.69%)

0.91 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.9

/spectractor/simulation/throughput.py

import os
import numpy as np
import matplotlib.pyplot as plt

from scipy.interpolate import interp1d

from spectractor.config import set_logger
import spectractor.parameters as parameters


def load_transmission(file_name):
    """Load the transmission files and crop in wavelength using LAMBDA_MIN and LAMBDA_MAX.

    The input file must have two or three columns:
        1. wavelengths in nm
        2. transmissions between 0 and 1.
        3. uncertainties on the transmissions (optional)

    Returns
    -------
    lambdas: array_like
        The wavelengths array in nm.
    transmissions: array_like
        The transmission array, values are between 0 and 1.
    uncertainties: array_like
        The uncertainty on the transmission array (0 if file does not contain the info).

    Examples
    --------
    >>> parameters.LAMBDA_MIN = 500
    >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "qecurve.txt"))
    >>> print(lambdas[:3])
    [500.81855389 508.18553888 519.23601637]
    >>> print(transmissions[:3])
    [0.74355972 0.75526932 0.76932084]
    >>> print(errors[:3])
    [0. 0. 0.]

    """
    if os.path.isabs(file_name):
        path = file_name
    else:
        mypath = os.path.dirname(os.path.dirname(__file__))
        path = os.path.join(mypath, file_name)
    data = np.loadtxt(path).T
    lambdas = data[0]
    sorted_indices = lambdas.argsort()
    lambdas = lambdas[sorted_indices]
    y = data[1][sorted_indices]
    err = np.zeros_like(y)
    if data.shape[0] == 3:
        err = data[2][sorted_indices]
    indexes = np.logical_and(lambdas > parameters.LAMBDA_MIN, lambdas < parameters.LAMBDA_MAX)
    return lambdas[indexes], y[indexes], err[indexes]


def plot_transmission_simple(ax, lambdas, transmissions,  uncertainties=None, label="", title="", lw=2):
    """Plot the transmission with respect to the wavelength.

    Parameters
    ----------
    ax: Axes
        An Axes instance.
    lambdas: array_like
        The wavelengths array in nm.
    transmissions: array_like
        The transmission array, values are between 0 and 1.
    uncertainties: array_like, optional
        The uncertainty on the transmission array (default: None).
    label: str, optional
        The label of the curve for the legend (default: "")
    title: str, optional
        The title of the plot (default: "")
    lw: int
        Line width (default: 2).

    Examples
    --------

    .. plot::
        :include-source:

        >>> from spectractor.simulation.atmosphere import plot_transmission_simple
        >>> from spectractor import parameters
        >>> fig = plt.figure()
        >>> ax = plt.gca()
        >>> parameters.LAMBDA_MIN = 500
        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "qecurve.txt"))
        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Quantum efficiency")
        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "lsst_mirrorthroughput.txt"))
        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Mirror 1")
        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "FGB37.txt"))
        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="FGB37")
        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "RG715.txt"))
        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="RG715")
        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, parameters.OBS_FULL_INSTRUMENT_TRANSMISSON))
        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Full instrument")
        >>> if parameters.DISPLAY: plt.show()

    """
    if uncertainties is None or np.all(np.isclose(uncertainties, np.zeros_like(transmissions))):
        ax.plot(lambdas, transmissions, "-", label=label, lw=lw)
    else:
        ax.errorbar(lambdas, transmissions, yerr=uncertainties, label=label, lw=lw)
    if title != "":
        ax.set_title(title)
    ax.set_xlabel(r"$\lambda$ [nm]")
    ax.set_ylabel("Transmission")
    ax.set_xlim(parameters.LAMBDA_MIN, parameters.LAMBDA_MAX)
    ax.grid()
    if label != "":
        ax.legend(loc="best")


class TelescopeTransmission:

    def __init__(self, filter_label=""):
        """Transmission of the telescope as product of the following transmissions:

        - mirrors
        - throughput
        - quantum efficiency
        - Filter

        Parameters
        ----------
        filter_label: str, optional
            The filter string name.

        Examples
        --------
        >>> t = TelescopeTransmission()
        >>> t.plot_transmission()
        """

        self.my_logger = set_logger(self.__class__.__name__)
        self.filter_label = filter_label
        self.transmission = None
        self.transmission_err = None
        self.load_transmission()

    def load_transmission(self):
        """Load the transmission files and make a function.

        Returns
        -------
        transmission: callable
            The transmission function using wavelengths in nm.

        Examples
        --------
        >>> t = TelescopeTransmission()
        >>> t.plot_transmission()

        >>> t2 = TelescopeTransmission(filter_label="RG715")
        >>> t2.plot_transmission()

        .. doctest:
            :hide:

            >>> assert t.transmission is not None
            >>> assert t.transmission_err is not None
            >>> assert t2.transmission is not None
            >>> assert t2.transmission_err is not None
            >>> assert np.sum(t.transmission(parameters.LAMBDAS)) > np.sum(t2.transmission(parameters.LAMBDAS))

        """
        wl, trm, err = load_transmission(os.path.join(parameters.THROUGHPUT_DIR,
                                                      parameters.OBS_FULL_INSTRUMENT_TRANSMISSON))
        to = interp1d(wl, trm, kind='linear', bounds_error=False, fill_value=0.)
        err = np.sqrt(err ** 2 + parameters.OBS_TRANSMISSION_SYSTEMATICS ** 2)
        to_err = interp1d(wl, err, kind='linear', bounds_error=False, fill_value=0.)

        TF = lambda x: 1
        TF_err = lambda x: 0
        if self.filter_label != "" and "empty" not in self.filter_label.lower():
            if ".txt" in self.filter_label:
                filter_filename = self.filter_label
            else:
                filter_filename = self.filter_label + ".txt"
            wl, trb, err = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, filter_filename))
            TF = interp1d(wl, trb, kind='linear', bounds_error=False, fill_value=0.)
            TF_err = interp1d(wl, err, kind='linear', bounds_error=False, fill_value=0.)

        # self.transmission=lambda x: self.qe(x)*self.to(x)*(self.tm(x)**2)*self.tf(x)
        self.transmission = lambda x: to(x) * TF(x)
        self.transmission_err = lambda x: np.sqrt(to_err(x)**2 + TF_err(x)**2)
        return self.transmission

    def plot_transmission(self):
        """Plot the transmission function and the associated uncertainties.

        Examples
        --------
        >>> t = TelescopeTransmission()
        >>> t.plot_transmission()

        """
        plt.figure()
        plot_transmission_simple(plt.gca(), parameters.LAMBDAS, self.transmission(parameters.LAMBDAS),
                                 uncertainties=self.transmission_err(parameters.LAMBDAS))
        if parameters.DISPLAY:
            plt.show()
        else:
            plt.close('all')

    def reset_lambda_range(self, transmission_threshold=1e-4):
        integral = np.cumsum(self.transmission(parameters.LAMBDAS))
        lambda_min = parameters.LAMBDAS[0]
        for k, tr in enumerate(integral):
            if tr > transmission_threshold:
                lambda_min = parameters.LAMBDAS[k]
                break
        lambda_max = parameters.LAMBDAS[0]
        for k, tr in enumerate(integral):
            if tr > integral[-1] - transmission_threshold:
                lambda_max = parameters.LAMBDAS[k]
                break
        parameters.LAMBDA_MIN = max(lambda_min, parameters.LAMBDA_MIN)
        parameters.LAMBDA_MAX = min(lambda_max, parameters.LAMBDA_MAX)
        parameters.LAMBDAS = np.arange(parameters.LAMBDA_MIN, parameters.LAMBDA_MAX, parameters.LAMBDA_STEP)
        self.my_logger.info(f"\n\tWith filter {self.filter_label}, set parameters.LAMBDA_MIN={parameters.LAMBDA_MIN} "
                            f"and parameters.LAMBDA_MAX={parameters.LAMBDA_MAX}.")


if __name__ == "__main__":
    import doctest

    doctest.testmod()

1	import os	1✔
2	import numpy as np	1✔
3	import matplotlib.pyplot as plt	1✔
4
5	from scipy.interpolate import interp1d	1✔
6
7	from spectractor.config import set_logger	1✔
8	import spectractor.parameters as parameters	1✔
9
10
11	def load_transmission(file_name):	1✔
12	"""Load the transmission files and crop in wavelength using LAMBDA_MIN and LAMBDA_MAX.
13
14	The input file must have two or three columns:
15	1. wavelengths in nm
16	2. transmissions between 0 and 1.
17	3. uncertainties on the transmissions (optional)
18
19	Returns
20	-------
21	lambdas: array_like
22	The wavelengths array in nm.
23	transmissions: array_like
24	The transmission array, values are between 0 and 1.
25	uncertainties: array_like
26	The uncertainty on the transmission array (0 if file does not contain the info).
27
28	Examples
29	--------
30	>>> parameters.LAMBDA_MIN = 500
31	>>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "qecurve.txt"))
32	>>> print(lambdas[:3])
33	[500.81855389 508.18553888 519.23601637]
34	>>> print(transmissions[:3])
35	[0.74355972 0.75526932 0.76932084]
36	>>> print(errors[:3])
37	[0. 0. 0.]
38
39	"""
40	if os.path.isabs(file_name):	1✔
41	path = file_name	×
42	else:
43	mypath = os.path.dirname(os.path.dirname(__file__))	1✔
44	path = os.path.join(mypath, file_name)	1✔
45	data = np.loadtxt(path).T	1✔
46	lambdas = data[0]	1✔
47	sorted_indices = lambdas.argsort()	1✔
48	lambdas = lambdas[sorted_indices]	1✔
49	y = data[1][sorted_indices]	1✔
50	err = np.zeros_like(y)	1✔
51	if data.shape[0] == 3:	1✔
52	err = data[2][sorted_indices]	1✔
53	indexes = np.logical_and(lambdas > parameters.LAMBDA_MIN, lambdas < parameters.LAMBDA_MAX)	1✔
54	return lambdas[indexes], y[indexes], err[indexes]	1✔
55
56
57	def plot_transmission_simple(ax, lambdas, transmissions, uncertainties=None, label="", title="", lw=2):	1✔
58	"""Plot the transmission with respect to the wavelength.
59
60	Parameters
61	----------
62	ax: Axes
63	An Axes instance.
64	lambdas: array_like
65	The wavelengths array in nm.
66	transmissions: array_like
67	The transmission array, values are between 0 and 1.
68	uncertainties: array_like, optional
69	The uncertainty on the transmission array (default: None).
70	label: str, optional
71	The label of the curve for the legend (default: "")
72	title: str, optional
73	The title of the plot (default: "")
74	lw: int
75	Line width (default: 2).
76
77	Examples
78	--------
79
80	.. plot::
81	:include-source:
82
83	>>> from spectractor.simulation.atmosphere import plot_transmission_simple
84	>>> from spectractor import parameters
85	>>> fig = plt.figure()
86	>>> ax = plt.gca()
87	>>> parameters.LAMBDA_MIN = 500
88	>>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "qecurve.txt"))
89	>>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Quantum efficiency")
90	>>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "lsst_mirrorthroughput.txt"))
91	>>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Mirror 1")
92	>>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "FGB37.txt"))
93	>>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="FGB37")
94	>>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "RG715.txt"))
95	>>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="RG715")
96	>>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, parameters.OBS_FULL_INSTRUMENT_TRANSMISSON))
97	>>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Full instrument")
98	>>> if parameters.DISPLAY: plt.show()
99
100	"""
101	if uncertainties is None or np.all(np.isclose(uncertainties, np.zeros_like(transmissions))):	1✔
102	ax.plot(lambdas, transmissions, "-", label=label, lw=lw)	1✔
103	else:
104	ax.errorbar(lambdas, transmissions, yerr=uncertainties, label=label, lw=lw)	1✔
105	if title != "":	1✔
106	ax.set_title(title)	1✔
107	ax.set_xlabel(r"$\lambda$ [nm]")	1✔
108	ax.set_ylabel("Transmission")	1✔
109	ax.set_xlim(parameters.LAMBDA_MIN, parameters.LAMBDA_MAX)	1✔
110	ax.grid()	1✔
111	if label != "":	1✔
112	ax.legend(loc="best")	1✔
113
114
115	class TelescopeTransmission:	1✔
116
117	def __init__(self, filter_label=""):	1✔
118	"""Transmission of the telescope as product of the following transmissions:
119
120	- mirrors
121	- throughput
122	- quantum efficiency
123	- Filter
124
125	Parameters
126	----------
127	filter_label: str, optional
128	The filter string name.
129
130	Examples
131	--------
132	>>> t = TelescopeTransmission()
133	>>> t.plot_transmission()
134	"""
135
136	self.my_logger = set_logger(self.__class__.__name__)	1✔
137	self.filter_label = filter_label	1✔
138	self.transmission = None	1✔
139	self.transmission_err = None	1✔
140	self.load_transmission()	1✔
141
142	def load_transmission(self):	1✔
143	"""Load the transmission files and make a function.
144
145	Returns
146	-------
147	transmission: callable
148	The transmission function using wavelengths in nm.
149
150	Examples
151	--------
152	>>> t = TelescopeTransmission()
153	>>> t.plot_transmission()
154
155	>>> t2 = TelescopeTransmission(filter_label="RG715")
156	>>> t2.plot_transmission()
157
158	.. doctest:
159	:hide:
160
161	>>> assert t.transmission is not None
162	>>> assert t.transmission_err is not None
163	>>> assert t2.transmission is not None
164	>>> assert t2.transmission_err is not None
165	>>> assert np.sum(t.transmission(parameters.LAMBDAS)) > np.sum(t2.transmission(parameters.LAMBDAS))
166
167	"""
168	wl, trm, err = load_transmission(os.path.join(parameters.THROUGHPUT_DIR,	1✔
169	parameters.OBS_FULL_INSTRUMENT_TRANSMISSON))
170	to = interp1d(wl, trm, kind='linear', bounds_error=False, fill_value=0.)	1✔
171	err = np.sqrt(err 2 + parameters.OBS_TRANSMISSION_SYSTEMATICS 2)	1✔
172	to_err = interp1d(wl, err, kind='linear', bounds_error=False, fill_value=0.)	1✔
173
174	TF = lambda x: 1	1✔
175	TF_err = lambda x: 0	1✔
176	if self.filter_label != "" and "empty" not in self.filter_label.lower():	1✔
177	if ".txt" in self.filter_label:	1✔
178	filter_filename = self.filter_label	×
179	else:
180	filter_filename = self.filter_label + ".txt"	1✔
181	wl, trb, err = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, filter_filename))	1✔
182	TF = interp1d(wl, trb, kind='linear', bounds_error=False, fill_value=0.)	1✔
183	TF_err = interp1d(wl, err, kind='linear', bounds_error=False, fill_value=0.)	1✔
184
185	# self.transmission=lambda x: self.qe(x)self.to(x)(self.tm(x)*2)self.tf(x)
186	self.transmission = lambda x: to(x) * TF(x)	1✔
187	self.transmission_err = lambda x: np.sqrt(to_err(x)2 + TF_err(x)2)	1✔
188	return self.transmission	1✔
189
190	def plot_transmission(self):	1✔
191	"""Plot the transmission function and the associated uncertainties.
192
193	Examples
194	--------
195	>>> t = TelescopeTransmission()
196	>>> t.plot_transmission()
197
198	"""
199	plt.figure()	1✔
200	plot_transmission_simple(plt.gca(), parameters.LAMBDAS, self.transmission(parameters.LAMBDAS),	1✔
201	uncertainties=self.transmission_err(parameters.LAMBDAS))
202	if parameters.DISPLAY:	1✔
203	plt.show()	×
204	else:
205	plt.close('all')	1✔
206
207	def reset_lambda_range(self, transmission_threshold=1e-4):	1✔
208	integral = np.cumsum(self.transmission(parameters.LAMBDAS))	1✔
209	lambda_min = parameters.LAMBDAS[0]	1✔
210	for k, tr in enumerate(integral):	1✔
211	if tr > transmission_threshold:	1✔
212	lambda_min = parameters.LAMBDAS[k]	1✔
213	break	1✔
214	lambda_max = parameters.LAMBDAS[0]	1✔
215	for k, tr in enumerate(integral):	1✔
216	if tr > integral[-1] - transmission_threshold:	1✔
217	lambda_max = parameters.LAMBDAS[k]	1✔
218	break	1✔
219	parameters.LAMBDA_MIN = max(lambda_min, parameters.LAMBDA_MIN)	1✔
220	parameters.LAMBDA_MAX = min(lambda_max, parameters.LAMBDA_MAX)	1✔
221	parameters.LAMBDAS = np.arange(parameters.LAMBDA_MIN, parameters.LAMBDA_MAX, parameters.LAMBDA_STEP)	1✔
222	self.my_logger.info(f"\n\tWith filter {self.filter_label}, set parameters.LAMBDA_MIN={parameters.LAMBDA_MIN} "	1✔
223	f"and parameters.LAMBDA_MAX={parameters.LAMBDA_MAX}.")
224
225
226	if __name__ == "__main__":	1✔
UNCOV 227	import doctest	×
228
UNCOV 229	doctest.testmod()	×

LSSTDESC / Spectractor / 7207821764

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous