Changeset 1116
 Timestamp:
 Oct 19, 2013 5:35:16 PM (8 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r1110 r1116 1342 1342 Hmax[2] = ((Hmax[2]+1)/4)*4 1343 1343 1344 def OmitMap(data,reflDict ):1344 def OmitMap(data,reflDict,pgbar=None): 1345 1345 '''default doc string 1346 1346 … … 1382 1382 h,k,l = hkl+Hmax 1383 1383 Fhkl[h,k,l] = F*phasem 1384 rho = fft.fftn(fft.fftshift(Fhkl))/cell[6] 1385 print 'NB: this is just an Fobs map for now  under development' 1384 rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6] 1385 M = np.mgrid[0:4,0:4,0:4] 1386 blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten())) 1387 iBeg = blkIds*rho0.shape/4 1388 iFin = (blkIds+1)*rho0.shape/4 1389 rho_omit = np.zeros_like(rho0) 1390 nBlk = 0 1391 for iB,iF in zip(iBeg,iFin): 1392 rho1 = np.copy(rho0) 1393 rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0. 1394 Fnew = fft.ifftshift(fft.ifftn(rho1)) 1395 Fnew = np.where(Fnew,Fnew,1.0) #avoid divide by zero 1396 phase = Fnew/np.absolute(Fnew) 1397 OFhkl = np.absolute(Fhkl)*phase 1398 rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j) 1399 rho_omit[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = np.copy(rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]]) 1400 nBlk += 1 1401 pgbar.Update(nBlk) 1402 mapData['rho'] = np.real(rho_omit)/cell[6] 1403 mapData['rhoMax'] = max(np.max(mapData['rho']),np.min(mapData['rho'])) 1386 1404 print 'Omit map time: %.4f'%(time.time()time0),'no. elements: %d'%(Fhkl.size) 1387 print rho.shape1388 mapData['rho'] = np.real(rho)1389 mapData['rhoMax'] = max(np.max(mapData['rho']),np.min(mapData['rho']))1390 1405 return mapData 1391 1406 
trunk/GSASIIphsGUI.py
r1110 r1116 4900 4900 reflData = G2frame.PatternTree.GetItemPyData(PatternId)[1] 4901 4901 if mapData['MapType'] == 'Omit': 4902 mapData.update(G2mth.OmitMap(data,reflData)) 4902 pgbar = wx.ProgressDialog('Omit map','Blocks done',65, 4903 style = wx.PD_ELAPSED_TIMEwx.PD_AUTO_HIDE) 4904 mapData.update(G2mth.OmitMap(data,reflData,pgbar)) 4905 pgbar.Destroy() 4903 4906 else: 4904 4907 mapData.update(G2mth.FourierMap(data,reflData))
Note: See TracChangeset
for help on using the changeset viewer.