source: trunk/GSASIIfpaGUI.py @ 4959

Last change on this file since 4959 was 4959, checked in by toby, 2 years ago

add NIST spectrum kludge for incident beam mono
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 35.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-06-15 18:27:03 +0000 (Tue, 15 Jun 2021) $
4# $Author: toby $
5# $Revision: 4959 $
6# $URL: trunk/GSASIIfpaGUI.py $
7# $Id: GSASIIfpaGUI.py 4959 2021-06-15 18:27:03Z toby $
8########### SVN repository information ###################
9'''
10*GSASIIfpaGUI: Fundamental Parameters Routines*
11===============================================
12
13This module contains routines for getting Fundamental Parameters
14Approach (FPA) input, setting up for running the NIST XRD Fundamental
15Parameters Code, plotting the convolutors and computing a set of peaks
16generated by that code.
17
18'''
19from __future__ import division, print_function
20import os.path
21import numpy as np
22import copy
23
24import wx
25import wx.lib.scrolledpanel as wxscroll
26
27import NIST_profile as FP
28
29import GSASIIpath
30import GSASIIctrlGUI as G2G
31import GSASIIdataGUI as G2gd
32import GSASIIplot as G2plt
33import GSASIImath as G2mth
34import GSASIIpwd as G2pwd
35WACV = wx.ALIGN_CENTER_VERTICAL
36
37simParms = {}
38'''Parameters to set range for pattern simulation
39'''
40
41parmDict = {}
42'''Parameter dict used for reading Topas-style values. These are
43converted to SI units and placed into :data:`NISTparms`
44'''
45
46NISTparms = {}
47'''Parameters in a nested dict, with an entry for each concolutor. Entries in
48those dicts have values in SI units (of course). NISTparms can be
49can be input directly or can be from created from :data:`parmDict`
50by :func:`XferFPAsettings`
51'''
52
53BraggBrentanoParms = [
54    ('divergence', 0.5, 'Bragg-Brentano divergence angle (degrees)'),
55    ('soller_angle', 2.0, 'Soller slit axial divergence (degrees)'),
56    ('Rs', 220, 'Diffractometer radius (mm)'),
57    ('filament_length', 12., 'X-ray tube line focus length (mm)'),
58    ('sample_length', 12., 'Illuminated sample length in axial direction (mm)'),
59    ('receiving_slit_length', 12., 'Length of receiving slit in axial direction (mm)'),
60    ('LAC_cm', 0.,'Linear absorption coef. adjusted for packing density (cm-1)'),
61    ('sample_thickness', 1., 'Depth of sample (mm)'),
62    ('convolution_steps', 8, 'Number of Fourier-space bins per two-theta step'),
63    ('source_width', 0.04,'Tube filament width, in projection at takeoff angle (mm)'),
64    ('tube-tails_L-tail', -1.,'Left-side tube tails width, in projection (mm)'),   
65    ('tube-tails_R-tail', 1.,'Right-side tube tails width, in projection (mm)'),
66    ('tube-tails_rel-I', 0.001,'Tube tails fractional intensity (no units)'),
67    ]
68'''FPA dict entries used in :func:`FillParmSizer`. Tuple contains
69a dict key, a default value and a description. These are the parameters
70needed for all Bragg Brentano instruments
71'''
72
73BBPointDetector = [
74    ('receiving_slit_width', 0.2, 'Width of receiving slit (mm)'),]
75'''Additional FPA dict entries used in :func:`FillParmSizer`
76needed for Bragg Brentano instruments with point detectors.
77'''
78
79BBPSDDetector = [
80    ('lpsd_th2_angular_range', 3.0, 'Angular range observed by PSD (degrees 2Theta)'),
81    ('lpsd_equitorial_divergence', 0.1, 'Equatorial divergence of the primary beam (degrees)'),]
82'''Additional FPA dict entries used in :func:`FillParmSizer`
83needed for Bragg Brentano instruments with linear (1-D) PSD detectors.
84'''
85
86IBmonoParms = [
87    ('src_mono_mm',119,'Distance from xray line source to monochromator crystal (mm)'),
88    ('focus_mono_mm',217,'Distance from monochromator crystal to focus slit (mm)'),
89    ('passband_mistune',-0.145,'Offset from the tuning of the IBM to the center of the reference line of the spectrum, in units of its bandwidth'),
90    ('mono_src_proj_mn',51,'Projection width of line-focus xray source on the monochromator along beam direction (microns), sets bandwidth'),
91    ('passband_shoulder',.087,'Width of transition region from flat top to tails, in units of the bandwidth'),
92    ('two_theta_mono',27.27,'The full diffraction angle of the IBM crystal, e.g. double 2theta-Bragg for the mono (deg)'),
93    ('mono_slit_attenuation',0.03,'Attenuation of Cu K alpha2 lines due to focal slit'),
94    ]
95'''Additional FPA dict entries used in :func:`FillParmSizer`, needed for Incident Beam Monochromator
96'''
97   
98Citation = '''MH Mendenhall, K Mullen && JP Cline (2015), J. Res. of NIST, 120, p223. DOI: 10.6028/jres.120.014
99
100For Incident Beam Mono model, also cite: MH Mendenhall, D. Black && JP Cline, J. Appl. Cryst. (2019), 52, p1087. DOI: 10.1107/S1600576719010951
101'''
102
103IBmono = False
104'''set to True if an incident beam monochromator is in use
105'''
106DetMode = 'BBpoint'
107'''The type of detector, either 'BBpoint' for Bragg-Brentano point detector or
108      or 'BBPSD' (linear) position sensitive detector
109'''
110
111def SetCu2Wave():
112    '''Set the parameters to the two-line Cu K alpha 1+2 spectrum
113    '''
114    parmDict['wave'] = {i:v for i,v in enumerate((1.5405925, 1.5443873))}
115    parmDict['int'] = {i:v for i,v in enumerate((0.653817, 0.346183))}
116    parmDict['lwidth'] = {i:v for i,v in enumerate((0.501844,0.626579))}
117
118def SetCu6wave():
119    '''Set the emission parameters to the NIST six-line Cu K alpha spectrum
120    '''
121    # values from Marcus Mendenhall from atan_windowed_FP_profile.py
122    parmDict['wave'] = {i:v for i,v in enumerate((1.5405925, 1.5443873, 1.5446782, 1.5410769, 1.53471, 1.53382, ))}
123    parmDict['int'] = {i:v for i,v in enumerate((0.58384351, 0.2284605 , 0.11258773, 0.07077796, 0.0043303, 0.00208613, ))}
124    parmDict['lwidth'] = {i:v for i,v in enumerate((0.436, 0.487, 0.63, 0.558, 2.93, 2.93,))}
125
126def SetMonoWave():
127    '''Eliminates the short-wavelength line from the six-line Cu K
128    alpha spectrum when incident beam mono; resets it to 6 if no mono
129    '''
130    if IBmono and len(parmDict['wave']) == 6:
131        for key in 'wave','int','lwidth':
132            if 5 in parmDict[key]: del parmDict[key][5]
133            if 4 in parmDict[key]: del parmDict[key][4]
134    if (not IBmono) and len(parmDict['wave']) == 4:
135        SetCu6wave()
136
137def writeNIST(filename):
138    '''Write the NIST FPA terms into a JSON-like file that can be reloaded
139    in _onReadFPA
140    ''' 
141    if not filename: return
142   
143    fp = open(filename,'w')
144    fp.write('# parameters to be used in the NIST XRD Fundamental Parameters program\n')
145    fp.write('{\n')
146    for key in sorted(NISTparms):
147        fp.write("   '"+key+"' : "+str(NISTparms[key])+",")
148        if not key: fp.write('  # global parameters')
149        fp.write('\n')
150    fp.write('}\n')
151    fp.close()
152       
153#SetCu2Wave() # use these as default
154SetCu6wave() # use these as default
155SetMonoWave()
156
157def FillParmSizer():
158    '''Create a list of input items for the parameter section of the
159    input window, sets default values when not set and displays them
160    in the scrolledpanel prmPnl.
161    '''
162    prmSizer = prmPnl.GetSizer()
163    prmSizer.Clear(True)
164    if IBmono:
165        itemList = [i for i in BraggBrentanoParms if not i[0].startswith('tube-tails')]
166    else:
167        itemList = copy.deepcopy(BraggBrentanoParms)
168    if DetMode == 'BBpoint':
169        itemList += BBPointDetector
170    elif DetMode == 'BBPSD':
171        itemList += BBPSDDetector
172    else:
173        raise Exception('Unknown DetMode in FillParmSizer: '+DetMode)
174    if IBmono:
175        itemList += IBmonoParms   
176    text = wx.StaticText(prmPnl,wx.ID_ANY,'label',style=wx.ALIGN_CENTER,
177                size=(170,-1)) # make just a bit bigger than largest item in column
178    text.SetBackgroundColour(wx.WHITE)
179    prmSizer.Add(text,0,wx.EXPAND)
180    text = wx.StaticText(prmPnl,wx.ID_ANY,'value',style=wx.ALIGN_CENTER)
181    text.SetBackgroundColour(wx.WHITE)
182    prmSizer.Add(text,0,wx.EXPAND)
183    text = wx.StaticText(prmPnl,wx.ID_ANY,'explanation',style=wx.ALIGN_CENTER)
184    text.SetBackgroundColour(wx.WHITE)
185    prmSizer.Add(text,0,wx.EXPAND)
186    for lbl,defVal,text in itemList:
187        prmSizer.Add(wx.StaticText(prmPnl,wx.ID_ANY,lbl),1,wx.ALIGN_RIGHT|WACV,1)
188        if lbl not in parmDict: parmDict[lbl] = defVal
189        ctrl = G2G.ValidatedTxtCtrl(prmPnl,parmDict,lbl,size=(70,-1))
190        prmSizer.Add(ctrl,1,wx.ALL|WACV,1)
191        txt = wx.StaticText(prmPnl,wx.ID_ANY,text,size=(400,-1))
192        txt.Wrap(380)
193        prmSizer.Add(txt)
194    px1,py = prmSizer.GetSize()
195    prmSizer.Layout()
196    FPdlg = prmPnl.GetParent()
197    FPdlg.SendSizeEvent()
198       
199def MakeTopasFPASizer(G2frame,FPdlg,SetButtonStatus):
200    '''Create a GUI with parameters for the NIST XRD Fundamental Parameters Code.
201    Parameter input is modeled after Topas input parameters.
202
203    :param wx.Frame G2frame: main GSAS-II window
204    :param wx.Window FPdlg: Frame or Dialog where GUI will appear
205    :param SetButtonStatus: a callback function to call to see what buttons in
206      this windows can be enabled. Called with done=True to trigger closing
207      the parent window as well.
208    :returns: a sizer with the GUI controls
209    '''
210    def _onOK(event):
211        XferFPAsettings(parmDict)
212        SetButtonStatus(done=True) # done=True triggers the simulation
213        FPdlg.Destroy()
214    def _onClose(event):
215        SetButtonStatus()
216        FPdlg.Destroy()
217    def _onAddWave(event):
218        newkey = max(parmDict['wave'].keys())+1
219        for key,defVal in zip(
220            ('wave','int','lwidth'),
221            (0.0,   1.0,   0.1),
222            ):
223            parmDict[key][newkey] = defVal
224        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
225    def _onRemWave(event):
226        lastkey = max(parmDict['wave'].keys())
227        for key in ('wave','int','lwidth'):
228            if lastkey in parmDict[key]:
229                del parmDict[key][lastkey]
230        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
231    def _onSetCu6wave(event):
232        SetCu6wave()
233        SetMonoWave()
234        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
235    def _onSetCu2Wave(event):
236        SetCu2Wave()
237        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
238    def _onSetDetBtn(event):
239        global DetMode
240        if detBtn1.GetValue():
241            DetMode = 'BBpoint'
242            wx.CallAfter(FillParmSizer)
243        else:
244            DetMode = 'BBPSD'
245            wx.CallAfter(FillParmSizer)
246    def _onSetMonoBtn(event):
247        global IBmono
248        IBmono = not monoBtn1.GetValue()
249        SetMonoWave()
250        #wx.CallAfter(FillParmSizer)
251        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
252    def PlotTopasFPA(event):
253        XferFPAsettings(parmDict)
254        ttArr = np.arange(max(0.5,
255                            simParms['plotpos']-simParms['calcwid']),
256                            simParms['plotpos']+simParms['calcwid'],
257                            simParms['step'])
258        intArr = np.zeros_like(ttArr)
259        NISTpk = setupFPAcalc()
260        try:
261            center_bin_idx,peakObj = doFPAcalc(
262                NISTpk,ttArr,simParms['plotpos'],simParms['calcwid'],
263                simParms['step'])
264        except Exception as err:
265            msg = "Error computing convolution, revise input"
266            print(msg)
267            print(err)
268            return
269        G2plt.PlotFPAconvolutors(G2frame,NISTpk)
270        pkPts = len(peakObj.peak)
271        pkMax = peakObj.peak.max()
272        startInd = center_bin_idx-(pkPts//2) #this should be the aligned start of the new data
273        # scale peak so max I=10,000 and add into intensity array
274        if startInd < 0:
275            intArr[:startInd+pkPts] += 10000 * peakObj.peak[-startInd:]/pkMax
276        elif startInd > len(intArr):
277            return
278        elif startInd+pkPts >= len(intArr):
279            offset = pkPts - len( intArr[startInd:] )
280            intArr[startInd:startInd+pkPts-offset] += 10000 * peakObj.peak[:-offset]/pkMax
281        else:
282            intArr[startInd:startInd+pkPts] += 10000 * peakObj.peak/pkMax
283        G2plt.PlotXY(G2frame, [(ttArr, intArr)],
284                     labelX=r'$2\theta, deg$',
285                     labelY=r'Intensity (arbitrary)',
286                     Title='FPA peak', newPlot=True, lines=True)
287    def _onSaveFPA(event):
288        XferFPAsettings(parmDict)
289        filename = G2G.askSaveFile(G2frame,'','.NISTfpa',
290                                       'dict of NIST FPA values',FPdlg)
291        writeNIST(filename)
292
293    if FPdlg.GetSizer(): FPdlg.GetSizer().Clear(True)
294    MainSizer = wx.BoxSizer(wx.VERTICAL)
295    MainSizer.Add((-1,5))
296    waveSizer = wx.FlexGridSizer(cols=len(parmDict['wave'])+1,hgap=3,vgap=5)
297    for lbl,prm,defVal in zip(
298            (u'Wavelength (\u212b)','Rel. Intensity',u'Lorentz Width\n(\u212b/1000)'),
299            ('wave','int','lwidth'),
300            (0.0,   1.0,   0.1),
301            ):
302        text = wx.StaticText(FPdlg,wx.ID_ANY,lbl,style=wx.ALIGN_CENTER)
303        text.SetBackgroundColour(wx.WHITE)
304        waveSizer.Add(text,0,wx.EXPAND)
305        if prm not in parmDict: parmDict[prm] = {}
306        for i in parmDict['wave'].keys():
307            if i not in parmDict[prm]: parmDict[prm][i] = defVal
308            if prm == 'wave':
309                ctrl = G2G.ValidatedTxtCtrl(FPdlg,parmDict[prm],i,size=(90,-1),nDig=(10,6))
310            else:
311                ctrl = G2G.ValidatedTxtCtrl(FPdlg,parmDict[prm],i,size=(90,-1),nDig=(10,3))
312            waveSizer.Add(ctrl,1,WACV,1)
313    MainSizer.Add(waveSizer)
314    MainSizer.Add((-1,5))
315    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
316    btn = wx.Button(FPdlg, wx.ID_ANY,'Add wave')
317    btnsizer.Add(btn)
318    btn.Bind(wx.EVT_BUTTON,_onAddWave)
319    btn = wx.Button(FPdlg, wx.ID_ANY,'Remove wave')
320    btnsizer.Add(btn)
321    btn.Bind(wx.EVT_BUTTON,_onRemWave)
322    btn = wx.Button(FPdlg, wx.ID_ANY,'CuKa1+2')
323    btnsizer.Add(btn)
324    btn.Bind(wx.EVT_BUTTON,_onSetCu2Wave)
325    btn = wx.Button(FPdlg, wx.ID_ANY,'NIST CuKa')
326    btnsizer.Add(btn)
327    btn.Bind(wx.EVT_BUTTON,_onSetCu6wave)
328    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
329    MainSizer.Add((-1,5))
330   
331    btnsizer = wx.GridBagSizer( 2, 5)
332    btnsizer.Add( wx.StaticText(FPdlg, wx.ID_ANY, 'Detector type'),
333                 (0,0), (2,1), wx.ALIGN_CENTER | wx.ALL, 5)   
334    detBtn1 = wx.RadioButton(FPdlg,wx.ID_ANY,'Point',style=wx.RB_GROUP)
335    detBtn1.SetValue(DetMode == 'BBpoint')
336    btnsizer.Add(detBtn1, (0,1))
337    detBtn1.Bind(wx.EVT_RADIOBUTTON,_onSetDetBtn)
338    detBtn2 = wx.RadioButton(FPdlg,wx.ID_ANY,'PSD')
339    detBtn2.SetValue(not DetMode == 'BBpoint')
340    btnsizer.Add(detBtn2, (1,1))
341    detBtn2.Bind(wx.EVT_RADIOBUTTON,_onSetDetBtn)
342    btnsizer.Add( (40,-1), (0,2), (1,1), wx.ALIGN_CENTER | wx.ALL, 5)   
343    btnsizer.Add( wx.StaticText(FPdlg, wx.ID_ANY, 'Incident Beam Mono'),
344                 (0,3), (2,1), wx.ALIGN_CENTER | wx.ALL, 5)   
345    monoBtn1 = wx.RadioButton(FPdlg,wx.ID_ANY,'No',style=wx.RB_GROUP)
346    monoBtn1.SetValue(not IBmono)
347    btnsizer.Add(monoBtn1, (0,4))
348    monoBtn1.Bind(wx.EVT_RADIOBUTTON,_onSetMonoBtn)
349    monoBtn2 = wx.RadioButton(FPdlg,wx.ID_ANY,'Yes')
350    monoBtn2.SetValue(IBmono)
351    btnsizer.Add(monoBtn2, (1,4))
352    monoBtn2.Bind(wx.EVT_RADIOBUTTON,_onSetMonoBtn)
353    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
354
355    global prmPnl
356    prmPnl = wxscroll.ScrolledPanel(FPdlg, wx.ID_ANY, #size=(200,200),
357        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
358    prmSizer = wx.FlexGridSizer(cols=3,hgap=3,vgap=5)
359    prmPnl.SetSizer(prmSizer)
360    FillParmSizer()   
361    MainSizer.Add(prmPnl,1,wx.EXPAND,1)
362    prmPnl.SetAutoLayout(1)
363    prmPnl.SetupScrolling()
364   
365    MainSizer.Add((-1,4))
366    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
367    btn = wx.Button(FPdlg, wx.ID_ANY, 'Plot peak')
368    btnsizer.Add(btn)
369    btn.Bind(wx.EVT_BUTTON,PlotTopasFPA)
370    btnsizer.Add(wx.StaticText(FPdlg,wx.ID_ANY,' at '))
371    if 'plotpos' not in simParms: simParms['plotpos'] =  simParms['minTT']
372    ctrl = G2G.ValidatedTxtCtrl(FPdlg,simParms,'plotpos',size=(70,-1))
373    btnsizer.Add(ctrl)
374    btnsizer.Add(wx.StaticText(FPdlg,wx.ID_ANY,' deg.  '))   
375    saveBtn = wx.Button(FPdlg, wx.ID_ANY,'Save FPA dict')
376    btnsizer.Add(saveBtn)
377    saveBtn.Bind(wx.EVT_BUTTON,_onSaveFPA)
378    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
379    MainSizer.Add((-1,4))
380    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
381    OKbtn = wx.Button(FPdlg, wx.ID_OK)
382    OKbtn.SetDefault()
383    btnsizer.Add(OKbtn)
384    Cbtn = wx.Button(FPdlg, wx.ID_CLOSE,"Cancel") 
385    btnsizer.Add(Cbtn)
386    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
387    MainSizer.Add((-1,4))
388    # bindings for close of window
389    OKbtn.Bind(wx.EVT_BUTTON,_onOK)
390    Cbtn.Bind(wx.EVT_BUTTON,_onClose)
391    FPdlg.SetSizer(MainSizer)
392    MainSizer.Layout()
393    MainSizer.Fit(FPdlg)
394    # control window size
395    px,py = prmSizer.GetSize()
396    dx,dy = FPdlg.GetSize()
397    FPdlg.SetMinSize((-1,-1))
398    FPdlg.SetMaxSize((-1,-1)) 
399    FPdlg.SetMinSize((dx,dy+200)) # leave a min of 200 points for scroll panel
400    FPdlg.SetMaxSize((max(dx,700),850)) 
401    FPdlg.SetSize((max(dx,px+20),min(750,dy+py+30))) # 20 for scroll bar, 30 for a bit of room at bottom
402   
403def XferFPAsettings(InpParms):
404    '''convert Topas-type parameters to SI units for NIST and place in a dict sorted
405    according to use in each convoluter
406
407    :param dict InpParms: a dict with Topas-like parameters, as set in
408      :func:`MakeTopasFPASizer`
409    :returns: a nested dict with global parameters and those for each convolution
410    '''
411    # cleanup old stuff
412    for key in "tube_tails","absorption","si_psd","displacement","receiver_slit":
413        if key in NISTparms:
414            del NISTparms[key]
415   
416    keys = list(InpParms['wave'].keys())
417    source_wavelengths_m = 1.e-10 * np.array([InpParms['wave'][i] for i in keys])
418    la = [InpParms['int'][i] for i in keys]
419 
420    if IBmono:  # kludge: apply mono_slit_attenuation since it is not part of NIST FPA code
421        norm = [InpParms['mono_slit_attenuation'] if
422                    (1.5443 < InpParms['wave'][i] < 1.5447) else 1.
423                for i in keys]
424        source_intensities = norm * np.array(la)/max(la)
425    else:
426        source_intensities = np.array(la)/max(la)
427    source_lor_widths_m = 1.e-10 * 1.e-3 * np.array([InpParms['lwidth'][i] for i in keys])
428    source_gauss_widths_m = 1.e-10 * 1.e-3 * np.array([0.001 for i in keys])
429   
430    NISTparms["emission"] = {'emiss_wavelengths' : source_wavelengths_m,
431                'emiss_intensities' : source_intensities,
432                'emiss_gauss_widths' : source_gauss_widths_m,
433                'emiss_lor_widths' : source_lor_widths_m,
434                'crystallite_size_gauss' : 1.e-9 * InpParms.get('Size_G',1e6),
435                'crystallite_size_lor' : 1.e-9 * InpParms.get('Size_L',1e6)}
436    if IBmono:
437        NISTparms["emission"]['a_mono'] = InpParms['src_mono_mm'] * 10**-3
438        NISTparms["emission"]['b_mono'] = InpParms['focus_mono_mm'] * 10**-3
439        NISTparms["emission"]['ibm_source_width'] = InpParms['mono_src_proj_mn'] * 10**-6
440        for i in  ('passband_mistune','passband_shoulder','two_theta_mono'):
441            NISTparms["emission"][i] = InpParms[i]
442    elif InpParms.get('source_width', 0) > 0 and InpParms.get(
443            'tube-tails_rel-I',0) > 0:
444        NISTparms["tube_tails"] = {
445            'main_width' : 1e-3 * InpParms.get('source_width', 0.),
446            'tail_left' : -1e-3 * InpParms.get('tube-tails_L-tail',0.),
447            'tail_right' : 1e-3 * InpParms.get('tube-tails_R-tail',0.),
448            'tail_intens' : InpParms.get('tube-tails_rel-I',0.),}
449   
450    if InpParms['filament_length'] == InpParms['receiving_slit_length']: # workaround:
451        InpParms['receiving_slit_length'] *= 1.00001 # avoid bug when slit lengths are identical
452    NISTparms["axial"]  = {
453            'axDiv':"full", 'slit_length_source' : 1e-3*InpParms['filament_length'],
454            'slit_length_target' : 1e-3*InpParms['receiving_slit_length'],
455            'length_sample' : 1e-3 * InpParms['sample_length'], 
456            'n_integral_points' : 10,
457            'angI_deg' : InpParms['soller_angle'],
458            'angD_deg': InpParms['soller_angle']
459            }
460    if InpParms.get('LAC_cm',0) > 0:
461        NISTparms["absorption"] = {
462            'absorption_coefficient': InpParms['LAC_cm']*100, #like LaB6, in m^(-1)
463            'sample_thickness': 1e-3 * InpParms['sample_thickness'],
464            }
465
466    if InpParms.get('lpsd_equitorial_divergence',0) > 0 and InpParms.get(
467            'lpsd_th2_angular_range',0) > 0 and DetMode == 'BBPSD':
468        PSDdetector_length_mm=np.arcsin(np.pi*InpParms['lpsd_th2_angular_range']/180.
469                                            )*InpParms['Rs'] # mm
470        NISTparms["si_psd"] = {
471            'equatorial_divergence_deg': InpParms['lpsd_equitorial_divergence'],
472            'si_psd_window_bounds': (0.,PSDdetector_length_mm/1000.)
473            }
474
475    if InpParms.get('Specimen_Displacement'):
476        NISTparms["displacement"] = {'specimen_displacement': 1e-3 * InpParms['Specimen_Displacement']}
477
478    if InpParms.get('receiving_slit_width'):
479        NISTparms["receiver_slit"] = {'slit_width':1e-3*InpParms['receiving_slit_width']}
480
481    # set Global parameters
482    max_wavelength = source_wavelengths_m[np.argmax(source_intensities)]
483    NISTparms[""] = {
484        'equatorial_divergence_deg' : InpParms['divergence'],
485        'dominant_wavelength' : max_wavelength,
486        'diffractometer_radius' : 1e-3* InpParms['Rs'],
487        'oversampling' : InpParms['convolution_steps'],
488        }
489
490def setupFPAcalc():
491    '''Create a peak profile object using the NIST XRD Fundamental
492    Parameters Code.
493   
494    :returns: a profile object that can provide information on
495      each convolution or compute the composite peak shape.
496    '''
497    if IBmono:
498        p=FP.FP_windowed(anglemode="twotheta",
499                    output_gaussian_smoother_bins_sigma=1.0,
500                    oversampling=NISTparms.get('oversampling',10))
501    else:
502        p=FP.FP_profile(anglemode="twotheta",
503                    output_gaussian_smoother_bins_sigma=1.0,
504                    oversampling=NISTparms.get('oversampling',10))
505       
506    p.debug_cache=False
507    #set parameters for each convolver
508    for key in NISTparms:
509        if key:
510            p.set_parameters(convolver=key,**NISTparms[key])
511        else:
512            p.set_parameters(**NISTparms[key])
513    return p
514       
515def doFPAcalc(NISTpk,ttArr,twotheta,calcwid,step):
516    '''Compute a single peak using a NIST profile object
517
518    :param object NISTpk: a peak profile computational object from the
519      NIST XRD Fundamental Parameters Code, typically established from
520      a call to :func:`SetupFPAcalc`
521    :param np.Array ttArr: an evenly-spaced grid of two-theta points (degrees)
522    :param float twotheta: nominal center of peak (degrees)
523    :param float calcwid: width to perform convolution (degrees)
524    :param float step: step size
525    '''
526    # find closest point to twotheta (may be outside limits of the array)
527    center_bin_idx=min(ttArr.searchsorted(twotheta),len(ttArr)-1)
528    NISTpk.set_optimized_window(twotheta_exact_bin_spacing_deg=step,
529                twotheta_window_center_deg=ttArr[center_bin_idx],
530                twotheta_approx_window_fullwidth_deg=calcwid,
531                )
532    NISTpk.set_parameters(twotheta0_deg=twotheta)
533    return center_bin_idx,NISTpk.compute_line_profile()
534
535def MakeSimSizer(G2frame, dlg):
536    '''Create a GUI to get simulation with parameters for Fundamental
537    Parameters fitting.
538
539    :param wx.Window dlg: Frame or Dialog where GUI will appear
540
541    :returns: a sizer with the GUI controls
542 
543    '''
544    def _onOK(event):
545        msg = ''
546        if simParms['minTT']-simParms['calcwid']/1.5 < 0.1:
547            msg += 'First peak minus half the calc width is too low'
548        if simParms['maxTT']+simParms['calcwid']/1.5 > 175:
549            if msg: msg += '\n'
550            msg += 'Last peak plus half the calc width is too high'
551        if simParms['npeaks'] < 8:
552            if msg: msg += '\n'
553            msg += 'At least 8 peaks are needed'
554        if msg:
555            G2G.G2MessageBox(dlg,msg,'Bad input, try again')
556            return
557        # compute "obs" pattern
558        ttArr = np.arange(max(0.5,
559                            simParms['minTT']-simParms['calcwid']/1.5),
560                            simParms['maxTT']+simParms['calcwid']/1.5,
561                            simParms['step'])
562        intArr = np.zeros_like(ttArr)
563        peaklist = np.linspace(simParms['minTT'],simParms['maxTT'],
564                               simParms['npeaks'],endpoint=True)
565        peakSpacing = (peaklist[-1]-peaklist[0])/(len(peaklist)-1)
566        NISTpk = setupFPAcalc()
567        minPtsHM = len(intArr)  # initialize points above half-max
568        maxPtsHM = 0
569        for num,twoth_peak in enumerate(peaklist):
570            try:
571                center_bin_idx,peakObj = doFPAcalc(
572                    NISTpk,ttArr,twoth_peak,simParms['calcwid'],
573                    simParms['step'])
574            except Exception as err:
575                if msg: msg += '\n'
576                msg = "Error computing convolution, revise input. Error =\n"+str(err)
577                continue
578            if num == 0: G2plt.PlotFPAconvolutors(G2frame,NISTpk)
579            pkMax = peakObj.peak.max()
580            pkPts = len(peakObj.peak)
581            minPtsHM = min(minPtsHM,sum(peakObj.peak >= 0.5*pkMax)) # points above half-max
582            maxPtsHM = max(maxPtsHM,sum(peakObj.peak >= 0.5*pkMax)) # points above half-max
583            startInd = center_bin_idx-(pkPts//2) #this should be the aligned start of the new data
584            # scale peak so max I=10,000 and add into intensity array
585            if startInd < 0:
586                intArr[:startInd+pkPts] += 10000 * peakObj.peak[-startInd:]/pkMax
587            elif startInd > len(intArr):
588                break
589            elif startInd+pkPts >= len(intArr):
590                offset = pkPts - len( intArr[startInd:] )
591                intArr[startInd:startInd+pkPts-offset] += 10000 * peakObj.peak[:-offset]/pkMax
592            else:
593                intArr[startInd:startInd+pkPts] += 10000 * peakObj.peak/pkMax
594        # check if peaks are too closely spaced
595        if maxPtsHM*simParms['step'] > peakSpacing/4:
596            if msg: msg += '\n'
597            msg += 'Maximum FWHM ({}) is too large compared to the peak spacing ({}). Decrease number of peaks or increase data range.'.format(
598                maxPtsHM*simParms['step'], peakSpacing)
599        # check if too few points across Hmax
600        if minPtsHM < 10:
601            if msg: msg += '\n'
602            msg += 'There are only {} points above the half-max. 10 are needed. Dropping step size.'.format(minPtsHM)
603            simParms['step'] *= 0.5
604        if msg:
605            G2G.G2MessageBox(dlg,msg,'Bad input, try again')
606            wx.CallAfter(MakeSimSizer,G2frame, dlg)
607            return
608        # pattern has been computed successfully
609        dlg.Destroy()
610        wx.CallAfter(FitFPApeaks,ttArr, intArr, peaklist, maxPtsHM) # do peakfit outside event callback
611
612    def FitFPApeaks(ttArr, intArr, peaklist, maxPtsHM):
613        '''Perform a peak fit to the FP simulated pattern
614        '''
615        pgbar = wx.ProgressDialog('FPA Simulation','Starting FPA simulation',100,
616            parent=G2frame,
617            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE #|wx.PD_CAN_ABORT
618            )
619        pgbar.Raise()
620        wx.BeginBusyCursor()
621        # pick out one or two most intense wavelengths
622        ints = list(NISTparms['emission']['emiss_intensities'])
623        Lam1 = NISTparms['emission']['emiss_wavelengths'][np.argmax(ints)]*1e10
624        if 'two_theta_mono' in NISTparms['emission']: # is there an IBM?
625            Lam2 = None   # Yes, ~monochromatic
626        else:  # get lambda 2
627            ints[np.argmax(ints)] = -1
628            Lam2 = NISTparms['emission']['emiss_wavelengths'][np.argmax(ints)]*1e10
629        histId = G2frame.AddSimulatedPowder(ttArr,intArr,
630                                       'NIST Fundamental Parameters simulation',
631                                       Lam1,Lam2)
632        controls = G2frame.GPXtree.GetItemPyData(
633            G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Controls'))
634        controldat = controls.get('data',
635                            {'deriv type':'analytic','min dM/M':0.001,})  #fil
636        Parms,Parms2 = G2frame.GPXtree.GetItemPyData(
637            G2gd.GetGPXtreeItemId(G2frame,histId,'Instrument Parameters'))
638        peakData = G2frame.GPXtree.GetItemPyData(
639            G2gd.GetGPXtreeItemId(G2frame,histId,'Peak List'))
640        # set background to 0 with one term = 0; disable refinement
641        bkg1,bkg2 = bkg = G2frame.GPXtree.GetItemPyData(
642            G2gd.GetGPXtreeItemId(G2frame,histId,'Background'))
643        bkg1[1]=False
644        bkg1[2]=0
645        bkg1[3]=0.0
646        limits = G2frame.GPXtree.GetItemPyData(
647            G2gd.GetGPXtreeItemId(G2frame,histId,'Limits'))
648        # approximate asym correction
649        try:
650            Parms['SH/L'][1] = 0.25 * (
651                NISTparms['axial']['length_sample']+
652                NISTparms['axial']['slit_length_source']
653                        ) / NISTparms['']['diffractometer_radius']
654        except:
655            pass
656        pgbar.Update(5,newmsg='Creating peak list')
657        pgbar.Raise()
658        for pos in peaklist:
659            i = ttArr.searchsorted(pos)
660            area = sum(intArr[max(0,i-maxPtsHM):min(len(intArr),i+maxPtsHM)])
661            peakData['peaks'].append(G2mth.setPeakparms(Parms,Parms2,pos,area))
662        pgbar.Update(10,newmsg='Refining peak positions')
663        histData = G2frame.GPXtree.GetItemPyData(histId)
664        # refine peak positions only
665        bxye = np.zeros(len(histData[1][1]))
666        peakData['sigDict'] = G2pwd.DoPeakFit('LSQ',peakData['peaks'],
667            bkg,limits[1],Parms,Parms2,histData[1],bxye,[],False,controldat)[0]
668        pgbar.Update(20,newmsg='Refining peak positions && areas')
669        # refine peak areas as well
670        for pk in peakData['peaks']:
671            pk[1] = True
672        peakData['sigDict'] = G2pwd.DoPeakFit('LSQ',peakData['peaks'],
673            bkg,limits[1],Parms,Parms2,histData[1],bxye,[],False,controldat)[0]
674        pgbar.Update(40,newmsg='Refining profile function')
675        # refine profile function
676        for p in ('U', 'V', 'W', 'X', 'Y'):
677            Parms[p][2] = True
678        peakData['sigDict'] = G2pwd.DoPeakFit('LSQ',peakData['peaks'],
679            bkg,limits[1],Parms,Parms2,histData[1],bxye,[],False,controldat)[0]
680        pgbar.Update(70,newmsg='Refining profile function && asymmetry')
681        # add in asymmetry
682        Parms['SH/L'][2] = True
683        peakData['sigDict'] = G2pwd.DoPeakFit('LSQ',peakData['peaks'],
684            bkg,limits[1],Parms,Parms2,histData[1],bxye,[],False,controldat)[0]
685        pgbar.Update(100,newmsg='Done')
686        # reset "initial" profile
687        for p in Parms:
688            if len(Parms[p]) == 3:
689                Parms[p][0] = Parms[p][1]
690                Parms[p][2] = False
691        pgbar.Destroy()
692        wx.EndBusyCursor()
693        # save Iparms
694        pth = G2G.GetExportPath(G2frame)
695        fldlg = wx.FileDialog(G2frame, 'Set name to save GSAS-II instrument parameters file', pth, '', 
696            'instrument parameter files (*.instprm)|*.instprm',wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
697        try:
698            if fldlg.ShowModal() == wx.ID_OK:
699                filename = fldlg.GetPath()
700                # make sure extension is .instprm
701                filename = os.path.splitext(filename)[0]+'.instprm'
702                File = open(filename,'w')
703                File.write("#GSAS-II instrument parameter file; do not add/delete items!\n")
704                for item in Parms:
705                    File.write(item+':'+str(Parms[item][1])+'\n')
706                File.close()
707                print ('Instrument parameters saved to: '+filename)
708        finally:
709            fldlg.Destroy()
710        #GSASIIpath.IPyBreak()
711       
712    def _onClose(event):
713        dlg.Destroy()
714    def SetButtonStatus(done=False):
715        OKbtn.Enable(bool(NISTparms))
716        #saveBtn.Enable(bool(NISTparms))
717        if done: _onOK(None)
718    def _onSetFPA(event):
719        # Create a non-modal dialog for Topas-style FP input.
720        FPdlg = wx.Dialog(dlg,wx.ID_ANY,'FPA parameters',
721                style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
722        MakeTopasFPASizer(G2frame,FPdlg,SetButtonStatus)
723        FPdlg.CenterOnParent()
724        FPdlg.Raise()
725        FPdlg.Show() 
726    def _onSaveFPA(event):
727        filename = G2G.askSaveFile(G2frame,'','.NISTfpa',
728                                       'dict of NIST FPA values',dlg)
729        writeNIST(filename)
730       
731    def _onReadFPA(event):
732        filename = G2G.GetImportFile(G2frame,
733                message='Read file with dict of values for NIST Fundamental Parameters',
734                parent=dlg,
735                wildcard='dict of NIST FPA values|*.NISTfpa')
736        if not filename: return
737        if not filename[0]: return
738        try:
739            txt = open(filename[0],'r').read()
740            NISTparms.clear()
741            d = eval(txt)
742            NISTparms.update(d)
743        except Exception as err:
744            G2G.G2MessageBox(dlg,
745                    u'Error reading file {}:{}\n'.format(filename,err),
746                    'Bad dict input')
747        #GSASIIpath.IPyBreak()
748        SetButtonStatus()
749
750    if dlg.GetSizer(): dlg.GetSizer().Clear(True)
751    MainSizer = wx.BoxSizer(wx.VERTICAL)
752    MainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
753            'Fit Profile Parameters to Peaks from Fundamental Parameters',
754            style=wx.ALIGN_CENTER),0,wx.EXPAND)
755    MainSizer.Add((-1,5))
756    prmSizer = wx.FlexGridSizer(cols=2,hgap=3,vgap=5)
757    text = wx.StaticText(dlg,wx.ID_ANY,'value',style=wx.ALIGN_CENTER)
758    text.SetBackgroundColour(wx.WHITE)
759    prmSizer.Add(text,0,wx.EXPAND)
760    text = wx.StaticText(dlg,wx.ID_ANY,'explanation',style=wx.ALIGN_CENTER)
761    text.SetBackgroundColour(wx.WHITE)
762    prmSizer.Add(text,0,wx.EXPAND)
763    for key,defVal,text in (
764            ('minTT',3.,'Location of first peak in 2theta (deg)'),
765            ('maxTT',123.,'Location of last peak in 2theta (deg)'),
766            ('step',0.01,'Pattern step size (deg 2theta)'),
767            ('npeaks',13,'Number of peaks'),
768            ('calcwid',2.,'Range to compute each peak (deg 2theta)'),
769            ):
770        if key not in simParms: simParms[key] = defVal
771        ctrl = G2G.ValidatedTxtCtrl(dlg,simParms,key,size=(70,-1))
772        prmSizer.Add(ctrl,1,wx.ALL|WACV,1)
773        txt = wx.StaticText(dlg,wx.ID_ANY,text,size=(300,-1))
774        txt.Wrap(280)
775        prmSizer.Add(txt)
776    MainSizer.Add(prmSizer)
777
778    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
779    btn = wx.Button(dlg, wx.ID_ANY,'Input FP vals')
780    btnsizer.Add(btn)
781    btn.Bind(wx.EVT_BUTTON,_onSetFPA)
782    #saveBtn = wx.Button(dlg, wx.ID_ANY,'Save FPA dict')
783    #btnsizer.Add(saveBtn)
784    #saveBtn.Bind(wx.EVT_BUTTON,_onSaveFPA)
785    readBtn = wx.Button(dlg, wx.ID_ANY,'Read FPA dict')
786    btnsizer.Add(readBtn)
787    readBtn.Bind(wx.EVT_BUTTON,_onReadFPA)
788    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
789    MainSizer.Add((-1,4),1,wx.EXPAND,1)
790    txt = wx.StaticText(dlg,wx.ID_ANY,'If you use this, please cite: '+Citation,size=(350,-1))
791    txt.Wrap(340)
792    MainSizer.Add(txt,0,wx.ALIGN_CENTER)
793    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
794    OKbtn = wx.Button(dlg, wx.ID_OK)
795    OKbtn.SetDefault()
796    btnsizer.Add(OKbtn)
797    Cbtn = wx.Button(dlg, wx.ID_CLOSE,"Cancel") 
798    btnsizer.Add(Cbtn)
799    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
800    MainSizer.Add((-1,4),1,wx.EXPAND,1)
801    # bindings for close of window
802    OKbtn.Bind(wx.EVT_BUTTON,_onOK)
803    Cbtn.Bind(wx.EVT_BUTTON,_onClose)
804    SetButtonStatus()
805    dlg.SetSizer(MainSizer)
806    MainSizer.Layout()
807    MainSizer.Fit(dlg)
808    dlg.SetMinSize(dlg.GetSize())
809    dlg.SendSizeEvent()
810    dlg.Raise()
811   
812def GetFPAInput(G2frame):
813    dlg = wx.Dialog(G2frame,wx.ID_ANY,'FPA input',
814            style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
815    MakeSimSizer(G2frame,dlg)
816    dlg.CenterOnParent()
817    dlg.Show()
818    return
819       
820if __name__ == "__main__":
821    app = wx.PySimpleApp()
822    GSASIIpath.InvokeDebugOpts()
823    frm = wx.Frame(None) # create a frame
824    frm.Show(True)
825    frm.TutorialImportDir = '/tmp'
826
827    GetFPAInput(frm)
828   
829    app.MainLoop()
Note: See TracBrowser for help on using the repository browser.