Changeset 4919 for trunk/GSASIIpwd.py
 Timestamp:
 Jun 3, 2021 10:12:41 AM (5 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIpwd.py
r4915 r4919 718 718 """) 719 719 720 721 #GSASII peak fitting routine: Finger, Cox & Jephcoat model722 723 720 724 721 class fcjde_gen(st.rv_continuous): … … 741 738 fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}1)  1/s]/cos(T) 742 739 743 * if x >= 0: fcj.pdf = 0 740 * if x >= 0: fcj.pdf = 0 741 744 742 """ 745 743 def _pdf(self,x,t,s,dx): … … 759 757 fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') 760 758 759 def getFCJVoigt(pos,intens,sig,gam,shl,xdata): 760 '''Compute the FingerCoxJepcoat modified Voigt function for a 761 CW powder peak by direct convolution. This version is not used. 762 ''' 763 DX = xdata[1]xdata[0] 764 widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl) 765 x = np.linspace(posfmin,pos+fmin,256) 766 dx = x[1]x[0] 767 Norm = norm.pdf(x,loc=pos,scale=widths[0]) 768 Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1]) 769 arg = [pos,shl/57.2958,dx,] 770 FCJ = fcjde.pdf(x,*arg,loc=pos) 771 if len(np.nonzero(FCJ)[0])>5: 772 z = np.column_stack([Norm,Cauchy,FCJ]).T 773 Z = fft.fft(z) 774 Df = fft.ifft(Z.prod(axis=0)).real 775 else: 776 z = np.column_stack([Norm,Cauchy]).T 777 Z = fft.fft(z) 778 Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real 779 Df /= np.sum(Df) 780 Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0) 781 return intens*Df(xdata)*DX/dx 782 783 #### GSASII peak fitting routine: Finger, Cox & Jephcoat model 784 761 785 def getWidthsCW(pos,sig,gam,shl): 762 786 '''Compute the peak widths used for computing the range of a peak 763 787 for constant wavelength data. On lowangle side, 50 FWHM are used, 764 on highangle side 75 are used, lowangle side extended for axial divergence788 on highangle side 75 are used, high angle side extended for axial divergence 765 789 (for peaks above 90 deg, these are reversed.) 790 791 param pos: peak position; 2theta in degrees 792 param sig: Gaussian peak variance in centideg^2 793 param gam: Lorentzian peak width in centidegrees 794 param shl: axial divergence parameter (S+H)/ 795 796 returns: widths; [Gaussian sigma, Lornetzian gamma] in degrees 797 returns: low angle, high angle ends of peak; 20FWHM & 50 FWHM from position 798 reverset for 2theta > 90 deg. 766 799 ''' 767 800 widths = [np.sqrt(sig)/100.,gam/100.] … … 777 810 for constant wavelength data. 50 FWHM are used on both sides each 778 811 extended by exponential coeff. 812 813 param pos: peak position; TOF in musec 814 param alp,bet: TOF peak exponential rise & decay parameters 815 param sig: Gaussian peak variance in musec^2 816 param gam: Lorentzian peak width in musec 817 818 returns: widths; [Gaussian sigma, Lornetzian gamma] in musec 819 returns: low TOF, high TOF ends of peak; 50FWHM from position 779 820 ''' 780 821 widths = [np.sqrt(sig),gam] … … 832 873 return gamFW(2.35482*s,g) #sqrt(8ln2)*sig = FWHM(G) 833 874 834 def getFCJVoigt(pos,intens,sig,gam,shl,xdata):835 '''Compute the FingerCoxJepcoat modified Voigt function for a836 CW powder peak by direct convolution. This version is not used.837 '''838 DX = xdata[1]xdata[0]839 widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl)840 x = np.linspace(posfmin,pos+fmin,256)841 dx = x[1]x[0]842 Norm = norm.pdf(x,loc=pos,scale=widths[0])843 Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])844 arg = [pos,shl/57.2958,dx,]845 FCJ = fcjde.pdf(x,*arg,loc=pos)846 if len(np.nonzero(FCJ)[0])>5:847 z = np.column_stack([Norm,Cauchy,FCJ]).T848 Z = fft.fft(z)849 Df = fft.ifft(Z.prod(axis=0)).real850 else:851 z = np.column_stack([Norm,Cauchy]).T852 Z = fft.fft(z)853 Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real854 Df /= np.sum(Df)855 Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)856 return intens*Df(xdata)*DX/dx857 858 875 def getBackground(pfx,parmDict,bakType,dataType,xdata,fixback=None): 859 876 '''Computes the background from vars pulled from gpx file or tree. … … 970 987 iFin = np.searchsorted(xdata,pkP+fmax) 971 988 if 'C' in dataType: 972 ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin]) 973 yb[iBeg:iFin] += ybi 989 ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])[0] 990 yb[iBeg:iFin] += ybi/cw[iBeg:iFin] 974 991 elif 'T' in dataType: 975 ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin]) 992 ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])[0] 976 993 yb[iBeg:iFin] += ybi 977 994 elif 'B' in dataType: 978 ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin]) 995 ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin])[0] 979 996 yb[iBeg:iFin] += ybi 980 997 sumBk[2] += np.sum(ybi) … … 1130 1147 '''Compute the FingerCoxJepcoat modified PseudoVoigt function for a 1131 1148 CW powder peak in external Fortran routine 1149 1150 param pos: peak position in degrees 1151 param sig: Gaussian variance in centideg^2 1152 param gam: Lorentzian width in centideg 1153 param shl: axial divergence parameter (S+H)/L 1154 param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv) 1155 1156 returns: array: calculated peak function at each xdata 1157 returns: integral of peak; nominally = 1.0 1132 1158 ''' 1133 Df = pyd.pypsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) 1134 Df /= np.sum(Df) 1135 return Df 1159 if len(xdata): 1160 cw = np.diff(xdata) 1161 cw = np.append(cw,cw[1]) 1162 Df = pyd.pypsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) 1163 return Df,np.sum(100.*Df*cw) 1164 else: 1165 return 0.,1. 1136 1166 1137 1167 def getdFCJVoigt3(pos,sig,gam,shl,xdata): 1138 1168 '''Compute analytic derivatives the FingerCoxJepcoat modified PseudoVoigt 1139 1169 function for a CW powder peak 1170 1171 param pos: peak position in degrees 1172 param sig: Gaussian variance in centideg^2 1173 param gam: Lorentzian width in centideg 1174 param shl: axial divergence parameter (S+H)/L 1175 param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv) 1176 1177 returns: arrays: function and derivatives of pos, sig, gam, & shl 1140 1178 ''' 1141 1179 Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdatapos,pos,sig,gam,shl) … … 1145 1183 '''Compute the simple PseudoVoigt function for a 1146 1184 small angle Bragg peak in external Fortran routine 1185 1186 param pos: peak position in degrees 1187 param sig: Gaussian variance in centideg^2 1188 param gam: Lorentzian width in centideg 1189 1190 returns: array: calculated peak function at each xdata 1191 returns: integral of peak; nominally = 1.0 1147 1192 ''' 1148 1193 1194 cw = np.diff(xdata) 1195 cw = np.append(cw,cw[1]) 1149 1196 Df = pyd.pypsvoigt(len(xdata),xdatapos,sig,gam) 1150 Df /= np.sum(Df) 1151 return Df 1197 return Df,np.sum(100.*Df*cw) 1152 1198 1153 1199 def getdPsVoigt(pos,sig,gam,xdata): 1154 1200 '''Compute the simple PseudoVoigt function derivatives for a 1155 1201 small angle Bragg peak peak in external Fortran routine 1202 1203 param pos: peak position in degrees 1204 param sig: Gaussian variance in centideg^2 1205 param gam: Lorentzian width in centideg 1206 1207 returns: arrays: function and derivatives of pos, sig, gam, & shl 1156 1208 ''' 1157 1209 … … 1164 1216 ''' 1165 1217 1218 cw = np.diff(xdata) 1219 cw = np.append(cw,cw[1]) 1166 1220 Df = pyd.pyepsvoigt(len(xdata),xdatapos,alp,bet,sig,gam) 1167 Df /= np.sum(Df) 1168 return Df 1221 return Df,np.sum(Df*cw) 1169 1222 1170 1223 def getdEpsVoigt(pos,alp,bet,sig,gam,xdata): … … 1316 1369 1317 1370 def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType): 1318 'Computes the profile for a powder pattern' 1371 '''Computes the profile for a powder pattern for single peak fitting 1372 NB: not used for Rietveld refinement 1373 ''' 1319 1374 1320 1375 yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0] 1321 1376 yc = np.zeros_like(yb) 1322 cw = np.diff(xdata)1323 cw = np.append(cw,cw[1])1377 # cw = np.diff(xdata) 1378 # cw = np.append(cw,cw[1]) 1324 1379 if 'C' in dataType: 1325 1380 shl = max(parmDict['SH/L'],0.002) … … 1355 1410 elif not iBegiFin: #peak above high limit 1356 1411 return yb+yc 1357 yc[iBeg:iFin] += intens*getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) 1412 fp = getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])[0] 1413 yc[iBeg:iFin] += intens*fp 1358 1414 if Ka2: 1359 1415 pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) … … 1361 1417 iFin = np.searchsorted(xdata,pos2+fmin) 1362 1418 if iBegiFin: 1363 yc[iBeg:iFin] += intens*kRatio*getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) 1419 fp2 = getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])[0] 1420 yc[iBeg:iFin] += intens*kRatio*fp2 1364 1421 iPeak += 1 1365 1422 except KeyError: #no more peaks to process … … 1405 1462 elif not iBegiFin: #peak above high limit 1406 1463 return yb+yc 1407 yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin]) 1464 yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])[0] 1408 1465 iPeak += 1 1409 1466 except KeyError: #no more peaks to process … … 1463 1520 elif not iBegiFin: #peak above high limit 1464 1521 return yb+yc 1465 yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin]) 1522 yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])[0] 1466 1523 iPeak += 1 1467 1524 except KeyError: #no more peaks to process … … 1469 1526 1470 1527 def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType): 1471 'needs a doc string' 1472 # needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order 1528 '''Computes the profile derivatives for a powder pattern for single peak fitting 1529 1530 return: np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order 1531 1532 NB: not used for Rietveld refinement 1533 ''' 1473 1534 dMdv = np.zeros(shape=(len(varyList),len(xdata))) 1474 1535 dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback) … … 1529 1590 dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) 1530 1591 for i in range(1,5): 1531 dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*dMdipk[i]1532 dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk[0]1592 dMdpk[i][iBeg:iFin] += intens*dMdipk[i] 1593 dMdpk[0][iBeg:iFin] += dMdipk[0] 1533 1594 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} 1534 1595 if Ka2: … … 1539 1600 dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) 1540 1601 for i in range(1,5): 1541 dMdpk[i][iBeg:iFin] += 100.*cw[iBeg:iFin]*intens*kRatio*dMdipk2[i]1542 dMdpk[0][iBeg:iFin] += 100.*cw[iBeg:iFin]*kRatio*dMdipk2[0]1543 dMdpk[5][iBeg:iFin] += 100.*cw[iBeg:iFin]*dMdipk2[0]1602 dMdpk[i][iBeg:iFin] += intens*kRatio*dMdipk2[i] 1603 dMdpk[0][iBeg:iFin] += kRatio*dMdipk2[0] 1604 dMdpk[5][iBeg:iFin] += dMdipk2[0] 1544 1605 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens} 1545 1606 for parmName in ['pos','int','sig','gam']: … … 2194 2255 yb *= 0. 2195 2256 yd *= 0. 2196 cw = x[1:]x[:1]2197 2257 xBeg = np.searchsorted(x,Limits[0]) 2198 2258 xFin = np.searchsorted(x,Limits[1])+1 … … 2266 2326 FWHM = getFWHM(peak[0],Inst) 2267 2327 try: 2268 binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])]) 2328 xpk = x.searchsorted(peak[0]) 2329 cw = x[xpk]x[xpk1] 2330 binsperFWHM.append(FWHM/cw) 2269 2331 except IndexError: 2270 2332 binsperFWHM.append(0.) … … 4244 4306 hplot = plotter.add('derivatives test for '+name).gca() 4245 4307 hplot.plot(xdata,100.*dx*getdFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[idx+1]) 4246 y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata) 4308 y0 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[0] 4247 4309 myDict[name] += delt 4248 y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata) 4310 y1 = getFCJVoigt3(myDict['pos'],myDict['sig'],myDict['gam'],myDict['shl'],xdata)[0] 4249 4311 hplot.plot(xdata,(y1y0)/delt,'r+') 4250 4312
Note: See TracChangeset
for help on using the changeset viewer.