source: trunk/GSASIIfpaGUI.py @ 4898

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

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