Changeset 1626 for trunk/GSASIImath.py


Ignore:
Timestamp:
Jan 9, 2015 4:00:01 PM (7 years ago)
Author:
vondreele
Message:

Changed "Transform atom " to "Transform draw atom" to fix menu problem
implement a 4D version of findOffset - works OK
implement facility to allow moving map on 4th dimension

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1625 r1626  
    17451745    DH = []
    17461746    Dphi = []
     1747    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     1748    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    17471749    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
    17481750    for F in Flist:
    17491751        hkl = np.unravel_index(Fdict[F],hklShape)
     1752        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
     1753            continue
    17501754        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
    17511755        Uniq = np.array(Uniq,dtype='i')
     
    17591763            dH = H-hkl
    17601764            dang = ang-ang0
    1761             if np.any(np.abs(dH)- Hmax > 0):    #keep low order DHs
    1762                 continue
    17631765            DH.append(dH)
    17641766            Dphi.append((dang+.5) % 1.0)
     
    18771879    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
    18781880    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
    1879     mapData['rollMap'] = [0,0,0]
    18801881    return mapData
     1882   
     1883def findSSOffset(SGData,SSGData,A,Fhklm):   
     1884    '''default doc string
     1885   
     1886    :param type name: description
     1887   
     1888    :returns: type name: description
     1889   
     1890    '''
     1891    if SGData['SpGrp'] == 'P 1':
     1892        return [0,0,0,0]   
     1893    hklmShape = Fhklm.shape
     1894    hklmHalf = np.array(hklmShape)/2
     1895    sortHKLM = np.argsort(Fhklm.flatten())
     1896    Fdict = {}
     1897    for hklm in sortHKLM:
     1898        HKLM = np.unravel_index(hklm,hklmShape)
     1899        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
     1900        if F == 0.:
     1901            break
     1902        Fdict['%.6f'%(np.absolute(F))] = hklm
     1903    Flist = np.flipud(np.sort(Fdict.keys()))
     1904    F = str(1.e6)
     1905    i = 0
     1906    DH = []
     1907    Dphi = []
     1908    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
     1909    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
     1910    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
     1911    for F in Flist:
     1912        hklm = np.unravel_index(Fdict[F],hklmShape)
     1913        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
     1914            continue
     1915        Uniq = np.inner(hklm-hklmHalf,SSGMT)
     1916        Phi = np.inner(hklm-hklmHalf,SSGT)
     1917        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
     1918        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
     1919        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
     1920        ang0 = np.angle(Fh0,deg=True)/360.
     1921        for H,phi in zip(Uniq,Phi)[1:]:
     1922            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
     1923            dH = H-hklm
     1924            dang = ang-ang0
     1925            DH.append(dH)
     1926            Dphi.append((dang+.5) % 1.0)
     1927        if i > 20 or len(DH) > 30:
     1928            break
     1929        i += 1
     1930    DH = np.array(DH)
     1931    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
     1932    Dphi = np.array(Dphi)
     1933    steps = np.array(hklmShape)
     1934    X,Y,Z,T = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2],0:1:1./steps[3]]
     1935    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
     1936    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
     1937    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
     1938    hist,bins = np.histogram(Mmap,bins=1000)
     1939#    for i,item in enumerate(hist[:10]):
     1940#        print item,bins[i]
     1941    chisq = np.min(Mmap)
     1942    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
     1943    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
     1944#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
     1945    return DX
    18811946   
    18821947def SSChargeFlip(data,reflDict,pgbar):
     
    19702035    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
    19712036    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
    1972     roll = findOffset(SGData,A,CEhkl[:,:,:,maxM+1])               #CEhkl needs to be just the observed set, not the full set!
     2037    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
    19732038
    19742039    mapData['Rcf'] = Rcf
    19752040    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
    19762041    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
    1977     mapData['rollMap'] = [0,0,0]
    19782042
    19792043    map4DData['Rcf'] = Rcf
    1980     map4DData['rho'] = np.real(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2))
     2044    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2),roll[3],axis=3))
    19812045    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
    19822046    return mapData,map4DData
Note: See TracChangeset for help on using the changeset viewer.