source: trunk/GSASIIfpaGUI.py @ 4403

Last change on this file since 4403 was 4403, checked in by toby, 3 years ago

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