Changeset 652

Timestamp:

Jun 22, 2012 1:53:54 PM (10 years ago)

Author:

vondreele

Message:

Add delete map peaks menu item
correct integration result for flat detector
remove print of map peak list
more work on findOffset

Location:

trunk

Files:

• GSASIIgrid.py (modified) (2 diffs)
• GSASIIimage.py (modified) (1 diff)
• GSASIImath.py (modified) (9 diffs)
• GSASIIphsGUI.py (modified) (8 diffs)

trunk/GSASIIgrid.py

@@ -38 +38 @@
 
 [ wxID_FOURCALC, wxID_FOURSEARCH, wxID_PEAKSMOVE, wxID_PEAKSCLEAR, wxID_CHARGEFLIP,
-    wxID_PEAKSUNIQUE,
-] = [wx.NewId() for item in range(6)]
+    wxID_PEAKSUNIQUE, wxID_PEAKSDELETE,
+] = [wx.NewId() for item in range(7)]
 
 [ wxID_PWDRADD, wxID_HKLFADD, wxID_DATADELETE,
@@ -568 +568 @@
         self.MapPeaksEdit.Append(id=wxID_PEAKSUNIQUE, kind=wx.ITEM_NORMAL,text='Unique peaks', 
             help='Reduce map peak list to unique set')
+        self.MapPeaksEdit.Append(id=wxID_PEAKSDELETE, kind=wx.ITEM_NORMAL,text='Delete peaks', 
+            help='Delete selected peaks')
         self.MapPeaksEdit.Append(id=wxID_PEAKSCLEAR, kind=wx.ITEM_NORMAL,text='Clear peaks', 
             help='Clear the map peak list')

trunk/GSASIIimage.py

@@ -722 +722 @@
         del NST
         if Dtth:
-            H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
+            H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
         else:
-            H2 = LUtth
+            H2 = np.array(LUtth)
         if Dazm:        
             H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
         else:
             H1 = LRazm
+        H0[0] /= npcosd(H2[:-1])
         Nup += 1
         dlg.Update(Nup)

trunk/GSASIImath.py

@@ -25 +25 @@
 import scipy.optimize as so
 import scipy.linalg as sl
+import numpy.fft as fft
 
 sind = lambda x: np.sin(x*np.pi/180.)
@@ -524 +525 @@
 def FourierMap(data,reflData):
     
-    import numpy.fft as fft
     generalData = data['General']
     if not generalData['Map']['MapType']:
@@ -583 +583 @@
     return mapData
     
+# map printing for testing purposes
+def printRho(SGLaue,rho,rhoMax):                          
+    dim = len(rho.shape)
+    if dim == 2:
+        ix,jy = rho.shape
+        for j in range(jy):
+            line = ''
+            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
+                line += (jy-j)*'  '
+            for i in range(ix):
+                r = int(100*rho[i,j]/rhoMax)
+                line += '%4d'%(r)
+            print line+'\n'
+    else:
+        ix,jy,kz = rho.shape
+        for k in range(kz):
+            print 'k = ',k
+            for j in range(jy):
+                line = ''
+                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
+                    line += (jy-j)*'  '
+                for i in range(ix):
+                    r = int(100*rho[i,j,k]/rhoMax)
+                    line += '%4d'%(r)
+                print line+'\n'
+## keep this
+                
 def findOffset(SGData,rho,Fhkl):
     
@@ -592 +619 @@
     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'
+# 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)
@@ -612 +639 @@
     DH = []
     Dphi = []
-    while i < 20:
+    while i < 20:       #use 20 strongest F's
         F = Flist[i]
         hkl = np.unravel_index(Fdict[F],hklShape)
@@ -629 +656 @@
                 continue
             DH.append(dH)
-            Dphi.append((dang+0.5) % 1.0)
+            Dphi.append((dang+.5) % 1.0)
         i += 1
     DH = np.array(DH)
@@ -636 +663 @@
 #        print item[0],'%.4f'%(item[1])
     DX = np.zeros(3)
-    X,Y,Z = np.mgrid[0:1:10j,0:1:10j,0:1:10j]
+    X,Y,Z = np.mgrid[0:12,0:12,0:12]
     XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
     Mmin = 1.e10
+    Ms = []
     
     for xyz in XYZ:             #do a global search for best roll
-        M = np.sum(calcPhase(xyz,DH,Dphi)**2)
+        M = np.sum(calcPhase(xyz/12.,DH,Dphi)**2)
+        Ms.append(M)
         if M < Mmin:
             DX = xyz
             Mmin = M
-    
+#    Ms = np.array(Ms)
+#    Ms = np.reshape(Ms,newshape=(12,12,12))
+#    printRho(SGData['SGLaue'],Ms,np.max(Ms))    
+    print ' Search result:',Mmin,DX
     result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
 #    for item in zip(DH,Dphi,result[2]['fvec']):
@@ -653 +685 @@
     print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
     return DX
+    
+def findOffset2(SGData,rho,Fhkl):
+    if SGData['SpGrp'] == 'P 1':
+        return [0,0,0]
+    mapShape = rho.shape
+    mapHalf = np.array(mapShape)/2
+    steps = np.array(mapShape)
+    hklShape = Fhkl.shape
+    hklHalf = np.array(hklShape)/2
+    for M,T in SGData['SGOps'][1:]:
+        FThkl = np.zeros(shape=hklShape,dtype='c16')
+        for ind,F in enumerate(Fhkl.flatten()):
+            if F:
+                HKL = np.unravel_index(ind,hklShape)-hklHalf
+                HKLT = np.inner(HKL,M.T)
+                phi = (np.angle(F,deg=True)/360.+np.vdot(T,HKL) % 1.)*360.
+                phasep = complex(cosd(phi),-sind(phi))      #want F*
+                h,k,l = HKLT+hklHalf
+                FThkl[h,k,l] = np.absolute(F)*phasep
+        Cd = np.real(fft.fftn(fft.fftshift(Fhkl*FThkl)))
+        CdMax = np.array(np.unravel_index(np.argmax(Cd),mapShape))
+        print CdMax,CdMax/2
+    DX = CdMax/2            
+    
+    print ' Test map offset : %d %d %d'%(DX[0],DX[1],DX[2])
+    return DX
+
 
 def ChargeFlip(data,reflData,pgbar):
-#    import scipy.fftpack as fft
-    import numpy.fft as fft
     generalData = data['General']
     mapData = generalData['Map']
@@ -728 +785 @@
     np.seterr(**old)
     print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
-    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
+    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
     CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
+#    roll = findOffset2(SGData,CErho,CEhkl) #doesn't work
     roll = findOffset(SGData,CErho,CEhkl)
-    
-#    for i,j,k in [[1,0,0],[0,1,0],[0,0,1]]:
-#        print i,j,k,CErho.shape,CEhkl.shape
-#    
-#        Rmap = np.roll(np.roll(np.roll(CErho,i,axis=0),j,axis=1),k,axis=2)
-#        Rhkl = fft.ifftshift(fft.ifftn(Rmap))
-#        R = findOffset(SGData,Rmap,Rhkl)
-    
+        
     mapData['Rcf'] = Rcf
     mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)

