Source code for tesfdmtools.methods.brems

#!/usr/bin/env python3
#
# explore sample bremstrahlung spectrum
#
'''
.. module:: brems
   :synopsis: provide a model Bremsstrahlung spectrum 
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>

Provide a model Bremsstrahlung spectrum based on a real MXS spectrum
measured at low resolution

'''

import pkg_resources

txtname=pkg_resources.resource_filename(__name__,'MXS-spectrum.txt')

# list of line positions, to be excluded from the Bremsstrahlung spectrum
elines=[(5179.,5641.),(5722.,6166.),(6247.,6600.),(7280.,7642.),(7760.,8339.),(8556.,9154.),(3527.,3847.)]
# range of the Bremsstrahlung spectrum to be considered
ebrange=[3000.,11340.]

import numpy
from numpy.polynomial import chebyshev as cheb

def rspectxt(fname):
   '''
   Read a low resolution MXS spectrum from the provided resource file

   Args:
      * `fname` = file name of the MXS spectrum file

   Returns:
      * `energy` = array of energy values
      * `counts` = array of corresponding counts values 
 
   '''
   energy=numpy.zeros(10000)
   counts=numpy.zeros(10000)
   ff=open(fname,'r')
   k=0
   for ll in ff:		# read the spectrum text file
      if ll[0] != '#': 
         its=ll.strip(' \n').split()   
         energy[k]=float(its[2])
         counts[k]=float(its[0])
         k=k+1
   ff.close()
   energy=energy[0:k]
   counts=counts[0:k]
   energy=energy*1000.0		# convert energies from keV to eV
   indx,=numpy.where( ( energy > 0.0 ) & ( energy < 12000.0 ) )  # only use this range
   energy=energy[indx]
   counts=counts[indx]
   return (energy,counts)

def brems(energies):
   '''
   Make model Bremsstrahlungs spectrum for given energies based on MXS spectrum

   Args:
      * `energies` = array of energies for which the Bremsstrahlung intensity will be given

   Returns:
      * `intensities` = array of normalized intensities of the Bremsstrahlung spectrum

   '''
   energy,counts=rspectxt(txtname)  # read MXS spectrum

   for ee in elines:                # delete known lines
      indx,=numpy.where( ( energy < ee[0] ) | ( energy > ee[1] ) )
      energy=energy[indx]
      counts=counts[indx]

   indx,=numpy.where( ( energy > ebrange[0] ) & ( energy < ebrange[1] ) )  # limit spectral range
   energy=energy[indx]
   counts=counts[indx] 

   xx=numpy.linspace(energy[0],energy[-1],500,dtype=float)
   cc=cheb.chebfit(energy,counts,6)					   # fit Chebyshev polynomial
   cf=cheb.chebval(xx,cc)
   indx,=numpy.where( ( xx < 3400.0 ) | ( cf < 0.0 ) )			   # use only positive values
   cf[indx]=0.0

   bb=numpy.interp(energies,xx,cf)			# interpolate for given energies
   bb=bb/numpy.sum(bb)					# normalize

   return bb

#============================================================================
if __name__ == "__main__":

   '''
   Test program
   '''

   import matplotlib.pyplot as plt

   ee=numpy.linspace(1000.0,12000.0,512)
   bb=brems(ee)   

   f,ax=plt.subplots()

   eee,sss=rspectxt(txtname)
   ax.plot(eee,sss,label='MXS spectrum')
   tt=numpy.sum(sss)
   bb=bb*tt*0.415

   ax.plot(ee,bb,label='Brems. model')
   ax.plot([ee[0],ee[-1]],[0.0,0.0],'k-')
   ax.set_xlabel('Energy (eV')
   ax.set_ylabel('Normalized Intensity')
   ax.set_title('Model Bremsstrahlungs spectrum')
   ax.legend()

   plt.show()
   plt.close('all')