eic-opticks/ana/bnd.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 ::
0023
0024     In [13]: m0.data.shape
0025     Out[13]: (38, 2, 39, 4)   # 38 materials, 2 groups, 39 wavelengths, 4 qwns in each group : 38 mats include 2 added ones: GlassSchottF2, MainH2OHale -> 36 standard ones
0026
0027         ## huh m0.names shows that GlassSchottF2 appears in the middle, not at the tail ???  MainH2OHale is last
0028
0029     In [14]: b0.data.shape
0030     Out[14]: (123, 4, 2, 39, 4)  # 123 bnds, 4 matsur, 2 groups, 39 wavelengths, 4 qwns in each group
0031
0032     In [15]: s0.data.shape       # 49 surs include 4 added "perfect" ones that dont appear in bnds, so 49-4 = 44 standard ones
0033     Out[15]: (48, 2, 39, 4)
0034
0035     In [17]: len(b0.bnd.surs)   # includes 1 blank, so 44 surs occuring in bnd
0036     Out[17]: 45
0037
0038     In [18]: len(b0.bnd.mats)
0039     Out[18]: 36
0040
0041
0042     In [55]: np.where( np.logical_or( rel < -0.02, rel > 0.02 ))   ## 221/219168 off by more than 2%
0043     Out[55]:
0044     (array([ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  1,
0045             1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  2,  2,  2,  2,
0046             2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,
0047             3,  4,  4,  4,  4,  5,  5,  5,  5,  5,  5,  5,  5,  6,  6,  6,  6,  6,  6,  6,  6,  7,  7,  7,  7,  7,  7,  7,  7,  8,  8,  8,  8,  8,  8,  8,  8, 13, 13, 13, 13, 27, 27, 27, 27, 27, 27, 27,
0048            27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27]),
0049      array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0050            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0051            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0052            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),   ## all group zero, ie not GROUPVEL
0053      array([141, 142, 143, 144, 145, 145, 146, 146, 147, 147, 148, 148, 149, 149, 150, 150, 151, 152, 153, 154, 155, 288, 323, 324, 325, 326, 327, 328, 329, 330, 341, 342, 343, 344, 345, 346, 347, 348,
0054            349, 350, 351, 653, 654, 655, 658, 659, 141, 142, 143, 144, 145, 145, 146, 146, 147, 147, 148, 148, 149, 149, 150, 150, 151, 152, 153, 154, 155, 288, 323, 324, 325, 326, 327, 328, 329, 330,
0055            341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 653, 654, 655, 658, 659, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 141, 142, 143, 144, 145, 145, 146,
0056            146, 147, 147, 148, 148, 149, 149, 150, 150, 151, 152, 153, 154, 155, 145, 146, 147, 148, 149, 150, 224, 225, 226, 228, 229, 263, 264, 265, 268, 269, 308, 733, 734, 738, 739, 124, 125, 126,
0057            127, 128, 129, 130, 131, 124, 125, 126, 127, 128, 129, 130, 131, 124, 125, 126, 127, 128, 129, 130, 131, 124, 125, 126, 127, 128, 129, 130, 131, 733, 734, 738, 739, 121, 122, 123, 124, 125,
0058            126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 141, 142, 143, 144, 145, 145, 146, 146, 147, 147, 148, 148, 149, 149, 150, 150, 151, 152, 153, 154, 155]),
0059      array([1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1,
0060            1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 1, 1,
0061            2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
0062            1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 1, 1]))
0063                    ## never 0, refractive index, mostly absorption and scattering lengths
0064
0065
0066 Material property min/max relative interpolation differences with 20nm domain step::
0067
0068                                                   RINDEX                   ABSLEN                 RAYLEIGH                 REEMPROB                 GROUPVEL
0069  0                      GdDopedLS      -0.0048     0.0053       -0.0096     0.0821        0.0000     0.0237       -0.0423     0.0032       -0.0125     0.0065
0070  1             LiquidScintillator      -0.0048     0.0053       -0.0100     0.0821        0.0000     0.0237       -0.0423     0.0032       -0.0125     0.0065
0071  2                        Acrylic      -0.0046     0.0053        0.0000     0.0968        0.0000     0.0237        0.0000     0.0000       -0.0123     0.0064
0072  3                     MineralOil      -0.0046     0.0053       -0.0083     0.0232        0.0000     0.0237        0.0000     0.0000       -0.0123     0.0063
0073  4                       Bialkali       0.0000     0.0000       -0.0396     0.0017        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0074  5                       IwsWater      -0.0001     0.0000       -0.0084     0.0254        0.0000     0.0000        0.0000     0.0000       -0.0006     0.0005
0075  6                          Water      -0.0001     0.0000       -0.0084     0.0254        0.0000     0.0000        0.0000     0.0000       -0.0006     0.0005
0076  7                      DeadWater      -0.0001     0.0000       -0.0084     0.0254        0.0000     0.0000        0.0000     0.0000       -0.0006     0.0005
0077  8                       OwsWater      -0.0001     0.0000       -0.0084     0.0254        0.0000     0.0000        0.0000     0.0000       -0.0006     0.0005
0078  9                            ESR       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0079 10                   OpaqueVacuum       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0080 11                           Rock       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0081 12                         Vacuum       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0082 13                          Pyrex       0.0000     0.0000       -0.0396     0.0017        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0083 14                            Air       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0084 15                  GlassSchottF2       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0085 16                            PPE       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0086 17                      Aluminium       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0087 18          ADTableStainlessSteel       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0088 19                           Foam       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0089 20                       Nitrogen      -0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000       -0.0000     0.0000
0090 21                    NitrogenGas       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0091 22                          Nylon       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0092 23                            PVC       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0093 24                          Tyvek       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0094 25                       Bakelite       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0095 26                         MixGas       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0096 27                           Iron      -0.0046     0.0053        0.0000     0.0968        0.0000     0.0237        0.0000     0.0000       -0.0123     0.0064
0097 28                         Teflon       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0098 29             UnstStainlessSteel       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0099 30                            BPE       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0100 31                          Ge_68       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0101 32                          Co_60       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0102 33                           C_13       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0103 34                         Silver       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0104 35                        RadRock       0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000        0.0000     0.0000
0105
0106
0107 ::
0108
0109     In [21]: dict(filter(lambda kv:kv[1]<299,zip(i1m.names,i1m.data[:,1,430-60,0])))
0110     Out[21]:
0111     {'Acrylic': 192.77956,
0112      'Bialkali': 205.61897,
0113      'DeadWater': 217.83527,
0114      'GdDopedLS': 194.5192,
0115      'IwsWater': 217.83527,
0116      'LiquidScintillator': 194.5192,
0117      'MineralOil': 197.13411,
0118      'OwsWater': 217.83527,
0119      'Pyrex': 205.61897,
0120      'Teflon': 192.77956,
0121      'Water': 217.83527}
0122
0123
0124 """
0125 import os, logging, numpy as np
0126 from collections import OrderedDict as odict
0127 log = logging.getLogger(__name__)
0128
0129 from opticks.ana.base import opticks_main
0130 from opticks.ana.proplib import PropLib, Bnd
0131
0132 idp_ = lambda _:os.path.expandvars("$IDPATH/%s" % _ )
0133
0134
0135
0136 if __name__ == '__main__':
0137     ok = opticks_main()
0138
0139     # from old geocache without groupvel setup
0140     m0 = PropLib("GMaterialLib")
0141     s0 = PropLib("GSurfaceLib")
0142     b0 = PropLib("GBndLib")
0143
0144     # persisted postcache modified bnd buffer including groupvel calc, identity wavelengths 20nm steps
0145     b1 = PropLib("GBndLib", data="$TMP/InterpolationTest/CInterpolationTest_identity.npy" )
0146     b2 = PropLib("GBndLib", data="$TMP/InterpolationTest/OInterpolationTest_identity.npy" )
0147     assert np.allclose( b1.data, b2.data ) # after unset alignment
0148
0149     # persisted postcache groupvel calc, interpol wavelengths 1nm steps
0150     i1 = PropLib("GBndLib", data="$TMP/InterpolationTest/CInterpolationTest_interpol.npy" )
0151     i2 = PropLib("GBndLib", data="$TMP/InterpolationTest/OInterpolationTest_interpol.npy" )
0152
0153     assert np.allclose( i1.data[:,:,:,::20,:], i2.data[:,:,:,::20,:] ) == True   # plucking identity wavelengths from the interpol ones
0154     assert np.allclose( i1.data[:,:,:,::20,:], b1.data  )
0155     assert np.allclose( i2.data[:,:,:,::20,:], b2.data  )
0156
0157
0158     # collapse material duplication for simpler comparisons
0159     order = m0.names
0160     b0m = b0.as_mlib(order=order)
0161     b1m = b1.as_mlib(order=order)
0162     b2m = b2.as_mlib(order=order)
0163     i1m = i1.as_mlib(order=order)
0164     i2m = i2.as_mlib(order=order)
0165
0166     assert np.all( i1m.names == i2m.names )
0167     matnames = i1m.names
0168     assert np.all( b0m.data == m0.data[~np.logical_or(m0.names == 'GlassSchottF2', m0.names == 'MainH2OHale')] )
0169     assert np.all( b1m.data == b2m.data )   ## identity match
0170
0171     assert np.all( i1m.data[:,:,::20,:] == i2m.data[:,:,::20,:] )   ## on the raster, where there is no interpolation to do, get perfect match
0172
0173
0174     avg = (i1m.data + i2m.data)/2.0
0175     dif = i1m.data - i2m.data                 # absolute difference
0176
0177     # signed relative difference of 36*2*761*4 = 219168 property values
0178     rel = np.where( np.logical_or(avg < 1e-6, dif == 0), 0, dif/avg )
0179
0180
0181     rd = np.zeros( (len(rel), 2, 4, 2), dtype=np.float32 )
0182     rd[:,:,:,0] = np.amin(rel, axis=2)
0183     rd[:,:,:,1] = np.amax(rel, axis=2)
0184
0185
0186     print "".join(["%2s %30s " % ("","")] + map(lambda plab:"  %21s  " % plab, PropLib.M_LABELS ))
0187     for i in range(len(rd)):
0188         labl = ["%2d %30s " % ( i, matnames[i] )]
0189         rnge =  map(lambda mimx:"  %10.4f %10.4f  " % ( float(mimx[0]), float(mimx[1])) , rd[i].reshape(-1,2)[:5] )
0190         print "".join(labl + rnge)
0191
0192
0193
0194     Gd,LS,Ac,MO = 0,1,2,3
0195     gvel  = i1m.data[(Gd,Ac,LS,Ac,MO),1,430-60,0]
0196     gvel2 = i1m.data[(Gd,LS,Ac,MO,Ac),1,430-60,0]
0197     gvel3 = i1m.data[(Gd,LS,Ac,LS,Ac),1,430-60,0]
0198     dist = np.array([0,3000-5,3000+5,4000-5,4000+5,5000-5], dtype=np.float32)   # tconcentric radii
0199     ddif = np.diff(dist)
0200
0201     tdif = ddif/gvel
0202     tabs = np.cumsum(tdif) + 0.1
0203
0204     tdif2 = ddif/gvel2
0205     tabs2 = np.cumsum(tdif2) + 0.1
0206
0207     tdif3 = ddif/gvel3
0208     tabs3 = np.cumsum(tdif3) + 0.1
0209
0210
0211     print "gvel: %r " %  gvel
0212     print "dist: %r " %  dist
0213     print "ddif: %r " %  ddif
0214     print "tdif: %r " %  tdif
0215     print "tabs: %r " %  tabs
0216
0217     print "gvel2: %r " %  gvel2
0218     print "tdif2: %r " %  tdif2
0219     print "tabs2: %r " %  tabs2
0220
0221     print "gvel3: %r " %  gvel3
0222     print "tdif3: %r " %  tdif3
0223     print "tabs3: %r " %  tabs3
0224
0225
0226
0227
0228
0229