Changeset 376
 Timestamp:
 Sep 20, 2011 1:42:17 PM (11 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r375 r376 128 128 129 129 def Polarization(Pola,Tth,Azm=0.0): 130 # Calculate angle dependent xray polarization correction (not scaled correctly!) 131 # Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized 132 # Azm: azimuthal angle e.g. 0.0 in plane of polarization 133 # Tth: 2theta scattering angle  can be numpy array 134 # which (if either) of these is "right"? 135 # return (Pola*npcosd(Azm)**2+(1.Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \ 136 # Pola*npsind(Azm)**2+(1.Pola)*npcosd(Azm)**2 130 """ Calculate angle dependent xray polarization correction (not scaled correctly!) 131 input: 132 Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized 133 Azm: azimuthal angle e.g. 0.0 in plane of polarization 134 Tth: 2theta scattering angle  can be numpy array 135 which (if either) of these is "right"? 136 return: 137 pola = (Pola*npcosd(Azm)**2+(1.Pola)*npsind(Azm)**2)*npcosd(Tth)**2+ \ 138 Pola*npsind(Azm)**2+(1.Pola)*npcosd(Azm)**2 139 dpdPola: derivative needed for least squares 140 """ 137 141 pola = ((1.0Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ 138 142 (1.0Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2 139 143 dpdPola = npsind(Tth)**2*(npsind(Azm)**2npcosd(Azm)**2) 140 return pola,dpdPola /pola144 return pola,dpdPola 141 145 142 146 def Oblique(ObCoeff,Tth): 
trunk/GSASIIstruct.py
r375 r376 169 169 Histograms = dictionary of histograms as {name:data,...} 170 170 Phases = dictionary of phases that use histograms 171 CovData = [varyList,covariance matrix]171 CovData = dictionary of refined variables, varyList, & covariance matrix 172 172 ''' 173 173 … … 205 205 data[iphase][1] = Phases[phaseName] 206 206 elif datum[0] == 'Covariance': 207 varyList,covMatrix = CovData 208 data[0][1] = {'varyList':varyList,'covariance':covMatrix} 207 data[0][1] = CovData 209 208 try: 210 209 histogram = Histograms[datum[0]] … … 1125 1124 parmdict.update(zip(varylist,values)) 1126 1125 1126 def GetPrefOri(refl,G,phfx,calcControls,parmDict): 1127 if calcControls[phfx+'poType'] == 'MD': 1128 MD = parmDict[phfx+'MD'] 1129 MDAxis = calcControls[phfx+'MDAxis'] 1130 sumMD = 0 1131 for H in refl[10]: 1132 cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) 1133 A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) 1134 sumMD += A**3 1135 POcorr = sumMD/len(refl[10]) 1136 else: #spherical harmonics 1137 POcorr = 1.0 1138 return POcorr 1139 1140 def GetPrefOriDerv(refl,G,phfx,calcControls,parmDict): 1141 POderv = {} 1142 if calcControls[phfx+'poType'] == 'MD': 1143 MD = parmDict[phfx+'MD'] 1144 MDAxis = calcControls[phfx+'MDAxis'] 1145 sumMD = 0 1146 sumdMD = 0 1147 for H in refl[10]: 1148 cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) 1149 A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) 1150 sumMD += A**3 1151 sumdMD = (1.5*A**5)*(2.0*MD*cosP**2(sinP/MD)**2) 1152 POcorr = sumMD/len(refl[10]) 1153 POderv[phfx+'MD'] = sumdMD/len(refl[10]) 1154 else: #spherical harmonics 1155 POcorr = 1.0 1156 return POcorr,POderv 1157 1158 def GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict): 1159 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity 1160 Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0] 1161 Icorr *= GetPrefOri(refl,G,phfx,calcControls,parmDict) 1162 return Icorr 1163 1164 def GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict): 1165 Icorr = GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict) 1166 pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth']) 1167 POcorr,dIdPO = GetPrefOriDerv(refl,G,phfx,calcControls,parmDict) 1168 dIdPola *= refl[3]/pola 1169 for iPO in dIdPO: 1170 dIdPO[iPO] *= refl[3]/POcorr 1171 dIdsh = refl[3]/parmDict[hfx+'Scale'] 1172 dIdsp = refl[3]/parmDict[phfx+'Scale'] 1173 return Icorr,dIdsh,dIdsp,dIdPola,dIdPO 1174 1127 1175 def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): 1128 1176 … … 1158 1206 gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum 1159 1207 return gam 1160 1161 def GetMarchDollase(refl,G,phfx,calcControls,parmDict):1162 MD = parmDict[phfx+'MD']1163 MDAxis = calcControls[phfx+'MDAxis']1164 sumMD = 01165 for H in refl[10]:1166 cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)1167 sumMD += 1/np.sqrt((MD*cosP)**2+sinP**2/MD)**31168 return sumMD/len(refl[10])1169 1170 def GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict):1171 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity1172 Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]1173 if calcControls[phfx+'poType'] == 'MD':1174 Icorr *= GetMarchDollase(refl,G,phfx,calcControls,parmDict)1175 1176 return Icorr1177 1208 1178 1209 def GetHStrainShift(refl,SGData,phfx,parmDict): … … 1282 1313 continue 1283 1314 elif not iBegiFin: #peak above high limit  done 1284 return yc,yb1315 break 1285 1316 yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here 1286 1317 if Ka2: … … 1294 1325 return yc,yb 1295 1326 yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here 1296 else: 1297 raise ValueError 1327 elif 'T' in calcControls[hfx+'histType']: 1328 print 'TOF Undefined at present' 1329 raise ValueError #no TOF yet 1298 1330 return yc,yb 1299 1331 … … 1341 1373 gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] 1342 1374 return gamDict 1343 1344 def GetMarchDollaseDerv(refl,G,phfx,calcControls,parmDict):1345 MD = parmDict[phfx+'MD']1346 MDAxis = calcControls[phfx+'MDAxis']1347 sumMD = 01348 sumdMD = 01349 for H in refl[10]:1350 cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)1351 A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)1352 sumMD += A**31353 sumdMD = (1.5*A**5)*(2.0*MD*cosP**2(sinP/MD)**2)1354 return sumMD/len(refl[10]),refl[8]*sumdMD1355 1356 def GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict):1357 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity1358 pola,dpdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])1359 dIdpola = refl[8]*Icorr*dpdPola/refl[3]1360 Icorr *= pola1361 MDcorr,dIdMD = GetMarchDollaseDerv(refl,G,phfx,calcControls,parmDict)1362 Icorr *= MDcorr1363 dIdsh = Icorr/parmDict[hfx+'Scale']1364 dIdsp = Icorr/parmDict[phfx+'Scale']1365 dIdMD *= Icorr1366 return Icorr,dIdsh,dIdsp,dIdpola,dIdMD1367 1375 1368 1376 def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict): … … 1476 1484 sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)]) 1477 1485 for refl in refList: 1478 if 'C' in calcControls[hfx+'histType']: 1486 if 'C' in calcControls[hfx+'histType']: #CW powder 1479 1487 h,k,l = refl[:3] 1480 1488 iref = pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)] 1481 Icorr,dIdsh,dIdsp,dIdpola,dIdMD = GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict) 1482 hkl = refl[:3] 1489 Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict) 1490 if 'Pawley' in Phase['General']['Type']: 1491 try: 1492 refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) 1493 except KeyError: 1494 # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) 1495 continue 1496 else: 1497 raise ValueError #wants strctrfacr deriv calc here 1498 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) 1499 iBeg = np.searchsorted(x,refl[5]fmin) 1500 iFin = np.searchsorted(x,refl[5]+fmax) 1501 if not iBeg+iFin: #peak below low limit  skip peak 1502 continue 1503 elif not iBegiFin: #peak above high limit  done 1504 break 1483 1505 pos = refl[5] 1484 1506 tanth = tand(pos/2.0) 1485 1507 costh = cosd(pos/2.0) 1486 dsdU = tanth**21487 dsdV = tanth1488 dsdW = 1.01489 dgdX = 1.0/costh1490 dgdY = tanth1491 Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)1492 iBeg = np.searchsorted(x,refl[5]fmin)1493 iFin = np.searchsorted(x,refl[5]+fmax)1494 1508 dMdpk = np.zeros(shape=(6,len(x))) 1495 1509 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) 1496 1510 for i in range(1,5): 1497 1511 dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i] 1498 dMdpk[0][iBeg:iFin] += 100.*dx*Icorr* dMdipk[0]1512 dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[0]/refl[3] 1499 1513 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} 1500 1514 if Ka2: 1501 pos2 = refl[5]+lamRatio*tan d(refl[5]/2.0)# + 360/pi * Dlam/lam * tan(th)1515 pos2 = refl[5]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) 1502 1516 kdelt = int((pos2refl[5])/dx) 1503 iBeg = min(lenX,iBeg+kdelt)1504 iFin = min(lenX,iFin+kdelt)1505 if iBeg iFin:1506 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg :iFin])1517 iBeg2 = min(lenX,iBeg+kdelt) 1518 iFin2 = min(lenX,iFin+kdelt) 1519 if iBeg2iFin2: 1520 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) 1507 1521 for i in range(1,5): 1508 dMdpk[i][iBeg :iFin] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]1509 dMdpk[0][iBeg :iFin] += 100.*dx*kRatio*dMdipk2[0]1510 dMdpk[5][iBeg :iFin] += 100.*dx*dMdipk2[0]1511 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]* Icorr*refl[8]}1522 dMdpk[i][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i] 1523 dMdpk[0][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[0]/refl[3] 1524 dMdpk[5][iBeg2:iFin2] += 100.*dx*Icorr*dMdipk2[0] 1525 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[8]} 1512 1526 try: 1513 1527 idx = varylist.index(pfx+'PWLref:'+str(iref)) … … 1517 1531 dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) 1518 1532 names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], 1519 hfx+'U':[ dsdU,'sig'],hfx+'V':[dsdV,'sig'],hfx+'W':[dsdW,'sig'],1520 hfx+'X':[ dgdX,'gam'],hfx+'Y':[dgdY,'gam'],hfx+'SH/L':[1.0,'shl'],1533 hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], 1534 hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], 1521 1535 hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], 1522 1536 hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], 1523 hfx+'DisplaceY':[dpdY,'pos'], phfx+'MD':[dIdMD,'int'],}1537 hfx+'DisplaceY':[dpdY,'pos'],} 1524 1538 for name in names: 1525 1539 if name in varylist: 1526 1540 item = names[name] 1527 1541 dMdv[varylist.index(name)] += item[0]*dervDict[item[1]] 1542 for iPO in dIdPO: 1543 if iPO in varylist: 1544 dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int'] 1528 1545 cellDervNames = cellVaryDerv(pfx,SGData,dpdA) 1529 1546 for name,dpdA in cellDervNames: … … 1538 1555 if name in varylist: 1539 1556 dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam'] 1540 else: 1541 raise ValueError 1557 elif 'T' in calcControls[hfx+'histType']: 1558 print 'TOF Undefined at present' 1559 raise ValueError #no TOF yet 1542 1560 1543 1561 return dMdv … … 1686 1704 # print 'indParmList: ',G2mv.indParmList 1687 1705 # print 'fixedDict: ',G2mv.fixedDict 1688 covFile = ospath.splitext(GPXfile)[0]+'.gpxcov'1689 file = open(covFile,'wb')1690 cPickle.dump(varyList,file,1)1691 cPickle.dump(result[0],file,1)1692 cPickle.dump(covMatrix,file,1)1693 file.close()1694 1706 sigDict = dict(zip(varyList,sig)) 1695 1707 SetPhaseData(parmDict,sigDict,Phases) 1696 1708 SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms) 1697 1709 SetHistogramData(parmDict,sigDict,Histograms) 1698 SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,[varyList,cov]) 1710 covData = {'variables':result[0],'varyList':varyList,'covariance':cov} 1711 SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData) 1699 1712 #for testing purposes!!! 1700 1713 # file = open('structTestdata.dat','wb')
Note: See TracChangeset
for help on using the changeset viewer.