Changeset 4519 for trunk/GSASIImath.py
 Timestamp:
 Jul 13, 2020 1:27:06 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r4514 r4519 1463 1463 psin = np.sin(twopi*phase) 1464 1464 pcos = np.cos(twopi*phase) 1465 MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3) #cos term1466 MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3) #sin term1465 MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)/2. #cos term 1466 MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)/2. #sin term 1467 1467 MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=1)*SGData['SpnFlp'][nxs,:,nxs,nxs] 1468 1468 MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=1)*SGData['SpnFlp'][nxs,:,nxs,nxs] 1469 1469 return MmodA,MmodB #Ntau,Nops,Natm,Mxyz; cos & sin parts; sum matches drawn atom moments 1470 1470 1471 def MagMod2(m, glTau,XYZ,modQ,MSSdata,SGData,SSGData):1471 def MagMod2(m,XYZ,modQ,MSSdata,SGData,SSGData): 1472 1472 ''' 1473 1473 this needs to make magnetic moment modulations & magnitudes as … … 1477 1477 Bm = np.array(MSSdata[:3]).T[:,0,:] #...sin pos mods 1478 1478 SGMT = np.array([ops[0] for ops in SGData['SGOps']]) #not .T!! 1479 SSGMT = np.array([ops[0] for ops in SSGData['SSGOps']]) #not .T!! 1479 1480 Sinv = np.array([nl.inv(ops[0]) for ops in SSGData['SSGOps']]) 1480 1481 SGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1481 1482 if SGData['SGInv']: 1482 1483 SGMT = np.vstack((SGMT,SGMT)) 1484 SSGMT = np.vstack((SSGMT,SSGMT)) 1483 1485 Sinv = np.vstack((Sinv,Sinv)) 1484 1486 SGT = np.vstack((SGT,SGT)) 1485 1487 SGMT = np.vstack([SGMT for cen in SGData['SGCen']]) 1488 SSGMT = np.vstack([SSGMT for cen in SGData['SGCen']]) 1486 1489 Sinv = np.vstack([Sinv for cen in SGData['SGCen']]) 1487 1490 SGT = np.vstack([SGT+cen for cen in SSGData['SSGCen']])%1. 1488 1491 if SGData['SGGray']: 1489 1492 SGMT = np.vstack((SGMT,SGMT)) 1493 SSGMT = np.vstack((SSGMT,SSGMT)) 1490 1494 Sinv = np.vstack((Sinv,Sinv)) 1491 1495 SGT = np.vstack((SGT,SGT+.5))%1. 1492 mst = Sinv[:,3,:3]1493 1496 epsinv = Sinv[:,3,3] 1494 phase = np.inner(XYZ,modQ).T+(np.inner(mst,modQ)epsinv)[:,nxs]+glTau 1495 1496 psin = np.sin(twopi*m*phase).T 1497 pcos = np.cos(twopi*m*phase).T 1498 MmodA = Am[nxs,nxs,:,:]*pcos[:,:,nxs,nxs] #cos term 1499 MmodB = Bm[nxs,nxs,:,:]*psin[:,:,nxs,nxs] #sin term 1500 MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=1)*SGData['SpnFlp'][nxs,:,nxs,nxs] 1501 MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=1)*SGData['SpnFlp'][nxs,:,nxs,nxs] 1502 return MmodA,MmodB #Nref,Ntau,Nops,Natm,Mxyz; cos & sin parts; sum matches drawn atom moments 1497 phi = np.inner(XYZ,modQ).T 1498 TA = phi+(epsinv*(np.inner(modQ,SGT[:,:3])SGT[:,3]))[:,nxs] #Nops,Natm 1499 phase = phi+(np.inner(modQ,SGT[:,:3])SGT[:,3])[:,nxs] 1500 1501 pcos = np.cos(twopi*m[:,nxs,nxs]*phase[nxs,:,:]) #Nref,Nops,Natm 1502 psin = np.sin(twopi*m[:,nxs,nxs]*phase[nxs,:,:]) 1503 MmodA = TA[nxs,:,:,nxs]*(Am[nxs,nxs,:,:]*pcos[:,:,:,nxs]Bm[nxs,nxs,:,:]*psin[:,:,:,nxs])/2. #Nref,Nops,Natm,Mxyz 1504 MmodB = TA[nxs,:,:,nxs]*(Am[nxs,nxs,:,:]*psin[:,:,:,nxs]+Bm[nxs,nxs,:,:]*pcos[:,:,:,nxs])/2. #Nref,Nops,Natm,Mxyz 1505 MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=1)*SGData['MagMom'][nxs,:,nxs,nxs] 1506 MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=1)*SGData['MagMom'][nxs,:,nxs,nxs] 1507 return MmodA,MmodB #Nref,Nops,Natm,Mxyz; cos & sin parts 1503 1508 1504 1509 def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt): … … 4163 4168 return 1./dsp 4164 4169 4170 def getPinkalpha(ins,tth): 4171 '''get TOF peak profile alpha 4172 4173 :param dict ins: instrument parameters with at least 'alpha' 4174 as values only 4175 :param float tth: 2theta of peak 4176 4177 :returns: flaot getPinkalpha: peak alpha 4178 4179 ''' 4180 return ins['alpha0']+ ins['alpha1']*tand(tth/2.) 4181 4182 def getPinkalphaDeriv(tth): 4183 '''get derivatives of TOF peak profile beta wrt alpha 4184 4185 :param float dsp: dspacing of peak 4186 4187 :returns: float getTOFalphaDeriv: d(alp)/d(alpha0), d(alp)/d(alpha1) 4188 4189 ''' 4190 return 1.0,tand(tth/2.) 4191 4192 def getPinkbeta(ins,tth): 4193 '''get TOF peak profile beta 4194 4195 :param dict ins: instrument parameters with at least 'beat0' & 'beta1' 4196 as values only 4197 :param float tth: 2theta of peak 4198 4199 :returns: float getaPinkbeta: peak beta 4200 4201 ''' 4202 return ins['beta0']+ins['beta1']*tand(tth/2.) 4203 4204 def getPinkbetaDeriv(tth): 4205 '''get derivatives of TOF peak profile beta wrt beta0 & beta1 4206 4207 :param float dsp: dspacing of peak 4208 4209 :returns: list getTOFbetaDeriv: d(beta)/d(beta0) & d(beta)/d(beta1) 4210 4211 ''' 4212 return 1.0,tand(tth/2.) 4213 4165 4214 def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False): 4166 4215 '''set starting peak parameters for single peak fits from plot selection or auto selection … … 4176 4225 for CW: [pos,0,mag,1,sig,0,gam,0] 4177 4226 for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] 4227 for Pink: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] 4178 4228 NB: mag refinement set by default, all others off 4179 4229 … … 4183 4233 ind = 1 4184 4234 ins = {} 4185 if 'C' in Parms['Type'][0]: #CW data  TOF later in an else 4186 for x in ['U','V','W','X','Y','Z']: 4187 ins[x] = Parms.get(x,[0.0,0.0])[ind] 4188 if ifQ: #qplot  convert back to 2theta 4189 pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi)) 4190 sig = getCWsig(ins,pos) 4191 gam = getCWgam(ins,pos) 4192 XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st 4193 else: 4235 if 'T' in Parms['Type'][0]: 4194 4236 if ifQ: 4195 4237 dsp = 2.*np.pi/pos … … 4211 4253 gam = getTOFgamma(ins,dsp) 4212 4254 XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] 4255 elif 'C' in Parms['Type'][0]: #CW data  TOF later in an else 4256 for x in ['U','V','W','X','Y','Z']: 4257 ins[x] = Parms.get(x,[0.0,0.0])[ind] 4258 if ifQ: #qplot  convert back to 2theta 4259 pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi)) 4260 sig = getCWsig(ins,pos) 4261 gam = getCWgam(ins,pos) 4262 XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st 4263 elif 'B' in Parms['Type'][0]: 4264 for x in ['U','V','W','X','Y','Z','alpha0','alpha1','beta0','beta1']: 4265 ins[x] = Parms.get(x,[0.0,0.0])[ind] 4266 if ifQ: #qplot  convert back to 2theta 4267 pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi)) 4268 alp = getPinkalpha(ins,pos) 4269 bet = getPinkbeta(ins,pos) 4270 sig = getCWsig(ins,pos) 4271 gam = getCWgam(ins,pos) 4272 XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0] #default refine intensity 1st 4213 4273 return XY 4214 4274
Note: See TracChangeset
for help on using the changeset viewer.