Merge pull request #649 from HEXRD/absorption-correction

Absorption correction

.codecov.yml

@@ -0,0 +1,3 @@
+coverage:
+  status:
+    patch: false

conda.recipe/meta.yaml

@@ -35,6 +35,7 @@ requirements:
 
   run:
     - appdirs
+    - chemparse
     - fabio >=0.11
     # fabio uses lxml, and lxml 4.9.1 has some crashing issue on Windows
     # so, make sure lxml >= 4.9.2 is being used

hexrd/constants.py

@@ -310,16 +310,17 @@ def is_writable_file(path):
 
 
 # some physical constants
-cAvogadro = 6.02214076E23
-cBoltzmann = 1.380649E-23
-cCharge = 1.602176634E-19
-cJ2eV = 1.602176565E-19
-cLight = 299792458.0
-cMoment = 9.2740100707E-24
-cPermea = 1.2566370616E-6
-cPermit = 8.8541878163E-12
-cPlanck = 6.62607015E-34
-cRestmass = 9.1093837090E-31
+cAvogadro = 6.02214076E23  # Avogadro's constant Na
+cBoltzmann = 1.380649E-23  # Boltzmann's constant, K
+cCharge = 1.602176634E-19  # charge of electron
+cJ2eV = 1.602176565E-19  # joule to ev, JperkeV*1e-3
+cLight = 299792458.0  # speed of light, same as c above but name is more descriptive
+cMoment = 9.2740100707E-24  # magnetic moment of electron
+cPermea = 1.2566370616E-6  # permeability of free space
+cPermit = 8.8541878163E-12  # permittivity of free space
+cPlanck = 6.62607015E-34  # same as h above but name is more descriptive
+cRestmass = 9.1093837090E-31  # rest mass of electron
+cClassicalelectronRad = 2.8179403e-6  # classical electron radius in nm
 
 '''
 adding another parametrization of the
@@ -708,6 +709,80 @@ def is_writable_file(path):
                          231.03588, 238.02891, 237.0, 244.0, 243.0,
                          247.0, 247.0, 251.0])
 
+"""
+dictionary of atomic weights
+"""
+ATOM_WEIGHTS_DICT = {'H': 1.00794,'He': 4.002602,'Li': 6.941,
+                'Be': 9.012182,'B': 10.811,'C': 12.0107,
+                'N': 14.0067,'O': 15.9994,'F': 18.9984032,
+                'Ne': 20.1797,'Na': 22.98976928,'Mg': 24.305,
+                'Al': 26.9815386,'Si': 28.0855,'P': 30.973762,
+                'S': 32.065,'Cl': 35.453,'Ar': 39.948,
+                'K': 39.0983,'Ca': 40.078,'Sc': 44.955912,
+                'Ti': 47.867,'V': 50.9415,'Cr': 51.9961,
+                'Mn': 54.938045,'Fe': 55.845,'Co': 58.933195,
+                'Ni': 58.6934,'Cu': 63.546,'Zn': 65.38,
+                'Ga': 69.723,'Ge': 72.64,'As': 74.9216,
+                'Se': 78.96,'Br': 79.904,'Kr': 83.798,
+                'Rb': 85.4678,'Sr': 87.62,'Y': 88.90585,
+                'Zr': 91.224,'Nb': 92.90638,'Mo': 95.96,
+                'Tc': 98.9062,'Ru': 101.07,'Rh': 102.9055,
+                'Pd': 106.42,'Ag': 107.8682,'Cd': 112.411,
+                'In': 114.818,'Sn': 118.71,'Sb': 121.76,
+                'Te': 127.6,'I': 126.90447,'Xe': 131.293,
+                'Cs': 132.9054519,'Ba': 137.327,'La': 138.90547,
+                'Ce': 140.116,'Pr': 140.90765,'Nd': 144.242,
+                'Pm': 145.0,'Sm': 150.36,'Eu': 151.964,
+                'Gd': 157.25,'Tb': 158.92535,'Dy': 162.5,
+                'Ho': 164.93032,'Er': 167.259,'Tm': 168.93421,
+                'Yb': 173.054,'Lu': 174.9668,'Hf': 178.49,
+                'Ta': 180.94788,'W': 183.84,'Re': 186.207,
+                'Os': 190.23,'Ir': 192.217,'Pt': 195.084,
+                'Au': 196.966569,'Hg': 200.59,'Tl': 204.3833,
+                'Pb': 207.2,'Bi': 208.9804,'Po': 209.0,
+                'At': 210.0,'Rn': 222.0,'Fr': 223.0,'Ra': 226.0,
+                'Ac': 227.0,'Th': 232.0377,'Pa': 231.03588,
+                'U': 238.02891,'Np': 237.0,'Pu': 244.0,
+                'Am': 243.0,'Cm': 247.0,'Bk': 247.0,
+                'Cf': 251.0}
+
+"""
+densities of elements in g/cc
+"""
+DENSITY = {'H': 8.99e-05,'He': 0.0001785,'Li': 0.535,'Be': 1.848,
+           'B': 2.46,'C': 2.26,'N': 0.001251,'O': 0.001429,
+           'F': 0.001696,'Ne': 0.0009,'Na': 0.968,'Mg': 1.738,
+           'Al': 2.7,'Si': 2.33,'P': 1.823,'S': 1.96,'Cl': 0.003214,
+           'Ar': 0.001784,'K': 0.856,'Ca': 1.55,'Sc': 2.985,
+           'Ti': 4.507,'V': 6.11,'Cr': 7.14,'Mn': 7.47,'Fe': 7.874,
+           'Co': 8.9,'Ni': 8.908,'Cu': 8.92,'Zn': 7.14,'Ga': 5.904,
+           'Ge': 5.323,'As': 5.727,'Se': 4.819,'Br': 3.12,'Kr': 0.00375,
+           'Rb': 1.532,'Sr': 2.63,'Y': 4.472,'Zr': 6.511,
+           'Nb': 8.57,'Mo': 10.28,'Tc': 11.5,'Ru': 12.37,
+           'Rh': 12.45,'Pd': 12.023,'Ag': 10.49,'Cd': 8.65,
+           'In': 7.31,'Sn': 7.31,'Sb': 6.697,'Te': 6.24,
+           'I': 4.94,'Xe': 0.0059,'Cs': 1.879,'Ba': 3.51,
+           'La': 6.146,'Ce': 6.689,'Pr': 6.64,'Nd': 7.01,
+           'Pm': 7.264,'Sm': 7.353,'Eu': 5.244,'Gd': 7.901,
+           'Tb': 8.219,'Dy': 8.551,'Ho': 8.795,'Er': 9.066,
+           'Tm': 9.321,'Yb': 6.57,'Lu': 9.841,'Hf': 13.31,
+           'Ta': 16.65,'W': 19.25,'Re': 21.02,'Os': 22.59,
+           'Ir': 22.56,'Pt': 21.09,'Au': 19.3,'Hg': 13.534,
+           'Tl': 11.85,'Pb': 11.34,'Bi': 9.78,'Po': 9.196,
+           'At': None,'Rn': 0.00973,'Fr': None,'Ra': 5.0,
+           'Ac': 10.07,'Th': 11.724,'Pa': 15.37,'U': 19.05,
+           'Np': 20.45,'Pu': 19.816,'Am': 13.67,'Cm': 13.51,
+           'Bk': 14.78,'Cf': 15.1}
+
+# some polymer densities commonly used in hexrd
+DENSITY_COMPOUNDS = {
+    'C10H8O4': 1.4,
+    'Ba2263F2263Br1923I339C741H1730N247O494': 3.3,
+    'LiF': 2.64,
+    'quartz': 2.65,
+    'diamond': 3.5,
+}
+
 '''
 dictionary of atomic numbers with element symbol as keys
 used in I/O from cif file

