Changeset 652


Ignore:
Timestamp:
Jun 22, 2012 1:53:54 PM (10 years ago)
Author:
vondreele
Message:

Add delete map peaks menu item
correct integration result for flat detector
remove print of map peak list
more work on findOffset

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r648 r652  
    3838
    3939[ wxID_FOURCALC, wxID_FOURSEARCH, wxID_PEAKSMOVE, wxID_PEAKSCLEAR, wxID_CHARGEFLIP,
    40     wxID_PEAKSUNIQUE,
    41 ] = [wx.NewId() for item in range(6)]
     40    wxID_PEAKSUNIQUE, wxID_PEAKSDELETE,
     41] = [wx.NewId() for item in range(7)]
    4242
    4343[ wxID_PWDRADD, wxID_HKLFADD, wxID_DATADELETE,
     
    568568        self.MapPeaksEdit.Append(id=wxID_PEAKSUNIQUE, kind=wx.ITEM_NORMAL,text='Unique peaks',
    569569            help='Reduce map peak list to unique set')
     570        self.MapPeaksEdit.Append(id=wxID_PEAKSDELETE, kind=wx.ITEM_NORMAL,text='Delete peaks',
     571            help='Delete selected peaks')
    570572        self.MapPeaksEdit.Append(id=wxID_PEAKSCLEAR, kind=wx.ITEM_NORMAL,text='Clear peaks',
    571573            help='Clear the map peak list')
  • trunk/GSASIIimage.py

    r584 r652  
    722722        del NST
    723723        if Dtth:
    724             H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
     724            H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
    725725        else:
    726             H2 = LUtth
     726            H2 = np.array(LUtth)
    727727        if Dazm:       
    728728            H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
    729729        else:
    730730            H1 = LRazm
     731        H0[0] /= npcosd(H2[:-1])
    731732        Nup += 1
    732733        dlg.Update(Nup)
  • trunk/GSASIImath.py

    r648 r652  
    2525import scipy.optimize as so
    2626import scipy.linalg as sl
     27import numpy.fft as fft
    2728
    2829sind = lambda x: np.sin(x*np.pi/180.)
     
    524525def FourierMap(data,reflData):
    525526   
    526     import numpy.fft as fft
    527527    generalData = data['General']
    528528    if not generalData['Map']['MapType']:
     
    583583    return mapData
    584584   
     585# map printing for testing purposes
     586def printRho(SGLaue,rho,rhoMax):                         
     587    dim = len(rho.shape)
     588    if dim == 2:
     589        ix,jy = rho.shape
     590        for j in range(jy):
     591            line = ''
     592            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     593                line += (jy-j)*'  '
     594            for i in range(ix):
     595                r = int(100*rho[i,j]/rhoMax)
     596                line += '%4d'%(r)
     597            print line+'\n'
     598    else:
     599        ix,jy,kz = rho.shape
     600        for k in range(kz):
     601            print 'k = ',k
     602            for j in range(jy):
     603                line = ''
     604                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     605                    line += (jy-j)*'  '
     606                for i in range(ix):
     607                    r = int(100*rho[i,j,k]/rhoMax)
     608                    line += '%4d'%(r)
     609                print line+'\n'
     610## keep this
     611               
    585612def findOffset(SGData,rho,Fhkl):
    586613   
     
    592619    if SGData['SpGrp'] == 'P 1':
    593620        return [0,0,0]   
    594 # will need to consider 'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111'
     621# will need to consider 'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111'?
    595622    mapShape = rho.shape
    596623    steps = np.array(mapShape)
     
    612639    DH = []
    613640    Dphi = []
    614     while i < 20:
     641    while i < 20:       #use 20 strongest F's
    615642        F = Flist[i]
    616643        hkl = np.unravel_index(Fdict[F],hklShape)
     
    629656                continue
    630657            DH.append(dH)
    631             Dphi.append((dang+0.5) % 1.0)
     658            Dphi.append((dang+.5) % 1.0)
    632659        i += 1
    633660    DH = np.array(DH)
     
    636663#        print item[0],'%.4f'%(item[1])
    637664    DX = np.zeros(3)
    638     X,Y,Z = np.mgrid[0:1:10j,0:1:10j,0:1:10j]
     665    X,Y,Z = np.mgrid[0:12,0:12,0:12]
    639666    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
    640667    Mmin = 1.e10
     668    Ms = []
    641669   
    642670    for xyz in XYZ:             #do a global search for best roll
    643         M = np.sum(calcPhase(xyz,DH,Dphi)**2)
     671        M = np.sum(calcPhase(xyz/12.,DH,Dphi)**2)
     672        Ms.append(M)
    644673        if M < Mmin:
    645674            DX = xyz
    646675            Mmin = M
    647    
     676#    Ms = np.array(Ms)
     677#    Ms = np.reshape(Ms,newshape=(12,12,12))
     678#    printRho(SGData['SGLaue'],Ms,np.max(Ms))   
     679    print ' Search result:',Mmin,DX
    648680    result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
    649681#    for item in zip(DH,Dphi,result[2]['fvec']):
     
    653685    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
    654686    return DX
     687   
     688def findOffset2(SGData,rho,Fhkl):
     689    if SGData['SpGrp'] == 'P 1':
     690        return [0,0,0]
     691    mapShape = rho.shape
     692    mapHalf = np.array(mapShape)/2
     693    steps = np.array(mapShape)
     694    hklShape = Fhkl.shape
     695    hklHalf = np.array(hklShape)/2
     696    for M,T in SGData['SGOps'][1:]:
     697        FThkl = np.zeros(shape=hklShape,dtype='c16')
     698        for ind,F in enumerate(Fhkl.flatten()):
     699            if F:
     700                HKL = np.unravel_index(ind,hklShape)-hklHalf
     701                HKLT = np.inner(HKL,M.T)
     702                phi = (np.angle(F,deg=True)/360.+np.vdot(T,HKL) % 1.)*360.
     703                phasep = complex(cosd(phi),-sind(phi))      #want F*
     704                h,k,l = HKLT+hklHalf
     705                FThkl[h,k,l] = np.absolute(F)*phasep
     706        Cd = np.real(fft.fftn(fft.fftshift(Fhkl*FThkl)))
     707        CdMax = np.array(np.unravel_index(np.argmax(Cd),mapShape))
     708        print CdMax,CdMax/2
     709    DX = CdMax/2           
     710   
     711    print ' Test map offset : %d %d %d'%(DX[0],DX[1],DX[2])
     712    return DX
     713
    655714
    656715def ChargeFlip(data,reflData,pgbar):
    657 #    import scipy.fftpack as fft
    658     import numpy.fft as fft
    659716    generalData = data['General']
    660717    mapData = generalData['Map']
     
    728785    np.seterr(**old)
    729786    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
    730     print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
     787    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
    731788    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     789#    roll = findOffset2(SGData,CErho,CEhkl) #doesn't work
    732790    roll = findOffset(SGData,CErho,CEhkl)
    733    
    734 #    for i,j,k in [[1,0,0],[0,1,0],[0,0,1]]:
    735 #        print i,j,k,CErho.shape,CEhkl.shape
    736 #   
    737 #        Rmap = np.roll(np.roll(np.roll(CErho,i,axis=0),j,axis=1),k,axis=2)
    738 #        Rhkl = fft.ifftshift(fft.ifftn(Rmap))
    739 #        R = findOffset(SGData,Rmap,Rhkl)
    740    
     791       
    741792    mapData['Rcf'] = Rcf
    742793    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
  • trunk/GSASIIphsGUI.py

    r650 r652  
    317317            generalData['POhkl'] = [0,0,1]
    318318        if 'Map' not in generalData:
    319             generalData['Map'] = {'MapType':'','RefList':'','Resolution':1.0,
     319            generalData['Map'] = {'MapType':'','RefList':'','Resolution':0.5,
    320320                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
    321321        if 'Flip' not in generalData:
    322             generalData['Flip'] = {'RefList':'','Resolution':1.0,'Norm element':'None',
    323                 'k-factor':1.1}
     322            generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
     323                'k-factor':0.1}
    324324        if 'doPawley' not in generalData:
    325325            generalData['doPawley'] = False
     
    12761276                cid = colLabels.index(parm)
    12771277            dlg.Destroy()
    1278             print parm,cid,indx
    12791278            if parm in ['Type']:
    12801279                dlg = G2elemGUI.PickElement(G2frame)
     
    12821281                    if dlg.Elem not in ['None']:
    12831282                        El = dlg.Elem.strip()
    1284                         print El,indx,cid
    12851283                        for r in indx:                       
    12861284                            atomData[r][cid] = El
     
    37373735        G2plt.PlotStructure(G2frame,data)
    37383736       
     3737    def OnPeaksDelete(event):
     3738        if 'Map Peaks' in data:
     3739            mapPeaks = data['Map Peaks']
     3740            Ind = MapPeaks.GetSelectedRows()
     3741            Ind.sort()
     3742            Ind.reverse()
     3743            print Ind
     3744            for ind in Ind:
     3745                mapPeaks = np.delete(mapPeaks,ind,0)
     3746            data['Map Peaks'] = mapPeaks
     3747        FillMapPeaksGrid()
     3748        G2plt.PlotStructure(G2frame,data)
     3749               
    37393750    def OnPeaksUnique(event):
    37403751        generalData = data['General']
     
    37563767                            mapPeaks[ind][1:] -= Dx/2.
    37573768        FillMapPeaksGrid()
     3769        G2plt.PlotStructure(G2frame,data)
    37583770   
    37593771    def OnFourierMaps(event):
     
    38223834            data['Map Peaks'] = np.concatenate((mags,peaks),axis=1)
    38233835           
    3824             print ' Map search peaks found:'
    3825             print '  No.    Mag.      x        y        z'
    3826             for j,i in enumerate(sortIdx[::-1]):
    3827                 print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
     3836#            print ' Map search peaks found:'
     3837#            print '  No.    Mag.      x        y        z'
     3838#            for j,i in enumerate(sortIdx[::-1]):
     3839#                print ' %3d %8.3f %8.5f %8.5f %8.5f'%(j,mags[i],peaks[i][0],peaks[i][1],peaks[i][2])
    38283840            print ' Map search finished, time = %.2fs'%(time.time()-time0)
    38293841        Page = G2frame.dataDisplay.FindPage('Map peaks')
     
    38323844        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
    38333845        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksUnique, id=G2gd.wxID_PEAKSUNIQUE)
     3846        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksDelete, id=G2gd.wxID_PEAKSDELETE)
    38343847        G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
    38353848        UpdateDrawAtoms()
     
    39493962            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksMove, id=G2gd.wxID_PEAKSMOVE)
    39503963            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksUnique, id=G2gd.wxID_PEAKSUNIQUE)
     3964            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksDelete, id=G2gd.wxID_PEAKSDELETE)
    39513965            G2frame.dataFrame.Bind(wx.EVT_MENU, OnPeaksClear, id=G2gd.wxID_PEAKSCLEAR)
    39523966            FillMapPeaksGrid()
Note: See TracChangeset for help on using the changeset viewer.