Changeset 658

Timestamp:

Jun 28, 2012 9:28:57 AM (10 years ago)

Author:

vondreele

Message:

use steps/4 or /6 for fourier & charge flipping
more mods to findOffset
remove a stray debug print

Location:

trunk

Files:

• GSASIImath.py (modified) (12 diffs)
• GSASIIphsGUI.py (modified) (1 diff)

trunk/GSASIImath.py

@@ -522 +522 @@
         else:
             return text
+            
+def adjHKLmax(SGData,Hmax):
+    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
+        Hmax[0] = ((Hmax[0]+3)/6)*6
+        Hmax[1] = ((Hmax[1]+3)/6)*6
+        Hmax[2] = ((Hmax[2]+1)/4)*4
+    else:
+        Hmax[0] = ((Hmax[0]+2)/4)*4
+        Hmax[1] = ((Hmax[1]+2)/4)*4
+        Hmax[2] = ((Hmax[2]+1)/4)*4
 
 def FourierMap(data,reflData):
@@ -535 +545 @@
     A = G2lat.cell2A(cell[:6])
     Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
+    adjHKLmax(SGData,Hmax)
     Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
 #    Fhkl[0,0,0] = generalData['F000X']
@@ -610 +621 @@
 ## keep this
                 
-def findOffset(SGData,rho,Fhkl):    
+def findOffset(SGData,Fhkl):    
     if SGData['SpGrp'] == 'P 1':
         return [0,0,0]    
-    mapShape = rho.shape
-    steps = np.array(mapShape)
     hklShape = Fhkl.shape
-    mapHalf = np.array(mapShape)/2
+    steps = np.array(hklShape)
     Fmax = np.max(np.absolute(Fhkl))
     hklHalf = np.array(hklShape)/2
@@ -646 +655 @@
             dH = H-hkl
             dang = ang-ang0
-            if np.any(np.abs(dH)-8 > 0):    #keep low order DHs
+            if np.any(np.abs(dH)-6 > 0):    #keep low order DHs
                 continue
             DH.append(dH)
@@ -658 +667 @@
     chisq = np.min(Mmap)
     DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
-    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
+    print ' map offset chi**2: %.3f, map offset: %d %d %d, no. terms: %d'%(chisq,DX[0],DX[1],DX[2],len(DH))
     return DX
     
@@ -678 +687 @@
     Vol = cell[6]
     Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
+    adjHKLmax(SGData,Hmax)
     Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
     time0 = time.time()
@@ -735 +745 @@
     print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
     CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
-    roll = findOffset(SGData,CErho,CEhkl)
+    roll = findOffset(SGData,CEhkl)
         
     mapData['Rcf'] = Rcf
@@ -747 +757 @@
     norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
     
-    def noDuplicate(xyz,peaks,SGData,incre):                  #be careful where this is used - it's slow
+    def noEquivalent(xyz,peaks,SGData):                  #be careful where this is used - it's slow
         equivs = G2spc.GenAtom(xyz,SGData)
         xyzs = [equiv[0] for equiv in equivs]
@@ -754 +764 @@
                 return False
         return True
-            
+        
+    def noDuplicate(xyz,peaks,incre):
+        if True in [np.allclose(xyz*incre,peak*incre,atol=2.0) for peak in peaks]:
+            return False
+        return True
+                        
     def findRoll(rhoMask,mapHalf):
         indx = np.array(ma.nonzero(rhoMask)).T
@@ -806 +821 @@
         mapHalf = np.array(rho.shape)/2
         res = mapData['Resolution']
-        incre = 1./np.array(rho.shape)
+        incre = np.array(rho.shape)
         step = max(1.0,1./res)+1
         steps = np.array(3*[step,])
@@ -829 +844 @@
         if np.any(x1 < 0):
             break
-        peak = (np.array(x1[1:4])-rMI)*incre
+        peak = (np.array(x1[1:4])-rMI)/incre
         if not len(peaks):
             peaks.append(peak)
@@ -835 +850 @@
         else:
             if keepDup:
-                peaks.append(peak)
-                mags.append(x1[0])
-            elif noDuplicate(peak,peaks,SGData,incre) and x1[0] > 0.:
+                if noDuplicate(peak,peaks,incre):
+                    peaks.append(peak)
+                    mags.append(x1[0])
+            elif noEquivalent(peak,peaks,SGData) and x1[0] > 0.:
                 peaks.append(peak)
                 mags.append(x1[0])

trunk/GSASIIphsGUI.py

@@ -3741 +3741 @@
             Ind.sort()
             Ind.reverse()
-            print Ind
             for ind in Ind:
                 mapPeaks = np.delete(mapPeaks,ind,0)