hexrd/instrument/constants.py

@@ -0,0 +1,64 @@
+from hexrd.constants import DENSITY, DENSITY_COMPOUNDS
+
+
+# default filter and coating materials
+class FILTER_DEFAULTS:
+    TARDIS = {
+        'material': 'Ge',
+        'density' : DENSITY['Ge'],
+        'thickness' : 10 # microns
+    }
+    PXRDIP = {
+        'material': 'Cu',
+        'density' : DENSITY['Cu'],
+        'thickness' : 10 # microns
+    }
+
+COATING_DEFAULT = {
+    'material': 'C10H8O4',
+    'density': DENSITY_COMPOUNDS['C10H8O4'],
+    'thickness': 9 # microns
+}
+
+PHOSPHOR_DEFAULT = {
+    'material': 'Ba2263F2263Br1923I339C741H1730N247O494',
+    'density': DENSITY_COMPOUNDS['Ba2263F2263Br1923I339C741H1730N247O494'], # g/cc
+    'thickness': 115, # microns
+    'readout_length': 222, #microns
+    'pre_U0': 0.695
+}
+
+
+class PHYSICS_PACKAGE_DEFAULTS:
+    # Default physics package for dynamic compression
+    HED = {
+        'sample_material': 'Fe',
+        'sample_density': DENSITY['Fe'],
+        'sample_thickness': 15, # in microns
+        'window_material': 'LiF',
+        'window_density': DENSITY_COMPOUNDS['LiF'],
+        'window_thickness': 150, # in microns
+    }
+    # # Template for HEDM type physics package
+    # HEDM = {
+    #     'sample_material': 'Fe',
+    #     'sample_density': DENSITY['Fe'],
+    #     'sample_thickness': 1000, # in microns
+    #     'sample_geometry': 'cylinder'
+    # }
+
+
+# Default pinhole area correction parameters
+class PINHOLE_DEFAULTS:
+    TARDIS = {
+        'pinhole_material' : 'Ta',
+        'pinhole_diameter' : 400, # in microns
+        'pinhole_thickness' : 100, # in microns
+        'pinhole_density' : 16.65, # g/cc
+    }
+    PXRDIP = {
+        'pinhole_material' : 'Ta',
+        'pinhole_diameter' : 130, # in microns
+        'pinhole_thickness' : 70, # in microns
+        'pinhole_density' : 16.65, # g/cc
+    }

