Changeset 1626

Timestamp:

Jan 9, 2015 4:00:01 PM (8 years ago)

Author:

vondreele

Message:

Changed "Transform atom " to "Transform draw atom" to fix menu problem
implement a 4D version of findOffset - works OK
implement facility to allow moving map on 4th dimension

Location:

trunk

Files:

• GSASIIgrid.py (modified) (1 diff)
• GSASIImath.py (modified) (4 diffs)
• GSASIIplot.py (modified) (3 diffs)

trunk/GSASIIgrid.py

@@ -2490 +2490 @@
         self.DrawAtomEdit.Append(id=wxID_DRAWADDEQUIV, kind=wx.ITEM_NORMAL,text='Add atoms',
             help='Add symmetry & cell equivalents to drawing set from selected atoms')
-        self.DrawAtomEdit.Append(id=wxID_DRAWTRANSFORM, kind=wx.ITEM_NORMAL,text='Transform atoms',
+        self.DrawAtomEdit.Append(id=wxID_DRAWTRANSFORM, kind=wx.ITEM_NORMAL,text='Transform draw atoms',
             help='Transform selected atoms by symmetry & cell translations')
         self.DrawAtomEdit.Append(id=wxID_DRAWFILLCOORD, kind=wx.ITEM_NORMAL,text='Fill CN-sphere',

trunk/GSASIImath.py

@@ -1745 +1745 @@
     DH = []
     Dphi = []
+    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
+    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
     for F in Flist:
         hkl = np.unravel_index(Fdict[F],hklShape)
+        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
+            continue
         iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
         Uniq = np.array(Uniq,dtype='i')
@@ -1759 +1763 @@
             dH = H-hkl
             dang = ang-ang0
-            if np.any(np.abs(dH)- Hmax > 0):    #keep low order DHs
-                continue
             DH.append(dH)
             Dphi.append((dang+.5) % 1.0)
@@ -1877 +1879 @@
     mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
-    mapData['rollMap'] = [0,0,0]
     return mapData
+    
+def findSSOffset(SGData,SSGData,A,Fhklm):    
+    '''default doc string
+    
+    :param type name: description
+    
+    :returns: type name: description
+    
+    '''
+    if SGData['SpGrp'] == 'P 1':
+        return [0,0,0,0]    
+    hklmShape = Fhklm.shape
+    hklmHalf = np.array(hklmShape)/2
+    sortHKLM = np.argsort(Fhklm.flatten())
+    Fdict = {}
+    for hklm in sortHKLM:
+        HKLM = np.unravel_index(hklm,hklmShape)
+        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
+        if F == 0.:
+            break
+        Fdict['%.6f'%(np.absolute(F))] = hklm
+    Flist = np.flipud(np.sort(Fdict.keys()))
+    F = str(1.e6)
+    i = 0
+    DH = []
+    Dphi = []
+    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
+    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
+    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
+    for F in Flist:
+        hklm = np.unravel_index(Fdict[F],hklmShape)
+        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
+            continue
+        Uniq = np.inner(hklm-hklmHalf,SSGMT)
+        Phi = np.inner(hklm-hklmHalf,SSGT)
+        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
+        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
+        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
+        ang0 = np.angle(Fh0,deg=True)/360.
+        for H,phi in zip(Uniq,Phi)[1:]:
+            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
+            dH = H-hklm
+            dang = ang-ang0
+            DH.append(dH)
+            Dphi.append((dang+.5) % 1.0)
+        if i > 20 or len(DH) > 30:
+            break
+        i += 1
+    DH = np.array(DH)
+    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
+    Dphi = np.array(Dphi)
+    steps = np.array(hklmShape)
+    X,Y,Z,T = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2],0:1:1./steps[3]]
+    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
+    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
+    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
+    hist,bins = np.histogram(Mmap,bins=1000)
+#    for i,item in enumerate(hist[:10]):
+#        print item,bins[i]
+    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 %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
+#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
+    return DX
     
 def SSChargeFlip(data,reflDict,pgbar):
@@ -1970 +2035 @@
     SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
-    roll = findOffset(SGData,A,CEhkl[:,:,:,maxM+1])               #CEhkl needs to be just the observed set, not the full set!
+    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
 
     mapData['Rcf'] = Rcf
     mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
-    mapData['rollMap'] = [0,0,0]
 
     map4DData['Rcf'] = Rcf
-    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2))
+    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2),roll[3],axis=3))
     map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
     return mapData,map4DData

trunk/GSASIIplot.py

@@ -2798 +2798 @@
             Off += 1
         elif event.key == '-':
-            Off -= 1 
+            Off -= 1
+        elif event.key in ['l','r',]:
+            roll = 1
+            if  event.key == 'l':
+                roll = -1
+            rho = Map['rho']
+            Map['rho'] = np.roll(rho,roll,axis=3)
         wx.CallAfter(ModulationPlot,G2frame,data,Atom,Ax,Off)
 
@@ -2815 +2821 @@
         Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
         
-    Page.Choice = ['+: shift up','-: shift down','0: reset']
+    Page.Choice = ['+: shift up','-: shift down','0: reset shift','l: move left','r: move right']
     Page.keyPress = OnPlotKeyPress
     Page.SetFocus()
@@ -2853 +2859 @@
     if Ax == 'x':
         slab = np.sum(np.sum(rho[:,ix[1]-ib:ix[1]+ib,ix[2]-ib:ix[2]+ib,:],axis=2),axis=1)
-        Plot.plot(wave[0],tau)
+        Plot.plot(tau,wave[0])
     elif Ax == 'y':
         slab = np.sum(np.sum(rho[ix[0]-ib:ix[0]+ib,:,ix[2]-ib:ix[2]+ib,:],axis=2),axis=0)
-        Plot.plot(wave[1],tau)
+        Plot.plot(tau,wave[1])
     elif Ax == 'z':
         slab = np.sum(np.sum(rho[ix[0]-ib:ix[0]+ib,ix[1]-ib:ix[1]+ib,:,:],axis=1),axis=0)
-        Plot.plot(wave[2],tau)
+        Plot.plot(tau,wave[2])
     Plot.set_title(Title)
     Plot.set_xlabel('t')
     Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax))
-    Slab = np.concatenate((slab,slab),axis=1).T
-    Plot.contour(Slab,20,extent=(-.5+Off*.005,.5+Off*.005,0.,2.))
-    
+    Slab = np.concatenate((slab,slab),axis=1)
+    Plot.contour(Slab,20,extent=(0.,2.,-.5+Off*.005,.5+Off*.005))
     Page.canvas.draw()