Changeset 1977
- Timestamp:
- Sep 23, 2015 3:28:12 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1976 r1977 34 34 import numpy.fft as fft 35 35 import scipy.optimize as so 36 import pypowder as p yd36 import pypowder as pwd 37 37 38 38 sind = lambda x: np.sin(x*np.pi/180.) … … 924 924 ################################################################################ 925 925 926 def Modulation(waveTypes,SSUniq,SGT,FSSdata,XSSdata,USSdata,Mast): 927 import pypowder as pwd 926 def Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast): 928 927 929 928 nxs = np.newaxis 930 929 glTau,glWt = pwd.pygauleg(0.,1.,32) #get Gauss-Legendre intervals & weights 931 930 932 def expModInt(H,Af,Bf,Ax,Bx,Au,Bu ,S):931 def expModInt(H,Af,Bf,Ax,Bx,Au,Bu): 933 932 ''' 934 933 H: array ops X hklt … … 971 970 else: 972 971 HbH = 1.0 973 D = H[:,3][:,nxs]* (glTau[nxs,:]+S[:,nxs])#m*tau; ops X 32972 D = H[:,3][:,nxs]*glTau[nxs,:] #m*tau; ops X 32 974 973 HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:] #atoms X 32 X ops 975 974 sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1) #atoms X ops; sum for G-L integration 976 975 cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1) #ditto 977 GSASIIpath.IPyBreak()976 # GSASIIpath.IPyBreak() 978 977 return cosHA.T,sinHA.T #ops X atoms 979 978 … … 984 983 Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 985 984 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)) 987 986 return GpA # SGops x atoms 988 987 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) 988 def 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 1005 1001 1006 1002 def posFourier(tau,psin,pcos,smul): -
trunk/GSASIIstrMath.py
r1976 r1977 983 983 waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 984 984 SStauM = list(GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata)) 985 SST = SSGT[:,3]986 985 if SGData['SGInv']: 987 986 SStauM[0] = np.hstack((SStauM[0],SStauM[0])) 988 987 SStauM[1] = np.hstack((SStauM[1],SStauM[1])) 989 SST = np.hstack((SST,-SST))990 988 modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) 991 989 FF = np.zeros(len(Tdata)) … … 1036 1034 SSPhi = np.hstack((SSPhi,-SSPhi)) 1037 1035 # 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]) 1040 1038 sinp = np.sin(phase) 1041 1039 cosp = np.cos(phase) … … 1103 1101 Phi = np.inner(H[:3],SGT) 1104 1102 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) 1107 1105 #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,:]) 1109 1107 sinp = np.sin(phase) 1110 1108 cosp = np.cos(phase) … … 1119 1117 fot = (FF+FP-Bab)*occ*Tcorr 1120 1118 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 1125 1121 1126 1122 fas = np.sum(np.sum(fa,axis=1),axis=1)
Note: See TracChangeset
for help on using the changeset viewer.