Changeset 1979
- Timestamp:
- Sep 28, 2015 8:59:29 AM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIplot.py
r1978 r1979 320 320 ''' 321 321 from matplotlib.patches import Circle,CirclePolygon 322 global HKL,HKLF,HKLref 322 323 HKLref = hklRef 323 global HKL,HKLF,HKLref324 324 325 325 def OnSCKeyPress(event): -
trunk/GSASIIstrMath.py
r1978 r1979 971 971 phfx = pfx.split(':')[0]+hfx 972 972 ast = np.sqrt(np.diag(G)) 973 Mast = twopisq*np.multiply.outer(ast,ast) 974 975 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 973 Mast = twopisq*np.multiply.outer(ast,ast) 974 # SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 976 975 SGInv = SGData['SGInv'] 977 SGT = np.array([ops[1] for ops in SGData['SGOps']])976 # SGT = np.array([ops[1] for ops in SGData['SGOps']]) 978 977 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 979 978 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 980 979 FFtables = calcControls['FFtables'] 981 980 BLtables = calcControls['BLtables'] 981 Flack = 1.0 982 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 983 Flack = 1.-2.*parmDict[phfx+'Flack'] 984 TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],]) 985 TwDict = refDict.get('TwDict',{}) 986 if 'S' in calcControls[hfx+'histType']: 987 NTL = calcControls[phfx+'NTL'] 988 NM = calcControls[phfx+'TwinNMN']+1 989 TwinLaw = calcControls[phfx+'TwinLaw'] 990 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 991 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 982 992 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 983 993 waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 984 SStauM = list(GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata))985 if SGData['SGInv']:986 SStauM[0] = np.hstack((SStauM[0],SStauM[0]))987 SStauM[1] = np.hstack((SStauM[1],SStauM[1]))988 994 modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) 989 995 FF = np.zeros(len(Tdata)) … … 995 1001 Uij = np.array(G2lat.U6toUij(Uijdata)) 996 1002 bij = Mast*Uij.T 1003 blkSize = 100 #no. of reflections in a block 1004 nRef = refDict['RefList'].shape[0] 997 1005 if not len(refDict['FF']): 998 1006 if 'N' in calcControls[hfx+'histType']: 999 1007 dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's 1008 refDict['FF']['El'] = dat.keys() 1009 refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values() 1000 1010 else: 1001 1011 dat = G2el.getFFvalues(FFtables,0.) 1002 refDict['FF']['El'] = dat.keys() 1003 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1012 refDict['FF']['El'] = dat.keys() 1013 refDict['FF']['FF'] = np.ones((nRef,len(dat))) 1014 for iref,ref in enumerate(refDict['RefList']): 1015 SQ = 1./(2.*ref[4])**2 1016 dat = G2el.getFFvalues(FFtables,SQ) 1017 refDict['FF']['FF'][iref] *= dat.values() 1004 1018 time0 = time.time() 1005 1019 nref = len(refDict['RefList'])/100 1006 for iref,refl in enumerate(refDict['RefList']): 1007 if 'NT' in calcControls[hfx+'histType']: 1008 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im]) 1009 fbs = np.array([0,0]) 1010 H = refl[:4] 1011 SQ = 1./(2.*refl[4+im])**2 1012 SQfactor = 4.0*SQ*twopisq 1013 Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor) 1014 if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor 1015 if 'N' in calcControls[hfx+'histType']: 1016 dat = G2el.getBLvalues(BLtables) 1017 refDict['FF']['FF'][iref] = dat.values() 1018 else: #'X' 1019 dat = G2el.getFFvalues(FFtables,SQ) 1020 refDict['FF']['FF'][iref] = dat.values() 1020 #reflection processing begins here - big arrays! 1021 iBeg = 0 1022 while iBeg < nRef: 1023 iFin = min(iBeg+blkSize,nRef) 1024 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1025 H = refl.T[:4] #array(blkSize,4) 1026 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,3) or (blkSize,3) 1027 TwMask = np.any(H,axis=-1) 1028 if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i] 1029 for ir in range(blkSize): 1030 iref = ir+iBeg 1031 if iref in TwDict: 1032 for i in TwDict[iref]: 1033 for n in range(NTL): 1034 H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1035 TwMask = np.any(H,axis=-1) 1036 SQ = 1./(2.*refl.T[5])**2 #array(blkSize) 1037 SQfactor = 4.0*SQ*twopisq #ditto prev. 1038 if 'T' in calcControls[hfx+'histType']: 1039 if 'P' in calcControls[hfx+'histType']: 1040 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14]) 1041 else: 1042 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 1043 FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0) 1044 FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0) 1045 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw)) 1021 1046 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1022 FF = refDict['FF']['FF'][iref][Tindx] 1023 HM = np.zeros(4) 1024 HM[:3] += H[3]*modQ 1025 HM += H 1026 Uniq = np.inner(HM[:3],SGMT) 1027 SSUniq = np.inner(H,SSGMT) 1028 Phi = np.inner(HM[:3],SGT) 1029 SSPhi = np.inner(H,SSGT) 1047 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0) 1048 Uniq = np.inner(H,SSGMT) 1049 Phi = np.inner(H,SSGT) 1030 1050 if SGInv: #if centro - expand HKL sets 1031 1051 Uniq = np.vstack((Uniq,-Uniq)) 1032 SSUniq = np.vstack((SSUniq,-SSUniq))1033 1052 Phi = np.hstack((Phi,-Phi)) 1034 SSPhi = np.hstack((SSPhi,-SSPhi))1035 1053 # GSASIIpath.IPyBreak() 1036 GfpuA = G2mth.Modulation(waveTypes, SSUniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms1037 phase = twopi*(np.inner(Uniq ,(dXdata.T+Xdata.T))+SSPhi[:,np.newaxis])1054 GfpuA = G2mth.Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms 1055 phase = twopi*(np.inner(Uniq[:,:3],(dXdata.T+Xdata.T))+Phi[:,np.newaxis]) 1038 1056 sinp = np.sin(phase) 1039 1057 cosp = np.cos(phase) 1040 1058 biso = -SQfactor*Uisodata 1041 1059 Tiso = np.where(biso<1.,np.exp(biso),1.0) 1042 HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq ])1060 HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq[:,:3]]) 1043 1061 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1044 1062 Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq) 1045 1063 fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr]) #2 x sym x atoms 1046 fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr]) #swapped around - better ?1064 fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr]) #swapped around - better 1047 1065 fa *= GfpuA 1048 1066 fb *= GfpuA
Note: See TracChangeset
for help on using the changeset viewer.