source: trunk/GSASIIfpaGUI.py

Last change on this file was 5082, checked in by toby, 11 months ago

fix scaling of ISODISTORT displacement modes so esds are correct; minor formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 35.5 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-11-14 16:35:31 +0000 (Sun, 14 Nov 2021) $
4# $Author: vondreele $
5# $Revision: 5082 $
6# $URL: trunk/GSASIIfpaGUI.py $
7# $Id: GSASIIfpaGUI.py 5082 2021-11-14 16:35:31Z vondreele $
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    ('SiPSD_th2_angular_range', 3.0, 'Angular range observed by PSD (degrees 2Theta)'),]
81'''Additional FPA dict entries used in :func:`FillParmSizer`
82needed for Bragg Brentano instruments with linear (1-D) Si PSD detectors.
83'''
84
85IBmonoParms = [
86    ('src_mono_mm',119,'Distance from xray line source to monochromator crystal (mm)'),
87    ('focus_mono_mm',217,'Distance from monochromator crystal to focus slit (mm)'),
88    ('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'),
89    ('mono_src_proj_mn',51,'Projection width of line-focus xray source on the monochromator along beam direction (microns), sets bandwidth'),
90    ('passband_shoulder',.087,'Width of transition region from flat top to tails, in units of the bandwidth'),
91    ('two_theta_mono',27.27,'The full diffraction angle of the IBM crystal, e.g. double 2theta-Bragg for the mono (deg)'),
92    ('mono_slit_attenuation',0.03,'Attenuation of Cu K alpha2 lines due to focal slit'),
93    ]
94'''Additional FPA dict entries used in :func:`FillParmSizer`, needed for Incident Beam Monochromator
95'''
96   
97Citation = '''MH Mendenhall, K Mullen && JP Cline (2015), J. Res. of NIST, 120, p223. DOI: 10.6028/jres.120.014
98
99For Incident Beam Mono model, also cite: MH Mendenhall, D Black && JP Cline (2019), J. Appl. Cryst., 52, p1087. DOI: 10.1107/S1600576719010951
100'''
101
102IBmono = False
103'''set to True if an incident beam monochromator is in use
104'''
105DetMode = 'BBpoint'
106'''The type of detector, either 'BBpoint' for Bragg-Brentano point detector or
107      or 'BBPSD' (linear) position sensitive detector
108'''
109
110def SetCu2Wave():
111    '''Set the parameters to the two-line Cu K alpha 1+2 spectrum
112    '''
113    parmDict['wave'] = {i:v for i,v in enumerate((1.5405925, 1.5443873))}
114    parmDict['int'] = {i:v for i,v in enumerate((0.653817, 0.346183))}
115    parmDict['lwidth'] = {i:v for i,v in enumerate((0.501844,0.626579))}
116
117def SetCu6wave():
118    '''Set the emission parameters to the NIST six-line 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    txtSizer = wx.BoxSizer(wx.HORIZONTAL)
183    text = wx.StaticText(prmPnl,wx.ID_ANY,'explanation',style=wx.ALIGN_CENTER)
184    text.SetBackgroundColour(wx.WHITE)
185    txtSizer.Add(text,1,wx.EXPAND)
186    txtSizer.Add(G2G.HelpButton(prmPnl,helpIndex='FPAinput'))
187    prmSizer.Add(txtSizer,0,wx.EXPAND)
188    for lbl,defVal,text in itemList:
189        prmSizer.Add(wx.StaticText(prmPnl,wx.ID_ANY,lbl),1,wx.ALIGN_RIGHT|WACV,1)
190        if lbl not in parmDict: parmDict[lbl] = defVal
191        ctrl = G2G.ValidatedTxtCtrl(prmPnl,parmDict,lbl,size=(70,-1))
192        prmSizer.Add(ctrl,1,wx.ALL|WACV,1)
193        txt = wx.StaticText(prmPnl,wx.ID_ANY,text,size=(400,-1))
194        txt.Wrap(380)
195        prmSizer.Add(txt)
196    px1,py = prmSizer.GetSize()
197    prmSizer.Layout()
198    FPdlg = prmPnl.GetParent()
199    FPdlg.SendSizeEvent()
200       
201def MakeTopasFPASizer(G2frame,FPdlg,SetButtonStatus):
202    '''Create a GUI with parameters for the NIST XRD Fundamental Parameters Code.
203    Parameter input is modeled after Topas input parameters.
204
205    :param wx.Frame G2frame: main GSAS-II window
206    :param wx.Window FPdlg: Frame or Dialog where GUI will appear
207    :param SetButtonStatus: a callback function to call to see what buttons in
208      this windows can be enabled. Called with done=True to trigger closing
209      the parent window as well.
210    :returns: a sizer with the GUI controls
211    '''
212    def _onOK(event):
213        XferFPAsettings(parmDict)
214        SetButtonStatus(done=True) # done=True triggers the simulation
215        FPdlg.Destroy()
216    def _onClose(event):
217        SetButtonStatus()
218        FPdlg.Destroy()
219    def _onAddWave(event):
220        newkey = max(parmDict['wave'].keys())+1
221        for key,defVal in zip(
222            ('wave','int','lwidth'),
223            (0.0,   1.0,   0.1),
224            ):
225            parmDict[key][newkey] = defVal
226        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
227    def _onRemWave(event):
228        lastkey = max(parmDict['wave'].keys())
229        for key in ('wave','int','lwidth'):
230            if lastkey in parmDict[key]:
231                del parmDict[key][lastkey]
232        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
233    def _onSetCu6wave(event):
234        SetCu6wave()
235        SetMonoWave()
236        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
237    def _onSetCu2Wave(event):
238        SetCu2Wave()
239        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
240    def _onSetDetBtn(event):
241        global DetMode
242        if detBtn1.GetValue():
243            DetMode = 'BBpoint'
244            wx.CallAfter(FillParmSizer)
245        else:
246            DetMode = 'BBPSD'
247            wx.CallAfter(FillParmSizer)
248    def _onSetMonoBtn(event):
249        global IBmono
250        IBmono = not monoBtn1.GetValue()
251        SetMonoWave()
252        #wx.CallAfter(FillParmSizer)
253        wx.CallAfter(MakeTopasFPASizer,G2frame,FPdlg,SetButtonStatus)
254    def PlotTopasFPA(event):
255        XferFPAsettings(parmDict)
256        ttArr = np.arange(max(0.5,
257                            simParms['plotpos']-simParms['calcwid']),
258                            simParms['plotpos']+simParms['calcwid'],
259                            simParms['step'])
260        intArr = np.zeros_like(ttArr)
261        NISTpk = setupFPAcalc()
262        try:
263            center_bin_idx,peakObj = doFPAcalc(
264                NISTpk,ttArr,simParms['plotpos'],simParms['calcwid'],
265                simParms['step'])
266        except Exception as err:
267            msg = "Error computing convolution, revise input"
268            print(msg)
269            print(err)
270            return
271        G2plt.PlotFPAconvolutors(G2frame,NISTpk)
272        pkPts = len(peakObj.peak)
273        pkMax = peakObj.peak.max()
274        startInd = center_bin_idx-(pkPts//2) #this should be the aligned start of the new data
275        # scale peak so max I=10,000 and add into intensity array
276        if startInd < 0:
277            intArr[:startInd+pkPts] += 10000 * peakObj.peak[-startInd:]/pkMax
278        elif startInd > len(intArr):
279            return
280        elif startInd+pkPts >= len(intArr):
281            offset = pkPts - len( intArr[startInd:] )
282            intArr[startInd:startInd+pkPts-offset] += 10000 * peakObj.peak[:-offset]/pkMax
283        else:
284            intArr[startInd:startInd+pkPts] += 10000 * peakObj.peak/pkMax
285        G2plt.PlotXY(G2frame, [(ttArr, intArr)],
286                     labelX=r'$2\theta, deg$',
287                     labelY=r'Intensity (arbitrary)',
288                     Title='FPA peak', newPlot=True, lines=True)
289    def _onSaveFPA(event):
290        XferFPAsettings(parmDict)
291        filename = G2G.askSaveFile(G2frame,'','.NISTfpa',
292                                       'dict of NIST FPA values',FPdlg)
293        writeNIST(filename)
294
295    if FPdlg.GetSizer(): FPdlg.GetSizer().Clear(True)
296    MainSizer = wx.BoxSizer(wx.VERTICAL)
297    MainSizer.Add((-1,5))
298    waveSizer = wx.FlexGridSizer(cols=len(parmDict['wave'])+1,hgap=3,vgap=5)
299    for lbl,prm,defVal in zip(
300            (u'Wavelength (\u212b)','Rel. Intensity',u'Lorentz Width\n(\u212b/1000)'),
301            ('wave','int','lwidth'),
302            (0.0,   1.0,   0.1),
303            ):
304        text = wx.StaticText(FPdlg,wx.ID_ANY,lbl,style=wx.ALIGN_CENTER)
305        text.SetBackgroundColour(wx.WHITE)
306        waveSizer.Add(text,0,wx.EXPAND)
307        if prm not in parmDict: parmDict[prm] = {}
308        for i in parmDict['wave'].keys():
309            if i not in parmDict[prm]: parmDict[prm][i] = defVal
310            if prm == 'wave':
311                ctrl = G2G.ValidatedTxtCtrl(FPdlg,parmDict[prm],i,size=(90,-1),nDig=(10,6))
312            else:
313                ctrl = G2G.ValidatedTxtCtrl(FPdlg,parmDict[prm],i,size=(90,-1),nDig=(10,3))
314            waveSizer.Add(ctrl,1,WACV,1)
315    MainSizer.Add(waveSizer)
316    MainSizer.Add((-1,5))
317    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
318    btn = wx.Button(FPdlg, wx.ID_ANY,'Add wave')
319    btnsizer.Add(btn)
320    btn.Bind(wx.EVT_BUTTON,_onAddWave)
321    btn = wx.Button(FPdlg, wx.ID_ANY,'Remove wave')
322    btnsizer.Add(btn)
323    btn.Bind(wx.EVT_BUTTON,_onRemWave)
324    btn = wx.Button(FPdlg, wx.ID_ANY,'CuKa1+2')
325    btnsizer.Add(btn)
326    btn.Bind(wx.EVT_BUTTON,_onSetCu2Wave)
327    btn = wx.Button(FPdlg, wx.ID_ANY,'NIST CuKa')
328    btnsizer.Add(btn)
329    btn.Bind(wx.EVT_BUTTON,_onSetCu6wave)
330    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
331    MainSizer.Add((-1,5))
332   
333    btnsizer = wx.GridBagSizer( 2, 5)
334    btnsizer.Add( wx.StaticText(FPdlg, wx.ID_ANY, 'Detector type'),
335                 (0,0), (2,1), wx.ALIGN_CENTER | wx.ALL, 5)   
336    detBtn1 = wx.RadioButton(FPdlg,wx.ID_ANY,'Point',style=wx.RB_GROUP)
337    detBtn1.SetValue(DetMode == 'BBpoint')
338    btnsizer.Add(detBtn1, (0,1))
339    detBtn1.Bind(wx.EVT_RADIOBUTTON,_onSetDetBtn)
340    detBtn2 = wx.RadioButton(FPdlg,wx.ID_ANY,'PSD')
341    detBtn2.SetValue(not DetMode == 'BBpoint')
342    btnsizer.Add(detBtn2, (1,1))
343    detBtn2.Bind(wx.EVT_RADIOBUTTON,_onSetDetBtn)
344    btnsizer.Add( (40,-1), (0,2), (1,1), wx.ALIGN_CENTER | wx.ALL, 5)   
345    btnsizer.Add( wx.StaticText(FPdlg, wx.ID_ANY, 'Incident Beam Mono'),
346                 (0,3), (2,1), wx.ALIGN_CENTER | wx.ALL, 5)   
347    monoBtn1 = wx.RadioButton(FPdlg,wx.ID_ANY,'No',style=wx.RB_GROUP)
348    monoBtn1.SetValue(not IBmono)
349    btnsizer.Add(monoBtn1, (0,4))
350    monoBtn1.Bind(wx.EVT_RADIOBUTTON,_onSetMonoBtn)
351    monoBtn2 = wx.RadioButton(FPdlg,wx.ID_ANY,'Yes')
352    monoBtn2.SetValue(IBmono)
353    btnsizer.Add(monoBtn2, (1,4))
354    monoBtn2.Bind(wx.EVT_RADIOBUTTON,_onSetMonoBtn)
355    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
356
357    global prmPnl
358    prmPnl = wxscroll.ScrolledPanel(FPdlg, wx.ID_ANY, #size=(200,200),
359        style = wx.TAB_TRAVERSAL|wx.SUNKEN_BORDER)
360    prmSizer = wx.FlexGridSizer(cols=3,hgap=3,vgap=5)
361    prmPnl.SetSizer(prmSizer)
362    FillParmSizer()   
363    MainSizer.Add(prmPnl,1,wx.EXPAND,1)
364    prmPnl.SetAutoLayout(1)
365    prmPnl.SetupScrolling()
366   
367    MainSizer.Add((-1,4))
368    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
369    btn = wx.Button(FPdlg, wx.ID_ANY, 'Plot peak')
370    btnsizer.Add(btn)
371    btn.Bind(wx.EVT_BUTTON,PlotTopasFPA)
372    btnsizer.Add(wx.StaticText(FPdlg,wx.ID_ANY,' at '))
373    if 'plotpos' not in simParms: simParms['plotpos'] =  simParms['minTT']
374    ctrl = G2G.ValidatedTxtCtrl(FPdlg,simParms,'plotpos',size=(70,-1))
375    btnsizer.Add(ctrl)
376    btnsizer.Add(wx.StaticText(FPdlg,wx.ID_ANY,' deg.  '))   
377    saveBtn = wx.Button(FPdlg, wx.ID_ANY,'Save FPA dict')
378    btnsizer.Add(saveBtn)
379    saveBtn.Bind(wx.EVT_BUTTON,_onSaveFPA)
380    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
381    MainSizer.Add((-1,4))
382    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
383    OKbtn = wx.Button(FPdlg, wx.ID_OK)
384    OKbtn.SetDefault()
385    btnsizer.Add(OKbtn)
386    Cbtn = wx.Button(FPdlg, wx.ID_CLOSE,"Cancel") 
387    btnsizer.Add(Cbtn)
388    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
389    MainSizer.Add((-1,4))
390    # bindings for close of window
391    OKbtn.Bind(wx.EVT_BUTTON,_onOK)
392    Cbtn.Bind(wx.EVT_BUTTON,_onClose)
393    FPdlg.SetSizer(MainSizer)
394    MainSizer.Layout()
395    MainSizer.Fit(FPdlg)
396    # control window size
397    px,py = prmSizer.GetSize()
398    dx,dy = FPdlg.GetSize()
399    FPdlg.SetMinSize((-1,-1))
400    FPdlg.SetMaxSize((-1,-1)) 
401    FPdlg.SetMinSize((dx,dy+200)) # leave a min of 200 points for scroll panel
402    FPdlg.SetMaxSize((max(dx,700),850)) 
403    FPdlg.SetSize((max(dx,px+20),min(750,dy+py+30))) # 20 for scroll bar, 30 for a bit of room at bottom
404   
405def XferFPAsettings(InpParms):
406    '''convert Topas-type parameters to SI units for NIST and place in a dict sorted
407    according to use in each convoluter
408
409    :param dict InpParms: a dict with Topas-like parameters, as set in
410      :func:`MakeTopasFPASizer`
411    :returns: a nested dict with global parameters and those for each convolution
412    '''
413    # cleanup old stuff
414    for key in "tube_tails","absorption","si_psd","displacement","receiver_slit":
415        if key in NISTparms:
416            del NISTparms[key]
417   
418    keys = list(InpParms['wave'].keys())
419    source_wavelengths_m = 1.e-10 * np.array([InpParms['wave'][i] for i in keys])
420    la = [InpParms['int'][i] for i in keys]
421 
422    if IBmono:  # kludge: apply mono_slit_attenuation since it is not part of NIST FPA code
423        norm = [InpParms['mono_slit_attenuation'] if
424                    (1.5443 < InpParms['wave'][i] < 1.5447) else 1.
425                for i in keys]
426        source_intensities = norm * np.array(la)/max(la)
427    else:
428        source_intensities = np.array(la)/max(la)
429    source_lor_widths_m = 1.e-10 * 1.e-3 * np.array([InpParms['lwidth'][i] for i in keys])
430    source_gauss_widths_m = 1.e-10 * 1.e-3 * np.array([0.001 for i in keys])
431   
432    NISTparms["emission"] = {'emiss_wavelengths' : source_wavelengths_m,
433                'emiss_intensities' : source_intensities,
434                'emiss_gauss_widths' : source_gauss_widths_m,
435                'emiss_lor_widths' : source_lor_widths_m,
436                'crystallite_size_gauss' : 1.e-9 * InpParms.get('Size_G',1e6),
437                'crystallite_size_lor' : 1.e-9 * InpParms.get('Size_L',1e6)}
438    if IBmono:
439        NISTparms["emission"]['a_mono'] = InpParms['src_mono_mm'] * 10**-3
440        NISTparms["emission"]['b_mono'] = InpParms['focus_mono_mm'] * 10**-3
441        NISTparms["emission"]['ibm_source_width'] = InpParms['mono_src_proj_mn'] * 10**-6
442        for i in  ('passband_mistune','passband_shoulder','two_theta_mono'):
443            NISTparms["emission"][i] = InpParms[i]
444    elif InpParms.get('source_width', 0) > 0 and InpParms.get(
445            'tube-tails_rel-I',0) > 0:
446        NISTparms["tube_tails"] = {
447            'main_width' : 1e-3 * InpParms.get('source_width', 0.),
448            'tail_left' : -1e-3 * InpParms.get('tube-tails_L-tail',0.),
449            'tail_right' : 1e-3 * InpParms.get('tube-tails_R-tail',0.),
450            'tail_intens' : InpParms.get('tube-tails_rel-I',0.),}
451   
452    if InpParms['filament_length'] == InpParms['receiving_slit_length']: # workaround:
453        InpParms['receiving_slit_length'] *= 1.00001 # avoid bug when slit lengths are identical
454    NISTparms["axial"]  = {
455            'axDiv':"full", 'slit_length_source' : 1e-3*InpParms['filament_length'],
456            'slit_length_target' : 1e-3*InpParms['receiving_slit_length'],
457            'length_sample' : 1e-3 * InpParms['sample_length'], 
458            'n_integral_points' : 10,
459            'angI_deg' : InpParms['soller_angle'],
460            'angD_deg': InpParms['soller_angle']
461            }
462    if InpParms.get('LAC_cm',0) > 0:
463        NISTparms["absorption"] = {
464            'absorption_coefficient': InpParms['LAC_cm']*100, #like LaB6, in m^(-1)
465            'sample_thickness': 1e-3 * InpParms['sample_thickness'],
466            }
467
468    if InpParms.get('SiPSD_th2_angular_range',0) > 0 and DetMode == 'BBPSD':
469        PSDdetector_length_mm=np.arcsin(np.pi*InpParms['SiPSD_th2_angular_range']/180.
470                                            )*InpParms['Rs'] # mm
471        NISTparms["si_psd"] = {
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    topSizer = wx.BoxSizer(wx.HORIZONTAL)
753    topSizer.Add(wx.StaticText(dlg,wx.ID_ANY,
754            'Fit Profile Parameters to Peaks from Fundamental Parameters',
755            style=wx.ALIGN_CENTER),0)
756    topSizer.Add((5,-1),1,wx.EXPAND)
757    topSizer.Add(G2G.HelpButton(dlg,helpIndex='FPA'))
758    MainSizer.Add(topSizer,0,wx.EXPAND)
759    G2G.HorizontalLine(MainSizer,dlg)
760    MainSizer.Add((5,5),0)   
761    prmSizer = wx.FlexGridSizer(cols=2,hgap=3,vgap=5)
762    text = wx.StaticText(dlg,wx.ID_ANY,'value',style=wx.ALIGN_CENTER)
763    text.SetBackgroundColour(wx.WHITE)
764    prmSizer.Add(text,0,wx.EXPAND)
765    text = wx.StaticText(dlg,wx.ID_ANY,'explanation',style=wx.ALIGN_CENTER)
766    text.SetBackgroundColour(wx.WHITE)
767    prmSizer.Add(text,0,wx.EXPAND)
768    for key,defVal,text in (
769            ('minTT',3.,'Location of first peak in 2theta (deg)'),
770            ('maxTT',123.,'Location of last peak in 2theta (deg)'),
771            ('step',0.01,'Pattern step size (deg 2theta)'),
772            ('npeaks',13,'Number of peaks'),
773            ('calcwid',2.,'Range to compute each peak (deg 2theta)'),
774            ):
775        if key not in simParms: simParms[key] = defVal
776        ctrl = G2G.ValidatedTxtCtrl(dlg,simParms,key,size=(70,-1))
777        prmSizer.Add(ctrl,1,wx.ALL|WACV,1)
778        txt = wx.StaticText(dlg,wx.ID_ANY,text,size=(300,-1))
779        txt.Wrap(280)
780        prmSizer.Add(txt)
781    MainSizer.Add(prmSizer)
782
783    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
784    btn = wx.Button(dlg, wx.ID_ANY,'Input FP vals')
785    btnsizer.Add(btn)
786    btn.Bind(wx.EVT_BUTTON,_onSetFPA)
787    #saveBtn = wx.Button(dlg, wx.ID_ANY,'Save FPA dict')
788    #btnsizer.Add(saveBtn)
789    #saveBtn.Bind(wx.EVT_BUTTON,_onSaveFPA)
790    readBtn = wx.Button(dlg, wx.ID_ANY,'Read FPA dict')
791    btnsizer.Add(readBtn)
792    readBtn.Bind(wx.EVT_BUTTON,_onReadFPA)
793    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
794    MainSizer.Add((-1,4),1,wx.EXPAND,1)
795    txt = wx.StaticText(dlg,wx.ID_ANY,'If you use this, please cite: '+Citation,size=(350,-1))
796    txt.Wrap(340)
797    MainSizer.Add(txt,0,wx.ALIGN_CENTER)
798    btnsizer = wx.BoxSizer(wx.HORIZONTAL)
799    OKbtn = wx.Button(dlg, wx.ID_OK)
800    OKbtn.SetDefault()
801    btnsizer.Add(OKbtn)
802    Cbtn = wx.Button(dlg, wx.ID_CLOSE,"Cancel") 
803    btnsizer.Add(Cbtn)
804    MainSizer.Add(btnsizer, 0, wx.ALIGN_CENTER, 0)
805    MainSizer.Add((-1,4),1,wx.EXPAND,1)
806    # bindings for close of window
807    OKbtn.Bind(wx.EVT_BUTTON,_onOK)
808    Cbtn.Bind(wx.EVT_BUTTON,_onClose)
809    SetButtonStatus()
810    dlg.SetSizer(MainSizer)
811    MainSizer.Layout()
812    MainSizer.Fit(dlg)
813    dlg.SetMinSize(dlg.GetSize())
814    dlg.SendSizeEvent()
815    dlg.Raise()
816   
817def GetFPAInput(G2frame):
818    dlg = wx.Dialog(G2frame,wx.ID_ANY,'FPA input',
819            style=wx.DEFAULT_DIALOG_STYLE|wx.RESIZE_BORDER)
820    MakeSimSizer(G2frame,dlg)
821    dlg.CenterOnParent()
822    dlg.Show()
823    return
824       
825if __name__ == "__main__":
826    app = wx.PySimpleApp()
827    GSASIIpath.InvokeDebugOpts()
828    frm = wx.Frame(None) # create a frame
829    frm.Show(True)
830    frm.TutorialImportDir = '/tmp'
831    size = wx.Size(700,600)               
832    frm.plotFrame = wx.Frame(None,-1,'GSASII Plots',size=size,
833                    style=wx.DEFAULT_FRAME_STYLE ^ wx.CLOSE_BOX)
834    frm.G2plotNB = G2plt.G2PlotNoteBook(frm.plotFrame,G2frame=frm)
835    frm.plotFrame.Show()
836
837    GetFPAInput(frm)
838   
839    app.MainLoop()
Note: See TracBrowser for help on using the repository browser.