Changeset 1981
- Timestamp:
- Sep 29, 2015 3:04:45 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1980 r1981 2249 2249 ''' 2250 2250 generalData = data['General'] 2251 mapData = generalData['4DmapData'] 2251 map4DData = generalData['4DmapData'] 2252 mapData = generalData['Map'] 2252 2253 dmin = mapData['Resolution'] 2253 2254 SGData = generalData['SGData'] … … 2266 2267 if ref[5] > dmin: 2267 2268 Fosq,Fcsq,ph = ref[9:12] 2269 Fosq = np.where(Fosq>0.,Fosq,0.) #can't use Fo^2 < 0 2268 2270 Uniq = np.inner(ref[:4],SSGMT) 2269 2271 Phi = np.inner(ref[:4],SSGT) … … 2276 2278 phasem = complex(a,-b) 2277 2279 if 'Fobs' in mapData['MapType']: 2278 F = np. where(Fosq>0.,np.sqrt(Fosq),0.)2280 F = np.sqrt(Fosq) 2279 2281 h,k,l,m = hkl+Hmax 2280 2282 Fhkl[h,k,l,m] = F*phasep … … 2282 2284 Fhkl[h,k,l,m] = F*phasem 2283 2285 elif 'Fcalc' in mapData['MapType']: 2284 F = np. where(Fcsq>0.,np.sqrt(Fcsq),0.)2286 F = np.sqrt(Fcsq) 2285 2287 h,k,l,m = hkl+Hmax 2286 2288 Fhkl[h,k,l,m] = F*phasep … … 2288 2290 Fhkl[h,k,l,m] = F*phasem 2289 2291 elif 'delt-F' in mapData['MapType']: 2290 dF = np. where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)2292 dF = np.sqrt(Fosq)-np.sqrt(Fcsq) 2291 2293 h,k,l,m = hkl+Hmax 2292 2294 Fhkl[h,k,l,m] = dF*phasep 2293 2295 h,k,l,m = -hkl+Hmax 2294 2296 Fhkl[h,k,l,m] = dF*phasem 2295 rho = fft.fftn(fft.fftshift(Fhkl))/cell[6] 2296 print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size) 2297 SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6] #4D map 2298 rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6] #3D map 2299 map4DData['rho'] = np.real(SSrho) 2300 map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho'])) 2301 map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])] 2302 map4DData['Type'] = reflDict['Type'] 2297 2303 mapData['Type'] = reflDict['Type'] 2298 2304 mapData['rho'] = np.real(rho) 2299 2305 mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho'])) 2300 2306 mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])] 2301 return mapData2307 print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size) 2302 2308 2303 2309 def FourierMap(data,reflDict): … … 2369 2375 mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho'])) 2370 2376 mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])] 2371 return mapData2372 2377 2373 2378 # map printing for testing purposes -
trunk/GSASIIphsGUI.py
r1978 r1981 2363 2363 atomSizer.Add(waveType,0,WACV) 2364 2364 axchoice = ['x','y','z'] 2365 if len( Map['rho']):2365 if len(D4Map['rho']): 2366 2366 atomSizer.Add(wx.StaticText(waveData,label=' Show contour map for axis:'),0,WACV) 2367 2367 mapSel = wx.ComboBox(waveData,value=' ',choices=axchoice, … … 2515 2515 cx,ct,cs,cia = generalData['AtomPtrs'] 2516 2516 atomData = data['Atoms'] 2517 Map = generalData['4DmapData'] 2517 Map = generalData['Map'] 2518 D4Map = generalData['4DmapData'] 2518 2519 if waveData.GetSizer(): 2519 2520 waveData.GetSizer().Clear(True) … … 2556 2557 ReflData = GetReflData(G2frame,phaseName,reflNames) 2557 2558 if ReflData == None: return 2558 mapData.update(G2mth.Fourier4DMap(data,ReflData)) 2559 G2mth.Fourier4DMap(data,ReflData) 2560 data['Drawing']['contourLevel'] = 1. 2561 data['Drawing']['mapSize'] = 10. 2559 2562 mapSig = np.std(mapData['rho']) 2560 print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f sigma = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']),mapSig)2563 print '4D '+mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f sigma = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']),mapSig) 2561 2564 wx.CallAfter(UpdateWavesData) 2562 2565 … … 6077 6080 pgbar.Destroy() 6078 6081 else: 6079 mapData.update(G2mth.FourierMap(data,ReflData)) 6082 if generalData['Type'] in ['modulated',]: 6083 dim = '4D ' 6084 G2mth.Fourier4DMap(data,ReflData) 6085 else: 6086 dim = '3D ' 6087 G2mth.FourierMap(data,ReflData) 6080 6088 mapData['Flip'] = False 6081 6089 mapSig = np.std(mapData['rho']) … … 6084 6092 data['Drawing']['contourLevel'] = 1. 6085 6093 data['Drawing']['mapSize'] = 10. 6086 print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f sigma = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']),mapSig)6094 print dim+mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f sigma = %.3f'%(np.max(mapData['rho']),np.min(mapData['rho']),mapSig) 6087 6095 UpdateDrawAtoms() 6088 6096 G2plt.PlotStructure(G2frame,data) -
trunk/GSASIIplot.py
r1979 r1981 64 64 npatan2d = lambda x,y: 180.*np.arctan2(x,y)/np.pi 65 65 GkDelta = unichr(0x0394) 66 Gkrho = unichr(0x03C1) 66 67 67 68 class G2PlotMpl(wx.Panel): … … 2991 2992 2992 2993 def ModulationPlot(G2frame,data,atom,ax,off=0): 2993 global Off,Atom,Ax 2994 global Off,Atom,Ax,Slab,Off 2994 2995 Off = off 2995 2996 Atom = atom 2996 2997 Ax = ax 2998 2997 2999 def OnMotion(event): 2998 3000 xpos = event.xdata 2999 3001 if xpos: #avoid out of frame mouse position 3000 3002 ypos = event.ydata 3003 ix = int(round(xpos*10)) 3004 iy = int(round((Slab.shape[0]-1)*(ypos+0.5-Off*0.005))) 3001 3005 Page.canvas.SetCursor(wx.CROSS_CURSOR) 3002 3006 try: 3003 G2frame.G2plotNB.status.SetStatusText('t =%9.3f %s =%9.3f'%(xpos,GkDelta+Ax,ypos),1) 3007 G2frame.G2plotNB.status.SetStatusText('t =%9.3f %s =%9.3f %s=%9.3f'%(xpos,GkDelta+Ax,ypos,Gkrho,Slab[iy,ix]/8.),1) 3008 # GSASIIpath.IPyBreak() 3004 3009 except TypeError: 3005 3010 G2frame.G2plotNB.status.SetStatusText('Select '+Title+' pattern first',1) … … 3046 3051 Page.keyPress = OnPlotKeyPress 3047 3052 Map = General['4DmapData'] 3048 MapType = Map['MapType']3053 MapType = mapData['MapType'] 3049 3054 rhoSize = np.array(Map['rho'].shape) 3050 3055 atxyz = np.array(atom[cx:cx+3]) … … 3067 3072 scof.append(spos[0][:3]) 3068 3073 ccof.append(spos[0][3:]) 3069 wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof),1) #hmm, why -1 work for Na2CO3?3074 wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof),1) 3070 3075 if mapData['Flip']: 3071 3076 Title = 'Charge flip' … … 3093 3098 Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax)) 3094 3099 Slab = np.hstack((slab,slab,slab)) 3095 Plot.contour(Slab[:,:21],20,extent=(0.,2.,-.5+Off*.005,.5+Off*.005)) 3100 acolor = mpl.cm.get_cmap('RdYlGn') 3101 if 'delt' in MapType: 3102 Plot.contour(Slab[:,:21],20,extent=(0.,2.,-.5+Off*.005,.5+Off*.005),cmap=acolor) 3103 else: 3104 Plot.contour(Slab[:,:21],20,extent=(0.,2.,-.5+Off*.005,.5+Off*.005)) 3096 3105 Page.canvas.draw() 3097 3106
Note: See TracChangeset
for help on using the changeset viewer.