Changeset 1966 for trunk/GSASIImath.py
 Timestamp:
 Aug 28, 2015 3:02:12 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r1958 r1966 924 924 ################################################################################ 925 925 926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM ):926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast): 927 927 import scipy.special as sp 928 928 import scipy.integrate as si 929 930 def cosMod(tau,H,A,B,S): 931 uTau = posFourier(tau,A,B,S) 932 SdotU = twopi*(np.inner(uTau.T,H[:3].T))+H[3]*tau 933 return list(np.cos(SdotU))[0] 934 935 def sinMod(tau,H,A,B,S): 936 uTau = posFourier(tau,A,B,S) 937 SdotU = twopi*(np.inner(uTau.T,H[:3]))+H[3]*tau 938 return list(np.sin(SdotU))[0] 939 940 def expModInt(H,A,B,S): 941 intReal = np.array([[si.romberg(cosMod,0,1,args=(h,a,b,s),rtol=1.e1) for a,b in zip(A,B)] for h,s in zip(H,S)]) 942 intImag = np.array([[si.romberg(sinMod,0,1,args=(h,a,b,s),rtol=1.e1) for a,b in zip(A,B)] for h,s in zip(H,S)]) 943 # return np.sqrt(intReal**2+intImag**2) 944 return intReal,intImag 945 929 946 Smult,TauT = SStauM # both atoms x SGops 930 947 m = SSUniq.T[3] 931 A = np.array(XSSdata[:3]) #sin mods x waves x atoms 932 B = np.array(XSSdata[3:]) #cos mods... 933 if XSSdata.ndim > 2: 934 nh = np.arange(XSSdata.shape[1]) #0..max harmonic order1 935 M = np.where(m>0,m+nh[:,np.newaxis],mnh[:,np.newaxis]) # waves x SGops 936 HdotA = (np.inner(A.T,SSUniq.T[:3].T)) #atoms X waves x SGops 937 HdotB = (np.inner(B.T,SSUniq.T[:3].T)) # ditto 938 Hdot = np.sqrt(HdotA**2+HdotB**2) 939 Hphi = np.sum(np.arctan2(HdotB,HdotA)*M,axis=1) 940 GpA = np.sum(sp.jn(M,twopi*Hdot),axis=1) # sum on waves 948 Ax = np.array(XSSdata[:3]).T #atoms x waves x sin pos mods 949 Bx = np.array(XSSdata[3:]).T #...cos pos mods 950 Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms 951 Bf = np.array(FSSdata[1]).T #cos frac mods... 952 Ab = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 953 Bb = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 954 # GSASIIpath.IPyBreak() 955 if Ax.ndim > 2: 956 GpA = np.array(expModInt(SSUniq,Ax,Bx,Smult)) 941 957 else: 942 HdotA = np.inner(A.T,SSUniq.T[:3].T) #atoms X SGops 943 HdotB = np.inner(B.T,SSUniq.T[:3].T) # ditto 944 Hdot = np.sqrt(HdotA**2+HdotB**2) 945 Hphi = np.arctan2(HdotB,HdotA) 946 GpA = sp.jn(m,twopi*Hdot) # atoms x SGops 947 948 # GSASIIpath.IPyBreak() 949 return GpA.T,Hphi.T # SGops x atoms 950 951 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM): 958 GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:],Smult)) 959 return GpA # SGops x atoms 960 961 def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast): 952 962 import scipy.special as sp 953 963 Smult,TauT = SStauM # both atoms x SGops … … 1063 1073 return drawAtoms,Fade 1064 1074 1075 def BessJn(nmax,x): 1076 ''' compute Bessel function J(n,x) from scipy routine & recurrance relation 1077 returns sequence of J(n,x) for n in range [nmax...0...nmax] 1078 1079 :param integer nmax: maximul order for Jn(x) 1080 :param float x: argument for Jn(x) 1081 1082 :returns numpy array: [J(nmax,x)...J(0,x)...J(nmax,x)] 1083 1084 ''' 1085 import scipy.special as sp 1086 bessJn = np.zeros(2*nmax+1) 1087 bessJn[nmax] = sp.j0(x) 1088 bessJn[nmax+1] = sp.j1(x) 1089 bessJn[nmax1] = bessJn[nmax+1] 1090 for i in range(2,nmax+1): 1091 bessJn[i+nmax] = 2*(i1)*bessJn[nmax+i1]/xbessJn[nmax+i2] 1092 bessJn[nmaxi] = bessJn[i+nmax]*(1)**i 1093 return bessJn 1094 1095 def BessIn(nmax,x): 1096 ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation 1097 returns sequence of I(n,x) for n in range [nmax...0...nmax] 1098 1099 :param integer nmax: maximul order for In(x) 1100 :param float x: argument for In(x) 1101 1102 :returns numpy array: [I(nmax,x)...I(0,x)...I(nmax,x)] 1103 1104 ''' 1105 import scipy.special as sp 1106 bessIn = np.zeros(2*nmax+1) 1107 bessIn[nmax] = sp.i0(x) 1108 bessIn[nmax+1] = sp.i1(x) 1109 bessIn[nmax1] = bessIn[nmax+1] 1110 for i in range(2,nmax+1): 1111 bessIn[i+nmax] = bessIn[nmax+i2]2*(i1)*bessIn[nmax+i1]/x 1112 bessIn[nmaxi] = bessIn[i+nmax] 1113 return bessIn 1114 1065 1115 1066 1116 ################################################################################
Note: See TracChangeset
for help on using the changeset viewer.