hexrd/instrument/cylindrical_detector.py

@@ -120,6 +120,74 @@ def cart_to_dvecs(self,
     def pixel_angles(self, origin=ct.zeros_3):
         return _pixel_angles(origin=origin, **self._pixel_angle_kwargs)
 
+    def local_normal(self):
+        """get the local normal of each pixel in the 
+        cylindrical detector
+
+        output will be of shape [nx*ny, 3]
+        """
+        y, x = self.pixel_coords
+        x = x.flatten()
+        y = y.flatten()
+        num = x.shape[0]
+        naxis = np.cross(self.paxis, self.caxis)
+
+        th = x/self.radius
+        xp = np.sin(th)
+        xn = -np.cos(th)
+
+        pcomp = np.tile(xp, [3, 1]).T * np.tile(self.paxis, [num, 1])
+        ncomp = np.tile(xn, [3, 1]).T * np.tile(naxis, [num, 1])
+
+        cart3d = pcomp + ncomp
+        return cart3d
+
+    def calc_filter_coating_transmission(self, energy):
+        """
+        calculate the transmission after x-ray beam interacts
+        with the filter and the mylar polymer coating.
+        Specifications of the polymer coating is taken from:
+
+        M. Stoeckl, A. A. Solodov
+        Readout models for BaFBr0.85I0.15:Eu image plates 
+        Rev. Sci. Instrum. 89, 063101 (2018)
+
+        Transmission Formulas are consistent with:
+
+        Rygg et al., X-ray diffraction at the National 
+        Ignition Facility, Rev. Sci. Instrum. 91, 043902 (2020)
+        """
+        al_f = self.filter.absorption_length(energy)
+        al_c = self.coating.absorption_length(energy)
+        al_p = self.phosphor.energy_absorption_length(energy)
+
+        t_f = self.filter.thickness
+        t_c = self.coating.thickness
+        t_p = self.phosphor.thickness
+        L   = self.phosphor.readout_length
+
+        det_normal = self.local_normal()
+
+        y, x = self.pixel_coords
+        xy_data = np.vstack((x.flatten(), y.flatten())).T
+        dvecs = self.cart_to_dvecs(xy_data)
+        dvecs = dvecs/np.tile(np.linalg.norm(dvecs, axis=1), [3, 1]).T
+
+        secb = (1./np.sum(dvecs*det_normal, axis=1)).reshape(self.shape)
+
+        transmission_filter  = self.calc_transmission_generic(secb, t_f, al_f)
+        transmission_coating = self.calc_transmission_generic(secb, t_c, al_c)
+        transmission_phosphor = (
+            self.phosphor.pre_U0 *
+            self.calc_transmission_phosphor(secb, t_p, al_p, L, energy))
+
+        transmission_filter  = transmission_filter.reshape(self.shape)
+        transmission_coating = transmission_coating.reshape(self.shape)
+        transmission_filter_coating = (
+            transmission_filter * transmission_coating)
+
+        return transmission_filter_coating, transmission_phosphor
+
     @property
     def _pixel_angle_kwargs(self):
         # kwargs used for pixel angles, pixel_tth_gradient,