tdsat/duet_filters.py at master · hpearnshaw/tdsat · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def make_o2_filter(in_wave,**kwargs):
    """
    Mocks up a filter function to remove the O2 line flux.
    and "rejection" out of band.

    Optional inputs (defaults):
    low_wave = lower wavelength of the band (180*ur.nm)
    high_wave = lower wavelength of the band (220*ur.nm)
    rejection = Out of band rejection level (1e-3)
    diag = Give diagnostic info (False)


    Syntax:
    filter = make_filter(wave)


    """
    import astropy.units as ur
    import numpy as np

    low_wave = kwargs.pop('low_wave', 2468*ur.AA)
    high_wave = kwargs.pop('high_wave', 2472*ur.AA)
    rejection = kwargs.pop('rejection', 1e-3)
    diag = kwargs.pop('diag', False)


    wave = in_wave.to(ur.Angstrom).value

    low_check = low_wave.to(ur.Angstrom).value
    high_check = high_wave.to(ur.Angstrom).value


    red_filter = np.zeros_like(wave)
    red_filter[(wave<low_check) | (wave>high_check)] = 1.0
    red_filter[(wave>=low_check) & (wave<=high_check)] += rejection

    return red_filter


def make_red_filter(in_wave, **kwargs):
    """
    Mocks up a filter function at all of the "wave" points with 1 "in band"
    and "rejection" out of band.

    Optional inputs (defaults):
    low_wave = lower wavelength of the band (180*ur.nm)
    high_wave = lower wavelength of the band (220*ur.nm)
    rejection = Out of band rejection level (1e-3)
    diag = Give diagnostic info (False)


    Syntax:
    filter = make_filter(wave)


    """
    import astropy.units as ur
    import numpy as np


    low_wave = kwargs.pop('low_wave', 180*ur.nm)
    high_wave = kwargs.pop('high_wave', 220*ur.nm)
    rejection = kwargs.pop('rejection', 1e-3)
    diag = kwargs.pop('diag', False)
    blue_filter = kwargs.pop('blue_filter', False)

    wave = in_wave.to(ur.Angstrom).value

    low_check = low_wave.to(ur.Angstrom).value
    high_check = high_wave.to(ur.Angstrom).value


    red_filter = np.zeros_like(wave)
    red_filter[(wave<low_check) | (wave>high_check)] += rejection
    red_filter[(wave>=low_check) & (wave<=high_check)] += 1.0

    if blue_filter:
        red_filter[(in_wave < low_wave)] = 1e-6


    if diag:
        print('Low wavelength: {}'.format(low_wave))
        print('High wavelength: {}'.format(high_wave))
        print('Blue filter?: {}'.format(blue_filter))
        print('Rejection level: {}'.format(rejection))


    return red_filter


def optimize_filter(low_wave, high_wave, **kwargs):
    """
    Optimizes out-of-band filters based on the input bandpass

    ---
    Inputs:
    low_wave = Lower side of bandpass (units consistent with length)
    high_wave = High side of the bandpass

    Option inputs:

    qe_band = Which QE file to use (Default is "1")
    target_ratio = Out-of-band to in-band counts (0.5)
    blue_filter = If there's an asymmetric filter, apply this to everything blue-ward
        of the low_wave

    """
    from tdsat_telescope import load_qe, load_reflectivity
    from apply_transmission import apply_trans
    from zodi import load_zodi
    import astropy.units as ur

    # Check if the inputs make sense
    assert low_wave.unit.is_equivalent(ur.m), "Low-side wavelength does not have unit of length"
    assert high_wave.unit.is_equivalent(ur.m), "High-side wavelength does not have unit of length"

    qe_band = kwargs.pop('qe_band', 1)
    target_ratio = kwargs.pop('target_ratio', 0.5)
    blue_filter = kwargs.pop('blue_filter', False)
    diag = kwargs.pop('diag', False)
    # Load zodiacal background. Note that the out-of-band Zodi dominates over the
    # atmospheric lines (which are present here). Using the lowest Zodi background
    # represents the "worst case".
    zodi = load_zodi()

    # Load reflectivity and QE curves:
    ref_wave, reflectivity = load_reflectivity()
    qe_wave, qe = load_qe(band=qe_band)

    # Apply these to the Zodi spectrum:
    ref_flux = apply_trans(zodi['wavelength'], zodi['flux'], ref_wave, reflectivity/100.)
    qe_flux = apply_trans(zodi['wavelength'], ref_flux, qe_wave, qe)

    # Make a "standard" red filter:
    rejection = 1.0
    red_filter = make_red_filter(zodi['wavelength'], low_wave = low_wave,
                                 high_wave=high_wave, rejection = rejection,
                                 blue_filter=blue_filter, diag=diag)

    band_flux = apply_trans(zodi['wavelength'], qe_flux, zodi['wavelength'], red_filter)

    # Get the in-band, out-of-band ratio:

    in_band = band_flux[(zodi['wavelength'] > low_wave) &
               (zodi['wavelength']<high_wave)].sum()

    out_of_band = band_flux[((zodi['wavelength'] < low_wave) |
                (zodi['wavelength']>high_wave)) &
                (zodi['wavelength'] < 1*ur.micron)].sum()

    # Comput ratio:
    ratio = out_of_band / in_band
    target_rejection = (rejection * (target_ratio / ratio)).value

    if diag:
        print()
        print('Optimize filter diagnostics:')
        print('Low wave:{}'.format(low_wave))
        print('High wave:{}'.format(high_wave))
        print('Blue filter? {}'.format(blue_filter))
        print('Target ratio: {}'.format(target_ratio))
        print('Target rejection: {}'.format(target_rejection))
        print()
    return target_rejection