source: trunk/GSASIIfpaGUI.py @ 4894

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

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