Changeset 658 for trunk/GSASIImath.py


Ignore:
Timestamp:
Jun 28, 2012 9:28:57 AM (10 years ago)
Author:
vondreele
Message:

use steps/4 or /6 for fourier & charge flipping
more mods to findOffset
remove a stray debug print

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r657 r658  
    522522        else:
    523523            return text
     524           
     525def adjHKLmax(SGData,Hmax):
     526    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
     527        Hmax[0] = ((Hmax[0]+3)/6)*6
     528        Hmax[1] = ((Hmax[1]+3)/6)*6
     529        Hmax[2] = ((Hmax[2]+1)/4)*4
     530    else:
     531        Hmax[0] = ((Hmax[0]+2)/4)*4
     532        Hmax[1] = ((Hmax[1]+2)/4)*4
     533        Hmax[2] = ((Hmax[2]+1)/4)*4
    524534
    525535def FourierMap(data,reflData):
     
    535545    A = G2lat.cell2A(cell[:6])
    536546    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     547    adjHKLmax(SGData,Hmax)
    537548    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
    538549#    Fhkl[0,0,0] = generalData['F000X']
     
    610621## keep this
    611622               
    612 def findOffset(SGData,rho,Fhkl):   
     623def findOffset(SGData,Fhkl):   
    613624    if SGData['SpGrp'] == 'P 1':
    614625        return [0,0,0]   
    615     mapShape = rho.shape
    616     steps = np.array(mapShape)
    617626    hklShape = Fhkl.shape
    618     mapHalf = np.array(mapShape)/2
     627    steps = np.array(hklShape)
    619628    Fmax = np.max(np.absolute(Fhkl))
    620629    hklHalf = np.array(hklShape)/2
     
    646655            dH = H-hkl
    647656            dang = ang-ang0
    648             if np.any(np.abs(dH)-8 > 0):    #keep low order DHs
     657            if np.any(np.abs(dH)-6 > 0):    #keep low order DHs
    649658                continue
    650659            DH.append(dH)
     
    658667    chisq = np.min(Mmap)
    659668    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
    660     print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
     669    print ' map offset chi**2: %.3f, map offset: %d %d %d, no. terms: %d'%(chisq,DX[0],DX[1],DX[2],len(DH))
    661670    return DX
    662671   
     
    678687    Vol = cell[6]
    679688    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     689    adjHKLmax(SGData,Hmax)
    680690    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
    681691    time0 = time.time()
     
    735745    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
    736746    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
    737     roll = findOffset(SGData,CErho,CEhkl)
     747    roll = findOffset(SGData,CEhkl)
    738748       
    739749    mapData['Rcf'] = Rcf
     
    747757    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
    748758   
    749     def noDuplicate(xyz,peaks,SGData,incre):                  #be careful where this is used - it's slow
     759    def noEquivalent(xyz,peaks,SGData):                  #be careful where this is used - it's slow
    750760        equivs = G2spc.GenAtom(xyz,SGData)
    751761        xyzs = [equiv[0] for equiv in equivs]
     
    754764                return False
    755765        return True
    756            
     766       
     767    def noDuplicate(xyz,peaks,incre):
     768        if True in [np.allclose(xyz*incre,peak*incre,atol=2.0) for peak in peaks]:
     769            return False
     770        return True
     771                       
    757772    def findRoll(rhoMask,mapHalf):
    758773        indx = np.array(ma.nonzero(rhoMask)).T
     
    806821        mapHalf = np.array(rho.shape)/2
    807822        res = mapData['Resolution']
    808         incre = 1./np.array(rho.shape)
     823        incre = np.array(rho.shape)
    809824        step = max(1.0,1./res)+1
    810825        steps = np.array(3*[step,])
     
    829844        if np.any(x1 < 0):
    830845            break
    831         peak = (np.array(x1[1:4])-rMI)*incre
     846        peak = (np.array(x1[1:4])-rMI)/incre
    832847        if not len(peaks):
    833848            peaks.append(peak)
     
    835850        else:
    836851            if keepDup:
    837                 peaks.append(peak)
    838                 mags.append(x1[0])
    839             elif noDuplicate(peak,peaks,SGData,incre) and x1[0] > 0.:
     852                if noDuplicate(peak,peaks,incre):
     853                    peaks.append(peak)
     854                    mags.append(x1[0])
     855            elif noEquivalent(peak,peaks,SGData) and x1[0] > 0.:
    840856                peaks.append(peak)
    841857                mags.append(x1[0])
Note: See TracChangeset for help on using the changeset viewer.