Changeset 2783 for trunk


Ignore:
Timestamp:
Apr 15, 2017 11:31:45 AM (8 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

Location:
trunk
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified trunk/GSASII.py

    r2778 r2783  
    20602060            Tmin = min(rd.reflectometrydata[0])
    20612061            Tmax = max(rd.reflectometrydata[0])
     2062            ifDQ = np.any(rd.reflectometrydata[5])
    20622063            valuesdict = {
    20632064                'wtFactor':1.0,
    … …  
    20652066                'ranId':ran.randint(0,sys.maxint),
    20662067                'Offset':[0.0,0.0],
     2068                'ifDQ':ifDQ
    20672069                }
    20682070            rd.Sample['ranId'] = valuesdict['ranId'] # this should be removed someday

  • TabularUnified 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'

  • TabularUnified trunk/GSASIIpwdGUI.py

    r2782 r2783  
    224224    return {'Layers':[{'Name':'vacuum','DenMul':[1.0,False],},                                  #top layer
    225225        {'Name':'vacuum','Rough':[0.,False],'Penetration':[0.,False],'DenMul':[1.0,False],}],   #bottom layer
    226         'Scale':[1.0,False],'FltBack':[0.0,False],'Zero':'Top',                                 #globals
    227         'Minimizer':'LMLS','Resolution':[0.,'Const dq/q'],'Recomb':0.5,'Toler':0.5,           #minimizer controls
     226        'Scale':[1.0,False],'FltBack':[0.0,False],'Zero':'Top','dQ type':'None',                #globals
     227        'Minimizer':'LMLS','Resolution':[0.,'Const dq/q'],'Recomb':0.5,'Toler':0.5,             #minimizer controls
    228228        'DualFitFiles':['',],'DualFltBacks':[[0.0,False],],'DualScales':[[1.0,False],]}         #optional stuff for multidat fits?
    229229       
    … …  
    48964896            XY = [[R[:2500],F[:2500]],]
    48974897            G2plt.PlotXY(G2frame,XY,labelX='thickness',labelY='F(R)',newPlot=True,
    4898                 Title='Fourier transform',lines=True) 
     4898                Title='Fourier transform',lines=True)
     4899           
     4900        def OndQSel(event):
     4901            data['dQ type'] = dQSel.GetStringSelection()
    48994902           
    49004903        controlSizer = wx.BoxSizer(wx.VERTICAL)
    49014904        resol = wx.BoxSizer(wx.HORIZONTAL)
    4902         resol.Add(wx.StaticText(G2frame.dataDisplay,label=' Instrument resolution (%'+GkDelta+'Q/Q): '),0,WACV)
    4903         resol.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data['Resolution'],0,nDig=(10,2),min=0.,max=0.2),0,WACV)
    4904         resol.Add(wx.StaticText(G2frame.dataDisplay,label=' Zero position location: '),0,WACV)
    4905         poslist = ['Top','Bottom']
    4906         refpos = wx.ComboBox(G2frame.dataDisplay,value=data['Zero'],choices=poslist,
    4907             style=wx.CB_READONLY|wx.CB_DROPDOWN)
    4908         refpos.Bind(wx.EVT_COMBOBOX, OnRefPos)
    4909         resol.Add(refpos,0,WACV)
     4905        choice = ['None','const '+GkDelta+'Q/Q',]
     4906        if ProfDict['ifDQ']:
     4907            choice += [GkDelta+'Q/Q in data']
     4908        dQSel = wx.RadioBox(G2frame.dataDisplay,wx.ID_ANY,'Instrument resolution type:',choices=choice,
     4909            majorDimension=0,style=wx.RA_SPECIFY_COLS)
     4910        dQSel.SetStringSelection(data['dQ type'])
     4911        dQSel.Bind(wx.EVT_RADIOBOX,OndQSel)
     4912        resol.Add(dQSel,0,WACV)
     4913        resol.Add(wx.StaticText(G2frame.dataDisplay,label=' (FWHM %): '),0,WACV)
     4914        resol.Add(G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data['Resolution'],0,nDig=(10,3),min=0.,max=5.),0,WACV)
    49104915        controlSizer.Add(resol,0,WACV)
    49114916        minimiz = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    49284933        sld.Bind(wx.EVT_CHECKBOX, OnSLDplot)
    49294934        plotSizer.Add(sld,0,WACV)
     4935        plotSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Zero position location: '),0,WACV)
     4936        poslist = ['Top','Bottom']
     4937        refpos = wx.ComboBox(G2frame.dataDisplay,value=data['Zero'],choices=poslist,
     4938            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     4939        refpos.Bind(wx.EVT_COMBOBOX, OnRefPos)
     4940        plotSizer.Add(refpos,0,WACV)
    49304941#        q4fft = wx.CheckBox(G2frame.dataDisplay,label='Plot fft?')
    49314942#        q4fft.Bind(wx.EVT_CHECKBOX, OnQ4fftplot)
    … …  
    51805191    Substances = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Substances'))['Substances']
    51815192    ProfDict,Profile,Name = G2frame.PatternTree.GetItemPyData(G2frame.PatternId)[:3]
     5193    if 'ifDQ' not in ProfDict:
     5194        ProfDict['ifDQ'] = np.any(Profile[5])
     5195        data['dQ type'] = 'None'
    51825196    Limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Limits'))
    51835197    Inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[0]

  • TabularUnified trunk/imports/G2pwd_xye.py

    r2780 r2783  
    2626            extensionlist=('.xye','.chi',),
    2727            strictExtension=False,
    28             formatName = 'Topas xye & Fit2D chi',
    29             longFormatName = 'Topas .xye & Fit2D .chi powder data file'
     28            formatName = 'Topas xye or 2th Fit2D chi',
     29            longFormatName = 'Topas .xye or 2th Fit2D .chi powder data file'
    3030            )
    3131
    … …  
    3939        if '.chi' in filepointer.name:
    4040            self.Chi = True
     41        if2theta = False
    4142        for i,S in enumerate(filepointer):
    4243            if not S:
    … …  
    4647                if self.Chi:
    4748                    if i < 4:
     49                        if  '2-theta' in S.lower():
     50                            if2theta = True
    4851                        continue
    4952                    else:
    5053                        begin = False
    5154                else:
     55                    if2theta = True
    5256                    if gotCcomment and S.find('*/') > -1:
    5357                        begin = False
    … …  
    6266                # valid line to read?
    6367            #vals = S.split()
     68            if not if2theta:
     69                self.errors = 'Not a 2-theta chi file'
     70                return False
    6471            vals = S.replace(',',' ').replace(';',' ').split()
    6572            if len(vals) == 2 or len(vals) == 3:
Note: See TracChangeset for help on using the changeset viewer.