source: trunk/GSASIIfpaGUI.py @ 4893

Last change on this file since 4893 was 4893, checked in by toby, 6 months ago

FPA: add extra Save button (perhaps original should be removed)

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