Changeset 1981 for trunk/GSASIImath.py


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

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

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.