Changeset 1928
- Timestamp:
- Jul 15, 2015 2:34:19 PM (8 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1927 r1928 490 490 Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm]) 491 491 Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ 492 return Hpos 492 Rhos = np.array([getRho(pos,mapData) for pos in Hpos]) 493 imax = np.argmax(Rhos) 494 return [Hpos[imax],] 493 495 return [] 494 496 … … 2049 2051 hkl = np.asarray(hkl,dtype='i') 2050 2052 dp = 360.*Phi[i] #and phi 2051 a = cosd(ph )2052 b = sind(ph )2053 phasep = complex(a,b) +dp2054 phasem = complex(a,-b) -dp2053 a = cosd(ph+dp) 2054 b = sind(ph+dp) 2055 phasep = complex(a,b) 2056 phasem = complex(a,-b) 2055 2057 if 'Fobs' in mapData['MapType']: 2056 2058 F = np.where(Fosq>0.,np.sqrt(Fosq),0.) … … 2456 2458 return mapData,map4DData 2457 2459 2460 def getRho(xyz,mapData): 2461 ''' get scattering density at a point by 8-point interpolation 2462 param xyz: coordinate to be probed 2463 param: mapData: dict of map data 2464 2465 :returns: density at xyz 2466 ''' 2467 rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2) 2468 if not len(mapData): 2469 return 0.0 2470 rho = copy.copy(mapData['rho']) #don't mess up original 2471 if not len(rho): 2472 return 0.0 2473 mapShape = np.array(rho.shape) 2474 mapStep = 1./mapShape 2475 X = np.array(xyz)%1. #get into unit cell 2476 I = np.array(X*mapShape,dtype='int') 2477 return rho[I[0],I[1],I[2]] 2478 D = X-I*mapStep #position inside map cell 2479 D12 = D[0]*D[1] 2480 D13 = D[0]*D[2] 2481 D23 = D[1]*D[2] 2482 D123 = np.prod(D) 2483 Rho = rollMap(rho,I) #shifts map so point is in corner 2484 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]+ \ 2485 Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+ \ 2486 Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)- \ 2487 Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123) 2488 return R 2489 2458 2490 def SearchMap(generalData,drawingData,Neg=False): 2459 2491 '''Does a search of a density map for peaks meeting the criterion of peak -
trunk/GSASIIphsGUI.py
r1927 r1928 1591 1591 neigh[1][1].append(nextneigh[1][1][0]) 1592 1592 neigh[2] = max(0,nH) #set expected no. H's needed 1593 AddHydIds.append(neigh[1][1]) 1594 Neigh.append(neigh) 1593 if len(neigh[1][0]): 1594 AddHydIds.append(neigh[1][1]) 1595 Neigh.append(neigh) 1595 1596 if Neigh: 1597 mapError = False 1596 1598 dlg = G2gd.AddHatomDialog(G2frame,Neigh,data) 1597 1599 if dlg.ShowModal() == wx.ID_OK: … … 1599 1601 Neigh = dlg.GetData() 1600 1602 mapData = generalData['Map'] 1601 mapError = False 1602 for ineigh,neigh in enumerate(Neigh): 1603 for ineigh,neigh in enumerate(Neigh): 1603 1604 AddHydIds[ineigh].append(neigh[2]) 1604 if 'O' in neigh[0] and not len(mapData['rho']) or not 'delt-F' in mapData['MapType']:1605 if 'O' in neigh[0] and (not len(mapData['rho']) or not 'delt-F' in mapData['MapType']): 1605 1606 mapError = True 1606 1607 continue -
trunk/GSASIIplot.py
r1909 r1928 4521 4521 SetTranslation(newxy) 4522 4522 Tx,Ty,Tz = drawingData['viewPoint'][0] 4523 G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1) 4523 rho = G2mth.getRho([Tx,Ty,Tz],mapData) 4524 G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f; density: %.4f'%(Tx,Ty,Tz,rho),1) 4524 4525 elif event.MiddleIsDown(): 4525 4526 SetRotationZ(newxy) … … 5229 5230 choice += ['+: increase tau','-: decrease tau','0: set tau = 0'] 5230 5231 5232 Tx,Ty,Tz = drawingData['viewPoint'][0] 5233 rho = G2mth.getRho([Tx,Ty,Tz],mapData) 5234 G2frame.G2plotNB.status.SetStatusText('View point: %.4f, %.4f, %.4f; density: %.4f'%(Tx,Ty,Tz,rho),1) 5235 5231 5236 cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice) 5232 5237 cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
Note: See TracChangeset
for help on using the changeset viewer.