eic-opticks/ana/cf5d.py

0001 #
0002 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0003 #
0004 # This file is part of Opticks
0005 # (see https://bitbucket.org/simoncblyth/opticks).
0006 #
0007 # Licensed under the Apache License, Version 2.0 (the "License");
0008 # you may not use this file except in compliance with the License.
0009 # You may obtain a copy of the License at
0010 #
0011 #   http://www.apache.org/licenses/LICENSE-2.0
0012 #
0013 # Unless required by applicable law or agreed to in writing, software
0014 # distributed under the License is distributed on an "AS IS" BASIS,
0015 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0016 # See the License for the specific language governing permissions and
0017 # limitations under the License.
0018 #
0019
0020 """
0021
0022 5d comparison (ie staying in blib form) is tedious, m
0023 instead moved to trasforming the blib into mlib form
0024 and comparing that over in bnd.py
0025
0026
0027 """
0028
0029
0030 def cf5d(b0,i1,i2):
0031
0032     KL = odict()
0033     KL["RINDEX"]   = (0,PropLib.M_REFRACTIVE_INDEX)
0034     KL["ABSLEN"]   = (0,PropLib.M_ABSORPTION_LENGTH)
0035     KL["RAYLEIGH"] = (0,PropLib.M_SCATTERING_LENGTH)
0036     KL["REEMPROB"] = (0,PropLib.M_REEMISSION_PROB)
0037     KL["GROUPVEL"] = (1,PropLib.L_GROUP_VELOCITY)
0038
0039     shape = b0.data.shape
0040
0041     ni = shape[0]   # bnd
0042     nj = shape[1]   # omat/osur/isur/imat
0043     nk = shape[2]   # g0,g1
0044     nl = shape[3]   # samples = 39 for identity, or 761 for interpolating
0045     nm = shape[4]   # props
0046
0047     cf = np.zeros( (ni, nj, nk, nm, 2), dtype=np.float32  )
0048
0049     i12 = i1.data - i2.data
0050     #np.where( i12 > 300)      difficult to interpret the 5-tuple of indices
0051
0052     for i in range(ni):
0053         for j in range(nj):
0054             for k in range(nk):
0055
0056                 b1.dat.ijk = i,j,k
0057                 b2.dat.ijk = i,j,k
0058                 i1.dat.ijk = i,j,k
0059                 i2.dat.ijk = i,j,k
0060
0061                 assert np.allclose( b1.dat.d, b1.dat.d )
0062
0063                 for m in range(nm):
0064                     i12 = i1.dat.d[:,m] - i2.dat.d[:,m]
0065                     cf[i,j,k,m,0] = i12.min()
0066                     cf[i,j,k,m,1] = i12.max()
0067                 pass
0068             pass
0069         pass
0070     pass
0071
0072     # big discreps in all Water flavors  absorption length
0073     for i in range(ni):
0074         for j in range(nj):
0075             if cf[i,j].min() < -1 or cf[i,j].max() > 1:
0076                  print i, j, b0.names[i],"\n", cf[i,j]
0077
0078
0079
0080     cfm = odict()
0081
0082     for mat in b0.bnd.mats:
0083         oms = b0.bnd.oms(mat)   # indices of omat
0084         ims = b0.bnd.ims(mat)
0085
0086         print mat
0087         for om in oms:
0088             cfom = cf[om,0]
0089             if mat in cfm:
0090                 assert np.all(cfm[mat] == cfom)
0091             else:
0092                 cfm[mat] = cfom
0093             pass
0094             #if cfom.min() < -1 or cfom.max() > 1:
0095             #     print "om", b0.names[om],"\n", cfom
0096
0097         for im in ims:
0098
0099             cfim = cf[im,3]
0100             if mat in cfm:
0101                 assert np.all(cfm[mat] == cfim)
0102             else:
0103                 cfm[mat] = cfim
0104             pass
0105             #if cfim.min() < -1 or cfim.max() > 1:
0106             #     print "im", b0.names[im],"\n", cfim
0107
0108
0109     # interpol discreps seem big in absolute terms
0110     # but maybe not in relative
0111
0112 if 0:
0113
0114     for mat in cfm.keys():
0115         if cfm[mat].min() < -0.001 or cfm[mat].max() > 0.001:
0116             print mat,"\n", cfm[mat].reshape(-1,2).T
0117
0118             oms = b0.bnd.oms(mat)
0119             if len(oms) > 0:
0120                 i1.dat.ijk = oms[0],0,0
0121                 i2.dat.ijk = oms[0],0,0
0122                 for m in range(nm):
0123                     print m,"\n",(i2.dat.d[:,m] - i1.dat.d[:,m]) / (i2.dat.d[:,m] + i1.dat.d[:,m])/2.
0124
0125                 i1.dat.ijk = oms[0],0,1
0126                 i2.dat.ijk = oms[0],0,1
0127                 for m in range(nm):
0128                     print m,"\n",(i2.dat.d[:,m] - i1.dat.d[:,m]) / (i2.dat.d[:,m] + i1.dat.d[:,m])/2.
0129
0130
0131
0132
0133