Changeset 1996
- Timestamp:
- Oct 9, 2015 2:57:18 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1995 r1996 965 965 else: 966 966 Fmod = 1.0 967 if Ax.shape[1]: 968 tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 969 XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32 970 XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto 971 Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X pos X 32; sum waves 967 tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 968 XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32 969 XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto 970 Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X pos X 32; sum waves 972 971 if Au.shape[1]: 973 972 tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau #Uwaves x 32 … … 982 981 sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1) #refBlk X ops x atoms; sum for G-L integration 983 982 cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1) #ditto 984 # GSASIIpath.IPyBreak() 985 983 # GSASIIpath.IPyBreak() 986 984 return np.array([cosHA,sinHA]) # 2 x refBlk x SGops x atoms 987 985 … … 994 992 Mast: array orthogonalization matrix for Uij 995 993 ''' 996 import scipy.misc as sm 997 994 995 nxs = np.newaxis 996 numeric = False 998 997 cosHA,sinHA = Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast) 999 deltx = 0.00001 1000 deltf = 0.0001 1001 deltu = 0.000001 1002 Mf = FSSdata.T.shape 998 Mf = [H.shape[0],]+list(FSSdata.T.shape) #ops x atoms x waves x 2 (sin+cos frac mods) 1003 999 dGdMfC = np.zeros(Mf) 1004 1000 dGdMfS = np.zeros(Mf) 1005 Mx = XSSdata.T.shape 1006 Mx = [H.shape[0],]+list(Mx) #atoms x waves x sin+cos pos mods 1001 Mx = [H.shape[0],]+list(XSSdata.T.shape) #ops x atoms x waves x 6 (sin+cos pos mods) 1007 1002 dGdMxC = np.zeros(Mx) 1008 1003 dGdMxS = np.zeros(Mx) 1009 for i in range(Mx[1]): #atoms 1010 for j in range(Mx[2]): #waves 1011 for k in range(Mx[3]): #coeff 1012 XSSdata[k,j,i] += deltx 1013 Cap,Sap = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)) 1014 XSSdata[k,j,i] -= 2*deltx 1015 Cam,Sam = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)) 1016 XSSdata[k,j,i] += deltx 1017 dGdMxC[:,:,j,k] += (Cap-Cam)/(2*deltx) 1018 dGdMxS[:,:,j,k] += (Sap-Sam)/(2*deltx) 1019 Mu = USSdata.T.shape 1004 Mu = [H.shape[0],]+list(USSdata.T.shape) #ops x atoms x waves x 12 (sin+cos Uij mods) 1020 1005 dGdMuC = np.zeros(Mu) 1021 1006 dGdMuS = np.zeros(Mu) 1022 # GSASIIpath.IPyBreak() 1023 1024 1025 # nxs = np.newaxis 1026 # glTau,glWt = pwd.pygauleg(0.,1.,32) #get Gauss-Legendre intervals & weights 1027 # Bx = np.array(XSSdata[3:]).T #...cos pos mods 1028 # dGdBx = np.zeros_like(Bx) 1029 # Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms 1030 # dGdAf = np.zeros_like(Af) 1031 # Bf = np.array(FSSdata[1]).T #cos frac mods... 1032 # dGdBf = np.zeros_like(Bf) 1033 # Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 1034 # dGdAu = np.zeros_like(Au) 1035 # Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 1036 # dGdBu = np.zeros_like(Bu) 1037 # 1038 # GSASIIpath.IPyBreak() 1039 # if 'Fourier' in waveTypes: 1040 # nf = 0 1041 # nx = 0 1042 # XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32)) 1043 # FmodZ = np.zeros((Af.shape[0],Af.shape[1],32)) 1044 # if 'Crenel' in waveTypes: 1045 # nC = np.where('Crenel' in waveTypes) 1046 # nf = 1 1047 # #FmodZ = 0 replace 1048 # else: 1049 # nx = 1 1050 # if 'Sawtooth' in waveTypes: 1051 # nS = np.where('Sawtooth' in waveTypes) 1052 # #XmodZ = 0 replace 1053 # if Ax.shape[1]: 1054 # tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 1055 # StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:] #also dXmod/dAx 1056 # CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:] #also dXmod/dBx 1057 # XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32 1058 # XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto 1059 # Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X pos X 32; sum waves 1060 # D = H[:,3][:,nxs]*glTau[nxs,:] #m*tau; ops X 32 1061 # HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:] #ops x atoms X 32 1062 # if Af.shape[1]: 1063 # tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau #Fwaves x 32 1064 # StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf 1065 # CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf 1066 # FmodA = Af[:,nf:,nxs]*StauF #atoms X Fwaves X 32 1067 # FmodB = Bf[:,nf:,nxs]*CtauF 1068 # Fmod = np.sum(FmodA+FmodB+FmodC,axis=1) #atoms X 32; sum waves 1069 # else: 1070 # Fmod = np.ones_like(HdotX) 1071 # if Au.shape[1]: 1072 # tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau #Uwaves x 32 1073 # StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:] #also dUmod/dAu 1074 # CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:] #also dUmod/dBu 1075 # UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32 1076 # UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto 1077 # Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3) #atoms x 3x3 x 32; sum waves 1078 # HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.? 1079 # else: 1080 # HbH = np.ones_like(HdotX) 1081 # HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:] #ops x atoms x waves x 32 x xyz 1082 # HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:] 1007 if numeric: #numeric derivatives - slow!! works for pos waves 1008 deltx = 0.00001 1009 deltf = 0.0001 1010 deltu = 0.000001 1011 for i in range(Mx[1]): #atoms 1012 for j in range(Mx[2]): #waves 1013 for k in range(Mx[3]): #coeff 1014 XSSdata[k,j,i] += deltx 1015 Cap,Sap = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)) 1016 XSSdata[k,j,i] -= 2*deltx 1017 Cam,Sam = np.squeeze(Modulation(waveTypes,np.array([H,]),FSSdata,XSSdata,USSdata,Mast)) 1018 XSSdata[k,j,i] += deltx 1019 dGdMxC[:,:,j,k] += (Cap-Cam)/(2*deltx) 1020 dGdMxS[:,:,j,k] += (Sap-Sam)/(2*deltx) 1021 else: #analytic derivatives 1022 glTau,glWt = pwd.pygauleg(0.,1.,32) #get Gauss-Legendre intervals & weights 1023 Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms 1024 Bf = np.array(FSSdata[1]).T #cos frac mods... 1025 Ax = np.array(XSSdata[:3]).T #...cos pos mods 1026 Bx = np.array(XSSdata[3:]).T #...cos pos mods 1027 Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 1028 Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 1029 1030 if 'Fourier' in waveTypes: 1031 nf = 0 1032 nx = 0 1033 XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32)) 1034 FmodZ = np.zeros((Af.shape[0],Af.shape[1],32)) 1035 if 'Crenel' in waveTypes: 1036 nC = np.where('Crenel' in waveTypes) 1037 nf = 1 1038 #FmodZ = 0 replace 1039 else: 1040 nx = 1 1041 if 'Sawtooth' in waveTypes: 1042 nS = np.where('Sawtooth' in waveTypes) 1043 #XmodZ = 0 replace 1044 tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 1045 StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:] #also dXmod/dAx 1046 CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:] #also dXmod/dBx 1047 XmodA = Ax[:,nx:,:,nxs]*StauX #atoms X waves X pos X 32 1048 XmodB = Bx[:,nx:,:,nxs]*CtauX #ditto 1049 Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X pos X 32; sum waves 1050 D = H[:,3][:,nxs]*glTau[nxs,:] #m*tau; ops X 32 1051 HdotX = np.inner(H[:,:3],np.swapaxes(Xmod,1,2))+D[:,nxs,:] #ops x atoms X 32 1052 if Af.shape[1]: 1053 tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau #Fwaves x 32 1054 StauF = np.ones_like(Af)[:,nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf 1055 CtauF = np.ones_like(Bf)[:,nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf 1056 FmodA = Af[:,nf:,nxs]*StauF #atoms X Fwaves X 32 1057 FmodB = Bf[:,nf:,nxs]*CtauF 1058 Fmod = np.sum(FmodA+FmodB+FmodC,axis=1) #atoms X 32; sum waves 1059 else: 1060 Fmod = np.ones_like(HdotX) 1061 if Au.shape[1]: 1062 tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau #Uwaves x 32 1063 StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:] #also dUmod/dAu 1064 CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:] #also dUmod/dBu 1065 UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x 32 1066 UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto 1067 Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3) #atoms x 3x3 x 32; sum waves 1068 HbH = np.exp(-np.sum(H[:,nxs,nxs,:3]*np.inner(H[:,:3],Umod),axis=-1)) # ops x atoms x 32 add Overhauser corr.? 1069 else: 1070 HbH = np.ones_like(HdotX) 1071 HdotXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodA,-1,-2)[nxs,:,:,:,:] #ops x atoms x waves x 32 x xyz 1072 HdotXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(XmodB,-1,-2)[nxs,:,:,:,:] 1073 sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1) #ops x atoms x waves x xyz; sum for G-L integration 1074 cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1) #ditto 1075 # GSASIIpath.IPyBreak() 1083 1076 # dHdXA = H[:,nxs,nxs,nxs,:3]*np.swapaxes(StauX,-1,-2)[nxs,:,:,:,:] #ops x atoms x waves x 32 x xyz 1084 1077 # dHdXB = H[:,nxs,nxs,nxs,:3]*np.swapaxes(CtauX,-1,-2)[nxs,:,:,:,:] #ditto 1085 # sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1) #ops x atoms x waves x xyz; sum for G-L integration1086 # cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1) #ditto1087 1078 # dGdAx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXA*np.sin(twopi*HdotXA)+np.cos(twopi*HdotXA))*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1088 1079 ## ops x atoms x waves x xyz - imaginary part -
trunk/GSASIIplot.py
r1993 r1996 5050 5050 glShadeModel(GL_SMOOTH) 5051 5051 5052 def Draw(caller='',Fade= None):5052 def Draw(caller='',Fade=[]): 5053 5053 #useful debug? 5054 5054 # if caller: … … 5062 5062 pageName = G2frame.dataDisplay.GetPageText(page) 5063 5063 rhoXYZ = [] 5064 rho = [] 5064 5065 if len(D4mapData.get('rho',[])): #preferentially select 4D map if there 5065 5066 rho = D4mapData['rho'][:,:,:,int(G2frame.tau*10)] #pick current tau 3D slice … … 5129 5130 # glEnable(GL_BLEND) 5130 5131 # glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA) 5131 if Fade == None:5132 if not len(Fade): 5132 5133 atmFade = np.ones(len(drawingData['Atoms'])) 5133 5134 else: -
trunk/GSASIIstrIO.py
r1978 r1996 2289 2289 else: 2290 2290 hapDict[pfx+'TwinFr:0'] = twin[1][0] 2291 controlDict[pfx+'TwinNMN'] = twin[1][ 2]2291 controlDict[pfx+'TwinNMN'] = twin[1][1] 2292 2292 if Twins[0][1][1]: 2293 2293 hapVary.append(pfx+'TwinFr:'+str(it)) -
trunk/GSASIIstrMath.py
r1993 r1996 754 754 else: 755 755 if len(TwinLaw) > 1: 756 refl.T[9] = np.sum(fas[:,:,0] **2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)#FcT from primary twin element756 refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2 #FcT from primary twin element 757 757 refl.T[7] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+ \ 758 758 np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1) #Fc sum over twins … … 1067 1067 fag = fa*GfpuA[0]-fb*GfpuA[1] 1068 1068 fbg = fb*GfpuA[0]+fa*GfpuA[1] 1069 # GSASIIpath.IPyBreak()1070 1069 fas = np.sum(np.sum(fag,axis=-1),axis=-1) 1071 1070 fbs = np.sum(np.sum(fbg,axis=-1),axis=-1) 1071 # GSASIIpath.IPyBreak() 1072 1072 if 'P' in calcControls[hfx+'histType']: 1073 1073 refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) … … 1075 1075 else: 1076 1076 if len(TwinLaw) > 1: 1077 refl.T[10] = np.sum(fas[:,:,0] **2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)#FcT from primary twin element1077 refl.T[10] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2 #FcT from primary twin element 1078 1078 refl.T[8] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+ \ 1079 1079 np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1) #Fc sum over twins
Note: See TracChangeset
for help on using the changeset viewer.