Changeset 2356


Ignore:
Timestamp:
Jul 1, 2016 11:52:49 AM (5 years ago)
Author:
vondreele
Message:

revision to RDF calculation/plot

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r2355 r2356  
    291291    XY = xydata['IofQ'][1]   
    292292    #convert to Q
    293     hc = 12.397639
    294     wave = G2mth.getWave(inst)
    295     keV = hc/wave
    296     minQ = npT2q(Tth[0],wave)
    297     maxQ = npT2q(Tth[-1],wave)   
    298     Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
    299     dq = Qpoints[1]-Qpoints[0]
    300     XY[0] = npT2q(XY[0],wave)   
    301 #    Qdata = np.nan_to_num(si.griddata(XY[0],XY[1],Qpoints,method='linear')) #only OK for scipy 0.9!
    302     T = si.interp1d(XY[0],XY[1],bounds_error=False,fill_value=XY[1][0])      #OK for scipy 0.8
    303     Qdata = T(Qpoints)
     293    if 'C' in inst['Type'][0]:
     294        hc = 12.397639
     295        wave = G2mth.getWave(inst)
     296        keV = hc/wave
     297        minQ = npT2q(Tth[0],wave)
     298        maxQ = npT2q(Tth[-1],wave)   
     299        Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
     300        dq = Qpoints[1]-Qpoints[0]
     301        XY[0] = npT2q(XY[0],wave)
     302    elif 'T' in inst['Type'][0]:
     303        difC = inst['difC'][1]
     304        minQ = 2.*np.pi*difC/Tth[-1]
     305        maxQ = 2.*np.pi*difC/Tth[0]
     306        Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
     307        dq = Qpoints[1]-Qpoints[0]
     308        XY[0] = 2.*np.pi*difC/XY[0]
     309    Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])
    304310    Qdata -= np.min(Qdata)*data['BackRatio']
    305311   
    … …  
    317323    FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2)  #these are <f^2>,<f>^2,Cf
    318324    Q = xydata['SofQ'][1][0]
     325    auxPlot.append([Q,np.copy(CF),'CF-unCorr'])
    319326    ruland = Ruland(data['Ruland'],wave,Q,CF)
    320 #    auxPlot.append([Q,ruland,'Ruland'])     
     327    auxPlot.append([Q,ruland,'Ruland'])     
    321328    CF *= ruland
    322 #    auxPlot.append([Q,CF,'CF-Corr'])
     329    auxPlot.append([Q,CF,'CF-Corr'])
    323330    scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
    324331    xydata['SofQ'][1][1] *= scale
    … …  
    327334    scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ])
    328335    xydata['SofQ'][1][1] *= scale
    329 #    if data.get('sinDamp',False):   #not the right thing but leave as place holder
    330 #        sinDamp = np.sin(np.pi*xydata['SofQ'][1][0]/qLimits[1])
    331 #        xydata['SofQ'][1][1] *= sinDamp
    332336    xydata['FofQ'] = copy.deepcopy(xydata['SofQ'])
    333337    xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0)
    … …  
    340344    return auxPlot
    341345   
    342 def MakeRDF(RDFcontrols,background,inst,pwddata,xydata):
    343 #    GSASIIpath.IPyBreak()
     346def MakeRDF(RDFcontrols,background,inst,pwddata):
     347    import scipy.fftpack as ft
    344348    auxPlot = []
    345349    if 'C' in inst['Type'][0]:
    … …  
    349353        keV = hc/wave
    350354        minQ = npT2q(Tth[0],wave)
    351         maxQ = npT2q(Tth[-1],wave)   
     355        maxQ = npT2q(Tth[-1],wave)
     356        powQ = npT2q(Tth,wave) 
    352357       
    353358
    … …  
    357362        minQ = 2.*np.pi*difC/TOF[-1]
    358363        maxQ = 2.*np.pi*difC/TOF[0]
     364        powQ = 2.*np.pi*difC/TOF
    359365       
    360     Qpoints = np.linspace(0.,maxQ,len(pwddata[0]),endpoint=True)
     366#UseObsCalc? = True for obs-calc; False for obs & calc separate plots
     367    piDQ = np.pi/(maxQ-minQ)
     368    Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True)
     369    if RDFcontrols['UseObsCalc']:
     370        Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.)
     371    else:
     372        Qdata = si.griddata(powQ,pwddata[1],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0])
     373    Qdata *= np.sin((Qpoints-minQ)*piDQ)/piDQ
     374#    auxPlot.append([Qpoints,Qdata,'interp1d:'+RDFcontrols['Smooth']])
     375#    GSASIIpath.IPyBreak()
    361376    dq = Qpoints[1]-Qpoints[0]
    362     T = si.interp1d(pwddata[0],pwddata[1],bounds_error=False,fill_value=0.0,kind=RDFcontrols['Smooth'])      #OK for scipy 0.8
    363     Qdata = T(Qpoints)
    364     auxPlot.append([Qpoints,Qdata,'interp1d:'+RDFcontrols['Smooth']])
     377    nR = len(Qdata)
     378    DofR = -dq*np.imag(ft.fft(Qdata,16*nR)[:nR])
     379    R = 0.5*np.pi*np.linspace(0,nR,nR)/(4.*maxQ)
     380    iFin = np.searchsorted(R,RDFcontrols['maxR'])
     381    auxPlot.append([R[:iFin],DofR[:iFin],'D(R)'])   
    365382    return auxPlot
    366     print 'make RDF'
    367383
    368384################################################################################       

  • trunk/GSASIIpwdGUI.py

    r2355 r2356  
    6666            pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
    6767        self.panel = wx.Panel(self)         #just a dummy - gets destroyed in Draw!
    68         self.result = {'UseObsCalc':False,'maxR':10.0,'Smooth':'linear'}
     68        self.result = {'UseObsCalc':True,'maxR':20.0,'Smooth':'linear'}
    6969       
    7070        self.Draw()
    … …  
    9393        mainSizer = wx.BoxSizer(wx.VERTICAL)
    9494        mainSizer.Add(wx.StaticText(self.panel,label='Background RDF controls:'),0,WACV)
    95         useOC = wx.CheckBox(self.panel,label=' Use obs && calc intensities?')
     95        useOC = wx.CheckBox(self.panel,label=' Use obs - calc intensities?')
    9696        useOC.SetValue(self.result['UseObsCalc'])
    9797        useOC.Bind(wx.EVT_CHECKBOX,OnUseOC)
    … …  
    9999        dataSizer = wx.BoxSizer(wx.HORIZONTAL)
    100100        dataSizer.Add(wx.StaticText(self.panel,label=' Smoothing type: '),0,WACV)
    101         smChoice = ['linear','nearest','zero','slinear',]
     101        smChoice = ['linear','nearest',]
    102102        smCombo = wx.ComboBox(self.panel,value=self.result['Smooth'],choices=smChoice,
    103103            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    … …  
    980980        finally:
    981981            dlg.Destroy()
    982         xydata = {}
    983982        PatternId = G2frame.PatternId       
    984983        background = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Background'))
    … …  
    986985        inst,inst2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
    987986        pwddata = G2frame.PatternTree.GetItemPyData(PatternId)[1]
    988         auxPlot = G2pwd.MakeRDF(RDFcontrols,background,inst,pwddata,xydata)
     987        auxPlot = G2pwd.MakeRDF(RDFcontrols,background,inst,pwddata)
    989988#        GSASIIpath.IPyBreak()
    990989        superMinusOne = unichr(0xaf)+unichr(0xb9)
    991990        for plot in auxPlot:
    992991            XY = np.array(plot[:2])
    993             G2plt.PlotXY(G2frame,[XY,],Title=plot[2],labelX=r'$Q,\AA$'+superMinusOne,labelY=r'$I(Q)$')     
    994         G2plt.PlotXY(G2frame,xydata,Title='D(r)') 
    995        
     992            if plot[2] == 'D(R)':
     993                xlabel = r'$R, \AA$'
     994                ylabel = r'$D(R), arb. units$'
     995            else:
     996                xlabel = r'$Q,\AA$'+superMinusOne
     997                ylabel = r'$I(Q)$'
     998            G2plt.PlotXY(G2frame,[XY,],Title=plot[2],labelX=xlabel,labelY=ylabel,lines=True)     
    996999       
    9971000    def BackSizer():
    … …  
    47954798        G2plt.PlotISFG(G2frame,newPlot=False)       
    47964799                       
    4797     def OnSinDamp(event):
    4798         data['sinDamp'] = sinDamp.GetValue()
    4799         auxPlot = ComputePDF(data)
    4800         G2plt.PlotISFG(G2frame,newPlot=False)       
    4801                        
    48024800    def OnPacking(event):
    48034801        try:
    … …  
    49544952        Status.SetStatusText('PDF computed')
    49554953        for plot in auxPlot:
    4956             G2plt.PlotXY(G2frame,plot[:2],Title=plot[2])
     4954            XY = np.array(plot[:2])
     4955            G2plt.PlotXY(G2frame,[XY,],Title=plot[2])
    49574956       
    49584957        G2plt.PlotISFG(G2frame,newPlot=True,type='I(Q)')
    … …  
    51125111    lorch.Bind(wx.EVT_CHECKBOX, OnLorch)
    51135112    sqBox.Add(lorch,0,WACV)
    5114 # this isn't the right thing but something is needed here - leave as place holder
    5115 #    sinDamp = wx.CheckBox(G2frame.dataDisplay,label='SinQ damping?')
    5116 #    sinDamp.SetValue(data['sinDamp'])
    5117 #    sinDamp.Bind(wx.EVT_CHECKBOX,OnSinDamp)
    5118 #    sqBox.Add(sinDamp,0,WACV)
    51195113    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' Scaling q-range: '),0,WACV)
    51205114    SQmin = wx.TextCtrl(G2frame.dataDisplay,value='%.1f'%(data['QScaleLim'][0]))
Note: See TracChangeset for help on using the changeset viewer.