Changeset 1984


Ignore:
Timestamp:
Oct 5, 2015 9:46:27 AM (6 years ago)
Author:
vondreele
Message:

some SS derv work

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1981 r1984  
    926926################################################################################
    927927   
    928 def Modulation2(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast):
     928def Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast):
    929929   
    930930    nxs = np.newaxis
     
    988988    return GpA             # 2 x refBlk x SGops x atoms
    989989   
    990 #def Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast):
    991 #   
    992 #    nxs = np.newaxis
    993 #    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    994 #   
    995 #    def expModInt(H,Af,Bf,Ax,Bx,Au,Bu):
    996 #        '''
    997 #        H: array ops X hklt
    998 #        Ax: array atoms X waves X xyz
    999 #        Bx: array atoms X waves X xyz
    1000 #        S: array ops
    1001 #        '''
    1002 #        if 'Fourier' in waveTypes:
    1003 #            nf = 0
    1004 #            nx = 0
    1005 #            XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
    1006 #            FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
    1007 #            if 'Crenel' in waveTypes:
    1008 #                nC = np.where('Crenel' in waveTypes)
    1009 #                nf = 1
    1010 #                #FmodZ = 0   replace
    1011 #        else:
    1012 #            nx = 1
    1013 #            if 'Sawtooth' in waveTypes:
    1014 #                nS = np.where('Sawtooth' in waveTypes)
    1015 #                #XmodZ = 0   replace
    1016 #        if Af.shape[1]:
    1017 #            tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
    1018 #            FmodA = Af[:,nf:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X 32
    1019 #            FmodB = Bf[:,nf:,nxs]*np.cos(twopi*tauF[nxs,:,:])
    1020 #            Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
    1021 #        else:
    1022 #            Fmod = 1.0           
    1023 #        if Ax.shape[1]:
    1024 #            tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
    1025 #            XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
    1026 #            XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
    1027 #            Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
    1028 #        if Au.shape[1]:
    1029 #            tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
    1030 #            UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x 32
    1031 #            UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
    1032 #            Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
    1033 #            HbH = np.exp(np.swapaxes(np.array([-np.inner(h,np.inner(Umod,h)) for h in H[:,:3]]),1,2)).T #Overhauser corr.?
    1034 #        else:
    1035 #            HbH = 1.0
    1036 #        D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
    1037 #        HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:]     #atoms X 32 X ops
    1038 #        sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #atoms X ops; sum for G-L integration
    1039 #        cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #ditto
    1040 ##        GSASIIpath.IPyBreak()                 
    1041 #        return cosHA.T,sinHA.T      #ops X atoms
    1042 #
    1043 #    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
    1044 #    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
    1045 #    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
    1046 #    Bf = np.array(FSSdata[1]).T    #cos frac mods...
    1047 #    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    1048 #    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    1049 #    GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu))
    1050 #    return GpA             # 2 x SGops x atoms
    1051 #   
    1052990def ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast):
    1053991   
     
    1055993    glTau,glWt = pwd.pygauleg(0.,1.,32)         #get Gauss-Legendre intervals & weights
    1056994   
     995    def expModInt(H,Af,Bf,Ax,Bx,Au,Bu):
     996        '''
     997        H: array ops X hklt
     998        Ax: array atoms X waves X xyz
     999        Bx: array atoms X waves X xyz
     1000        S: array ops
     1001        '''
     1002        if 'Fourier' in waveTypes:
     1003            nf = 0
     1004            nx = 0
     1005            XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32))
     1006            FmodZ = np.zeros((Af.shape[0],Af.shape[1],32))
     1007            if 'Crenel' in waveTypes:
     1008                nC = np.where('Crenel' in waveTypes)
     1009                nf = 1
     1010                #FmodZ = 0   replace
     1011        else:
     1012            nx = 1
     1013            if 'Sawtooth' in waveTypes:
     1014                nS = np.where('Sawtooth' in waveTypes)
     1015                #XmodZ = 0   replace
     1016        if Af.shape[1]:
     1017            tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau  #Fwaves x 32
     1018            FmodA = Af[:,nf:,nxs]*np.sin(twopi*tauF[nxs,:,:])   #atoms X Fwaves X 32
     1019            FmodB = Bf[:,nf:,nxs]*np.cos(twopi*tauF[nxs,:,:])
     1020            Fmod = np.sum(FmodA+FmodB+FmodC,axis=1)             #atoms X 32; sum waves
     1021        else:
     1022            Fmod = 1.0           
     1023        if Ax.shape[1]:
     1024            tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau  #Xwaves x 32
     1025            XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32
     1026            XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto
     1027            Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1)                #atoms X pos X 32; sum waves
     1028        if Au.shape[1]:
     1029            tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau     #Uwaves x 32
     1030            UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x 32
     1031            UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto
     1032            Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3)      #atoms x 3x3 x 32; sum waves
     1033            HbH = np.exp(np.swapaxes(np.array([-np.inner(h,np.inner(Umod,h)) for h in H[:,:3]]),1,2)).T #Overhauser corr.?
     1034        else:
     1035            HbH = 1.0
     1036        D = H[:,3][:,nxs]*glTau[nxs,:]              #m*tau; ops X 32
     1037        HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:]     #atoms X 32 X ops
     1038        sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #atoms X ops; sum for G-L integration
     1039        cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1)       #ditto
     1040#        GSASIIpath.IPyBreak()                 
     1041        return cosHA.T,sinHA.T      #ops X atoms
     1042
    10571043    Ax = np.array(XSSdata[:3]).T   #atoms x waves x sin pos mods
    10581044    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
     
    10611047    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
    10621048    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
    1063 
    1064     return dGpdk
     1049    GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu))
     1050    return GpA             # 2 x SGops x atoms
    10651051   
    10661052def posFourier(tau,psin,pcos,smul):
  • trunk/GSASIIstrMath.py

    r1980 r1984  
    10651065            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    10661066            fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr])
    1067         GfpuA = G2mth.Modulation2(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
     1067        GfpuA = G2mth.Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms
    10681068        fa *= GfpuA
    10691069        fb *= GfpuA       
     
    12211221        Phi = np.inner(H[:3],SGT)
    12221222        SSPhi = np.inner(H,SSGT)
    1223         GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast)
    1224         dGAdk = G2mth.ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast)
     1223        GfpuA,dGAdk = G2mth.ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast)
    12251224        #need ModulationDerv here dGAdXsin, etc 
    12261225        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+SSPhi[np.newaxis,:])
Note: See TracChangeset for help on using the changeset viewer.