Changeset 4919 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Jun 3, 2021 10:12:41 AM (5 months ago)
Author:
vondreele
Message:

change powder peak shape functions to include effect of varying step size.
Impacts scale factors & peak intensities.
New notification system invoked to let users know of change.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r4915 r4919  
    718718""")
    719719   
    720    
    721 #GSASII peak fitting routine: Finger, Cox & Jephcoat model       
    722 
    723720
    724721class fcjde_gen(st.rv_continuous):
     
    741738        fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)|
    742739
    743      * if x >= 0: fcj.pdf = 0   
     740     * if x >= 0: fcj.pdf = 0   
     741     
    744742    """
    745743    def _pdf(self,x,t,s,dx):
     
    759757fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
    760758               
     759def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
     760    '''Compute the Finger-Cox-Jepcoat 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(pos-fmin,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
    761785def getWidthsCW(pos,sig,gam,shl):
    762786    '''Compute the peak widths used for computing the range of a peak
    763787    for constant wavelength data. On low-angle side, 50 FWHM are used,
    764     on high-angle side 75 are used, low angle side extended for axial divergence
     788    on high-angle side 75 are used, high angle side extended for axial divergence
    765789    (for peaks above 90 deg, these are reversed.)
     790   
     791    param pos: peak position; 2-theta 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 2-theta > 90 deg.
    766799    '''
    767800    widths = [np.sqrt(sig)/100.,gam/100.]
     
    777810    for constant wavelength data. 50 FWHM are used on both sides each
    778811    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
    779820    '''
    780821    widths = [np.sqrt(sig),gam]
     
    832873    return gamFW(2.35482*s,g)   #sqrt(8ln2)*sig = FWHM(G)
    833874               
    834 def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
    835     '''Compute the Finger-Cox-Jepcoat modified Voigt function for a
    836     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(pos-fmin,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]).T
    848         Z = fft.fft(z)
    849         Df = fft.ifft(Z.prod(axis=0)).real
    850     else:
    851         z = np.column_stack([Norm,Cauchy]).T
    852         Z = fft.fft(z)
    853         Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real
    854     Df /= np.sum(Df)
    855     Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
    856     return intens*Df(xdata)*DX/dx
    857 
    858875def getBackground(pfx,parmDict,bakType,dataType,xdata,fixback=None):
    859876    '''Computes the background from vars pulled from gpx file or tree.
     
    970987                iFin = np.searchsorted(xdata,pkP+fmax)
    971988            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]
    974991            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]
    976993                yb[iBeg:iFin] += ybi
    977994            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]
    979996                yb[iBeg:iFin] += ybi
    980997            sumBk[2] += np.sum(ybi)
     
    11301147    '''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a
    11311148    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
    11321158    '''
    1133     Df = pyd.pypsvfcj(len(xdata),xdata-pos,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),xdata-pos,pos,sig,gam,shl)
     1163        return Df,np.sum(100.*Df*cw)
     1164    else:
     1165        return 0.,1.
    11361166
    11371167def getdFCJVoigt3(pos,sig,gam,shl,xdata):
    11381168    '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt
    11391169    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
    11401178    '''
    11411179    Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
     
    11451183    '''Compute the simple Pseudo-Voigt function for a
    11461184    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
    11471192    '''
    11481193   
     1194    cw = np.diff(xdata)
     1195    cw = np.append(cw,cw[-1])
    11491196    Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam)
    1150     Df /= np.sum(Df)
    1151     return Df
     1197    return Df,np.sum(100.*Df*cw)
    11521198
    11531199def getdPsVoigt(pos,sig,gam,xdata):
    11541200    '''Compute the simple Pseudo-Voigt function derivatives for a
    11551201    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
    11561208    '''
    11571209   
     
    11641216    '''
    11651217   
     1218    cw = np.diff(xdata)
     1219    cw = np.append(cw,cw[-1])
    11661220    Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam)
    1167     Df /= np.sum(Df)
    1168     return Df 
     1221    return Df,np.sum(Df*cw) 
    11691222   
    11701223def getdEpsVoigt(pos,alp,bet,sig,gam,xdata):
     
    13161369
    13171370def 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    '''
    13191374   
    13201375    yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0]
    13211376    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])
    13241379    if 'C' in dataType:
    13251380        shl = max(parmDict['SH/L'],0.002)
     
    13551410                elif not iBeg-iFin:     #peak above high limit
    13561411                    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
    13581414                if Ka2:
    13591415                    pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     
    13611417                    iFin = np.searchsorted(xdata,pos2+fmin)
    13621418                    if iBeg-iFin:
    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
    13641421                iPeak += 1
    13651422            except KeyError:        #no more peaks to process
     
    14051462                elif not iBeg-iFin:     #peak above high limit
    14061463                    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]
    14081465                iPeak += 1
    14091466            except KeyError:        #no more peaks to process
     
    14631520                elif not iBeg-iFin:     #peak above high limit
    14641521                    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]
    14661523                iPeak += 1
    14671524            except KeyError:        #no more peaks to process
     
    14691526           
    14701527def 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    '''
    14731534    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
    14741535    dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback)
     
    15291590                dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])
    15301591                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]
    15331594                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
    15341595                if Ka2:
     
    15391600                        dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])
    15401601                        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]
    15441605                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens}
    15451606                for parmName in ['pos','int','sig','gam']:
     
    21942255    yb *= 0.
    21952256    yd *= 0.
    2196     cw = x[1:]-x[:-1]
    21972257    xBeg = np.searchsorted(x,Limits[0])
    21982258    xFin = np.searchsorted(x,Limits[1])+1
     
    22662326        FWHM = getFWHM(peak[0],Inst)
    22672327        try:
    2268             binsperFWHM.append(FWHM/cw[x.searchsorted(peak[0])])
     2328            xpk = x.searchsorted(peak[0])
     2329            cw = x[xpk]-x[xpk-1]
     2330            binsperFWHM.append(FWHM/cw)
    22692331        except IndexError:
    22702332            binsperFWHM.append(0.)
     
    42444306    hplot = plotter.add('derivatives test for '+name).gca()
    42454307    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]
    42474309    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]
    42494311    hplot.plot(xdata,(y1-y0)/delt,'r+')
    42504312
Note: See TracChangeset for help on using the changeset viewer.