Changeset 1928

Timestamp:

Jul 15, 2015 2:34:19 PM (8 years ago)

Author:

vondreele

Message:

fix error in Fourier calc.
new getRho routine
new add OH hydrogens looks in delt-F map
put xyz & rho of view point on status line

Location:

trunk

Files:

• GSASIImath.py (modified) (3 diffs)
• GSASIIphsGUI.py (modified) (2 diffs)
• GSASIIplot.py (modified) (2 diffs)

trunk/GSASIImath.py

@@ -490 +490 @@
             Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm])
             Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
-            return Hpos
+            Rhos = np.array([getRho(pos,mapData) for pos in Hpos])
+            imax = np.argmax(Rhos)
+            return [Hpos[imax],]
     return []
         
@@ -2049 +2051 @@
                 hkl = np.asarray(hkl,dtype='i')
                 dp = 360.*Phi[i]                #and phi
-                a = cosd(ph)
-                b = sind(ph)
-                phasep = complex(a,b)+dp
-                phasem = complex(a,-b)-dp
+                a = cosd(ph+dp)
+                b = sind(ph+dp)
+                phasep = complex(a,b)
+                phasem = complex(a,-b)
                 if 'Fobs' in mapData['MapType']:
                     F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
@@ -2456 +2458 @@
     return mapData,map4DData
     
+def getRho(xyz,mapData):
+    ''' get scattering density at a point by 8-point interpolation
+    param xyz: coordinate to be probed
+    param: mapData: dict of map data
+    
+    :returns: density at xyz
+    '''
+    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
+    if not len(mapData):
+        return 0.0
+    rho = copy.copy(mapData['rho'])     #don't mess up original
+    if not len(rho):
+        return 0.0
+    mapShape = np.array(rho.shape)
+    mapStep = 1./mapShape
+    X = np.array(xyz)%1.    #get into unit cell
+    I = np.array(X*mapShape,dtype='int')
+    return rho[I[0],I[1],I[2]]
+    D = X-I*mapStep         #position inside map cell
+    D12 = D[0]*D[1]
+    D13 = D[0]*D[2]
+    D23 = D[1]*D[2]
+    D123 = np.prod(D)
+    Rho = rollMap(rho,I)    #shifts map so point is in corner
+    R = Rho[0,0,0]*(1.-np.sum(D))+Rho[1,0,0]*D[0]+Rho[0,1,0]*D[1]+Rho[0,0,1]*D[2]+  \
+        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
+        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
+        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)    
+    return R
+       
 def SearchMap(generalData,drawingData,Neg=False):
     '''Does a search of a density map for peaks meeting the criterion of peak

trunk/GSASIIphsGUI.py

@@ -1591 +1591 @@
                     neigh[1][1].append(nextneigh[1][1][0])
                 neigh[2] = max(0,nH)  #set expected no. H's needed
-                AddHydIds.append(neigh[1][1])
-                Neigh.append(neigh)
+                if len(neigh[1][0]):
+                    AddHydIds.append(neigh[1][1])
+                    Neigh.append(neigh)
             if Neigh:
+                mapError = False
                 dlg = G2gd.AddHatomDialog(G2frame,Neigh,data)
                 if dlg.ShowModal() == wx.ID_OK:
@@ -1599 +1601 @@
                     Neigh = dlg.GetData()
                     mapData = generalData['Map']
-                    mapError = False
-                    for ineigh,neigh in enumerate(Neigh):                        
+                    for ineigh,neigh in enumerate(Neigh):
                         AddHydIds[ineigh].append(neigh[2])
-                        if 'O' in neigh[0] and not len(mapData['rho']) or not 'delt-F' in mapData['MapType']:
+                        if 'O' in neigh[0] and (not len(mapData['rho']) or not 'delt-F' in mapData['MapType']):
                             mapError = True
                             continue                            

trunk/GSASIIplot.py

@@ -4521 +4521 @@
                     SetTranslation(newxy)
                     Tx,Ty,Tz = drawingData['viewPoint'][0]
-                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
+                    rho = G2mth.getRho([Tx,Ty,Tz],mapData)
+                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f; density: %.4f'%(Tx,Ty,Tz,rho),1)
                 elif event.MiddleIsDown():
                     SetRotationZ(newxy)
@@ -5229 +5230 @@
         choice += ['+: increase tau','-: decrease tau','0: set tau = 0']
 
+    Tx,Ty,Tz = drawingData['viewPoint'][0]
+    rho = G2mth.getRho([Tx,Ty,Tz],mapData)
+    G2frame.G2plotNB.status.SetStatusText('View point: %.4f, %.4f, %.4f; density: %.4f'%(Tx,Ty,Tz,rho),1)
+
     cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice)
     cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)