geometry fixes for deformation calcs. Now correct!

Author: vondreele

GSASII/GSASIIlattice.py

@@ -2938,7 +2938,7 @@ def H2ThPh2(H,Bmat):
     :returns array Ph: HKL polar angles
     '''
     Hcart = np.inner(H,Bmat)
-    R = np.sqrt(np.sum(np.square(Hcart),axis=1))
+    R = nl.norm(Hcart,axis=1)
     Hcart /= R[:,nxs]
     Pl = acosd(Hcart[:,2])
     Az = atan2d(Hcart[:,1],Hcart[:,0])
@@ -2991,7 +2991,7 @@ def SHarmcal(SytSym,SHFln,psi,gam):
             SHVal += (SHFln[term][0]*Ksl)
     return SHVal
 
-def KslCalc(trm,psi,gam):
+def KslCalc(trm,psi,gam,RI=False):
     '''Compute one angular part term in spherical harmonics
 
     :param str trm:sp. harm term name in the form of 'C(l,m)' or 'C(l,m)c' for cubic
@@ -3002,9 +3002,21 @@ def KslCalc(trm,psi,gam):
     '''
     l,m = eval(trm.strip('C').strip('c'))
     if 'c' in trm:
-        return CubicSHarm(l,m,psi,gam)
+        if not RI:
+            return CubicSHarm(l,m,psi,gam)
+        else:
+            return CubicSHarm(l,m,psi,gam),0.0
+
     else:
-        return SphHarmAng(l,m,1.0,psi,gam)
+        if not RI:
+            return SphHarmAng(l,m,1.0,psi,gam)
+        else:
+            try:
+            #### TODO: this will be deprecated in scipy 1.17.0
+                ylmp = SQ2*spsp.sph_harm(m,l,rpd*psi,rpd*gam)*(-1)**m   #wants radians; order then degree
+            except AttributeError: #new one sph_harm_y in scipy 1.15.1 but buggy?
+                ylmp = SQ2*spsp.sph_harm_y(l,m,rpd*psi,rpd*gam)*(-1)**m #order L,M makes more sense
+            return np.real(ylmp),np.imag(ylmp)
 
 def SphHarmAng(L,M,P,Th,Ph):
     ''' Compute spherical harmonics values using scipy.special.sph_harm

GSASII/GSASIIphsGUI2.py

@@ -243,7 +243,7 @@ def MakeUVmat(defData,U,V):
             MY /= nl.norm(MY)
             MX = np.cross(MY,MZ)
             MX /= nl.norm(MX)
-        return np.array([MX,MY,MZ]).T
+        return np.vstack((MX,MY,MZ)).T
 
     def OnDeformRef(event):
         Obj = event.GetEventObject()

GSASII/GSASIIplot.py

@@ -7359,8 +7359,8 @@ def Draw(caller='',Fade=[],NPkey=False):
                             Npsi,Ngam = 60,30 
                             PSI,GAM = np.mgrid[0:Npsi,0:Ngam]   #[azm,pol]
                             PSI = PSI.flatten()*360./Npsi  #azimuth 0-360 incl
-                            PSI += 90.
                             GAM = GAM.flatten()*180./Ngam  #polar 0-180 incl
+                            PSI += 90.
                             Rp,PSIp,GAMp = G2mth.RotPolbyM(np.ones_like(PSI),PSI,GAM,UVMat)
                             P = G2lat.SHarmcal(SytSym,SHC,PSIp,GAMp).reshape((Npsi,Ngam))
                             if np.min(P) < np.max(P):

GSASII/GSASIIstrMath.py

@@ -450,10 +450,11 @@ def MakePolar(Orient,QB):
         QA = G2mth.invQ(Orient)       #rotates about chosen axis
         Q = G2mth.prodQQ(QA,QB)     #might be switched? QB,QA is order for plotting
         M = np.inner(G2mth.Q2Mat(Q),Bmat)
-        return G2lat.H2ThPh2(np.reshape(HKL,(-1,3)),M)[1:]
+        return G2lat.H2ThPh2(hkl,M)[1:]
 
     dFFdS = {}
     atFlg = []
+    hkl = np.reshape(HKL,(-1,3))
     SQR = np.repeat(SQ,HKL.shape[1])
     for iAt,Atype in enumerate(Tdata):
         if 'Q' in Atype:        #spinning RB
@@ -551,9 +552,9 @@ def MakePolar(Orient,QB):
             orKeys = [item for item in ORBtables[Atype] if item not in ['Slater','ZSlater','NSlater','SZE','popCore','popVal']]
             orKeys = [item for item in orKeys if 'Sl' in item]
             orbs = SHCdict[iAt]['1']
-            UVmat = np.inner(SHCdict[-iAt]['UVmat'],Bmat)
-            R,Th,Ph = G2lat.H2ThPh2(np.reshape(HKL,(-1,3)),UVmat) #radius,azimuth,polar
-#            R = 1/R     # correct dspacings - not used
+            UVmat = np.inner(SHCdict[-iAt]['UVmat'],Bmat) #OK
+            R,Th,Ph = G2lat.H2ThPh2(hkl,UVmat) #radius,azimuth,polar - OK
+            R = 1./R  # d-spacing for diagnostics - OK
             atFlg.append(1.0)
             orbTable = ORBtables[Atype]  # should point at Sl core
             ffOrb = orbTable['Sl core']