Changeset 2783 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Apr 15, 2017 11:31:45 AM (5 years ago)
Author:
vondreele
Message:

include imports/G2rfd_xye.py
implement instrument smearing of reflectometry calculated patterns
this can use selected dq/q value or dq values from data
add check for 2-theta version of chi file import

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r2777 r2783  
    6161npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave)    #=2pi*d*
    6262ateln2 = 8.0*math.log(2.0)
     63sateln2 = np.sqrt(ateln2)
    6364nxs = np.newaxis
    6465
     
    19381939        values = []
    19391940        bounds = []
    1940         parmDict['Res'] = data['Resolution'][0]/(100.*np.sqrt(ateln2))     #% FWHM-->decimal sig
     1941        parmDict['dQ type'] = data['dQ type']
     1942        parmDict['Res'] = data['Resolution'][0]/(100.*sateln2)     #% FWHM-->decimal sig
    19411943        for parm in ['Scale','FltBack']:
    19421944            parmDict[parm] = data[parm][0]
     
    19941996            print line2
    19951997   
    1996     def calcREFD(values,Q,Io,wt,parmDict,varyList):
     1998    def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
    19971999        parmDict.update(zip(varyList,values))
    1998         M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io)
     2000        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
    19992001        return M
    20002002   
    2001     def sumREFD(values,Q,Io,wt,parmDict,varyList):
     2003    def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList):
    20022004        parmDict.update(zip(varyList,values))
    2003         M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io)
     2005        M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io)
    20042006        return np.sum(M**2)
    20052007   
    2006     def getREFD(Q,parmDict):
     2008    def getREFD(Q,Qsig,parmDict):
    20072009        Ic = np.ones_like(Q)*parmDict['FltBack']
    20082010        Scale = parmDict['Scale']
     
    20252027            if cid+'Mag SLD' in parmDict:
    20262028                rho[ilay] += parmDict[cid+'Mag SLD']
    2027         A,B = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
    2028         Ic += (A**2+B**2)*Scale
    2029         #TODO: smear Ic by instrument resolution Qsig
     2029        if parmDict['dQ type'] == 'None':
     2030            AB = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
     2031        elif 'const' in parmDict['dQ type']:
     2032            AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:])
     2033        else:       #dQ/Q in data
     2034            AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:])
     2035        Ic += AB*Scale
    20302036        return Ic
    2031        
    2032         def Smear(f,w,z,dq):
    2033             y = f(w,z)
    2034             s = dq/ateln2
    2035             y += 0.1354*(f(w,z+2*s)+f(w,z-2*s))
    2036             y += 0.24935*(f(w,z-1.666667*s)+f(w,z+1.666667*s))
    2037             y += 0.4111*(f(w,z-1.333333*s)+f(w,z+1.333333*s))
    2038             y += 0.60653*(f(w,z-s) +f(w,z+s))
    2039             y += 0.80074*(f(w,z-0.6666667*s)+f(w,z+0.6666667*s))
    2040             y += 0.94596*(f(w,z-0.3333333*s)+f(w,z+0.3333333*s))
    2041             y *= 0.137023
    2042             return y
    20432037       
    20442038    def estimateT0(takestep):
     
    20472041        for i in range(100):
    20482042            x0 = takestep(values)
    2049             M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
     2043            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
    20502044            Mmin = min(M,Mmin)
    20512045            MMax = max(M,Mmax)
     
    20692063        if data['Minimizer'] == 'LMLS':
    20702064            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
    2071                 args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
     2065                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
    20722066            parmDict.update(zip(varyList,result[0]))
    20732067            chisq = np.sum(result[2]['fvec']**2)
     
    20822076            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
    20832077                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
    2084                 'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)})
     2078                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)})
    20852079            chisq = result.fun
    20862080            ncalc = result.nfev
     
    20892083        elif data['Minimizer'] == 'MC/SA Anneal':
    20902084            xyrng = np.array(bounds).T
    2091             result = G2mth.anneal(sumREFD, values, args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList),
     2085            result = G2mth.anneal(sumREFD, values,
     2086                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList),
    20922087                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
    20932088                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
     
    21012096        elif data['Minimizer'] == 'L-BFGS-B':
    21022097            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
    2103                 args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
     2098                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList))
    21042099            parmDict.update(zip(varyList,result['x']))
    21052100            chisq = result.fun
     
    21082103            covM = []
    21092104    else:   #nothing varied
    2110         M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
     2105        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)
    21112106        chisq = np.sum(M**2)
    21122107        ncalc = 0
     
    21182113    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
    21192114    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
    2120     Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],parmDict)
     2115    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict)
    21212116    Ib[Ibeg:Ifin] = parmDict['FltBack']
    21222117    try:
     
    22192214        if 'Mag SLD' in layer:
    22202215            rho[ilay] += layer['Mag SLD'][0]
    2221     A,B = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
    2222     Ic[iBeg:iFin] = (A**2+B**2)*Scale+Ib[iBeg:iFin]
     2216    if data['dQ type'] == 'None':
     2217        AB = abeles(0.5*Q[iBeg:iFin],depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
     2218    elif 'const' in data['dQ type']:
     2219        res = data['Resolution'][0]/(100.*sateln2)
     2220        AB = SmearAbeles(0.5*Q[iBeg:iFin],res*Q[iBeg:iFin],depth,rho,irho,sigma[1:])
     2221    else:       #dQ/Q in data
     2222        AB = SmearAbeles(0.5*Q[iBeg:iFin],Qsig[iBeg:iFin],depth,rho,irho,sigma[1:])
     2223    Ic[iBeg:iFin] = AB*Scale+Ib[iBeg:iFin]
    22232224
    22242225def abeles(kz, depth, rho, irho=0, sigma=0):
     
    22892290   
    22902291        r = B12/B11
    2291         return np.real(r),np.imag(r)
     2292        return np.absolute(r)**2
    22922293
    22932294    if np.isscalar(kz): kz = np.asarray([kz], 'd')
     
    23082309    return calc(kz, depth, rho, irho, sigma)
    23092310   
     2311def SmearAbeles(kz,dq, depth, rho, irho=0, sigma=0):
     2312    y = abeles(kz, depth, rho, irho, sigma)
     2313    s = dq/2.
     2314    y += 0.1354*(abeles(kz+2*s, depth, rho, irho, sigma)+abeles(kz-2*s, depth, rho, irho, sigma))
     2315    y += 0.24935*(abeles(kz-5*s/3., depth, rho, irho, sigma)+abeles(kz+5*s/3., depth, rho, irho, sigma))
     2316    y += 0.4111*(abeles(kz-4*s/3., depth, rho, irho, sigma)+abeles(kz+4*s/3., depth, rho, irho, sigma))
     2317    y += 0.60653*(abeles(kz-s, depth, rho, irho, sigma) +abeles(kz+s, depth, rho, irho, sigma))
     2318    y += 0.80074*(abeles(kz-2*s/3., depth, rho, irho, sigma)+abeles(kz-2*s/3., depth, rho, irho, sigma))
     2319    y += 0.94596*(abeles(kz-s/3., depth, rho, irho, sigma)+abeles(kz-s/3., depth, rho, irho, sigma))
     2320    y *= 0.137023
     2321    return y
     2322       
    23102323def makeRefdFFT(Limits,Profile):
    23112324    print 'make fft'
Note: See TracChangeset for help on using the changeset viewer.