deploy v2.1.0

build/lib/powerxrd/main.py

@@ -1,10 +1,12 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import scipy.optimize as optimize
+import pandas
 
 
-'''interplanar spacing "d_hkl" from Braggs law'''
 def braggs(twotheta,lmda=1.54):
+    '''interplanar spacing "d_hkl" from Braggs law'''
+
     'lambda in Angstroms'
     twothet_rad=twotheta*np.pi/180
 
@@ -36,45 +38,70 @@ def braggs_s(twotheta,lmda=1.54):
         dhkl = lmda /(2*np.sin(twothet_rad/2))
         dhkl = np.round(dhkl,2)
 
-
     return dhkl
 
 
-'''Scherrer equation'''
 def scherrer(K,lmda,beta,theta):
-
+    '''Scherrer equation'''
     print('Scherrer Width == K*lmda / (FWHM*cos(theta))')
     return K*lmda / (beta*np.cos(theta))    #tau
 
 
-'''Gaussian fit for FWHM'''
 def funcgauss(x,y0,a,mean,sigma):
-
+    '''Gaussian equation'''
     return y0+(a/(sigma*np.sqrt(2*np.pi)))*np.exp(-(x-mean)**2/(2*sigma*sigma))
 
-#def funcgauss(x,y0,a,mean,fwhm):
 
-#    return y0 + (a/(fwhm*np.sqrt(np.pi/(4*np.log(2)) )))*np.exp(-(4*np.log(2))*(x-mean)**2/(fwhm*fwhm))
+class Data:
+    def __init__(self,file):
+        '''Data structure.
+        Parameters
+        ----------
+        file : str
+            file name and/or path for XRD file in .xy format
+        '''
+        self.file = file
+
+
+    def importfile(self):
 
+        df = pandas.read_csv(self.file, sep='\t', header=None)   #'https://www.statology.org/pandas-read-text-file/'
+        x,y = np.array(df).T
 
+        return x,y
 
 class Chart:
 
     def __init__(self,x,y):
-        '''Model structure. Calls `3-D` array to process into 3-D model.
+        '''Chart structure. Constructs x-y XRD data to manipulate and analyze. 
 
         Parameters
         ----------
-        array : np.array(int)
-            array of the third-order populated with discrete, non-zero integers which may represent a voxel block type
-        hashblocks : dict[int]{str, float }
-            a dictionary for which the keys are the integer values on the discrete arrays (above) an the values are the color (str) for the specific key and the alpha for the voxel object (float)
+        x : np.array(float)
+            array with x-data 2-theta values
+        y : np.array(float)
+            array with y-data peak intensity values
+        K : float
+            dimensionless shape factor for Scherrer equation (default 0.9)
+        lambdaKa : float
+            X-ray wavelength of \alpha radiation
+        lambdaKi : float
+            X-ray wavelength of "i" radiation (\beta, \gamma, other)
         '''
-        self.x = x          # x values
-        self.y = y          # y values
+        self.x          = x          # x values
+        self.y          = y          # y values
+        self.K          = 0.9       
+        self.lambdaKa   = 0.15406
+        self.lambdaKi   = 0.139
 
-    '''Local maxima finder'''
     def local_max(self,xrange=[12,13]):
+        '''Local maxima finder
+
+        Parameters
+        ----------
+        xrange_Ka : [](float)
+            range of x to find local max
+        '''
         x1,x2=xrange
         xsearch_index=[]
         for n in self.x:
@@ -90,14 +117,22 @@ def local_max(self,xrange=[12,13]):
 
         return max_x, max_y
 
-    '''Emission lines arising from different types of radiation i.e. K_beta radiation
-    wavelength of K_beta == 0.139 nm'''
-    def emission_lines(self, show = True, twothet_range_Ka=[10,20], lmda_Ka = 0.154,lmda_Ki=0.139):
 
-        twothet_Ka_deg, int_Ka = Chart(self.x, self.y).local_max(xrange=twothet_range_Ka)
+    def emission_lines(self, show = True, xrange_Ka=[10,20]):
+        '''Emission lines arising from different types of radiation i.e. K_beta radiation
+        wavelength of K_beta == 0.139 nm
+        
+        Parameters
+        ----------
+        show: bool
+            show plot of XRD chart
+        xrange_Ka : [](float)
+            range of x-axis (2-theta) for K_alpha radiation
+        '''
+        twothet_Ka_deg, int_Ka = Chart(self.x, self.y).local_max(xrange=xrange_Ka)
         twothet_Ka=twothet_Ka_deg*np.pi/180
 
-        twothet_Ki = 2*np.arcsin((lmda_Ki/lmda_Ka)*np.sin(twothet_Ka/2))
+        twothet_Ki = 2*np.arcsin((self.lambdaKi/self.lambdaKa)*np.sin(twothet_Ka/2))
         twothet_Ki_deg = twothet_Ki*180/np.pi
 
         # return twothet_Ka_deg, int_Ka, twothet_Ki_deg
@@ -114,60 +149,120 @@ def emission_lines(self, show = True, twothet_range_Ka=[10,20], lmda_Ka = 0.154,
 
 
     def gaussfit(self):
+        '''Fit of a Gaussian curve ("bell curve") to raw x-y data'''
         meanest = self.x[list(self.y).index(max(self.y))]
         sigest = meanest - min(self.x)
-    #    print('estimates',meanest,sigest)
         popt, pcov = optimize.curve_fit(funcgauss,self.x,self.y,p0 = [min(self.y),max(self.y),meanest,sigest])
-        print('-Gaussian fit results-')
-    #    print('amplitude {}\nmean {}\nsigma {}'.format(*popt))
+        print('\n-Gaussian fit results-')
         print('y-shift {}\namplitude {}\nmean {}\nsigma {}'.format(*popt))
-
         print('covariance matrix \n{}'.format(pcov))
-    #    print('pcov',pcov)
         return popt
+
+
+    def SchPeak(self,show=True,xrange=[12,13]):
+        '''Scherrer width calculation for peak within a specified range
         
