Changeset 1977


Ignore:
Timestamp:
Sep 23, 2015 3:28:12 PM (6 years ago)
Author:
vondreele
Message:

More modulation structure factor stuff, begin SSF derivatives

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1976 r1977  
    3434import numpy.fft as fft
    3535import scipy.optimize as so
    36 import pypowder as pyd
     36import pypowder as pwd
    3737
    3838sind = lambda x: np.sin(x*np.pi/180.)
     
    924924################################################################################
    925925   
    926 def Modulation(waveTypes,SSUniq,SGT,FSSdata,XSSdata,USSdata,Mast):
    927     import pypowder as pwd
     926def Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast):
    928927   
    929928    nxs = np.newaxis
    930929    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    931930   
    932     def expModInt(H,Af,Bf,Ax,Bx,Au,Bu,S):
     931    def expModInt(H,Af,Bf,Ax,Bx,Au,Bu):
    933932        '''
    934933        H: array ops X hklt
     
    971970        else:
    972971            HbH = 1.0
    973         D = H[:,3][:,nxs]*(glTau[nxs,:]+S[:,nxs])              #m*tau; ops X 32
     972        D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
    974973        HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:]     #atoms X 32 X ops
    975974        sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #atoms X ops; sum for G-L integration
    976975        cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #ditto
    977         GSASIIpath.IPyBreak()                 
     976#        GSASIIpath.IPyBreak()                 
    978977        return cosHA.T,sinHA.T      #ops X atoms
    979978
     
    984983    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    985984    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    986     GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu,SGT))
     985    GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu))
    987986    return GpA             # SGops x atoms
    988987   
    989 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast):
    990     import scipy.special as sp
    991     Smult,TauT = SStauM             # both atoms x SGops
    992     nh = np.zeros(1)
    993     if XSSdata.ndim > 2:
    994         nh = np.arange(XSSdata.shape[1])       
    995     M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
    996     A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])])
    997     HdotA = twopi*(np.inner(A.T,SSUniq.T[:3].T)+SSPhi+TauT[:,np.newaxis,:])
    998     Gpm = sp.jn(M-1,HdotA)
    999     Gpp = sp.jn(M+1,HdotA)
    1000     if Gpm.ndim > 3: #sum over multiple harmonics
    1001         Gpm = np.sum(Gpm,axis=0)
    1002         Gpp = np.sum(Gpp,axis=0)
    1003     dGpdk = 0.5*(Gpm+Gpp)
    1004     return np.real(dGpdk),np.imag(dGpdk)
     988def ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast):
     989   
     990    nxs = np.newaxis
     991    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
     992   
     993    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
     994    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
     995    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
     996    Bf = np.array(FSSdata[1]).T    #cos frac mods...
     997    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
     998    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
     999
     1000    return dGpdk
    10051001   
    10061002def posFourier(tau,psin,pcos,smul):
  • trunk/GSASIIstrMath.py

    r1976 r1977  
    983983    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    984984    SStauM = list(GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata))
    985     SST = SSGT[:,3]
    986985    if SGData['SGInv']:
    987986        SStauM[0] = np.hstack((SStauM[0],SStauM[0]))
    988987        SStauM[1] = np.hstack((SStauM[1],SStauM[1]))
    989         SST = np.hstack((SST,-SST))
    990988    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    991989    FF = np.zeros(len(Tdata))
     
    10361034            SSPhi = np.hstack((SSPhi,-SSPhi))
    10371035#        GSASIIpath.IPyBreak()
    1038         GfpuA = G2mth.Modulation(waveTypes,SSUniq,SST,FSSdata,XSSdata,USSdata,Mast)
    1039         phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
     1036        GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast)
     1037        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+SSPhi[:,np.newaxis])
    10401038        sinp = np.sin(phase)
    10411039        cosp = np.cos(phase)
     
    11031101        Phi = np.inner(H[:3],SGT)
    11041102        SSPhi = np.inner(H,SSGT)
    1105         GfpuA,GfpuB = G2mth.Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast)
    1106         dGAdk,dGBdk = G2mth.ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast)
     1103        GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast)
     1104        dGAdk = G2mth.ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast)
    11071105        #need ModulationDerv here dGAdXsin, etc 
    1108         phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
     1106        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+SSPhi[np.newaxis,:])
    11091107        sinp = np.sin(phase)
    11101108        cosp = np.cos(phase)
     
    11191117        fot = (FF+FP-Bab)*occ*Tcorr
    11201118        fotp = FPP*occ*Tcorr
    1121         fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
    1122         fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
    1123         fa = fa*GfpuA[:,:,np.newaxis]-fb*GfpuB[:,:,np.newaxis]
    1124         fb = fb*GfpuA[:,:,np.newaxis]+fa*GfpuB[:,:,np.newaxis]
     1119        fa *= GfpuA
     1120        fb *= GfpuA       
    11251121       
    11261122        fas = np.sum(np.sum(fa,axis=1),axis=1)
Note: See TracChangeset for help on using the changeset viewer.