Changeset 646

Timestamp:

Jun 4, 2012 2:43:15 PM (13 years ago)

Author:

vondreele

Message:

eliminate Friedel argument in genHKLf - didn't do anything
work on molecule rotation in structure drawing stuff
a working version of findOffset, but slow

Location:

trunk

Files:

• GSASIImath.py (modified) (2 diffs)
• GSASIIphsGUI.py (modified) (1 diff)
• GSASIIplot.py (modified) (1 diff)
• GSASIIspc.py (modified) (4 diffs)
• GSASIIstruct.py (modified) (2 diffs)

trunk/GSASIImath.py ¶

@@ -584 +584 @@
     
 def findOffset(SGData,rho,Fhkl):
+    
+    def calcPhase(DX,DH,Dphi):
+        H,K,L = DH.T
+        Phi = (H*DX[0]+K*DX[1]+L*DX[2]+0.5) % 1.
+        return Dphi-Phi
+    
+    if SGData['SpGrp'] == 'P 1':
+        return [0,0,0]    
+# will need to consider 'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111'
     mapShape = rho.shape
+    steps = np.array(mapShape)
     hklShape = Fhkl.shape
     mapHalf = np.array(mapShape)/2
@@ -600 +610 @@
     F = str(1.e6)
     i = 0
+    DH = []
+    Dphi = []
     while float(F) > 0.5*Fmax:
         F = Flist[i]
         hkl = np.unravel_index(Fdict[F],hklShape)
-        iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData,Friedel=True)
-        Uniq = np.array(Uniq,dtype='i')+hklHalf
-        print hkl-hklHalf
+        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
+        Uniq = np.array(Uniq,dtype='i')
+        Phi = np.array(Phi)
+        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
+        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
         Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
-        ph0 = np.angle(Fh0,deg=True)/360.
-        for j,H in enumerate(Uniq):
-            Fh = Fhkl[H[0],H[1],H[2]]
-            h,k,l = H-hklHalf
-            ang = (np.angle(Fh,deg=True)/360.-phi[j]) % 1.
-            print '(%3d,%3d,%3d) %9.5f'%(h,k,l,ang)        
+        ang0 = np.angle(Fh0,deg=True)/360.
+        for j,H in enumerate(Uniq[1:]):
+            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1])
+            dH = H-hkl
+            dang = ang-ang0
+            DH.append(dH)
+            Dphi.append((dang+0.5) % 1.0)
         i += 1
-        
-        
-    
-    return [0,0,0]
+#    for item in zip(DH,Dphi):
+#        print item
+    DH = np.array(DH)
+    Dphi = np.array(Dphi)
+    DX = np.zeros(3)
+    X,Y,Z = np.mgrid[0:1:steps[0]*(0+1j),0:1:steps[1]*(0+1j),0:1:steps[2]*(0+1j)]
+    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
+    Mmin = 1.e10
+    
+    for xyz in XYZ:
+        M = np.sum(calcPhase(xyz,DH,Dphi)**2)
+        if M < Mmin:
+            DX = xyz
+            Mmin = M
+    
+    print DX,Mmin
+    result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
+    chisq = np.sum(result[2]['fvec']**2)
+    DX = np.array(np.rint(-result[0]*steps),dtype='i')
+    print chisq,DX
+    print result[0],result[0]*steps
+    return DX
 
 def ChargeFlip(data,reflData,pgbar):

trunk/GSASIIphsGUI.py ¶

@@ -3497 +3497 @@
             H = ref[:3]
             ref[4] = np.sqrt(1./G2lat.calc_rDsq2(H,G))
-            iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(H,SGData,Friedel=False)
+            iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(H,SGData)
             ref[3] = mulp/2             #convert from powder mulp.
             ref[11] = Uniq

trunk/GSASIIplot.py ¶

@@ -2863 +2863 @@
         glMatrixMode(GL_MODELVIEW)
         glLoadIdentity()
-        glRotate(anglez,0,0,1)
-        glRotate(angley,sind(anglez),cosd(anglez),0)
-        glRotate(anglex,cosd(anglez),-sind(anglez),0)
+        matX = G2lat.rotdMat(anglex,axis=0)
+        matY = G2lat.rotdMat(angley,axis=1)
+        matZ = G2lat.rotdMat(anglez,axis=2)
+        matRot = np.inner(matX,np.inner(matY,matZ))
+        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
+        glMultMatrixf(matRot.T)
         glMultMatrixf(A4mat.T)
         glTranslate(-Tx,-Ty,-Tz)

trunk/GSASIIspc.py ¶

@@ -35 +35 @@
              'SGOps': symmetry operations as [M,T] so that M*x+T = x'
              'SGSys': one of 'triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic'
+             'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111' for arbitrary axes
        '''
     LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
@@ -312 +313 @@
         return zip(XYZEquiv,Idup,Cell)
 
-def GenHKLf(HKL,SGData,Friedel=False):
+def GenHKLf(HKL,SGData):
     '''
     Uses old GSAS Fortran routine genhkl.for
@@ -318 +319 @@
         HKL - [h,k,l]
         SGData - space group data obtained from SpcGroup
-        Friedel = True to retain Friedel pairs for centrosymmetric case
     returns:
         iabsnt = True is reflection is forbidden by symmetry
-        mulp = reflection multiplicity including Fridel pairs
+        mulp = reflection multiplicity including Friedel pairs
         Uniq = numpy array of equivalent hkl in descending order of h,k,l
     '''
@@ -335 +335 @@
     Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
     phi = f[:Nuniq]
-    Uniq = np.array(Uniq)
-    
-    return iabsnt,2*mulp,Uniq,phi       #include Friedel pairs in powder mulp
+    
+    return iabsnt,mulp,Uniq,phi
                                   
 def GetOprPtrName(key):            

trunk/GSASIIstruct.py ¶

@@ -1211 +1211 @@
                 for h,k,l,d in HKLd:
                     ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
+                    mul *= 2      # for powder overlap of Friedel pairs
                     if ext:
                         continue
@@ -1983 +1984 @@
         biso = -SQfactor*Uisodata
         Tiso = np.where(biso<1.,np.exp(biso),1.0)
-#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
         HbH = -np.inner(H,np.inner(bij,H))
         Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])