Changeset 2561


Ignore:
Timestamp:
Dec 6, 2016 10:43:26 AM (7 years ago)
Author:
vondreele
Message:

add Rmax to PDF Controls; default = 100.; will not be exact
PDF always generates 5000 points for 0-Rmax; independent of chosen Qmax or image binning.
Trap attempt to calculate PDF without chemical formula - now get error message.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r2550 r2561  
    34553455                            'Geometry':'Cylinder','Diam':1.0,'Pack':0.50,'Form Vol':10.0,
    34563456                            'DetType':'Image plate','ObliqCoeff':0.2,'Ruland':0.025,'QScaleLim':[0,100],
    3457                             'Lorch':True,'BackRatio':0.0}
     3457                            'Lorch':True,'BackRatio':0.0,'Rmax':100.}
    34583458                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='PDF Controls'),Data)
    34593459                        self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='I(Q)'+PWDRname),[])       

  • trunk/GSASIIIO.py

    r2555 r2561  
    968968def PDFSave(G2frame,exports):
    969969    'Save a PDF G(r) and F(Q), S(Q) in column formats'
     970    import scipy.interpolate as scintp
    970971    for export in exports:
    971972        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
    … …  
    979980        fqId = G2gd.GetPatternTreeItemId(G2frame, PickId, FQname)
    980981        grId = G2gd.GetPatternTreeItemId(G2frame, PickId, GRname)
    981         sqdata = np.array(G2frame.PatternTree.GetItemPyData(sqId)[1][:2]).T
    982         fqdata = np.array(G2frame.PatternTree.GetItemPyData(fqId)[1][:2]).T
    983         grdata = np.array(G2frame.PatternTree.GetItemPyData(grId)[1][:2]).T
     982        sqdata = np.array(G2frame.PatternTree.GetItemPyData(sqId)[1][:2])
     983        sqfxn = scintp.interp1d(sqdata[0],sqdata[1],kind='linear')
     984        fqdata = np.array(G2frame.PatternTree.GetItemPyData(fqId)[1][:2])
     985        fqfxn = scintp.interp1d(fqdata[0],fqdata[1],kind='linear')
     986        grdata = np.array(G2frame.PatternTree.GetItemPyData(grId)[1][:2])
     987        grfxn = scintp.interp1d(grdata[0],grdata[1],kind='linear')
    984988        sqfile = open(sqfilename,'w')
    985989        fqfile = open(sqfilename,'w')
    … …  
    989993        grfile.write('#T G(R) %s\n'%(export))
    990994        sqfile.write('#L Q     S(Q)\n')
    991         sqfile.write('#L Q     F(Q)\n')
     995        fqfile.write('#L Q     F(Q)\n')
    992996        grfile.write('#L R     G(R)\n')
    993         for q,sq in sqdata:
     997        qnew = np.arange(sqdata[0][0],sqdata[0][-1],0.005)
     998        rnew = np.arange(grdata[0][0],grdata[0][-1],0.010)
     999        sqnew = zip(qnew,sqfxn(qnew))
     1000        for q,sq in sqnew:
    9941001            sqfile.write("%15.6g %15.6g\n" % (q,sq))
    9951002        sqfile.close()
    9961003        print ' S(Q) saved to: ',sqfilename
    997         for q,fq in fqdata:
     1004        fqnew = zip(qnew,fqfxn(qnew))
     1005        for q,fq in fqnew:
    9981006            fqfile.write("%15.6g %15.6g\n" % (q,fq))
    9991007        fqfile.close()
    10001008        print ' F(Q) saved to: ',fqfilename
    1001         for r,gr in grdata:
     1009        grnew = zip(rnew,grfxn(rnew))
     1010        for r,gr in grnew:
    10021011            grfile.write("%15.6g %15.6g\n" % (r,gr))
    10031012        grfile.close()

  • trunk/GSASIIpwd.py

    r2546 r2561  
    288288    #subtract backgrounds - if any & use PWDR limits
    289289#    GSASIIpath.IPyBreak()
    290     xydata['IofQ'] = copy.deepcopy(xydata['Sample'])
    291     xydata['IofQ'][1] = np.array(xydata['IofQ'][1])[:,Ibeg:Ifin]
     290    IofQ = copy.deepcopy(xydata['Sample'])
     291    IofQ[1] = np.array(IofQ[1])[:,Ibeg:Ifin]
    292292    if data['Sample Bkg.']['Name']:
    293         xydata['IofQ'][1][1] += (xydata['Sample Bkg.'][1][1][Ibeg:Ifin]+
     293        IofQ[1][1] += (xydata['Sample Bkg.'][1][1][Ibeg:Ifin]+
    294294            data['Sample Bkg.']['Add'])*data['Sample Bkg.']['Mult']
    295295    if data['Container']['Name']:
    … …  
    298298            xycontainer += (xydata['Container Bkg.'][1][1][Ibeg:Ifin]+
    299299                data['Container Bkg.']['Add'])*data['Container Bkg.']['Mult']
    300         xydata['IofQ'][1][1] += xycontainer[Ibeg:Ifin]
     300        IofQ[1][1] += xycontainer[Ibeg:Ifin]
    301301    #get element data & absorption coeff.
    302302    ElList = data['ElList']
    303303    Abs = G2lat.CellAbsorption(ElList,data['Form Vol'])
    304304    #Apply angle dependent corrections
    305     Tth = xydata['IofQ'][1][0]
     305    Tth = IofQ[1][0]
    306306    MuR = Abs*data['Diam']/20.0
    307     xydata['IofQ'][1][1] /= Absorb(data['Geometry'],MuR,Tth)
     307    IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth)
    308308    if 'X' in inst['Type'][0]:
    309         xydata['IofQ'][1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
     309        IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0]
    310310    if data['DetType'] == 'Image plate':
    311         xydata['IofQ'][1][1] *= Oblique(data['ObliqCoeff'],Tth)
    312     XY = xydata['IofQ'][1]   
     311        IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth)
     312    XY = IofQ[1]   
    313313    #convert to Q
     314#    nQpoints = len(XY[0])     #points for Q interpolation
     315    nQpoints = 5000
    314316    if 'C' in inst['Type'][0]:
    315317        wave = G2mth.getWave(inst)
    316318        minQ = npT2q(Tth[0],wave)
    317319        maxQ = npT2q(Tth[-1],wave)   
    318         Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
     320        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
    319321        dq = Qpoints[1]-Qpoints[0]
    320322        XY[0] = npT2q(XY[0],wave)
    … …  
    323325        minQ = 2.*np.pi*difC/Tth[-1]
    324326        maxQ = 2.*np.pi*difC/Tth[0]
    325         Qpoints = np.linspace(0.,maxQ,len(XY[0]),endpoint=True)
     327        Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True)
    326328        dq = Qpoints[1]-Qpoints[0]
    327329        XY[0] = 2.*np.pi*difC/XY[0]
    328     Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])
     330    Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0])    #interpolate I(Q)
    329331    Qdata -= np.min(Qdata)*data['BackRatio']
    330332   
    … …  
    333335    maxQ = np.searchsorted(Qpoints,qLimits[1])
    334336    newdata = []
    335     xydata['IofQ'][1][0] = Qpoints
    336     xydata['IofQ'][1][1] = Qdata
     337    xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]]
    337338    for item in xydata['IofQ'][1]:
    338339        newdata.append(item[:maxQ])
    … …  
    359360    xydata['GofR'] = copy.deepcopy(xydata['FofQ'])
    360361    nR = len(xydata['GofR'][1][1])
    361     xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],8*nR)[:nR])
    362     xydata['GofR'][1][0] = 0.25*np.pi*np.linspace(0,nR,nR)/qLimits[1]
     362#    mul = 12
     363    mul = int(round(2.*np.pi*nR/(data['Rmax']*qLimits[1])))
     364    xydata['GofR'][1][0] = 2.*np.pi*np.linspace(0,nR,nR)/(mul*qLimits[1])
     365    xydata['GofR'][1][1] = -dq*np.imag(ft.fft(xydata['FofQ'][1][1],mul*nR)[:nR])
    363366    if data.get('noRing',True):
    364367        xydata['GofR'][1][1] = np.where(xydata['GofR'][1][0]<0.5,0.,xydata['GofR'][1][1])

  • trunk/GSASIIpwdGUI.py

    r2560 r2561  
    47684768    if 'noRing' not in data:
    47694769        data['noRing'] = False
     4770    if 'Rmax' not in data:
     4771        data['Rmax'] = 100.
    47704772   
    47714773    def FillFileSizer(fileSizer,key):
    … …  
    50165018        wx.CallAfter(OnComputePDF,None)
    50175019       
     5020    def OnRmax(event):
     5021        event.Skip()
     5022        try:
     5023            value = float(rmax.GetValue())
     5024            if value > 200. or value < 10.:
     5025                raise ValueError
     5026        except ValueError:
     5027            value = data['Rmax']
     5028        data['Rmax'] = value
     5029        rmax.SetValue('%.1f'%(value))
     5030        wx.CallAfter(OnComputePDF,None)
     5031       
    50185032    def OnResetQ(event):
    50195033        resetQ.SetValue(False)
    … …  
    51575171    def OnComputePDF(event):
    51585172#        print 'Calculating PDF:'
     5173        if not data['ElList']:
     5174            G2frame.ErrorDialog('PDF error','Chemical formula not defined')
     5175            return
    51595176        auxPlot = ComputePDF(data)
    51605177#        print 'Done calculating PDF:'
    … …  
    51815198                if 'PDF' in Name.split()[0]:
    51825199                    Data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,id,'PDF Controls'))
     5200                    if not Data['ElList']:
     5201                        G2frame.ErrorDialog('PDF error','Chemical formula not defined for \n'+Name)
     5202                        return
    51835203                    ComputePDF(Data)                   
    51845204                id, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
    … …  
    53295349   
    53305350    sqBox = wx.BoxSizer(wx.HORIZONTAL)
     5351    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' Scaling q-range: '),0,WACV)
     5352    SQmin = wx.TextCtrl(G2frame.dataDisplay,value='%.1f'%(data['QScaleLim'][0]),size=wx.Size(50,20))
     5353    SQmin.Bind(wx.EVT_KILL_FOCUS,OnSQmin)   
     5354    SQmin.Bind(wx.EVT_TEXT_ENTER,OnSQmin)       
     5355    sqBox.Add(SQmin,0,WACV)
     5356    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' to Qmax '),0,WACV)
     5357    SQmax = wx.TextCtrl(G2frame.dataDisplay,value='%.1f'%(data['QScaleLim'][1]),size=wx.Size(50,20))
     5358    SQmax.Bind(wx.EVT_KILL_FOCUS,OnSQmax)
     5359    SQmax.Bind(wx.EVT_TEXT_ENTER,OnSQmax)       
     5360    sqBox.Add(SQmax,0,WACV)
     5361    resetQ = wx.CheckBox(parent=G2frame.dataDisplay,label='Reset?')
     5362    sqBox.Add(resetQ,0,WACV)
     5363    resetQ.Bind(wx.EVT_CHECKBOX, OnResetQ)
     5364    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' Rmax: '),0,WACV)
     5365    rmax = wx.TextCtrl(G2frame.dataDisplay,value='%.1f'%(data['Rmax']),size=wx.Size(50,20))
     5366    rmax.Bind(wx.EVT_KILL_FOCUS,OnRmax)
     5367    rmax.Bind(wx.EVT_TEXT_ENTER,OnRmax)       
     5368    sqBox.Add(rmax,0,WACV)
    53315369    lorch = wx.CheckBox(parent=G2frame.dataDisplay,label='Lorch damping?')
    53325370    lorch.SetValue(data['Lorch'])
    53335371    lorch.Bind(wx.EVT_CHECKBOX, OnLorch)
    53345372    sqBox.Add(lorch,0,WACV)
    5335     sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' Scaling q-range: '),0,WACV)
    5336     SQmin = wx.TextCtrl(G2frame.dataDisplay,value='%.1f'%(data['QScaleLim'][0]),size=wx.Size(50,20))
    5337     SQmin.Bind(wx.EVT_TEXT_ENTER,OnSQmin)       
    5338     SQmin.Bind(wx.EVT_KILL_FOCUS,OnSQmin)   
    5339     sqBox.Add(SQmin,0)
    5340     sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' to Qmax '),0,WACV)
    5341     SQmax = wx.TextCtrl(G2frame.dataDisplay,value='%.1f'%(data['QScaleLim'][1]),size=wx.Size(50,20))
    5342     SQmax.Bind(wx.EVT_TEXT_ENTER,OnSQmax)       
    5343     SQmax.Bind(wx.EVT_KILL_FOCUS,OnSQmax)
    5344     sqBox.Add(SQmax,0)
    5345     resetQ = wx.CheckBox(parent=G2frame.dataDisplay,label='Reset?')
    5346     sqBox.Add(resetQ,0)
    5347     resetQ.Bind(wx.EVT_CHECKBOX, OnResetQ)
    53485373    noRing = wx.CheckBox(parent=G2frame.dataDisplay,label='Suppress G(0) ringing?')
    53495374    noRing.SetValue(data['noRing'])
Note: See TracChangeset for help on using the changeset viewer.