Changeset 1970 for trunk/GSASIImath.py
 Timestamp:
 Sep 2, 2015 4:01:41 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r1966 r1970 925 925 926 926 def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata,SStauM,Mast): 927 import scipy.special as sp928 import scipy.integrate as si929 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]*tau933 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]*tau938 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,intImag927 import pypowder as pwd 928 929 nxs = np.newaxis 930 glTau,glWt = pwd.pygauleg(0.,1.,32) 931 932 def expModInt(H,A,B): 933 tau = np.arange(1.,A.shape[1]+1)[:,nxs]*glTau #waves x 32 934 XmodA = np.swapaxes(A,1,2)[:,:,nxs]*np.sin(twopi*tau.T) #atoms X pos X 32 X waves 935 XmodB = np.swapaxes(B,1,2)[:,:,nxs]*np.cos(twopi*tau.T) 936 Xmod = np.sum(XmodA+XmodB,axis=1) #atoms X pos X 32; sum waves 937 #Xmod is ParSup in J2K  values match after Calcm2 with klic=1 938 D = H[:,3][:,nxs]*tau[0,nxs] 939 HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T #atoms X 32 X ops 940 sinHA = np.sum(np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1) 941 cosHA = np.sum(np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1) 942 943 # GSASIIpath.IPyBreak() 944 return cosHA.T,sinHA.T 945 945 946 946 Smult,TauT = SStauM # both atoms x SGops … … 952 952 Ab = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 953 953 Bb = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 954 # GSASIIpath.IPyBreak()955 954 if Ax.ndim > 2: 956 GpA = np.array(expModInt(SSUniq,Ax,Bx ,Smult))955 GpA = np.array(expModInt(SSUniq,Ax,Bx)) 957 956 else: 958 GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:] ,Smult))957 GpA = np.array(expModInt(SSUniq,Ax[:,np.newaxis,:],Bx[:,np.newaxis,:])) 959 958 return GpA # SGops x atoms 960 959 … … 1072 1071 drawatom[dcx:dcx+3] = X 1073 1072 return drawAtoms,Fade 1073 1074 # gauleg.py Gauss Legendre numerical quadrature, x and w computation 1075 # integrate from a to b using n evaluations of the function f(x) 1076 # usage: from gauleg import gaulegf 1077 # x,w = gaulegf( a, b, n) 1078 # area = 0.0 1079 # for i in range(1,n+1): # yes, 1..n 1080 # area += w[i]*f(x[i]) 1081 1082 import math 1083 def gaulegf(a, b, n): 1084 x = range(n+1) # x[0] unused 1085 w = range(n+1) # w[0] unused 1086 eps = 3.0E14 1087 m = (n+1)/2 1088 xm = 0.5*(b+a) 1089 xl = 0.5*(ba) 1090 for i in range(1,m+1): 1091 z = math.cos(3.141592654*(i0.25)/(n+0.5)) 1092 while True: 1093 p1 = 1.0 1094 p2 = 0.0 1095 for j in range(1,n+1): 1096 p3 = p2 1097 p2 = p1 1098 p1 = ((2.0*j1.0)*z*p2(j1.0)*p3)/j 1099 1100 pp = n*(z*p1p2)/(z*z1.0) 1101 z1 = z 1102 z = z1  p1/pp 1103 if abs(zz1) <= eps: 1104 break 1105 1106 x[i] = xm  xl*z 1107 x[n+1i] = xm + xl*z 1108 w[i] = 2.0*xl/((1.0z*z)*pp*pp) 1109 w[n+1i] = w[i] 1110 return np.array(x), np.array(w) 1111 # end gaulegf 1112 1074 1113 1075 1114 def BessJn(nmax,x):
Note: See TracChangeset
for help on using the changeset viewer.