Changeset 655


Ignore:
Timestamp:
Jun 25, 2012 3:57:48 PM (9 years ago)
Author:
vondreele
Message:

implement single peak fits for background

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r643 r655  
    20202020        parmDict.update(histDict)
    20212021        for parm in parmDict:
    2022             if parm.split(':')[-1] in ['Azimuth','Gonio. radius','Lam1','Lam2','Omega','Chi','Phi']:
     2022            if parm.split(':')[-1] in ['Azimuth','Gonio. radius','Lam1','Lam2',
     2023                'Omega','Chi','Phi','nDebye','nPeaks']:
    20232024                parmDict[parm] = [parmDict[parm],' ']
    20242025            elif parm.split(':')[-2] in ['Ax','Ay','Az','SHmodel','SHord']:
  • trunk/GSASIImath.py

    r654 r655  
    674674   
    675675    result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
    676 #    for item in zip(DH,Dphi,result[2]['fvec']):
    677 #        print item[0],'%.4f %.4f'%(item[1],item[2])
     676    for item in zip(DH,Dphi,result[2]['fvec']):
     677        print item[0],'%.4f %.4f'%(item[1],item[2])
    678678    chisq = np.sum(result[2]['fvec']**2)
    679679    DX = np.array(np.fix(-result[0]*steps),dtype='i')
  • trunk/GSASIIpwd.py

    r650 r655  
    671671            iD += 1       
    672672        except KeyError:
    673             break   
     673            break
     674    iD = 0
     675    while True:
     676        try:
     677            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
     678            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
     679            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
     680            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
     681            shl = 0.002
     682            Wd,fmin,fmax = getWidths(pkP,pkS,pkG,shl)
     683            iBeg = np.searchsorted(xdata,pkP-fmin)
     684            iFin = np.searchsorted(xdata,pkP+fmax)
     685            yb[iBeg:iFin] += pkI*getFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
     686            iD += 1       
     687        except KeyError:
     688            break       
    674689    return yb
    675690   
     
    684699    dydb = np.zeros(shape=(nBak,len(xdata)))
    685700    dyddb = np.zeros(shape=(3*parmDict[pfx+'nDebye'],len(xdata)))
     701    dydpk = np.zeros(shape=(4*parmDict[pfx+'nPeaks'],len(xdata)))
     702    dx = xdata[1]-xdata[0]
    686703
    687704    if bakType in ['chebyschev','cosine']:
     
    735752            cqr = np.cos(q*dbR)
    736753            temp = np.exp(-dbU*q**2)
    737             dyddb[3*iD] = ff*sqr*temp
    738             dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/dbR
    739             dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2
     754            dyddb[3*iD] = 100.*dx*ff*sqr*temp
     755            dyddb[3*iD+1] = 100.*dx*ff*dbA*temp*(cqr-sqr)/dbR
     756            dyddb[3*iD+2] = -100.*dx*ff*dbA*sqr*temp*q**2
    740757            iD += 1       
    741758        except KeyError:
    742759            break
    743     return dydb,dyddb
     760    iD = 0
     761    while True:
     762        try:
     763            pkP = parmDict[pfx+'BkPkpos:'+str(iD)]
     764            pkI = parmDict[pfx+'BkPkint:'+str(iD)]
     765            pkS = parmDict[pfx+'BkPksig:'+str(iD)]
     766            pkG = parmDict[pfx+'BkPkgam:'+str(iD)]
     767            shl = 0.002
     768            Wd,fmin,fmax = getWidths(pkP,pkS,pkG,shl)
     769            iBeg = np.searchsorted(xdata,pkP-fmin)
     770            iFin = np.searchsorted(xdata,pkP+fmax)
     771            Df,dFdp,dFds,dFdg,dFdsh = getdFCJVoigt3(pkP,pkS,pkG,shl,xdata[iBeg:iFin])
     772            dydpk[4*iD][iBeg:iFin] += 100.*dx*pkI*dFdp
     773            dydpk[4*iD+1][iBeg:iFin] += 100.*dx*Df
     774            dydpk[4*iD+2][iBeg:iFin] += 100.*dx*pkI*dFds
     775            dydpk[4*iD+3][iBeg:iFin] += 100.*dx*pkI*dFdg
     776            iD += 1       
     777        except KeyError:
     778            break
     779    return dydb,dyddb,dydpk
    744780
    745781#use old fortran routine
     
    839875# needs to return np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order
    840876    dMdv = np.zeros(shape=(len(varyList),len(xdata)))
    841     dMdb,dMddb = getBackgroundDerv('',parmDict,bakType,xdata)
     877    dMdb,dMddb,dMdpk = getBackgroundDerv('',parmDict,bakType,xdata)
    842878    if 'Back:0' in varyList:            #background derivs are in front if present
    843879        dMdv[0:len(dMdb)] = dMdb
     
    848884            ip = names.index(parm)
    849885            dMdv[varyList.index(name)] = dMddb[3*int(id)+ip]
     886    names = ['BkPkpos','BkPkint','BkPksig','BkPkgam']
     887    for name in varyList:
     888        if 'BkPk' in name:
     889            parm,id = name.split(':')
     890            ip = names.index(parm)
     891            dMdv[varyList.index(name)] = dMdpk[4*int(id)+ip]
    850892    dx = xdata[1]-xdata[0]
    851893    U = parmDict['U']
     
    9801022        backVary += debyeVary
    9811023   
     1024        backDict['nPeaks'] = Debye['nPeaks']
    9821025        peaksDict = {}
    9831026        peaksList = []
  • trunk/GSASIIstruct.py

    r650 r655  
    14951495    def GetBackgroundParms(hId,Background):
    14961496        Back = Background[0]
    1497         Debye = Background[1]
     1497        DebyePeaks = Background[1]
    14981498        bakType,bakFlag = Back[:2]
    14991499        backVals = Back[3:]
     
    15031503        if bakFlag:
    15041504            backVary = backNames
    1505         backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
     1505        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
     1506        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
    15061507        debyeDict = {}
    15071508        debyeList = []
    1508         for i in range(Debye['nDebye']):
     1509        for i in range(DebyePeaks['nDebye']):
    15091510            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
    1510             debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
    1511             debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
     1511            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
     1512            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
    15121513        debyeVary = []
    15131514        for item in debyeList:
     
    15151516                debyeVary.append(item[0])
    15161517        backDict.update(debyeDict)
    1517         backVary += debyeVary   
     1518        backVary += debyeVary
     1519        peakDict = {}
     1520        peakList = []
     1521        for i in range(DebyePeaks['nPeaks']):
     1522            peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
     1523                ':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
     1524            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
     1525            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
     1526        peakVary = []
     1527        for item in peakList:
     1528            if item[1]:
     1529                peakVary.append(item[0])
     1530        backDict.update(peakDict)
     1531        backVary += peakVary
    15181532        return bakType,backDict,backVary           
    15191533       
     
    15541568    def PrintBackground(Background):
    15551569        Back = Background[0]
    1556         Debye = Background[1]
     1570        DebyePeaks = Background[1]
    15571571        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
    15581572        line = ' Coefficients: '
     
    15621576                line += '\n'+15*' '
    15631577        print line
    1564         if Debye['nDebye']:
     1578        if DebyePeaks['nDebye']:
    15651579            print '\n Debye diffuse scattering coefficients'
    15661580            parms = ['DebyeA','DebyeR','DebyeU']
    1567             line = ' names :'
     1581            line = ' names :  '
    15681582            for parm in parms:
    1569                 line += '%16s'%(parm)
     1583                line += '%8s refine?'%(parm)
    15701584            print line
    15711585            for j,term in enumerate(Debye['debyeTerms']):
     
    15731587                for i in range(3):
    15741588                    line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1]))                   
     1589                print line
     1590        if DebyePeaks['nPeaks']:
     1591            print '\n Single peak coefficients'
     1592            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
     1593            line = ' names :  '
     1594            for parm in parms:
     1595                line += '%8s refine?'%(parm)
     1596            print line
     1597            for j,term in enumerate(DebyePeaks['peaksList']):
     1598                line = ' peak'+'%2d'%(j)+':'
     1599                for i in range(4):
     1600                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
    15751601                print line
    15761602       
     
    16751701    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
    16761702        Back = Background[0]
    1677         Debye = Background[1]
     1703        DebyePeaks = Background[1]
    16781704        lenBack = len(Back[3:])
    1679         backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
     1705        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
    16801706        for i in range(lenBack):
    16811707            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
    16821708            if pfx+'Back:'+str(i) in sigDict:
    16831709                backSig[i] = sigDict[pfx+'Back:'+str(i)]
    1684         if Debye['nDebye']:
    1685             for i in range(Debye['nDebye']):
     1710        if DebyePeaks['nDebye']:
     1711            for i in range(DebyePeaks['nDebye']):
    16861712                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
    16871713                for j,name in enumerate(names):
    1688                     Debye['debyeTerms'][i][2*j] = parmDict[name]
     1714                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
    16891715                    if name in sigDict:
    16901716                        backSig[lenBack+3*i+j] = sigDict[name]           
     1717        if DebyePeaks['nPeaks']:
     1718            for i in range(DebyePeaks['nPeaks']):
     1719                names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
     1720                    pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
     1721                for j,name in enumerate(names):
     1722                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
     1723                    if name in sigDict:
     1724                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
    16911725        return backSig
    16921726       
     
    17181752    def PrintBackgroundSig(Background,backSig):
    17191753        Back = Background[0]
    1720         Debye = Background[1]
     1754        DebyePeaks = Background[1]
    17211755        lenBack = len(Back[3:])
    17221756        valstr = ' value : '
     
    17341768            print valstr
    17351769            print sigstr
    1736         if Debye['nDebye']:
     1770        if DebyePeaks['nDebye']:
    17371771            ifAny = False
    17381772            ptfmt = "%12.5f"
     
    17481782            if ifAny:
    17491783                print '\n Debye diffuse scattering coefficients'
     1784                print names
     1785                print ptstr
     1786                print sigstr
     1787        if DebyePeaks['nPeaks']:
     1788            ifAny = False
     1789            ptfmt = "%14.3f"
     1790            names =  ' names :'
     1791            ptstr =  ' values:'
     1792            sigstr = ' esds  :'
     1793            for item in sigDict:
     1794                if 'BkPk' in item:
     1795                    ifAny = True
     1796                    names += '%14s'%(item)
     1797                    ptstr += ptfmt%(parmDict[item])
     1798                    sigstr += ptfmt%(sigDict[item])
     1799            if ifAny:
     1800                print '\n Single peak coefficients'
    17501801                print names
    17511802                print ptstr
     
    26382689    bakType = calcControls[hfx+'bakType']
    26392690    dMdv = np.zeros(shape=(len(varylist),len(x)))
    2640     dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
     2691    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
    26412692    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
    26422693        bBpos =varylist.index(hfx+'Back:0')
     
    26492700            ip = names.index(parm)
    26502701            dMdv[varylist.index(name)] = dMddb[3*id+ip]
     2702    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
     2703    for name in varylist:
     2704        if 'BkPk' in name:
     2705            id = int(name.split(':')[-1])
     2706            parm = name[:int(name.rindex(':'))]
     2707            ip = names.index(parm)
     2708            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
    26512709    if 'C' in calcControls[hfx+'histType']:   
    26522710        dx = x[1]-x[0]
Note: See TracChangeset for help on using the changeset viewer.