trunk/GSASIIphsGUI.py

@@ -317 +317 @@
             generalData['POhkl'] = [0,0,1]
         if 'Map' not in generalData:
-            generalData['Map'] = {'MapType':'','RefList':'','Resolution':1.0,
+            generalData['Map'] = {'MapType':'','RefList':'','Resolution':0.5,
                 'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
         if 'Flip' not in generalData:
-            generalData['Flip'] = {'RefList':'','Resolution':1.0,'Norm element':'None',
-                'k-factor':1.1}
+            generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
+                'k-factor':0.1}
         if 'doPawley' not in generalData:
             generalData['doPawley'] = False
@@ -1276 +1276 @@
                 cid = colLabels.index(parm)
             dlg.Destroy()
-            print parm,cid,indx
             if parm in ['Type']:
                 dlg = G2elemGUI.PickElement(G2frame)
@@ -1282 +1281 @@
                     if dlg.Elem not in ['None']:
                         El = dlg.Elem.strip()
-                        print El,indx,cid
                         for r in indx:                        
                             atomData[r][cid] = El
@@ -3737 +3735 @@
         G2plt.PlotStructure(G2frame,data)
         
+    def OnPeaksDelete(event):
+        if 'Map Peaks' in data:
+            mapPeaks = data['Map Peaks']
+            Ind = MapPeaks.GetSelectedRows()
+            Ind.sort()
+            Ind.reverse()
+            print Ind
+            for ind in Ind:
+                mapPeaks = np.delete(mapPeaks,ind,0)
+            data['Map Peaks'] = mapPeaks
+        FillMapPeaksGrid()
+        G2plt.PlotStructure(G2frame,data)
+                
     def OnPeaksUnique(event):
         generalData = data['General']
@@ -3756 +3767 @@
                             mapPeaks[ind][1:] -= Dx/2.
         FillMapPeaksGrid()
+        G2plt.PlotStructure(G2frame,data)
     
     def OnFourierMaps(event):
@@ -3822 +3834 @@
             data['Map Peaks'] = np.concatenate((mags,peaks),axis=1)
             
-            print ' Map search peaks found:'
-            print '  No.    Mag.      x        y        z'
-            for j,i in enumerate(sortIdx[::-1]):
-                print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
+#            print ' Map search peaks found:'
+#            print '  No.    Mag.      x        y        z'
+#            for j,i in enumerate(sortIdx[::-1]):
+#                print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
             print ' Map search finished, time = %.2fs'%(time.time()-time0)
         Page = G2frame.dataDisplay.FindPage('Map peaks')
@@ -3832 +3844 @@
         G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
         G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksUnique, id=G2gd.wxID_PEAKSUNIQUE)
+        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksDelete, id=G2gd.wxID_PEAKSDELETE)
         G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
         UpdateDrawAtoms()
@@ -3949 +3962 @@
             G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
             G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksUnique, id=G2gd.wxID_PEAKSUNIQUE)
+            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksDelete, id=G2gd.wxID_PEAKSDELETE)
             G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
             FillMapPeaksGrid()