Changeset 2356 for trunk/GSASIIpwd.py


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

revision to RDF calculation/plot

File:
1 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################################################################################       
Note: See TracChangeset for help on using the changeset viewer.