Changeset 1981


Ignore:
Timestamp:
Sep 29, 2015 3:04:45 PM (6 years ago)
Author:
vondreele
Message:

fixes to modulation plots & 3D & 4D Fourier map processing

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIImath.py

    r1980 r1981  
    22492249    '''
    22502250    generalData = data['General']
    2251     mapData = generalData['4DmapData']
     2251    map4DData = generalData['4DmapData']
     2252    mapData = generalData['Map']
    22522253    dmin = mapData['Resolution']
    22532254    SGData = generalData['SGData']
    … …  
    22662267        if ref[5] > dmin:
    22672268            Fosq,Fcsq,ph = ref[9:12]
     2269            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
    22682270            Uniq = np.inner(ref[:4],SSGMT)
    22692271            Phi = np.inner(ref[:4],SSGT)
    … …  
    22762278                phasem = complex(a,-b)
    22772279                if 'Fobs' in mapData['MapType']:
    2278                     F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
     2280                    F = np.sqrt(Fosq)
    22792281                    h,k,l,m = hkl+Hmax
    22802282                    Fhkl[h,k,l,m] = F*phasep
    … …  
    22822284                    Fhkl[h,k,l,m] = F*phasem
    22832285                elif 'Fcalc' in mapData['MapType']:
    2284                     F = np.where(Fcsq>0.,np.sqrt(Fcsq),0.)
     2286                    F = np.sqrt(Fcsq)
    22852287                    h,k,l,m = hkl+Hmax
    22862288                    Fhkl[h,k,l,m] = F*phasep
    … …  
    22882290                    Fhkl[h,k,l,m] = F*phasem                   
    22892291                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)
    22912293                    h,k,l,m = hkl+Hmax
    22922294                    Fhkl[h,k,l,m] = dF*phasep
    22932295                    h,k,l,m = -hkl+Hmax
    22942296                    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']
    22972303    mapData['Type'] = reflDict['Type']
    22982304    mapData['rho'] = np.real(rho)
    22992305    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
    23002306    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
    2301     return mapData
     2307    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
    23022308
    23032309def FourierMap(data,reflDict):
    … …  
    23692375    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
    23702376    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
    2371     return mapData
    23722377   
    23732378# map printing for testing purposes

  • trunk/GSASIIphsGUI.py

    r1978 r1981  
    23632363            atomSizer.Add(waveType,0,WACV)
    23642364            axchoice = ['x','y','z']
    2365             if len(Map['rho']):
     2365            if len(D4Map['rho']):
    23662366                atomSizer.Add(wx.StaticText(waveData,label=' Show contour map for axis:'),0,WACV)
    23672367                mapSel = wx.ComboBox(waveData,value=' ',choices=axchoice,
    … …  
    25152515        cx,ct,cs,cia = generalData['AtomPtrs']
    25162516        atomData = data['Atoms']
    2517         Map = generalData['4DmapData']
     2517        Map = generalData['Map']
     2518        D4Map = generalData['4DmapData']
    25182519        if waveData.GetSizer():
    25192520            waveData.GetSizer().Clear(True)
    … …  
    25562557        ReflData = GetReflData(G2frame,phaseName,reflNames)
    25572558        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.
    25592562        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)
    25612564        wx.CallAfter(UpdateWavesData)
    25622565           
    … …  
    60776080            pgbar.Destroy()
    60786081        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)
    60806088        mapData['Flip'] = False
    60816089        mapSig = np.std(mapData['rho'])
    … …  
    60846092        data['Drawing']['contourLevel'] = 1.
    60856093        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)
    60876095        UpdateDrawAtoms()
    60886096        G2plt.PlotStructure(G2frame,data)

  • trunk/GSASIIplot.py

    r1979 r1981  
    6464npatan2d = lambda x,y: 180.*np.arctan2(x,y)/np.pi
    6565GkDelta = unichr(0x0394)
     66Gkrho = unichr(0x03C1)
    6667   
    6768class G2PlotMpl(wx.Panel):   
    … …  
    29912992
    29922993def ModulationPlot(G2frame,data,atom,ax,off=0):
    2993     global Off,Atom,Ax
     2994    global Off,Atom,Ax,Slab,Off
    29942995    Off = off
    29952996    Atom = atom
    29962997    Ax = ax
     2998   
    29972999    def OnMotion(event):
    29983000        xpos = event.xdata
    29993001        if xpos:                                        #avoid out of frame mouse position
    30003002            ypos = event.ydata
     3003            ix = int(round(xpos*10))
     3004            iy = int(round((Slab.shape[0]-1)*(ypos+0.5-Off*0.005)))
    30013005            Page.canvas.SetCursor(wx.CROSS_CURSOR)
    30023006            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()                 
    30043009            except TypeError:
    30053010                G2frame.G2plotNB.status.SetStatusText('Select '+Title+' pattern first',1)
    … …  
    30463051    Page.keyPress = OnPlotKeyPress
    30473052    Map = General['4DmapData']
    3048     MapType = Map['MapType']
     3053    MapType = mapData['MapType']
    30493054    rhoSize = np.array(Map['rho'].shape)
    30503055    atxyz = np.array(atom[cx:cx+3])
    … …  
    30673072                scof.append(spos[0][:3])
    30683073                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)
    30703075    if mapData['Flip']:
    30713076        Title = 'Charge flip'
    … …  
    30933098    Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax))
    30943099    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))
    30963105    Page.canvas.draw()
    30973106   
Note: See TracChangeset for help on using the changeset viewer.