-    def SchPeak(self,show=True,xrange=[12,13],K=0.9,lambdaKa=0.15406):
+        Parameters
+        ----------
+        xrange : [](float)
+            range of x-axis (2-theta) where peak to be calculated is found
+        show: bool
+            show plot of XRD chart
+        '''
+
+        print('\nSchPeak: Scherrer width calc. for peak in range of [{},{}]'.format(*xrange))
 
-        x1,x2=xrange
         'xseg and yseg:x and y segments of data in selected xrange'
         xseg,yseg = [],[]
-        for n in self.x:
-            if n >= x1 and  n <= x2:
+        for n, j in zip(self.x,self.y):
+            if n >= xrange[0] and n <= xrange[1]:
                 xseg.append(n)
-                yseg.append(self.y[list(self.x).index(n)]) 
-        
+                yseg.append(j) 
+
 
-        y0,a,mean,sigma = Chart(self.x, self.y).gaussfit()
+        y0,a,mean,sigma = Chart(xseg,yseg).gaussfit()
         ysegfit = funcgauss(np.array(xseg),y0,a,mean,sigma)
-        
+
         'FULL WIDTH AT HALF MAXIMUM'
         FWHM_deg = sigma*2*np.sqrt(2*np.log(2))
         FWHM = FWHM_deg*np.pi/180
         print('\nFWHM == sigma*2*sqrt(2*ln(2)): {} degrees'.format(FWHM_deg))
 
+        HWMIN = sigma*np.sqrt(2*np.log((50)))
+        print('\nHalf-width Minimum (HWMIN) (1/50 max) == sigma*sqrt(2*ln(50)): {} degrees'.\
+            format(HWMIN))
+
         'scherrer width peak calculations'
         max_twotheta = xseg[list(yseg).index(max(yseg))]
 
         theta=max_twotheta/2
         theta=theta*np.pi/180
 
         print('K (shape factor): {}\nK-alpha: {} nm \nmax 2-theta: {} degrees'.\
-            format(K,lambdaKa,max_twotheta))
+            format(self.K,self.lambdaKa,max_twotheta))
 
-        Sch=scherrer(K,lambdaKa,FWHM,theta)
+        Sch=scherrer(self.K,self.lambdaKa,FWHM,theta)
         X,Y = xseg,ysegfit
 
         print('\nSCHERRER WIDTH: {} nm'.format(Sch))
 
         if show:
-            plt.plot(xseg,yseg,color='m')
+            plt.plot(X,Y,'c--')             # gauss fit 
+            plt.plot(xseg,yseg,color='m')   # fitted segment
 
-        return Sch,X,Y
+        left = mean - HWMIN
+        right = mean + HWMIN 
 
+        # return Sch,X,Y
+        return max_twotheta, Sch, left,right
 
-    '''Function for an "n" point moving average: '''
-    def mav(self,n=1,show=False):
+    def allpeaks_recur(self,show = True, left=0, right=1, tol=0.2,schpeaks=[]):
+        '''recursion component function for main allpeaks function below'''
+        maxx, maxy = Chart(self.x, self.y).local_max(xrange=[left,right])
+
+        Sch_x, Sch, l,r = Chart(self.x, self.y).SchPeak(show=show,xrange=[maxx-0.5,maxx+0.5])
 
+        # self.schpeaks.append(Sch)
+        schpeaks.append([Sch_x,Sch])
+        peak_max = maxy
+
+        if peak_max > tol*max(self.y):
+            Chart(self.x, self.y).allpeaks_recur(True, r, right,tol,schpeaks)
+            Chart(self.x, self.y).allpeaks_recur(True, left, l,tol,schpeaks)
+
+
+    def allpeaks(self, tol=0.2, show = True):
+        '''Driver code for allpeaks recursion : Automated Scherrer width calculation of all peaks
+        
+        Parameters
+        ----------
+        tol : float
+            tolerance for recursion - tol relates to the minimum peak height to be 
+            calculated as a percent of maximum peak in chart i.e. peak_max > tol*max(self.y)
+        show: bool
+            show plot of XRD chart
+        '''
+
+        #init xrange [left, right]
+        left = min(self.x)
+        right = max(self.x)
+        schpeaks_ = []
+
+        Chart(self.x, self.y).allpeaks_recur(True, left, right, tol, schpeaks_)
+
+
+        print('\nallpeaks : Automated Scherrer width calculation of all peaks'+\
+             ' [within a certain tolerance]\nSUMMARY:')
+
+        sortidcs = np.argsort(np.array(schpeaks_).T[0])
+        # print(sortidcs)
+        for i in sortidcs:
+            print('2-theta: {} deg - Sch width: {} nm'.format(*schpeaks_[i]))
+
+
+
+    def mav(self,n=1,show=False):
+        '''Function for an "n" point moving average. '''
         L=int(len(self.x)//n)
         newy=np.zeros(L)
         for i in range(L):
@@ -191,23 +286,21 @@ def mav(self,n=1,show=False):
         return newx,newy
 
 
-    '''Calculate relative peak intensity (i.e. comparing one peak to another)'''
     def XRD_int_ratio(self,xR1=[8.88,9.6],xR2=[10.81,11.52]):
-        'XRD b/t two intensities ratio'
+        '''Calculate relative peak intensity (i.e. comparing one peak to another)'''
+        # 'XRD b/t two intensities ratio'
         return Chart(self.x, self.y).local_max(xR2)[1]/Chart(self.x, self.y).local_max(xR1)[1]
 
 
 
     def backsub(self,tol=1,show=False):
         '''Background subtraction operation
         inputs:
-            x - x-data (e.g. 2Theta values)
-            y - y-data (e.g. Intensity)
             tol - tolerance (see below)
         outputs: 
             x
             y
-        
+
         This function is a running conditional statement 
         which evaluates whether a small increase 
         in the x-direction will increase the magnitude of the 
@@ -233,6 +326,4 @@ def backsub(self,tol=1,show=False):
         if show:
             plt.plot(self.x,self.y)
 
-        return self.x,backsub_y
-
-
+        return self.x,backsub_y

setup.py

@@ -16,7 +16,7 @@
 # This call to setup() does all the work
 setup(
     name="powerxrd",
-    version="2.0.0",
+    version="2.1.0",
     description="Simple tools to handle powder XRD (and XRD) data with Python.",
     long_description=long_description,
     long_description_content_type="text/markdown",