eic-opticks/ana/prism_spectrum.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 prism_spectrum.py : compare ggv-newton evts with PrismExpected
0023 ==================================================================
0024
0025 Scatter plot of wavelength against the appropriate prism coordinate
0026 shows expected behaviour, except at the endpoints with outlier at high wavelength
0027 and a cliff drop off at low wavelenth.
0028 Perhaps these are due to restricted range of refractive index values or
0029 problem with spectral edges of the black body light source.
0030
0031 TODO: check undispersed spot, extend refractive index values
0032
0033 """
0034
0035 import os, sys, logging, numpy as np
0036 log = logging.getLogger(__name__)
0037
0038 import matplotlib.pyplot as plt
0039 from matplotlib.colors import LogNorm
0040 from mpl_toolkits.mplot3d import Axes3D
0041
0042 from opticks.ana.base import opticks_main
0043 from opticks.ana.evt import Evt
0044 from opticks.ana.boundary import Boundary
0045 from opticks.ana.tprism import Prism, PrismCheck, PrismExpected
0046 from opticks.ana.cie  import cie_hist1d, cie_hist2d
0047
0048 deg = np.pi/180.
0049
0050
0051 def spectrum_plot(w, d, db, ntile=50):
0052
0053     hRGB_raw,hXYZ_raw, bx = cie_hist1d(w, d, db, colorspace="sRGB/D65", norm=2)
0054
0055     hRGB_1 = np.clip(hRGB_raw, 0, 1)
0056
0057     hRGB_2 = np.tile(hRGB_1, ntile ).reshape(-1,ntile,3)
0058
0059     extent = [0,ntile,bx[-1],bx[0]]
0060
0061     ax.imshow(hRGB_2, origin="upper", extent=extent, alpha=1, vmin=0, vmax=1, aspect='auto')
0062
0063
0064 def deviation_plot(w, d, db, ntile=50):
0065
0066     h, hx = np.histogram(d, bins=db)
0067
0068     extent = [0,1,hx[-1],hx[0]]
0069
0070     ht = np.repeat(h,ntile).reshape(-1, ntile)
0071
0072     im = ax.matshow(ht, origin="upper", extent=extent, alpha=1, aspect='auto')
0073
0074     fig.colorbar(im)
0075
0076     return ht
0077
0078
0079 def uv_deviation_spike(d):
0080     """
0081     Deviation spike at 61 degrees
0082     coming from UV wavelengths all way from 60 to 340
0083
0084     This is due to the refractive index high plateau
0085     for GlassSchottF2 over that range...
0086
0087     This invalidates the range, so cannot believe anything
0088     below about 350nm
0089
0090     ::
0091
0092         In [57]: n[d>60].min()
0093         Out[57]: 1.6774575260175819
0094
0095         In [58]: n[d>60].max()
0096         Out[58]: 1.6846611499786377
0097
0098     """
0099     plt.hist(d, bins=100)
0100
0101
0102
0103
0104 if __name__ == '__main__':
0105     logging.basicConfig(level=logging.INFO)
0106     args = opticks_main(det="newton", tag="1", src="torch")
0107
0108     plt.ion()
0109
0110     try:
0111         sel = Evt(tag=args.tag, det=args.det, src=args.src, seqs=["TO BT BT SA"], args=args)  # newton, uses single incident angle
0112     except IOError as err:
0113         log.fatal(err)
0114         sys.exit(args.mrc)
0115
0116     log.info("loaded %s " % repr(sel))
0117
0118     boundary = Boundary("Vacuum///GlassSchottF2")
0119
0120     prism = Prism("60.,300,300,0", boundary)
0121
0122     n = boundary.imat.refractive_index(sel.wl)
0123
0124     xprism = PrismExpected(prism.a, n)
0125
0126     pc = PrismCheck(prism, xprism, sel )
0127
0128     w = pc.wx
0129     x = pc.p3[:,0]
0130     y = pc.p3[:,1]
0131     z = pc.p3[:,2]
0132     d = pc.dv/deg
0133
0134     off = x != 700
0135     #assert np.all(x == 700.)
0136
0137     fig = plt.figure()
0138
0139     bands = [400,800]
0140     #bands = [w.min(),w.max()]
0141
0142     nc = 2
0143     nr = len(bands) - 1
0144
0145     nb = 100
0146
0147
0148
0149     for ib in range(len(bands)-1):
0150
0151
0152         b = np.logical_and( w > bands[ib], w < bands[ib+1] )
0153
0154         bins = np.linspace(d[b].min(), d[b].max(), nb )
0155
0156
0157         ax= fig.add_subplot(nr,nc,ib*2+1)
0158
0159         spectrum_plot(w[b], d[b], bins)
0160
0161
0162         ax= fig.add_subplot(nr,nc,ib*2+2)
0163
0164         deviation_plot(w[b], d[b], bins)
0165
0166