diff --git a/demo/xrayphysics_example.py b/demo/xrayphysics_example.py index bde7c71..a7c9fdf 100644 --- a/demo/xrayphysics_example.py +++ b/demo/xrayphysics_example.py @@ -8,21 +8,22 @@ # The next line prints out the version number and some other information about this package #physics.about() -whichPlot = 2 +whichPlot = 10 ######################################################################################### # Example 1: Getting linear attenuation coefficients for a compound ######################################################################################### # Set a numpy array of the x-ray energies -gammas = np.array(range(20,100),dtype=np.float32)+1.0 +gammas = np.array(range(20,100))+1.0 # Set the LAC (cm^-1) of water and its components chemForm = 'H2O' -mu_water = physics.mu(chemForm,gammas,1.0) -mu_water_PE = physics.muPE(chemForm,gammas,1.0) -mu_water_CS = physics.muCS(chemForm,gammas,1.0) -mu_water_RS = physics.muRS(chemForm,gammas,1.0) +density = 1.0 +mu_water = physics.mu(chemForm,gammas,density) +mu_water_PE = physics.muPE(chemForm,gammas,density) +mu_water_CS = physics.muCS(chemForm,gammas,density) +mu_water_RS = physics.muRS(chemForm,gammas,density) # Plot results if whichPlot == 1: @@ -38,7 +39,7 @@ # Example 2: Calculate cross section of bone by mass fractions ######################################################################################### if whichPlot == 2: - gammas = np.array(range(20,100),dtype=np.float32)+1.0 + gammas = np.array(range(20,100))+1.0 # Mass Density of Bone (g/cm^3) massDensity = 1.85 @@ -178,3 +179,25 @@ plt.ylabel('normalized differential cross section') plt.xlim((0.0,180.0)) plt.show() + + +######################################################################################### +# Example 10: Material Library +######################################################################################### +# All functions demonstrated above that take a chemical formula as an input can also take +# the material name of a material in the material list that is defined at the end of xrayphysics.py +# One can retrieve the chemical formula and mass density by material names +# To see the whole list of items in the material library, use the following: +#physics.print_material_library() + +# Now let's demonstrate how this works +if whichPlot == 10: + materialName = 'lung' + #materialName = 'bone_compact' + gammas = np.array(range(20,100))+1.0 + + plt.plot(gammas, physics.mu(materialName, gammas)) + plt.title(str(materialName) + ' LAC') + plt.xlabel('x-ray energy (keV)') + plt.ylabel('linear attenuation coefficient (cm^-1)') + plt.show() diff --git a/src/xrayphysics.py b/src/xrayphysics.py index ba51168..e948584 100644 --- a/src/xrayphysics.py +++ b/src/xrayphysics.py @@ -154,17 +154,45 @@ def fixArray(self, x): x = np.ascontiguousarray(x, dtype=np.float32) return x + def elementSymbolToAtomicNumber(self, elementStr): + if sys.version_info[0] == 3: + elementStr = bytes(str(elementStr), 'ascii') + self.libxrayphysics.elementSymbolToAtomicNumber.restype = ctypes.c_int + self.libxrayphysics.elementSymbolToAtomicNumber.argtypes = [ctypes.c_char_p] + return self.libxrayphysics.elementSymbolToAtomicNumber(elementStr) + def atomicMass(self, Z): """Returns the atomic mass of the given atomic number""" self.libxrayphysics.atomicMass.restype = ctypes.c_float self.libxrayphysics.atomicMass.argtypes = [ctypes.c_int] + if isinstance(Z, str): + Z = self.elementSymbolToAtomicNumber(Z) return self.libxrayphysics.atomicMass(Z) + + def density(self, Z): + """Returns the mass density of the given atomic number or element""" + return self.massDensity(Z) + + def massDensity(self, Z): + """Returns the mass density of the given atomic number or element""" + if isinstance(Z, str): + atomicNumber = self.elementSymbolToAtomicNumber(Z) + if atomicNumber < 1: + if Z in materialDensities.keys(): + return materialDensities[Z] * self.density_scalar() + else: + return 0.0 + else: + Z = atomicNumber + self.libxrayphysics.massDensity.restype = ctypes.c_float + self.libxrayphysics.massDensity.argtypes = [ctypes.c_int] + return self.libxrayphysics.massDensity(Z) * self.density_scalar() - def mu(self, Z, gamma, massDensity): + def mu(self, Z, gamma, massDensity=None): """Returns the Linear Attenuation Coefficient (LAC) of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the LAC massDensity (scalar): the mass density (g/cm^3) of the material @@ -172,13 +200,15 @@ def mu(self, Z, gamma, massDensity): The LAC (cm^-1) of the material at the specified energy or energies """ + if massDensity is None: + massDensity = self.massDensity(Z) return massDensity * self.sigma(Z, gamma) - def muPE(self, Z, gamma, massDensity): + def muPE(self, Z, gamma, massDensity=None): """Returns the Photoelectric component of the Linear Attenuation Coefficient (LAC) of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the LAC massDensity (scalar): the mass density (g/cm^3) of the material @@ -186,13 +216,15 @@ def muPE(self, Z, gamma, massDensity): The Photoelectric component of the LAC (cm^-1) of the material at the specified energy or energies """ + if massDensity is None: + massDensity = self.massDensity(Z) return massDensity * self.sigmaPE(Z, gamma) - def muCS(self, Z, gamma, massDensity): + def muCS(self, Z, gamma, massDensity=None): """Returns the Compton Scatter component of the Linear Attenuation Coefficient (LAC) of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the LAC massDensity (scalar): the mass density (g/cm^3) of the material @@ -200,13 +232,15 @@ def muCS(self, Z, gamma, massDensity): The Compton Scatter component of the LAC (cm^-1) of the material at the specified energy or energies """ + if massDensity is None: + massDensity = self.massDensity(Z) return massDensity * self.sigmaCS(Z, gamma) - def muRS(self, Z, gamma, massDensity): + def muRS(self, Z, gamma, massDensity=None): """Returns the Rayleigh Scatter component of the Linear Attenuation Coefficient (LAC) of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the LAC massDensity (scalar): the mass density (g/cm^3) of the material @@ -214,13 +248,15 @@ def muRS(self, Z, gamma, massDensity): The Rayleigh Scatter component of the LAC (cm^-1) of the material at the specified energy or energies """ + if massDensity is None: + massDensity = self.massDensity(Z) return massDensity * self.sigmaRS(Z, gamma) - def muPP(self, Z, gamma, massDensity): + def muPP(self, Z, gamma, massDensity=None): """Returns the Pair Production component of the Linear Attenuation Coefficient (LAC) of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the LAC massDensity (scalar): the mass density (g/cm^3) of the material @@ -228,13 +264,15 @@ def muPP(self, Z, gamma, massDensity): The Pair Production component of the LAC (cm^-1) of the material at the specified energy or energies """ + if massDensity is None: + massDensity = self.massDensity(Z) return massDensity * self.sigmaPP(Z, gamma) - def muTP(self, Z, gamma, massDensity): + def muTP(self, Z, gamma, massDensity=None): """Returns the Triplet Production component of the Linear Attenuation Coefficient (LAC) of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the LAC massDensity (scalar): the mass density (g/cm^3) of the material @@ -242,13 +280,15 @@ def muTP(self, Z, gamma, massDensity): The Triplet Production component of the LAC (cm^-1) of the material at the specified energy or energies """ + if massDensity is None: + massDensity = self.massDensity(Z) return massDensity * self.sigmaTP(Z, gamma) def sigma(self, Z, gamma): """Returns the mass cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -257,8 +297,11 @@ def sigma(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaCompound.restype = ctypes.c_float self.libxrayphysics.sigmaCompound.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -283,7 +326,7 @@ def sigmaPE(self, Z, gamma): """Returns the Photoelectric component of the mass cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -292,8 +335,11 @@ def sigmaPE(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaCompoundPE.restype = ctypes.c_float self.libxrayphysics.sigmaCompoundPE.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -318,7 +364,7 @@ def sigmaCS(self, Z, gamma): """Returns the Compton Scatter component of the mass cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -327,8 +373,11 @@ def sigmaCS(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaCompoundCS.restype = ctypes.c_float self.libxrayphysics.sigmaCompoundCS.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -353,7 +402,7 @@ def sigmaRS(self, Z, gamma): """Returns the Rayleigh Scatter component of the mass cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -362,8 +411,11 @@ def sigmaRS(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaCompoundRS.restype = ctypes.c_float self.libxrayphysics.sigmaCompoundRS.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -388,7 +440,7 @@ def sigmaPP(self, Z, gamma): """Returns the Pair Production component of the mass cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -397,8 +449,11 @@ def sigmaPP(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaCompoundPP.restype = ctypes.c_float self.libxrayphysics.sigmaCompoundPP.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -423,7 +478,7 @@ def sigmaTP(self, Z, gamma): """Returns the Triplet Production component of the mass cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -432,8 +487,11 @@ def sigmaTP(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaCompoundTP.restype = ctypes.c_float self.libxrayphysics.sigmaCompoundTP.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -486,7 +544,7 @@ def sigma_e(self, Z, gamma): """Returns the electron cross section of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar or numpy array): the energy or energies (keV) for which to calculate the mass cross section Returns: @@ -495,8 +553,11 @@ def sigma_e(self, Z, gamma): """ if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.sigmaeCompound.restype = ctypes.c_float self.libxrayphysics.sigmaeCompound.argtypes = [ctypes.c_char_p, ctypes.c_float] if type(gamma) is np.ndarray: @@ -529,7 +590,7 @@ def incoherentScatterDistribution(self, Z, gamma, theta, doNormalize=False): """Returns the Incoherent Scatter distribution of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar): the energy (keV) for which to calculate the Scatter distribution theta (scalar or numpy array): the angle(s) (degrees) for which to calculate the Scatter distribution doNormalize (bool): if True, normalizes the values so that they can express as a probability density function @@ -541,8 +602,11 @@ def incoherentScatterDistribution(self, Z, gamma, theta, doNormalize=False): #float incoherentScatterDistribution(int Z, float gamma, float theta) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.incoherentScatterDistributionCompound.restype = ctypes.c_float self.libxrayphysics.incoherentScatterDistributionCompound.argtypes = [ctypes.c_char_p, ctypes.c_float, ctypes.c_float] if type(theta) is np.ndarray: @@ -577,7 +641,7 @@ def coherentScatterDistribution(self, Z, gamma, theta, doNormalize=False): """Returns the Coherent Scatter distribution of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar): the energy (keV) for which to calculate the Scatter distribution theta (scalar or numpy array): the angle(s) (degrees) for which to calculate the Scatter distribution doNormalize (bool): if True, normalizes the values so that they can express as a probability density function @@ -589,8 +653,11 @@ def coherentScatterDistribution(self, Z, gamma, theta, doNormalize=False): #float coherentScatterDistribution(int Z, float gamma, float theta) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.coherentScatterDistributionCompound.restype = ctypes.c_float self.libxrayphysics.coherentScatterDistributionCompound.argtypes = [ctypes.c_char_p, ctypes.c_float, ctypes.c_float] if type(theta) is np.ndarray: @@ -626,7 +693,7 @@ def incoherentScatterDistribution_normalizationFactor(self, Z, gamma): """Returns the Incoherent Scatter distribution normalization factor of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar): the energy (keV) for which to calculate the normalization factor Returns: @@ -636,8 +703,11 @@ def incoherentScatterDistribution_normalizationFactor(self, Z, gamma): #float incoherentScatterDistribution(int Z, float gamma, float theta) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.incoherentScatterDistributionCompound_normalizationFactor.restype = ctypes.c_float self.libxrayphysics.incoherentScatterDistributionCompound_normalizationFactor.argtypes = [ctypes.c_char_p, ctypes.c_float] return self.libxrayphysics.incoherentScatterDistributionCompound_normalizationFactor(chemicalFormula, gamma) * self.cross_section_scalar() @@ -658,7 +728,7 @@ def coherentScatterDistribution_normalizationFactor(self, Z, gamma): """Returns the Coherent Scatter distribution normalization factor of the given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library gamma (scalar): the energy (keV) for which to calculate the normalization factor Returns: @@ -668,8 +738,11 @@ def coherentScatterDistribution_normalizationFactor(self, Z, gamma): #float coherentScatterDistribution(int Z, float gamma, float theta) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.coherentScatterDistributionCompound_normalizationFactor.restype = ctypes.c_float self.libxrayphysics.coherentScatterDistributionCompound_normalizationFactor.argtypes = [ctypes.c_char_p, ctypes.c_float] return self.libxrayphysics.coherentScatterDistributionCompound_normalizationFactor(chemicalFormula, gamma) * self.cross_section_scalar() @@ -702,6 +775,8 @@ def simulateSpectra(self, kV, takeOffAngle=11.0, Z=74, gammas=None): N = int(np.ceil(float(maxEnergy - minEnergy) / float(T_E))) minEnergy = maxEnergy - N*T_E gammas = np.ascontiguousarray(np.array(range(minEnergy,maxEnergy+1,T_E)), dtype=np.float32) + else: + gammas = self.fixArray(gammas) if Z is None: Z = 74 @@ -728,6 +803,8 @@ def changeTakeOffAngle(self, kV, takeOffAngle_cur, takeOffAngle_new, Z, gammas, """ #bool changeTakeOffAngle(float kV, float takeOffAngle_cur, float takeOffAngle_new, int Z, float* gammas, int N, float* s) + gammas = self.fixArray(gammas) + spectrum_cur = self.fixArray(spectrum_cur) self.libxrayphysics.changeTakeOffAngle.restype = ctypes.c_bool self.libxrayphysics.changeTakeOffAngle.argtypes = [ctypes.c_float, ctypes.c_float, ctypes.c_float, ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS")] spectrum_new = spectrum_cur.copy() @@ -738,7 +815,7 @@ def detectorResponse(self, chemicalFormula, density, thickness, gammas): """Analytic detector response model Args: - chemicalFormula (string): chemical formula, or mixture of compounds with mass fractions of the scintillator material + chemicalFormula (string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library density (scalar): the mass density (g/cm^3) of the scintillator thickness (scalar): the thickness (cm) of the scintillator gammas (numpy array): the energy samples (keV) at which to calculate the detector response @@ -746,6 +823,7 @@ def detectorResponse(self, chemicalFormula, density, thickness, gammas): Returns: A model of the detector response """ + gammas = self.fixArray(gammas) detResp = np.ascontiguousarray(np.zeros(gammas.size), dtype=np.float32) for n in range(gammas.size): detResp[n] = gammas[n]*(1.0-np.exp(-self.mu(chemicalFormula, gammas[n], density)*thickness)) @@ -755,7 +833,7 @@ def filterResponse(self, chemicalFormula, density, thickness, gammas): """X-ray Filter response model Args: - chemicalFormula (string): chemical formula, or mixture of compounds with mass fractions of the filter material + chemicalFormula (string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library density (scalar): the mass density (g/cm^3) of the filter material thickness (scalar): the thickness (cm) of the filter material gammas (numpy array): the energy samples (keV) at which to calculate the filter response @@ -764,6 +842,7 @@ def filterResponse(self, chemicalFormula, density, thickness, gammas): A model of the response to the given x-ray filter """ + gammas = self.fixArray(gammas) filtResp = np.ascontiguousarray(np.zeros(gammas.size), dtype=np.float32) for n in range(gammas.size): filtResp[n] = np.exp(-self.mu(chemicalFormula, gammas[n], density)*thickness) @@ -780,6 +859,8 @@ def meanEnergy(self, spectralResponse, gammas): The mean energy (keV) of the given spectra """ + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) self.libxrayphysics.meanEnergy.restype = ctypes.c_float self.libxrayphysics.meanEnergy.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int] return self.libxrayphysics.meanEnergy(spectralResponse, gammas, gammas.size) @@ -795,6 +876,8 @@ def normalizeSpectrum(self, spectralResponse, gammas): The normalized spectra """ + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) self.libxrayphysics.normalizeSpectrum.restype = ctypes.c_bool self.libxrayphysics.normalizeSpectrum.argtypes = [ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int] return self.libxrayphysics.normalizeSpectrum(spectralResponse, gammas, gammas.size) @@ -803,7 +886,7 @@ def effectiveZ(self, chemicalFormula, min_energy=10.0, max_energy=100.0, arealDe """Calculate the effective atomic number of a material Args: - chemicalFormula (string): chemical formula, or mixture of compounds with mass fractions of the material + chemicalFormula (string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library min_energy (scalar): the minimum of the energy in which to perform the calculation max_energy (scalar): the maximum of the energy in which to perform the calculation arealDensity (scalar): the areal density (g/cm^2) to use for the calculation @@ -813,8 +896,11 @@ def effectiveZ(self, chemicalFormula, min_energy=10.0, max_energy=100.0, arealDe """ if isinstance(chemicalFormula, str): + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + #float effectiveZ(const char* chemForm, float min_energy, float max_energy, float arealDensity); self.libxrayphysics.effectiveZ.restype = ctypes.c_float self.libxrayphysics.effectiveZ.argtypes = [ctypes.c_char_p, ctypes.c_float, ctypes.c_float, ctypes.c_float] @@ -827,7 +913,7 @@ def effectiveAttenuation(self, Z, density, thickness, spectralResponse, gammas): """Calculate the effective attenuation of a material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library density (scalar): mass density (g/cm^3) of the material thickness (scalar): the thickness (cm) of the material spectralResponse (numpy array): spectra model @@ -838,10 +924,15 @@ def effectiveAttenuation(self, Z, density, thickness, spectralResponse, gammas): """ #float effectiveAttenuation(float Z, float density, float thickness, float* spectralResponse, float* gammas, int N); + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.effectiveAttenuation_compound.restype = ctypes.c_float self.libxrayphysics.effectiveAttenuation_compound.argtypes = [ctypes.c_char_p, ctypes.c_float, ctypes.c_float, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int] return self.libxrayphysics.effectiveAttenuation_compound(chemicalFormula, density/self.density_scalar(), thickness/self.thickness_scalar(), spectralResponse, gammas, gammas.size) @@ -854,7 +945,7 @@ def effectiveEnergy(self, Z, density, thickness, spectralResponse, gammas): """Calculate the effective energy of a material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library density (scalar): mass density (g/cm^3) of the material thickness (scalar): the thickness (cm) of the material spectralResponse (numpy array): spectra model @@ -865,10 +956,15 @@ def effectiveEnergy(self, Z, density, thickness, spectralResponse, gammas): """ #float effectiveEnergy(float Z, float density, float thickness, float* spectralResponse, float* gammas, int N); + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.effectiveEnergy_compound.restype = ctypes.c_float self.libxrayphysics.effectiveEnergy_compound.argtypes = [ctypes.c_char_p, ctypes.c_float, ctypes.c_float, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int] return self.libxrayphysics.effectiveEnergy_compound(chemicalFormula, density/self.density_scalar(), thickness/self.thickness_scalar(), spectralResponse, gammas, gammas.size) @@ -881,7 +977,7 @@ def transmission(self, Z, density, thickness, spectralResponse, gammas): """Calculate the transmission through a given material Args: - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library density (scalar): mass density (g/cm^3) of the material thickness (scalar): the thickness (cm) of the material spectralResponse (numpy array): spectra model @@ -892,10 +988,15 @@ def transmission(self, Z, density, thickness, spectralResponse, gammas): """ #float transmission(float Z, float density, float thickness, float* spectralResponse, float* gammas, int N); + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.transmission_compound.restype = ctypes.c_float self.libxrayphysics.transmission_compound.argtypes = [ctypes.c_char_p, ctypes.c_float, ctypes.c_float, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int] return self.libxrayphysics.transmission_compound(chemicalFormula, density/self.density_scalar(), thickness/self.thickness_scalar(), spectralResponse, gammas, gammas.size) @@ -910,7 +1011,7 @@ def setBHlookupTable(self, spectralResponse, gammas, Z, referenceEnergy=0.0, T_a Args: spectralResponse (numpy array): spectra model gammas (numpy array): energies at which the spectra model is defined - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library referenceEnergy (scalar): the energy (keV) of the monochromatic attenuation (if not specified uses the mean energy of the spectra) T_atten (scalar): the sampling rate of the monochromatic attenuation N_atten (int): the number of samples of the monochromatic attenuation @@ -929,6 +1030,8 @@ def setBHlookupTable(self, spectralResponse, gammas, Z, referenceEnergy=0.0, T_a N_atten = int(np.ceil(max_lac / T_atten)) + 1 max_lac = float(N_atten-1)*T_atten + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) N_gamma = gammas.size if referenceEnergy <= 0.0: #referenceEnergy = self.meanEnergy(spectralResponse, gammas) @@ -937,8 +1040,11 @@ def setBHlookupTable(self, spectralResponse, gammas, Z, referenceEnergy=0.0, T_a LUT = np.zeros(N_atten, dtype=np.float32) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.setBHlookupTable_compound.restype = ctypes.c_bool self.libxrayphysics.setBHlookupTable_compound.argtypes = [ctypes.c_char_p, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_float, ctypes.c_int, ctypes.c_float] self.libxrayphysics.setBHlookupTable_compound(chemicalFormula, spectralResponse, gammas, N_gamma, LUT, T_atten, N_atten, referenceEnergy) @@ -954,7 +1060,7 @@ def setBHClookupTable(self, spectralResponse, gammas, Z, referenceEnergy=0.0, T_ Args: spectralResponse (numpy array): spectra model gammas (numpy array): energies at which the spectra model is defined - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library referenceEnergy (scalar): the energy (keV) of the monochromatic attenuation (if not specified uses the mean energy of the spectra) T_atten (scalar): the sampling rate of the polychromatic attenuation N_atten (int): the number of samples of the polychromatic attenuation @@ -975,6 +1081,8 @@ def setBHClookupTable(self, spectralResponse, gammas, Z, referenceEnergy=0.0, T_ N_atten = int(np.ceil(max_atten / T_atten)) + 1 max_atten = float(N_atten-1)*T_atten + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) N_gamma = gammas.size if referenceEnergy <= 0.0: #referenceEnergy = self.meanEnergy(spectralResponse, gammas) @@ -983,8 +1091,11 @@ def setBHClookupTable(self, spectralResponse, gammas, Z, referenceEnergy=0.0, T_ LUT = np.zeros(N_atten, dtype=np.float32) if isinstance(Z, str): chemicalFormula = Z + if chemicalFormula in materialFormulas.keys(): + chemicalFormula = materialFormulas[chemicalFormula] if sys.version_info[0] == 3: chemicalFormula = bytes(str(chemicalFormula), 'ascii') + self.libxrayphysics.setBHClookupTable_compound.restype = ctypes.c_bool self.libxrayphysics.setBHClookupTable_compound.argtypes = [ctypes.c_char_p, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_int, ndpointer(ctypes.c_float, flags="C_CONTIGUOUS"), ctypes.c_float, ctypes.c_int, ctypes.c_float] self.libxrayphysics.setBHClookupTable_compound(chemicalFormula, spectralResponse, gammas, N_gamma, LUT, T_atten, N_atten, referenceEnergy) @@ -1000,7 +1111,7 @@ def polynomialBHC(self, spectralResponse, gammas, Z, density, referenceEnergy=0. Args: spectralResponse (numpy array): spectra model gammas (numpy array): energies at which the spectra model is defined - Z (scalar or string): atomic number, chemical formula, or mixture of compounds with mass fractions + Z (scalar or string): atomic number, chemical formula, mixture of compounds with mass fractions, or member of the material library density (scalar): the mass density (g/cm^3) of the matial referenceEnergy (scalar): the energy (keV) of the monochromatic attenuation (if not specified uses the mean energy of the spectra) maxThickness (scalar): the maximum thickness (cm) that one expects to penetrate of this material (defines the range of values for the polynomial fit) @@ -1010,6 +1121,8 @@ def polynomialBHC(self, spectralResponse, gammas, Z, density, referenceEnergy=0. A numpy array of the coefficients of the BHC polynomial """ + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) order = max(1, min(order, 10)) if maxThickness <= 0.0: return None @@ -1053,6 +1166,10 @@ def setTwoMaterialBHClookupTable(self, spectralResponse, gammas, sigma_1, sigma_ A 2D numpy array that maps polychromatic attenuation to monochromatic attenuation """ + gammas = self.fixArray(gammas) + spectralResponse = self.fixArray(spectralResponse) + sigma_1 = self.fixArray(sigma_1) + sigma_2 = self.fixArray(sigma_2) N = gammas.size if spectralResponse.size != N or sigma_1.size != N or sigma_2.size != N: print('Error: Input arrays must all be the same size!') @@ -1096,6 +1213,11 @@ def setDEDlookupTable(self, spectralResponse_L, spectralResponse_H, gammas, basi """ #bool generateDEDlookUpTables(float* spectralResponses, float* gammas, int N_gamma, float* referenceEnergies, float* basisFunctions, float* LUT, float T_lac, int N_lac); + gammas = self.fixArray(gammas) + spectralResponse_L = self.fixArray(spectralResponse_L) + spectralResponse_H = self.fixArray(spectralResponse_H) + basisFunction_1 = self.fixArray(basisFunction_1) + basisFunction_2 = self.fixArray(basisFunction_2) N = gammas.size if spectralResponse_L.size != N or spectralResponse_H.size != N or basisFunction_1.size != N or basisFunction_2.size != N: print('Error: Input arrays must all be the same size!') @@ -1133,6 +1255,7 @@ def setDEDlookupTable(self, spectralResponse_L, spectralResponse_H, gammas, basi def PhotoelectricBasis(self, gammas): """Returns the photoelectric basis function""" + gammas = self.fixArray(gammas) basisFcn = gammas.copy() basisFcn[:] = gammas[:]**-3.0 basisFcn *= self.cross_section_scalar() @@ -1140,6 +1263,7 @@ def PhotoelectricBasis(self, gammas): def ComptonBasis(self, gammas): """Returns the Compton basis function""" + gammas = self.fixArray(gammas) ELECTRON_REST_MASS_ENERGY = 510.975 CLASSICAL_ELECTRON_RADIUS = 2.8179403267e-13 AVOGANDROS_NUMBER = 6.0221414107e23 @@ -1156,7 +1280,11 @@ def ComptonBasis(self, gammas): def PCAbases(self, listOfMaterials, gammas): """Returns the two basis functions derived from the PCA of several material cross sections""" + gammas = self.fixArray(gammas) N = len(listOfMaterials) + if N < 2: + print('Error: must give at least two materials') + return None, None X = np.zeros((N,gammas.size)) for n in range(N): X[n,:] = self.sigma(listOfMaterials[n], gammas) @@ -1171,3 +1299,142 @@ def PCAbases(self, listOfMaterials, gammas): b_2 = b_2/np.sqrt(np.sum(b_2**2)) return b_1, b_2 + + def print_material_library(self): + """Prints all materials in the material library""" + for key in materialFormulas: + print(key + ': ' + str(materialFormulas[key]) + ', ' + str(materialDensities[key]) + ' g/cm^3') + +# Define the material library +materialFormulas = dict([('vacuum', '78.084*N2 + 20.946*O2 + 0.9340*Ar + 0.0397*CO2'), +('air', '78.084*N2 + 20.946*O2 + 0.9340*Ar + 0.0397*CO2'), +('foam', 'C8H8'), +('water', 'H2O'), +('bone_cortical', '0.047234*H + 0.144330*C + 0.041990*N + 0.446096*O + 0.002200*Mg + 0.104970*P + 0.003150*S + 0.209930*Ca + 0.000100*Zn'), +('bone_compact','0.063984*H+0.278000*C+0.027000*N+0.410016*O+0.002000*Mg+0.070000*P+0.002000*S+0.147000*Ca'), +('adipose_tissue', '0.119477*H + 0.637240*C + 0.007970*N + 0.232333*O + 0.000500*Na + 0.000020*Mg + 0.000160*P + 0.000730*S + 0.001190*Cl + 0.000320*K + 0.000020*Ca + 0.000020*Fe + 0.000020*Zn'), +('lung', '0.10132562304283013*H + 0.10235810831090612*C + 0.028663471831760927*N + 0.75742799115584314*O + 0.0018408652066471242*Na + 0.00073034326133282642*Mg + 0.00080037617680309744*P + 0.0022510579972587113*S + 0.0026612507878702989*Cl + 0.0019409122287475113*K'), +('ethanol', 'C2H6O'), +('HP', 'H2O2'), +('saltWater', '96.5*H2O + 3.5*NaCl'), +('graphite', 'C'), +('diamond', 'C'), +('paper', 'C6H10O5'), +('rubber', 'C5H8'), +('wood', 'C6H12O6'), +('delrin', 'CH2O'), +('teflon', 'C2F4'), +('PMMA', 'C5O2H8'), +('acrylic', 'C5O2H8'), +('Polyethylene', 'C2H4'), +('PE', 'C2H4'), +('LDPE', 'C2H4'), +('HDPE', 'C2H4'), +('Polypropylene', 'C3H6'), +('PP', 'C3H6'), +('PET', 'C10H8O4'), +('Polycarbonate', 'C15H16O2'), +('PVC', 'C2H3Cl'), +('PVDF', 'C2H2F2'), +('melamine', 'C3H6N6'), +('PETG', 'C14H20O5S'), +('PPSU', 'C24H16O4S'), +('PEEK', 'C19H12O3'), +('silicone', '17.57539*C+3.666654*H+37.47571*O+41.28225*Si'), +('glass', '72*SiO2 + 14.2*Na2O + 10*CaO + 2.5*MgO + 0.6*Al2O3'), +('eglass', '54*SiO2 + 14*Al2O3 + 22*CaO + 10*B2O3'), +('aluminumOxide', 'Al2O3'), +('titaniumDioxide', 'TiO2'), +('brass', 'Cu3Zn2'), +('Inconel', '52.5*Ni + 19*Cr + 18*Fe + 5*Nb + 3*Mo + 1*Ti + 0.5*Al + 1*Co'), +('stainlessSteel', '68.5*Fe+17*Cr+12*Ni+2.5*Mo'), +('Ti64', '6*Al+4*V+90*Ti'), +('Ti5553', '5*Al+5*V+5*Mo+3*Cr+82*Ti'), +('304L', '70*Fe+20*Cr+10*Ni'), +('316L', '68.5*Fe+17*Cr+12*Ni+2.5*Mo'), +('Al6061', '97.35*Al+1.0*Mg+0.6*Si+0.35*Fe+0.3*Cu+0.1*Cr+0.1*Zn+0.1*Ti+0.1*Mn'), +('aermet100', '0.23*C+13.4*Co+11.1*Ni+3.1*Cr+1.2*Mo+70.97*Fe'), +('NM', 'CH3NO2'), +('PTM', 'C5H8N4O12'), +('RDX', '94*C3H6N6O6 + 6*C10F6Cl6H9'), +('PBXN', 'C50H80N57O57'), +('PBXN301', '0.8*C5H8N4O12+0.2*C2H6O1Si1'), +('TATB', 'C6H6N6O6'), +('HMX', 'C4H8N8O8'), +('LX17', '92.5*C6H6N6O6+7.5*C2ClF3'), +('TNT', 'C7H5N3O6'), +('GOS', 'Gd2O2S'), +('CsI', 'CsI'), +('LSO', 'Lu2SiO5'), +('LYSO', 'Lu1.8Y0.2SiO5'), +('NaI', 'NaI'), +('BGO', 'Bi4Ge3O12'), +('GSO', 'Gd2SiO5'), +('CZT', 'Cd0.9Zn0.1Te')]) + +# g/cm^3 +materialDensities = dict([('vacuum', 1.0e-16), +('air', 1.2e-3), +('foam', 0.035), +('water', 1.0), +('bone_cortical', 1.85), +('bone_compact',1.85), +('adipose_tissue', 0.92), +('lung', 0.25), +('ethanol', 0.789), +('HP', 1.450), +('saltWater', 1.025), +('graphite', 1.804), +('diamond', 3.51), +('paper', 0.8), +('rubber', 1.1), +('wood', 0.53), +('delrin', 1.3962), +('teflon', 2.1609), +('PMMA', 1.18), +('acrylic', 1.18), +('Polyethylene', 0.94), +('PE', 0.94), +('LDPE', 0.92), +('HDPE', 0.95), +('Polypropylene', 0.946), +('PP', 0.946), +('PET', 1.38), +('Polycarbonate', 1.21), +('PVC', 1.35), +('PVDF', 1.78), +('melamine', 1.57), +('PETG', 1.4), +('PPSU', 1.28), +('PEEK', 1.32), +('silicone', 1.159604139), +('glass', 2.53), +('eglass', 2.6), +('aluminumOxide', 4.0), +('titaniumDioxide', 4.23), +('brass', 8.73), +('Inconel', 8.1933), +('stainlessSteel', 8.0), +('Ti64', 4.429), +('Ti5553', 4.65), +('304L', 8.0), +('316L', 8.0), +('Al6061', 2.7), +('aermet100', 7.888), +('NM', 1.1371), +('PTM', 0.9), +('RDX', 1.63), +('PBXN', 1.65), +('PBXN301', 1.5615), +('TATB', 1.93), +('HMX', 1.91), +('LX17', 1.9), +('TNT', 1.65), +('GOS', 7.32), +('CsI', 4.51), +('LSO', 7.4), +('LYSO', 7.1), +('NaI', 3.67), +('BGO', 7.13), +('GSO', 6.71), +('CZT', 5.78)]) diff --git a/src/xrayphysics_c_interface.cpp b/src/xrayphysics_c_interface.cpp index 611d3f3..142cbc8 100644 --- a/src/xrayphysics_c_interface.cpp +++ b/src/xrayphysics_c_interface.cpp @@ -20,6 +20,17 @@ float atomicMass(int Z) return physics.atomicMass(Z); } +float massDensity(int Z) +{ + return physics.xsecTables.getMassDensity(Z); +} + +int elementSymbolToAtomicNumber(const char* chemForm) +{ + string chemForm_str = chemForm; + return physics.xsecTables.elementStringToAtomicNumber(chemForm_str); +} + bool simulateSpectra(float kVp, float takeOffAngle, int Z, float* gammas, int N, float* output) { return physics.simulateSpectra(kVp, takeOffAngle, Z, gammas, N, output); diff --git a/src/xrayphysics_c_interface.h b/src/xrayphysics_c_interface.h index 8d32b87..df17bef 100644 --- a/src/xrayphysics_c_interface.h +++ b/src/xrayphysics_c_interface.h @@ -13,6 +13,8 @@ extern "C" XRAYPHYSICS_API void about(); extern "C" XRAYPHYSICS_API float atomicMass(int Z); +extern "C" XRAYPHYSICS_API float massDensity(int Z); +extern "C" XRAYPHYSICS_API int elementSymbolToAtomicNumber(const char* chemForm); extern "C" XRAYPHYSICS_API bool simulateSpectra(float kVp, float takeOffAngle, int Z, float* gammas, int N, float* output); extern "C" XRAYPHYSICS_API bool changeTakeOffAngle(float kVp, float takeOffAngle_cur, float takeOffAngle_new, int Z, float* gammas, int N, float* s); diff --git a/src/xsec.h b/src/xsec.h index bedd969..196b194 100644 --- a/src/xsec.h +++ b/src/xsec.h @@ -96,10 +96,10 @@ class xsec Pm=61,Sm=62,Eu=63,Gd=64,Tb=65,Dy=66,Ho=67,Er=68,Tm=69,Yb=70,Lu=71,Hf=72,Ta=73,W=74,Re=75,Os=76,Ir=77,Pt=78,Au=79,Hg=80, Tl=81,Pb=82,Bi=83,Po=84,At=85,Rn=86,Fr=87,Ra=88,Ac=89,Th=90,Pa=91,U=92,Np=93,Pu=94,Am=95,Cm=96,Bk=97,Cf=98,Es=99,Fm=100}; -private: - int elementStringToAtomicNumber(string); +private: + size_t findNextCapitalLetter(string str, int pos = 0); size_t findNextPlusSign(string str, int pos = 0); vector splitIntoElements(string str);