changes to Absorbcorrect matrix math for deformation

Author: vondreele

GSASII/Absorb.py

@@ -89,6 +89,7 @@ class Absorb(wx.Frame):
     Pack = 0.50
     Radius = 0.4
     muT = 0.0
+    muR = 0.0
     def _init_coll_ABOUT_Items(self, parent):
 
         parent.Append(wxID_ABSORBABOUT,'About')
@@ -465,8 +466,9 @@ def SetWaveEnergy(self,Wave):
             self.muT = muT
 
         if self.Volume:
+            self.muR = self.Radius*self.Pack*muT/(10.0*self.Volume)
             Text += "%s %s%10.4g %s" % ("Total",' '+Gkmu+' =',self.Pack*muT/self.Volume,'cm'+Pwrm1+', ')
-            Text += "%s%10.4g%s" % ('Total '+Gkmu+'R =',self.Radius*self.Pack*muT/(10.0*self.Volume),', ')
+            Text += "%s%10.4g%s" % ('Total '+Gkmu+'R =',self.muR,', ')
             Text += "%s%10.4f%s\n" % ('Transmission exp(-2'+Gkmu+'R) =', \
                 100.0*math.exp(-2.*self.Radius*self.Pack*muT/(10.0*self.Volume)),'%')
             if muT > 0.:
@@ -475,6 +477,9 @@ def SetWaveEnergy(self,Wave):
                     Text += "%s %10.4g mm\n"%("1/e (~27.2%) penetration depth = ",pene)
                 else:
                     Text += "%s %10.4g %sm\n"%("1/e (~27.2%) penetration depth = ",1000.0*pene,Gkmu)
+                Tth = np.array([0.0,140.0])
+                Xabs = 1./G2pwd.Absorb('Cylinder',self.muR,Tth)
+                Text += "%s %10.4g to %10.4g\n"%("Cylinder absorption correction extremes:",Xabs[0],Xabs[1])
             self.Results.SetValue(Text)
             den = Mass/(0.602*self.Volume)
             if self.ifVol:
@@ -602,8 +607,7 @@ def UpDateAbsPlot(self,Wave,rePlot=True):
         self.bx.set_title(r"Cylinder abs. corr. for $\lambda= %.4f\AA$"%self.Wave,x=0,ha='left')
         self.bx.set_xlabel(r"$2\theta$",fontsize=14)
         Tth = np.arange(0.,140.0,0.1)
-        pmut = self.Radius*self.Pack*self.muT/(10.0*self.Volume)
-        Xabs = 1./G2pwd.Absorb('Cylinder',pmut,Tth)
+        Xabs = 1./G2pwd.Absorb('Cylinder',self.muR,Tth)
         self.bx.plot(Tth,Xabs)
         Ymin = 0.0
         Ymax = 0.0

GSASII/GSASIIlattice.py

@@ -2887,6 +2887,20 @@ def GetKclKsl(L,N,SGLaue,psi,phi,beta):
             Kcl = pcrs*(cosd(N*beta)+sind(N*beta))
     return Kcl*Ksl,Lnorm(L)
 
+def H2ThPh2(H,Bmat):
+    '''Convert HKL to spherical polar & azimuth angles
+
+    :param array H: array of hkl as [n,3]
+    :param [3,3] array Bmat: inv crystal to Cartesian transformation
+    :returns array Th: HKL azimuth angles
+    :returns array Ph: HKL polar angles
+    '''
+    Hcart = np.inner(H,Bmat)
+    R = np.sqrt(np.sum(np.square(Hcart),axis=1))
+    Pl = acosd(Hcart[:,2]/R)
+    Az = atan2d(Hcart[:,1],Hcart[:,0])
+    return R,Az,Pl
+
 def H2ThPh(H,Bmat,Q):
     '''Convert HKL to spherical polar & azimuth angles
 

GSASII/GSASIIstrMath.py

@@ -428,7 +428,7 @@ def extendChanges(prms):
         RBsu[k] = np.sqrt(np.inner(Avec.T,np.inner(covMatrix,Avec)))
     return RBsu
 
-def MakeSpHarmFF(HKL,Bmat,SHCdict,Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ,ifDeriv=False):
+def MakeSpHarmFF(HKL,Amat,Bmat,SHCdict,Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ,ifDeriv=False):
     ''' Computes hkl dependent form factors & derivatives from spinning rigid bodies
     :param array HKL: reflection hkl set to be considered
     :param array Bmat: inv crystal to Cartesian transfomation matrix
@@ -555,7 +555,8 @@ def MakePolar(Orient,QB):
                 orKeys = [item for item in orKeys if 'Sl' in item]
             orbs = SHCdict[iAt]
             UVmat = np.inner(nl.inv(SHCdict[-iAt]['UVmat']),Bmat)
-            Th,Ph = G2lat.H2ThPh(np.reshape(HKL,(-1,3)),UVmat,[1.,0.,0.,1.])
+            R,Th,Ph = G2lat.H2ThPh2(np.reshape(HKL,(-1,3)),UVmat)
+            R = 1/R     # correct dspacings
             atFlg.append(1.0)
             orbTable = ORBtables[Atype][orKeys[0]]  # should point at either Sl core or a Bessel core
             ffOrb = {item:orbTable[item] for item in orbTable if item not in ['Slater','ZSlater','NSlater','SZE','popCore','popVal']}
@@ -1214,7 +1215,7 @@ def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
         #FF has to have the Bessel*Sph.Har.*atm form factor for each refletion in Uniq for Q atoms; otherwise just normal FF
         #this must be done here. NB: same place for non-spherical atoms; same math except no Bessel part.
         if pfx in SHCdict:
-            MakeSpHarmFF(Uniq,Bmat,SHCdict[pfx],Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ)     #Not Amat!
+            MakeSpHarmFF(Uniq,Amat,Bmat,SHCdict[pfx],Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ)     #Not Amat!
         Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
         if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
             fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
@@ -1337,7 +1338,7 @@ def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
         Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])      #Nref*Nops,3,3
         Hij = np.reshape(np.array([G2lat.UijtoU6(uij) for uij in Hij]),(-1,len(SGT),6))     #Nref,Nops,6
         if pfx in SHCdict:
-            dffdsh,atFlg = MakeSpHarmFF(Uniq,Bmat,SHCdict[pfx],Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ,True)
+            dffdsh,atFlg = MakeSpHarmFF(Uniq,Amat,Bmat,SHCdict[pfx],Tdata,hType,FFtables,ORBtables,BLtables,FF,SQ,True)
             if len(dffdSH):
                 for item in dffdSH:
                     dffdSH[item] = np.concatenate((dffdSH[item],dffdsh[item]))