source: trunk/GSASIIplot.py @ 1571

Last change on this file since 1571 was 1571, checked in by vondreele, 8 years ago

fix typos in small angle tutorials
add menu with copy & reset to instrument parameters for SASD data
copy masks now copies the lower threshold
work on indexing incommensurate powder patterns, plots

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 211.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIplot: plotting routines*
4===============================
5
6'''
7########### SVN repository information ###################
8# $Date: 2014-11-15 18:01:30 +0000 (Sat, 15 Nov 2014) $
9# $Author: vondreele $
10# $Revision: 1571 $
11# $URL: trunk/GSASIIplot.py $
12# $Id: GSASIIplot.py 1571 2014-11-15 18:01:30Z vondreele $
13########### SVN repository information ###################
14import math
15import time
16import copy
17import sys
18import os.path
19import numpy as np
20import numpy.ma as ma
21import numpy.linalg as nl
22import wx
23import wx.aui
24import wx.glcanvas
25import matplotlib as mpl
26import mpl_toolkits.mplot3d.axes3d as mp3d
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 1571 $")
29import GSASIIgrid as G2gd
30import GSASIIimage as G2img
31import GSASIIpwd as G2pwd
32import GSASIIIO as G2IO
33import GSASIIpwdGUI as G2pdG
34import GSASIIimgGUI as G2imG
35import GSASIIphsGUI as G2phG
36import GSASIIlattice as G2lat
37import GSASIIspc as G2spc
38import GSASIImath as G2mth
39import pytexture as ptx
40from  OpenGL.GL import *
41from OpenGL.GLU import *
42from OpenGL.GLE import *
43import gltext
44from matplotlib.backends.backend_wx import _load_bitmap
45from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as Canvas
46from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
47
48# useful degree trig functions
49sind = lambda x: math.sin(x*math.pi/180.)
50cosd = lambda x: math.cos(x*math.pi/180.)
51tand = lambda x: math.tan(x*math.pi/180.)
52asind = lambda x: 180.*math.asin(x)/math.pi
53acosd = lambda x: 180.*math.acos(x)/math.pi
54atan2d = lambda x,y: 180.*math.atan2(y,x)/math.pi
55atand = lambda x: 180.*math.atan(x)/math.pi
56# numpy versions
57npsind = lambda x: np.sin(x*np.pi/180.)
58npcosd = lambda x: np.cos(x*np.pi/180.)
59nptand = lambda x: np.tan(x*np.pi/180.)
60npacosd = lambda x: 180.*np.arccos(x)/np.pi
61npasind = lambda x: 180.*np.arcsin(x)/np.pi
62npatand = lambda x: 180.*np.arctan(x)/np.pi
63npatan2d = lambda x,y: 180.*np.arctan2(x,y)/np.pi
64   
65class G2PlotMpl(wx.Panel):   
66    'needs a doc string'
67    def __init__(self,parent,id=-1,dpi=None,**kwargs):
68        wx.Panel.__init__(self,parent,id=id,**kwargs)
69        mpl.rcParams['legend.fontsize'] = 10
70        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(5,6))
71        self.canvas = Canvas(self,-1,self.figure)
72        self.toolbar = GSASIItoolbar(self.canvas)
73
74        self.toolbar.Realize()
75       
76        sizer=wx.BoxSizer(wx.VERTICAL)
77        sizer.Add(self.canvas,1,wx.EXPAND)
78        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
79        self.SetSizer(sizer)
80       
81class G2PlotOgl(wx.Panel):
82    'needs a doc string'
83    def __init__(self,parent,id=-1,dpi=None,**kwargs):
84        self.figure = wx.Panel.__init__(self,parent,id=id,**kwargs)
85        if 'win' in sys.platform:           #Windows (& Mac) already double buffered
86            self.canvas = wx.glcanvas.GLCanvas(self,-1,**kwargs)
87        else:                               #fix from Jim Hester for X systems
88            attribs = (wx.glcanvas.WX_GL_DOUBLEBUFFER,)         
89            self.canvas = wx.glcanvas.GLCanvas(self,-1,attribList=attribs,**kwargs)
90        # create GL context for wx > 2.8
91        i,j= wx.__version__.split('.')[0:2]
92        if int(i)+int(j)/10. > 2.8:
93            self.context = wx.glcanvas.GLContext(self.canvas)
94            self.canvas.SetCurrent(self.context)
95        else:
96            self.context = None
97        self.camera = {}
98        sizer=wx.BoxSizer(wx.VERTICAL)
99        sizer.Add(self.canvas,1,wx.EXPAND)
100        self.SetSizer(sizer)
101       
102class G2Plot3D(wx.Panel):
103    'needs a doc string'
104    def __init__(self,parent,id=-1,dpi=None,**kwargs):
105        wx.Panel.__init__(self,parent,id=id,**kwargs)
106        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(6,6))
107        self.canvas = Canvas(self,-1,self.figure)
108        self.toolbar = GSASIItoolbar(self.canvas)
109
110        self.toolbar.Realize()
111       
112        sizer=wx.BoxSizer(wx.VERTICAL)
113        sizer.Add(self.canvas,1,wx.EXPAND)
114        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
115        self.SetSizer(sizer)
116                             
117class G2PlotNoteBook(wx.Panel):
118    'create a tabbed window for plotting'
119    def __init__(self,parent,id=-1):
120        wx.Panel.__init__(self,parent,id=id)
121        #so one can't delete a plot page!!
122        self.nb = wx.aui.AuiNotebook(self, \
123            style=wx.aui.AUI_NB_DEFAULT_STYLE ^ wx.aui.AUI_NB_CLOSE_ON_ACTIVE_TAB)
124        sizer = wx.BoxSizer()
125        sizer.Add(self.nb,1,wx.EXPAND)
126        self.SetSizer(sizer)
127        self.status = parent.CreateStatusBar()
128        self.status.SetFieldsCount(2)
129        self.status.SetStatusWidths([150,-1])
130        self.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, self.OnPageChanged)
131        self.nb.Bind(wx.EVT_KEY_UP,self.OnNotebookKey)
132       
133        self.plotList = []
134           
135    def OnNotebookKey(self,event):
136        '''Called when a keystroke event gets picked up by the notebook window
137        rather the child. This is not expected, but somehow it does sometimes
138        on the Mac and perhaps Linux.
139
140        Assume that the page associated with the currently displayed tab
141        has a child, .canvas; give that child the focus and pass it the event.
142        '''
143        try:
144            Page = self.nb.GetPage(self.nb.GetSelection())
145        except ValueError: # occurs with no plot tabs
146            return
147        try:
148            Page.canvas.SetFocus()
149            wx.PostEvent(Page.canvas,event)
150        except AttributeError:
151            pass
152
153    def addMpl(self,name=""):
154        'Add a tabbed page with a matplotlib plot'
155        page = G2PlotMpl(self.nb)
156        self.nb.AddPage(page,name)
157       
158        self.plotList.append(name)
159       
160        return page.figure
161       
162    def add3D(self,name=""):
163        'Add a tabbed page with a 3D plot'
164        page = G2Plot3D(self.nb)
165        self.nb.AddPage(page,name)
166       
167        self.plotList.append(name)
168       
169        return page.figure
170       
171    def addOgl(self,name=""):
172        'Add a tabbed page with an openGL plot'
173        page = G2PlotOgl(self.nb)
174        self.nb.AddPage(page,name)
175       
176        self.plotList.append(name)
177       
178        return page.figure
179       
180    def Delete(self,name):
181        'delete a tabbed page'
182        try:
183            item = self.plotList.index(name)
184            del self.plotList[item]
185            self.nb.DeletePage(item)
186        except ValueError:          #no plot of this name - do nothing
187            return     
188               
189    def clear(self):
190        'clear all pages from plot window'
191        while self.nb.GetPageCount():
192            self.nb.DeletePage(0)
193        self.plotList = []
194        self.status.DestroyChildren()
195       
196    def Rename(self,oldName,newName):
197        'rename a tab'
198        try:
199            item = self.plotList.index(oldName)
200            self.plotList[item] = newName
201            self.nb.SetPageText(item,newName)
202        except ValueError:          #no plot of this name - do nothing
203            return     
204       
205    def OnPageChanged(self,event):
206        'respond to someone pressing a tab on the plot window'
207        if self.plotList:
208            self.status.SetStatusText('Better to select this from GSAS-II data tree',1)
209        self.status.DestroyChildren()                           #get rid of special stuff on status bar
210       
211class GSASIItoolbar(Toolbar):
212    'needs a doc string'
213    ON_MPL_HELP = wx.NewId()
214    ON_MPL_KEY = wx.NewId()
215    def __init__(self,plotCanvas):
216        Toolbar.__init__(self,plotCanvas)
217        POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6
218        self.DeleteToolByPos(POSITION_OF_CONFIGURE_SUBPLOTS_BTN)
219        parent = self.GetParent()
220        key = os.path.join(os.path.split(__file__)[0],'key.ico')
221        self.AddSimpleTool(self.ON_MPL_KEY,_load_bitmap(key),'Key press','Select key press')
222        wx.EVT_TOOL(self,self.ON_MPL_KEY,self.OnKey)
223        help = os.path.join(os.path.split(__file__)[0],'help.ico')
224        self.AddSimpleTool(self.ON_MPL_HELP,_load_bitmap(help),'Help on','Show help on')
225        wx.EVT_TOOL(self,self.ON_MPL_HELP,self.OnHelp)
226    def OnHelp(self,event):
227        'needs a doc string'
228        Page = self.GetParent().GetParent()
229        pageNo = Page.GetSelection()
230        bookmark = Page.GetPageText(pageNo)
231        bookmark = bookmark.strip(')').replace('(','_')
232        G2gd.ShowHelp(bookmark,self.TopLevelParent)
233    def OnKey(self,event):
234        'needs a doc string'
235        parent = self.GetParent()
236        if parent.Choice:
237            dlg = wx.SingleChoiceDialog(parent,'Select','Key press',list(parent.Choice))
238            if dlg.ShowModal() == wx.ID_OK:
239                sel = dlg.GetSelection()
240                event.key = parent.Choice[sel][0]
241                parent.keyPress(event)
242            dlg.Destroy()
243           
244################################################################################
245##### PlotSngl
246################################################################################
247           
248def PlotSngl(G2frame,newPlot=False,Data=None,hklRef=None,Title=''):
249    '''Structure factor plotting package - displays zone of reflections as rings proportional
250        to F, F**2, etc. as requested
251    '''
252    from matplotlib.patches import Circle,CirclePolygon
253    global HKL,HKLF
254   
255    def OnSCKeyPress(event):
256        i = zones.index(Data['Zone'])
257        newPlot = False
258        pwdrChoice = {'f':'Fo','s':'Fosq','u':'Unit Fc'}
259        hklfChoice = {'1':'|DFsq|>sig','3':'|DFsq|>3sig','w':'|DFsq|/sig','f':'Fo','s':'Fosq','i':'Unit Fc'}
260        if event.key == 'h':
261            Data['Zone'] = '100'
262            newPlot = True
263        elif event.key == 'k':
264            Data['Zone'] = '010'
265            newPlot = True
266        elif event.key == 'l':
267            Data['Zone'] = '001'
268            newPlot = True
269        elif event.key == 'i':
270            Data['Scale'] *= 1.1
271        elif event.key == 'd':
272            Data['Scale'] /= 1.1
273        elif event.key == '+':
274            Data['Layer'] = min(Data['Layer']+1,HKLmax[i])
275        elif event.key == '-':
276            Data['Layer'] = max(Data['Layer']-1,HKLmin[i])
277        elif event.key == '0':
278            Data['Layer'] = 0
279            Data['Scale'] = 1.0
280        elif event.key in hklfChoice and 'HKLF' in Name:
281            Data['Type'] = hklfChoice[event.key]           
282            newPlot = True
283        elif event.key in pwdrChoice and 'PWDR' in Name:
284            Data['Type'] = pwdrChoice[event.key]           
285            newPlot = True       
286        PlotSngl(G2frame,newPlot,Data,hklRef,Title)
287
288    def OnSCMotion(event):
289        xpos = event.xdata
290        if xpos:
291            xpos = round(xpos)                                        #avoid out of frame mouse position
292            ypos = round(event.ydata)
293            zpos = Data['Layer']
294            if '100' in Data['Zone']:
295                HKLtxt = '(%3d,%3d,%3d)'%(zpos,xpos,ypos)
296            elif '010' in Data['Zone']:
297                HKLtxt = '(%3d,%3d,%3d)'%(xpos,zpos,ypos)
298            elif '001' in Data['Zone']:
299                HKLtxt = '(%3d,%3d,%3d)'%(xpos,ypos,zpos)
300            Page.canvas.SetToolTipString(HKLtxt)
301            G2frame.G2plotNB.status.SetStatusText('HKL = '+HKLtxt,0)
302            G2frame.G2plotNB.status.SetStatusText('Use K-box to set plot controls',1)
303               
304    def OnSCPress(event):
305        zpos = Data['Layer']
306        xpos = event.xdata
307        if xpos:
308            pos = int(round(event.xdata)),int(round(event.ydata))
309            if '100' in Data['Zone']:
310                Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(zpos,pos[0],pos[1]))
311                hkl = np.array([zpos,pos[0],pos[1]])
312            elif '010' in Data['Zone']:
313                Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],zpos,pos[1]))
314                hkl = np.array([pos[0],zpos,pos[1]])
315            elif '001' in Data['Zone']:
316                Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],pos[1],zpos))
317                hkl = np.array([pos[0],pos[1],zpos])
318            h,k,l = hkl
319            hklf = HKLF[np.where(np.all(HKL-hkl == [0,0,0],axis=1))]
320            if len(hklf):
321                Fosq,sig,Fcsq = hklf[0]
322                HKLtxt = '( %.2f %.3f %.2f %.2f)'%(Fosq,sig,Fcsq,(Fosq-Fcsq)/(scale*sig))
323                G2frame.G2plotNB.status.SetStatusText('Fosq, sig, Fcsq, delFsq/sig = '+HKLtxt,1)
324                                 
325    Name = G2frame.PatternTree.GetItemText(G2frame.PatternId)
326    if not Title:
327        Title = Name
328    try:
329        plotNum = G2frame.G2plotNB.plotList.index('Structure Factors')
330        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
331        if not newPlot:
332            Plot = Page.figure.gca()          #get previous powder plot & get limits
333            xylim = Plot.get_xlim(),Plot.get_ylim()
334        Page.figure.clf()
335        Plot = Page.figure.gca()          #get a fresh plot after clf()
336    except ValueError:
337        Plot = G2frame.G2plotNB.addMpl('Structure Factors').gca()
338        plotNum = G2frame.G2plotNB.plotList.index('Structure Factors')
339        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
340        Page.canvas.mpl_connect('button_press_event', OnSCPress)
341        Page.canvas.mpl_connect('motion_notify_event', OnSCMotion)
342        Page.canvas.mpl_connect('key_press_event', OnSCKeyPress)
343        Page.keyPress = OnSCKeyPress
344        Page.Choice = (' key press','i: increase scale','d: decrease scale',
345            'h: select 100 zone','k: select 010 zone','l: select 001 zone',
346            'f: select Fo','s: select Fosq','u: select unit Fc',
347            '+: increase index','-: decrease index','0: zero layer',)
348        if 'HKLF' in Name:
349            Page.Choice += ('w: select |DFsq|/sig','1: select |DFsq|>sig','3: select |DFsq|>3sig',)
350    Page.SetFocus()
351   
352    G2frame.G2plotNB.status.SetStatusText('Use K-box to set plot controls',1)
353    Plot.set_aspect(aspect='equal')
354   
355    Type = Data['Type']           
356    scale = Data['Scale']
357    HKLmax = Data['HKLmax']
358    HKLmin = Data['HKLmin']
359    FosqMax = Data['FoMax']
360    Super = Data['Super']
361    SuperVec = []
362    for i in range(Super):
363        SuperVec.append(Data['SuperVec'][i][0])
364        SuperVec = np.array(SuperVec)
365    FoMax = math.sqrt(FosqMax)
366    xlabel = ['k, h=','h, k=','h, l=']
367    ylabel = ['l','l','k']
368    zones = ['100','010','001']
369    pzone = [[1,2],[0,2],[0,1]]
370    izone = zones.index(Data['Zone'])
371    Plot.set_title(Data['Type']+' for '+Title)
372    HKL = []
373    HKLF = []
374    time0 = time.time()
375    for refl in hklRef:
376        H = refl[:3]
377        if 'HKLF' in Name:
378            Fosq,sig,Fcsq = refl[5+Super:8+Super]
379        else:
380            Fosq,sig,Fcsq = refl[8+Super],1.0,refl[9+Super]
381        HKL.append(H+np.sum(SuperVec*refl[3:3+Super],axis=0))
382        HKLF.append([Fosq,sig,Fcsq])
383        if H[izone] == Data['Layer']:
384            A = 0
385            B = 0
386            if Type == 'Fosq':
387                A = scale*Fosq/FosqMax
388                B = scale*Fcsq/FosqMax
389                C = abs(A-B)
390            elif Type == 'Fo':
391                A = scale*math.sqrt(max(0,Fosq))/FoMax
392                B = scale*math.sqrt(max(0,Fcsq))/FoMax
393                C = abs(A-B)
394            elif Type == 'Unit Fc':
395                A = scale/2
396                B = scale/2
397                C = 0.0
398                if Fcsq and Fosq > 0:
399                    A *= min(1.0,Fosq/Fcsq)
400                    C = abs(A-B)
401            elif Type == '|DFsq|/sig':
402                if sig > 0.:
403                    A = (Fosq-Fcsq)/(3*sig)
404                B = 0
405            elif Type == '|DFsq|>sig':
406                if sig > 0.:
407                    A = (Fosq-Fcsq)/(3*sig)
408                if abs(A) < 1.0: A = 0
409                B = 0                   
410            elif Type == '|DFsq|>3sig':
411                if sig > 0.:
412                    A = (Fosq-Fcsq)/(3*sig)
413                if abs(A) < 3.0: A = 0
414                B = 0 
415            h = H+np.sum(SuperVec*refl[3:3+Super],axis=0)                   
416            xy = (h[pzone[izone][0]],h[pzone[izone][1]])
417            if Type in ['|DFsq|/sig','|DFsq|>sig','|DFsq|>3sig']:
418                if A > 0.0:
419                    Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w'))
420                else:
421                    Plot.add_artist(Circle(xy,radius=-A,ec='r',fc='w'))
422            else:
423                if A > 0.0 and A > B:
424                    Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w'))
425                if B:
426                    Plot.add_artist(Circle(xy,radius=B,ec='b',fc='w'))
427                    if A < B:
428                        Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w'))
429                    radius = C
430                    if radius > 0:
431                        if A > B:
432                            Plot.add_artist(Circle(xy,radius=radius,ec='g',fc='g'))
433                        else:                   
434                            Plot.add_artist(Circle(xy,radius=radius,ec='r',fc='r'))
435#    print 'plot time: %.3f'%(time.time()-time0)
436    HKL = np.array(HKL)
437    HKLF = np.array(HKLF)
438    Plot.set_xlabel(xlabel[izone]+str(Data['Layer']),fontsize=12)
439    Plot.set_ylabel(ylabel[izone],fontsize=12)
440    if not newPlot:
441        Page.toolbar.push_current()
442        Plot.set_xlim(xylim[0])
443        Plot.set_ylim(xylim[1])
444#        xylim = []
445        Page.toolbar.push_current()
446        Page.toolbar.draw()
447    else:
448        Plot.set_xlim((HKLmin[pzone[izone][0]],HKLmax[pzone[izone][0]]))
449        Plot.set_ylim((HKLmin[pzone[izone][1]],HKLmax[pzone[izone][1]]))
450        Page.canvas.draw()
451       
452################################################################################
453##### Plot3DSngl
454################################################################################
455
456def Plot3DSngl(G2frame,newPlot=False,Data=None,hklRef=None,Title=''):
457    '''3D Structure factor plotting package - displays reflections as rings proportional
458        to F, F**2, etc. as requested as 3D array
459    '''
460
461    def OnKeyBox(event):
462        mode = cb.GetValue()
463        if mode in ['jpeg','bmp','tiff',]:
464            try:
465                import Image as Im
466            except ImportError:
467                try:
468                    from PIL import Image as Im
469                except ImportError:
470                    print "PIL/pillow Image module not present. Cannot save images without this"
471                    raise Exception("PIL/pillow Image module not found")
472            Fname = os.path.join(Mydir,generalData['Name']+'.'+mode)
473            size = Page.canvas.GetSize()
474            glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
475            if mode in ['jpeg',]:
476                Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
477                im = Image.new("RGBA", (size[0],size[1]))
478            else:
479                Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
480                im = Image.new("RGB", (size[0],size[1]))
481            im.fromstring(Pix)
482            im.save(Fname,mode)
483            cb.SetValue(' save as/key:')
484            G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
485        else:
486            event.key = cb.GetValue()[0]
487            cb.SetValue(' save as/key:')
488            wx.CallAfter(OnKey,event)
489        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
490       
491    def OnKey(event):           #on key UP!!
492        Choice = {'F':'Fo','S':'Fosq','U':'Unit','D':'dFsq','W':'dFsq/sig'}
493        try:
494            keyCode = event.GetKeyCode()
495            if keyCode > 255:
496                keyCode = 0
497            key = chr(keyCode)
498        except AttributeError:       #if from OnKeyBox above
499            key = str(event.key).upper()
500        if key in ['C']:
501            drawingData['viewPoint'][0] = drawingData['default']
502            drawingData['viewDir'] = [0,0,1]
503            drawingData['oldxy'] = []
504            V0 = np.array([0,0,1])
505            V = np.inner(Amat,V0)
506            V /= np.sqrt(np.sum(V**2))
507            A = np.arccos(np.sum(V*V0))
508            Q = G2mth.AV2Q(A,[0,1,0])
509            drawingData['Quaternion'] = Q
510            Q = drawingData['Quaternion']
511        elif key == '+':
512            Data['Scale'] *= 1.25
513        elif key == '-':
514            Data['Scale'] /= 1.25
515        elif key == '0':
516            Data['Scale'] = 1.0
517        elif key == 'I':
518            Data['Iscale'] = not Data['Iscale']
519        elif key in Choice:
520            Data['Type'] = Choice[key]
521        Draw('key')
522           
523    Name = G2frame.PatternTree.GetItemText(G2frame.PatternId)
524    if Title != 'no phase':
525        generalData = G2frame.GetPhaseData()[Title]['General']
526        cell = generalData['Cell'][1:7]
527    else:
528        cell = [10,10,10,90,90,90]
529    drawingData = Data['Drawing']
530    Super = Data['Super']
531    SuperVec = []
532    for i in range(Super):
533        SuperVec.append(Data['SuperVec'][i][0])
534        SuperVec = np.array(SuperVec)
535    defaultViewPt = copy.copy(drawingData['viewPoint'])
536    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
537    Gmat,gmat = G2lat.cell2Gmat(cell)
538    invcell = G2lat.Gmat2cell(Gmat)
539    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
540    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
541    drawingData['Quaternion'] = G2mth.AV2Q(2*np.pi,np.inner(Bmat,[0,0,1]))
542    Wt = np.array([255,255,255])
543    Rd = np.array([255,0,0])
544    Gr = np.array([0,255,0])
545    wxGreen = wx.Colour(0,255,0)
546    Bl = np.array([0,0,255])
547    Or = np.array([255,128,0])
548    wxOrange = wx.Colour(255,128,0)
549    uBox = np.array([[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]])
550    uEdges = np.array([
551        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]]])
552    uColors = [Rd,Gr,Bl]
553   
554    def FillHKLRC():
555        R = np.zeros(len(hklRef))
556        C = []
557        HKL = []
558        RC = []
559        for i,refl in enumerate(hklRef):
560            H = refl[:3]
561            if 'HKLF' in Name:
562                Fosq,sig,Fcsq = refl[5+Super:8+Super]
563            else:
564                Fosq,sig,Fcsq = refl[8+Super],1.0,refl[9+Super]
565            HKL.append(H+np.sum(SuperVec*refl[3:3+Super],axis=0))
566            if Data['Type'] == 'Unit':
567                R[i] = 0.1
568                C.append(Gr)
569            elif Data['Type'] == 'Fosq':
570                if Fosq > 0:
571                    R[i] = Fosq
572                    C.append(Gr)
573                else:
574                    R[i] = -Fosq
575                    C.append(Rd)
576            elif Data['Type'] == 'Fo':
577                if Fosq > 0:
578                    R[i] = np.sqrt(Fosq)
579                    C.append(Gr)
580                else:
581                    R[i] = np.sqrt(-Fosq)
582                    C.append(Rd)
583            elif Data['Type'] == 'dFsq/sig':
584                dFsig = (Fosq-Fcsq)/sig
585                if dFsig > 0:
586                    R[i] = dFsig
587                    C.append(Gr)
588                else:
589                    R[i] = -dFsig
590                    C.append(Rd)
591            elif Data['Type'] == 'dFsq':
592                dF = Fosq-Fcsq
593                if dF > 0:
594                    R[i] = dF
595                    C.append(Gr)
596                else:
597                    R[i] = -dF
598                    C.append(Rd)
599        R /= np.max(R)
600        R *= Data['Scale']
601        R = np.where(R<1.e-5,1.e-5,R)
602        if Data['Iscale']:
603            R = np.where(R<=1.,R,1.)
604            C = np.array(C)
605            C = (C.T*R).T
606            R = np.ones_like(R)*0.05     
607        return HKL,zip(list(R),C)
608
609    def SetTranslation(newxy):
610#first get translation vector in screen coords.       
611        oldxy = drawingData['oldxy']
612        if not len(oldxy): oldxy = list(newxy)
613        dxy = newxy-oldxy
614        drawingData['oldxy'] = list(newxy)
615        V = np.array([-dxy[0],dxy[1],0.])
616#then transform to rotated crystal coordinates & apply to view point       
617        Q = drawingData['Quaternion']
618        V = np.inner(Bmat,G2mth.prodQVQ(G2mth.invQ(Q),V))
619        Tx,Ty,Tz = drawingData['viewPoint'][0]
620        Tx += V[0]*0.1
621        Ty += V[1]*0.1
622        Tz += V[2]*0.1
623        drawingData['viewPoint'][0] =  Tx,Ty,Tz
624       
625    def SetRotation(newxy):
626        'Perform a rotation in x-y space due to a left-mouse drag'
627    #first get rotation vector in screen coords. & angle increment       
628        oldxy = drawingData['oldxy']
629        if not len(oldxy): oldxy = list(newxy)
630        dxy = newxy-oldxy
631        drawingData['oldxy'] = list(newxy)
632        V = np.array([dxy[1],dxy[0],0.])
633        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
634        if not A: return # nothing changed, nothing to do
635    # next transform vector back to xtal coordinates via inverse quaternion
636    # & make new quaternion
637        Q = drawingData['Quaternion']
638        V = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,V))
639        DQ = G2mth.AVdeg2Q(A,V)
640        Q = G2mth.prodQQ(Q,DQ)
641        drawingData['Quaternion'] = Q
642    # finally get new view vector - last row of rotation matrix
643        VD = np.inner(Bmat,G2mth.Q2Mat(Q)[2])
644        VD /= np.sqrt(np.sum(VD**2))
645        drawingData['viewDir'] = VD
646       
647    def SetRotationZ(newxy):                       
648#first get rotation vector (= view vector) in screen coords. & angle increment       
649        View = glGetIntegerv(GL_VIEWPORT)
650        cent = [View[2]/2,View[3]/2]
651        oldxy = drawingData['oldxy']
652        if not len(oldxy): oldxy = list(newxy)
653        dxy = newxy-oldxy
654        drawingData['oldxy'] = list(newxy)
655        V = drawingData['viewDir']
656        A = [0,0]
657        A[0] = dxy[1]*.25
658        A[1] = dxy[0]*.25
659        if newxy[0] > cent[0]:
660            A[0] *= -1
661        if newxy[1] < cent[1]:
662            A[1] *= -1       
663# next transform vector back to xtal coordinates & make new quaternion
664        Q = drawingData['Quaternion']
665        V = np.inner(Amat,V)
666        Qx = G2mth.AVdeg2Q(A[0],V)
667        Qy = G2mth.AVdeg2Q(A[1],V)
668        Q = G2mth.prodQQ(Q,Qx)
669        Q = G2mth.prodQQ(Q,Qy)
670        drawingData['Quaternion'] = Q
671
672    def OnMouseDown(event):
673        xy = event.GetPosition()
674        drawingData['oldxy'] = list(xy)
675       
676    def OnMouseMove(event):
677        if event.ShiftDown():           #don't want any inadvertant moves when picking
678            return
679        newxy = event.GetPosition()
680                               
681        if event.Dragging():
682            if event.LeftIsDown():
683                SetRotation(newxy)
684                Q = drawingData['Quaternion']
685            elif event.RightIsDown():
686                SetTranslation(newxy)
687                Tx,Ty,Tz = drawingData['viewPoint'][0]
688            elif event.MiddleIsDown():
689                SetRotationZ(newxy)
690                Q = drawingData['Quaternion']
691            Draw('move')
692       
693    def OnMouseWheel(event):
694        if event.ShiftDown():
695            return
696        drawingData['cameraPos'] += event.GetWheelRotation()/120.
697        drawingData['cameraPos'] = max(0.1,min(20.00,drawingData['cameraPos']))
698        Draw('wheel')
699       
700    def SetBackground():
701        R,G,B,A = Page.camera['backColor']
702        glClearColor(R,G,B,A)
703        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
704       
705    def SetLights():
706        glEnable(GL_DEPTH_TEST)
707        glShadeModel(GL_SMOOTH)
708        glEnable(GL_LIGHTING)
709        glEnable(GL_LIGHT0)
710        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
711        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
712        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
713       
714    def RenderUnitVectors(x,y,z):
715        xyz = np.array([x,y,z])
716        glEnable(GL_COLOR_MATERIAL)
717        glLineWidth(1)
718        glPushMatrix()
719        glTranslate(x,y,z)
720        glBegin(GL_LINES)
721        for line,color in zip(uEdges,uColors):
722            glColor3ubv(color)
723            glVertex3fv(-line[1])
724            glVertex3fv(line[1])
725        glEnd()
726        glPopMatrix()
727        glColor4ubv([0,0,0,0])
728        glDisable(GL_COLOR_MATERIAL)
729               
730    def RenderDots(XYZ,RC):
731        glEnable(GL_COLOR_MATERIAL)
732        XYZ = np.array(XYZ)
733        glPushMatrix()
734        for xyz,rc in zip(XYZ,RC):
735            x,y,z = xyz
736            r,c = rc
737            glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,c)
738            glPointSize(r*50)
739            glBegin(GL_POINTS)
740            glVertex3fv(xyz)
741            glEnd()
742        glPopMatrix()
743        glColor4ubv([0,0,0,0])
744        glDisable(GL_COLOR_MATERIAL)
745       
746    def Draw(caller=''):
747#useful debug?       
748#        if caller:
749#            print caller
750# end of useful debug
751        G2frame.G2plotNB.status.SetStatusText('Plot type = %s for %s'%(Data['Type'],Name),1)
752        VS = np.array(Page.canvas.GetSize())
753        aspect = float(VS[0])/float(VS[1])
754        cPos = drawingData['cameraPos']
755        Zclip = drawingData['Zclip']*cPos/20.
756        Q = drawingData['Quaternion']
757        Tx,Ty,Tz = drawingData['viewPoint'][0]
758        G,g = G2lat.cell2Gmat(cell)
759        GS = G
760        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
761        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
762        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
763       
764        HKL,RC = FillHKLRC()
765       
766        SetBackground()
767        glInitNames()
768        glPushName(0)
769       
770        glMatrixMode(GL_PROJECTION)
771        glLoadIdentity()
772        glViewport(0,0,VS[0],VS[1])
773        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
774        gluLookAt(0,0,cPos,0,0,0,0,1,0)
775        SetLights()           
776           
777        glMatrixMode(GL_MODELVIEW)
778        glLoadIdentity()
779        matRot = G2mth.Q2Mat(Q)
780        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
781        glMultMatrixf(matRot.T)
782        glMultMatrixf(B4mat.T)
783        glTranslate(-Tx,-Ty,-Tz)
784        x,y,z = drawingData['viewPoint'][0]
785        RenderUnitVectors(x,y,z)
786        RenderUnitVectors(0,0,0)
787        RenderDots(HKL,RC)
788        time0 = time.time()
789        if Page.context: Page.canvas.SetCurrent(Page.context)    # wx 2.9 fix
790        Page.canvas.SwapBuffers()
791
792    # PlotStructure execution starts here (N.B. initialization above)
793    try:
794        plotNum = G2frame.G2plotNB.plotList.index('3D Structure Factors')
795        Page = G2frame.G2plotNB.nb.GetPage(plotNum)       
796    except ValueError:
797        Plot = G2frame.G2plotNB.addOgl('3D Structure Factors')
798        plotNum = G2frame.G2plotNB.plotList.index('3D Structure Factors')
799        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
800        Page.views = False
801        view = False
802        altDown = False
803    Font = Page.GetFont()
804    Page.SetFocus()
805    Page.Choice = None
806    choice = [' save as/key:','jpeg','tiff','bmp','c: recenter to default','+: increase scale',
807    '-: decrease scale','f: Fobs','s: Fobs**2','u: unit','d: Fo-Fc','w: DF/sig','i: toggle intensity scaling']
808    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice)
809    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
810    cb.SetValue(' save as/key:')
811    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
812    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
813    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
814    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
815    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
816    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
817#    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
818    Page.camera['position'] = drawingData['cameraPos']
819    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
820    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
821    try:
822        Page.canvas.SetCurrent()
823    except:
824        pass
825    Draw('main')
826#    if firstCall: Draw('main') # draw twice the first time that graphics are displayed
827
828       
829################################################################################
830##### PlotPatterns
831################################################################################
832           
833def PlotPatterns(G2frame,newPlot=False,plotType='PWDR'):
834    '''Powder pattern plotting package - displays single or multiple powder patterns as intensity vs
835    2-theta, q or TOF. Can display multiple patterns as "waterfall plots" or contour plots. Log I
836    plotting available.
837    '''
838    global HKL
839    global exclLines
840    global DifLine
841    global Ymax
842    plottype = plotType
843   
844    def OnPlotKeyPress(event):
845        try:        #one way to check if key stroke will work on plot
846            Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
847        except TypeError:
848            G2frame.G2plotNB.status.SetStatusText('Select '+plottype+' pattern first',1)
849            return
850        newPlot = False
851        if event.key == 'w':
852            G2frame.Weight = not G2frame.Weight
853            if not G2frame.Weight and 'PWDR' in plottype:
854                G2frame.SinglePlot = True
855            newPlot = True
856        elif event.key == 'e' and 'SASD' in plottype:
857            G2frame.ErrorBars = not G2frame.ErrorBars
858        elif event.key == 'b':
859            G2frame.SubBack = not G2frame.SubBack
860            if not G2frame.SubBack:
861                G2frame.SinglePlot = True               
862        elif event.key == 'n':
863            if G2frame.Contour:
864                pass
865            else:
866                G2frame.logPlot = not G2frame.logPlot
867                if not G2frame.logPlot:
868                    G2frame.Offset[0] = 0
869                newPlot = True
870        elif event.key == 's' and 'PWDR' in plottype:
871            if G2frame.Contour:
872                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
873                choice.sort()
874                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
875                if dlg.ShowModal() == wx.ID_OK:
876                    sel = dlg.GetSelection()
877                    G2frame.ContourColor = choice[sel]
878                else:
879                    G2frame.ContourColor = 'Paired'
880                dlg.Destroy()
881            elif G2frame.SinglePlot:
882                G2frame.SqrtPlot = not G2frame.SqrtPlot
883                if G2frame.SqrtPlot:
884                    G2frame.delOffset = .002
885                    G2frame.refOffset = -1.0
886                    G2frame.refDelt = .001
887                else:
888                    G2frame.delOffset = .02
889                    G2frame.refOffset = -1.0
890                    G2frame.refDelt = .01
891            newPlot = True
892        elif event.key == 'u' and (G2frame.Contour or not G2frame.SinglePlot):
893            if G2frame.Contour:
894                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
895            elif G2frame.Offset[0] < 100.:
896                G2frame.Offset[0] += 1.
897        elif event.key == 'd' and (G2frame.Contour or not G2frame.SinglePlot):
898            if G2frame.Contour:
899                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
900            elif G2frame.Offset[0] > 0.:
901                G2frame.Offset[0] -= 1.
902        elif event.key == 'l' and not G2frame.SinglePlot:
903            G2frame.Offset[1] -= 1.
904        elif event.key == 'r' and not G2frame.SinglePlot:
905            G2frame.Offset[1] += 1.
906        elif event.key == 'o' and not G2frame.SinglePlot:
907            G2frame.Cmax = 1.0
908            G2frame.Offset = [0,0]
909        elif event.key == 'c' and 'PWDR' in plottype:
910            newPlot = True
911            if not G2frame.Contour:
912                G2frame.SinglePlot = False
913                G2frame.Offset = [0.,0.]
914            else:
915                G2frame.SinglePlot = True               
916            G2frame.Contour = not G2frame.Contour
917        elif event.key == 'q': 
918            if 'PWDR' in plottype:
919                newPlot = True
920                G2frame.qPlot = not G2frame.qPlot
921                G2frame.dPlot = False
922            elif 'SASD' in plottype:
923                newPlot = True
924                G2frame.sqPlot = not G2frame.sqPlot
925        elif event.key == 't' and 'PWDR' in plottype:
926            G2frame.dPlot = not G2frame.dPlot
927            G2frame.qPlot = False
928            newPlot = True     
929        elif event.key == 'm':
930            G2frame.SqrtPlot = False
931            G2frame.SinglePlot = not G2frame.SinglePlot               
932            newPlot = True
933        elif event.key == '+':
934            if G2frame.PickId:
935                G2frame.PickId = False
936        elif event.key == 'i' and G2frame.Contour:                  #for smoothing contour plot
937            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
938               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
939               'mitchell','sinc','lanczos']
940            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
941            if dlg.ShowModal() == wx.ID_OK:
942                sel = dlg.GetSelection()
943                G2frame.Interpolate = choice[sel]
944            else:
945                G2frame.Interpolate = 'nearest'
946            dlg.Destroy()
947           
948        wx.CallAfter(PlotPatterns,G2frame,newPlot=newPlot,plotType=plottype)
949       
950    def OnMotion(event):
951        xpos = event.xdata
952        if xpos:                                        #avoid out of frame mouse position
953            ypos = event.ydata
954            Page.canvas.SetCursor(wx.CROSS_CURSOR)
955            try:
956                Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
957                if G2frame.qPlot and 'PWDR' in plottype:
958                    q = xpos
959                    dsp = 2.*np.pi/q
960                    try:
961                        xpos = G2lat.Dsp2pos(Parms,2.0*np.pi/xpos)
962                    except ValueError:      #avoid bad value in asin beyond upper limit
963                        pass
964                elif 'SASD' in plottype:
965                    q = xpos
966                    dsp = 2.*np.pi/q
967                elif G2frame.dPlot:
968                    dsp = xpos
969                    q = 2.*np.pi/dsp
970                    xpos = G2lat.Dsp2pos(Parms,xpos)
971                elif G2frame.Contour and 'T' in Parms['Type'][0]:
972                    xpos = X[xpos]                   
973                    dsp = G2lat.Pos2dsp(Parms,xpos)
974                    q = 2.*np.pi/dsp
975                else:
976                    dsp = G2lat.Pos2dsp(Parms,xpos)
977                    q = 2.*np.pi/dsp
978                if G2frame.Contour: #PWDR only
979                    if 'C' in Parms['Type'][0]:
980                        G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f q = %9.5f pattern ID =%5d'%(xpos,dsp,q,int(ypos)),1)
981                    else:
982                        G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q = %9.5f pattern ID =%5d'%(xpos,dsp,q,int(ypos)),1)
983                else:
984                    if 'C' in Parms['Type'][0]:
985                        if 'PWDR' in plottype:
986                            if G2frame.SqrtPlot:
987                                G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f q = %9.5f sqrt(Intensity) =%9.2f'%(xpos,dsp,q,ypos),1)
988                            else:
989                                G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f q = %9.5f Intensity =%9.2f'%(xpos,dsp,q,ypos),1)
990                        elif 'SASD' in plottype:
991                            G2frame.G2plotNB.status.SetStatusText('q =%12.5g Intensity =%12.5g d =%9.1f'%(q,ypos,dsp),1)
992                    else:
993                        if G2frame.SqrtPlot:
994                            G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q =%9.5f sqrt(Intensity) =%9.2f'%(xpos,dsp,q,ypos),1)
995                        else:
996                            G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f q =%9.5f Intensity =%9.2f'%(xpos,dsp,q,ypos),1)
997                if G2frame.itemPicked:
998                    Page.canvas.SetToolTipString('%9.5f'%(xpos))
999                if G2frame.PickId:
1000                    found = []
1001                    if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists'] or \
1002                        'PWDR' in G2frame.PatternTree.GetItemText(PickId):
1003                        if len(HKL):
1004                            view = Page.toolbar._views.forward()[0][:2]
1005                            wid = view[1]-view[0]
1006                            found = HKL[np.where(np.fabs(HKL.T[-1]-xpos) < 0.002*wid)]
1007                        if len(found):
1008                            if len(found[0]) > 6:   #SS reflections
1009                                h,k,l,m = found[0][:4]
1010                                Page.canvas.SetToolTipString('%d,%d,%d,%d'%(int(h),int(k),int(l),int(m)))
1011                            else:
1012                                h,k,l = found[0][:3] 
1013                                Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
1014                        else:
1015                            Page.canvas.SetToolTipString('')
1016
1017            except TypeError:
1018                G2frame.G2plotNB.status.SetStatusText('Select '+plottype+' pattern first',1)
1019               
1020    def OnPress(event): #ugh - this removes a matplotlib error for mouse clicks in log plots                 
1021        olderr = np.seterr(invalid='ignore')
1022                                                   
1023    def OnPick(event):
1024        if G2frame.itemPicked is not None: return
1025        PatternId = G2frame.PatternId
1026        try:
1027            Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
1028        except TypeError:
1029            return
1030        PickId = G2frame.PickId
1031        pick = event.artist
1032        mouse = event.mouseevent
1033        xpos = pick.get_xdata()
1034        ypos = pick.get_ydata()
1035        ind = event.ind
1036        xy = list(zip(np.take(xpos,ind),np.take(ypos,ind))[0])
1037        if G2frame.PatternTree.GetItemText(PickId) == 'Peak List':
1038            if ind.all() != [0] and ObsLine[0].get_label() in str(pick):                                    #picked a data point
1039                data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
1040                if G2frame.qPlot:                              #qplot - convert back to 2-theta
1041                    xy[0] = G2lat.Dsp2pos(Parms,2*np.pi/xy[0])
1042                elif G2frame.dPlot:                            #dplot - convert back to 2-theta
1043                    xy[0] = G2lat.Dsp2pos(Parms,xy[0])
1044                XY = G2mth.setPeakparms(Parms,Parms2,xy[0],xy[1])
1045                data['peaks'].append(XY)
1046                data['sigDict'] = {}    #now invalid
1047                G2pdG.UpdatePeakGrid(G2frame,data)
1048                PlotPatterns(G2frame,plotType=plottype)
1049            else:                                                   #picked a peak list line
1050                G2frame.itemPicked = pick
1051        elif G2frame.PatternTree.GetItemText(PickId) == 'Limits':
1052            if ind.all() != [0]:                                    #picked a data point
1053                LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
1054                data = G2frame.PatternTree.GetItemPyData(LimitId)
1055                if G2frame.qPlot:                              #qplot - convert back to 2-theta
1056                    xy[0] = G2lat.Dsp2pos(Parms,2*np.pi/xy[0])
1057                elif G2frame.dPlot:                            #dplot - convert back to 2-theta
1058                    xy[0] = G2lat.Dsp2pos(Parms,xy[0])
1059                if G2frame.ifGetExclude:
1060                    excl = [0,0]
1061                    excl[0] = max(data[1][0],min(xy[0],data[1][1]))
1062                    excl[1] = excl[0]+0.1
1063                    data.append(excl)
1064                    G2frame.ifGetExclude = False
1065                else:
1066                    if mouse.button==1:
1067                        data[1][0] = min(xy[0],data[1][1])
1068                    if mouse.button==3:
1069                        data[1][1] = max(xy[0],data[1][0])
1070                G2frame.PatternTree.SetItemPyData(LimitId,data)
1071                G2pdG.UpdateLimitsGrid(G2frame,data,plottype)
1072                wx.CallAfter(PlotPatterns,G2frame,plotType=plottype)
1073            else:                                                   #picked a limit line
1074                G2frame.itemPicked = pick
1075        elif G2frame.PatternTree.GetItemText(PickId) == 'Models':
1076            if ind.all() != [0]:                                    #picked a data point
1077                LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
1078                data = G2frame.PatternTree.GetItemPyData(LimitId)
1079                if mouse.button==1:
1080                    data[1][0] = min(xy[0],data[1][1])
1081                if mouse.button==3:
1082                    data[1][1] = max(xy[0],data[1][0])
1083                G2frame.PatternTree.SetItemPyData(LimitId,data)
1084                wx.CallAfter(PlotPatterns,G2frame,plotType=plottype)
1085            else:                                                   #picked a limit line
1086                G2frame.itemPicked = pick
1087        elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \
1088            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
1089            G2frame.itemPicked = pick
1090            pick = str(pick)
1091       
1092    def OnRelease(event):
1093        if G2frame.itemPicked is None: return
1094        if str(DifLine[0]) == str(G2frame.itemPicked):
1095            ypos = event.ydata
1096            G2frame.delOffset = -ypos/Ymax
1097            G2frame.itemPicked = None
1098            PlotPatterns(G2frame,plotType=plottype)
1099            return
1100        Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
1101        xpos = event.xdata
1102        PickId = G2frame.PickId
1103        if G2frame.PatternTree.GetItemText(PickId) in ['Peak List','Limits'] and xpos:
1104            lines = []
1105            for line in G2frame.Lines: 
1106                lines.append(line.get_xdata()[0])
1107            try:
1108                lineNo = lines.index(G2frame.itemPicked.get_xdata()[0])
1109            except ValueError:
1110                lineNo = -1
1111            if  lineNo in [0,1] or lineNo in exclLines:
1112                LimitId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Limits')
1113                data = G2frame.PatternTree.GetItemPyData(LimitId)
1114                id = lineNo/2+1
1115                id2 = lineNo%2
1116                if G2frame.qPlot and 'PWDR' in plottype:
1117                    data[id][id2] = G2lat.Dsp2pos(Parms,2.*np.pi/xpos)
1118                elif G2frame.dPlot and 'PWDR' in plottype:
1119                    data[id][id2] = G2lat.Dsp2pos(Parms,xpos)
1120                else:
1121                    data[id][id2] = xpos
1122                if id > 1 and data[id][0] > data[id][1]:
1123                        data[id].reverse()
1124                data[1][0] = min(max(data[0][0],data[1][0]),data[1][1])
1125                data[1][1] = max(min(data[0][1],data[1][1]),data[1][0])
1126                G2frame.PatternTree.SetItemPyData(LimitId,data)
1127                if G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Limits':
1128                    G2pdG.UpdateLimitsGrid(G2frame,data,plottype)
1129            elif lineNo > 1:
1130                PeakId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Peak List')
1131                data = G2frame.PatternTree.GetItemPyData(PeakId)
1132                if event.button == 3:
1133                    del data['peaks'][lineNo-2]
1134                else:
1135                    if G2frame.qPlot:
1136                        data['peaks'][lineNo-2][0] = G2lat.Dsp2pos(Parms,2.*np.pi/xpos)
1137                    elif G2frame.dPlot:
1138                        data['peaks'][lineNo-2][0] = G2lat.Dsp2pos(Parms,xpos)
1139                    else:
1140                        data['peaks'][lineNo-2][0] = xpos
1141                    data['sigDict'] = {}        #no longer valid
1142                G2frame.PatternTree.SetItemPyData(PeakId,data)
1143                G2pdG.UpdatePeakGrid(G2frame,data)
1144        elif G2frame.PatternTree.GetItemText(PickId) in ['Models',] and xpos:
1145            lines = []
1146            for line in G2frame.Lines: 
1147                lines.append(line.get_xdata()[0])
1148            try:
1149                lineNo = lines.index(G2frame.itemPicked.get_xdata()[0])
1150            except ValueError:
1151                lineNo = -1
1152            if  lineNo in [0,1]:
1153                LimitId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Limits')
1154                data = G2frame.PatternTree.GetItemPyData(LimitId)
1155                data[1][lineNo] = xpos
1156                data[1][0] = min(max(data[0][0],data[1][0]),data[1][1])
1157                data[1][1] = max(min(data[0][1],data[1][1]),data[1][0])
1158                G2frame.PatternTree.SetItemPyData(LimitId,data)       
1159        elif (G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \
1160            'PWDR' in G2frame.PatternTree.GetItemText(PickId)) and xpos:
1161            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
1162            pick = str(G2frame.itemPicked).split('(')[1].strip(')')
1163            if 'line' not in pick:       #avoid data points, etc.
1164                num = Phases.keys().index(pick)
1165                if num:
1166                    G2frame.refDelt = -(event.ydata-G2frame.refOffset)/(num*Ymax)
1167                else:       #1st row of refl ticks
1168                    G2frame.refOffset = event.ydata
1169        PlotPatterns(G2frame,plotType=plottype)
1170        G2frame.itemPicked = None   
1171
1172    try:
1173        plotNum = G2frame.G2plotNB.plotList.index('Powder Patterns')
1174        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1175        if not newPlot:
1176            Plot = Page.figure.gca()          #get previous powder plot & get limits
1177            xylim = Plot.get_xlim(),Plot.get_ylim()
1178        Page.figure.clf()
1179        Plot = Page.figure.gca()          #get a fresh plot after clf()
1180    except ValueError:
1181        if plottype == 'SASD':
1182            G2frame.logPlot = True
1183            G2frame.ErrorBars = True
1184        newPlot = True
1185        G2frame.Cmax = 1.0
1186        Plot = G2frame.G2plotNB.addMpl('Powder Patterns').gca()
1187        plotNum = G2frame.G2plotNB.plotList.index('Powder Patterns')
1188        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1189        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1190        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1191        Page.canvas.mpl_connect('pick_event', OnPick)
1192        Page.canvas.mpl_connect('button_release_event', OnRelease)
1193        Page.canvas.mpl_connect('button_press_event',OnPress)
1194    if plottype == 'PWDR':  # avoids a very nasty clash with KILL_FOCUS in SASD TextCtrl?
1195        Page.SetFocus()
1196    G2frame.G2plotNB.status.DestroyChildren()
1197    if G2frame.Contour:
1198        Page.Choice = (' key press','d: lower contour max','u: raise contour max','o: reset contour max',
1199            'i: interpolation method','s: color scheme','c: contour off')
1200    else:
1201        if G2frame.logPlot:
1202            if 'PWDR' in plottype:
1203                if G2frame.SinglePlot:
1204                    Page.Choice = (' key press','n: log(I) off',
1205                        'c: contour on','q: toggle q plot','t: toggle d-spacing plot',
1206                            'm: toggle multidata plot','w: toggle divide by sig','+: no selection')
1207                else:
1208                    Page.Choice = (' key press','n: log(I) off',
1209                        'd: offset down','l: offset left','r: offset right','u: offset up','o: reset offset',
1210                        'c: contour on','q: toggle q plot','t: toggle d-spacing plot',
1211                        'm: toggle multidata plot','w: toggle divide by sig','+: no selection')
1212            elif 'SASD' in plottype:
1213                if G2frame.SinglePlot:
1214                    Page.Choice = (' key press','b: toggle subtract background file','n: semilog on',
1215                        'q: toggle S(q) plot','m: toggle multidata plot','w: toggle (Io-Ic)/sig plot','+: no selection')
1216                else:
1217                    Page.Choice = (' key press','b: toggle subtract background file','n: semilog on',
1218                        'd: offset down','l: offset left','r: offset right','u: offset up','o: reset offset',
1219                        'q: toggle S(q) plot','m: toggle multidata plot','w: toggle (Io-Ic)/sig plot','+: no selection')
1220        else:
1221            if 'PWDR' in plottype:
1222                if G2frame.SinglePlot:
1223                    Page.Choice = (' key press',
1224                        'b: toggle subtract background','n: log(I) on','s: toggle sqrt plot','c: contour on',
1225                        'q: toggle q plot','t: toggle d-spacing plot','m: toggle multidata plot',
1226                        'w: toggle divide by sig','+: no selection')
1227                else:
1228                    Page.Choice = (' key press','l: offset left','r: offset right','d: offset down',
1229                        'u: offset up','o: reset offset','b: toggle subtract background','n: log(I) on','c: contour on',
1230                        'q: toggle q plot','t: toggle d-spacing plot','m: toggle multidata plot',
1231                        'w: toggle divide by sig','+: no selection')
1232            elif 'SASD' in plottype:
1233                if G2frame.SinglePlot:
1234                    Page.Choice = (' key press','b: toggle subtract background file','n: loglog on','e: toggle error bars',
1235                        'q: toggle S(q) plot','m: toggle multidata plot','w: toggle (Io-Ic)/sig plot','+: no selection')
1236                else:
1237                    Page.Choice = (' key press','b: toggle subtract background file','n: loglog on','e: toggle error bars',
1238                        'd: offset down','l: offset left','r: offset right','u: offset up','o: reset offset',
1239                        'q: toggle S(q) plot','m: toggle multidata plot','w: toggle (Io-Ic)/sig plot','+: no selection')
1240
1241    Page.keyPress = OnPlotKeyPress   
1242    PickId = G2frame.PickId
1243    PatternId = G2frame.PatternId
1244    colors=['b','g','r','c','m','k']
1245    Lines = []
1246    exclLines = []
1247    if G2frame.SinglePlot:
1248        Pattern = G2frame.PatternTree.GetItemPyData(PatternId)
1249        Pattern.append(G2frame.PatternTree.GetItemText(PatternId))
1250        PlotList = [Pattern,]
1251        Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
1252            G2frame.PatternId, 'Instrument Parameters'))
1253        Sample = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Sample Parameters'))
1254        ParmList = [Parms,]
1255        SampleList = [Sample,]
1256        Title = Pattern[-1]
1257    else:       
1258        Title = os.path.split(G2frame.GSASprojectfile)[1]
1259        PlotList = []
1260        ParmList = []
1261        SampleList = []
1262        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
1263        while item:
1264            if plottype in G2frame.PatternTree.GetItemText(item):
1265                Pattern = G2frame.PatternTree.GetItemPyData(item)
1266                if len(Pattern) < 3:                    # put name on end if needed
1267                    Pattern.append(G2frame.PatternTree.GetItemText(item))
1268                PlotList.append(Pattern)
1269                ParmList.append(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
1270                    item,'Instrument Parameters'))[0])
1271                SampleList.append(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
1272                    item, 'Sample Parameters')))
1273            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)               
1274    lenX = 0
1275    if PickId:
1276        if G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
1277            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
1278            HKL = []
1279            if Phases:
1280                try:
1281                    for peak in Phases[G2frame.RefList]['RefList']:
1282                        HKL.append(peak[:6])
1283                except TypeError:
1284                    for peak in Phases[G2frame.RefList]:
1285                        HKL.append(peak[:6])                   
1286                HKL = np.array(HKL)
1287        else:
1288            HKL = np.array(G2frame.HKL)
1289    Ymax = None
1290    for Pattern in PlotList:
1291        xye = Pattern[1]
1292        if xye[1] is None: continue
1293        if Ymax is None: Ymax = max(xye[1])
1294        Ymax = max(Ymax,max(xye[1]))
1295    if Ymax is None: return # nothing to plot
1296    offset = G2frame.Offset[0]*Ymax/100.0
1297    if G2frame.logPlot:
1298        Title = 'log('+Title+')'
1299    Plot.set_title(Title)
1300    if G2frame.qPlot or 'SASD' in plottype and not G2frame.Contour:
1301        Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=16)
1302    elif G2frame.dPlot and 'PWDR' in plottype and not G2frame.Contour:
1303        Plot.set_xlabel(r'$d, \AA$',fontsize=16)
1304    else:
1305        if 'C' in ParmList[0]['Type'][0]:       
1306            Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=16)
1307        else:
1308            if G2frame.Contour:
1309                Plot.set_xlabel(r'Channel no.',fontsize=16)           
1310            else:
1311                Plot.set_xlabel(r'$TOF, \mathsf{\mu}$s',fontsize=16)           
1312    if G2frame.Weight:
1313        if 'PWDR' in plottype:
1314            Plot.set_ylabel(r'$\mathsf{I/\sigma(I)}$',fontsize=16)
1315        elif 'SASD' in plottype:
1316            Plot.set_ylabel(r'$\mathsf{\Delta(I)/\sigma(I)}$',fontsize=16)
1317    else:
1318        if 'C' in ParmList[0]['Type'][0]:
1319            if 'PWDR' in plottype:
1320                if G2frame.SqrtPlot:
1321                    Plot.set_ylabel(r'$\sqrt{Intensity}$',fontsize=16)
1322                else:
1323                    Plot.set_ylabel(r'$Intensity$',fontsize=16)
1324            elif 'SASD' in plottype:
1325                if G2frame.sqPlot:
1326                    Plot.set_ylabel(r'$S(Q)=I*Q^{4}$',fontsize=16)
1327                else:
1328                    Plot.set_ylabel(r'$Intensity, cm^{-1}$',fontsize=16)
1329        else:
1330            if G2frame.SqrtPlot:
1331                Plot.set_ylabel(r'$\sqrt{Normalized\ intensity}$',fontsize=16)
1332            else:
1333                Plot.set_ylabel(r'$Normalized\ intensity$',fontsize=16)
1334    if G2frame.Contour:
1335        ContourZ = []
1336        ContourY = []
1337        Nseq = 0
1338    for N,Pattern in enumerate(PlotList):
1339        Parms = ParmList[N]
1340        Sample = SampleList[N]
1341        if 'C' in Parms['Type'][0]:
1342            wave = G2mth.getWave(Parms)
1343        else:
1344            difC = Parms['difC'][1]
1345        ifpicked = False
1346        LimitId = 0
1347        if Pattern[1] is None: continue # skip over uncomputed simulations
1348        xye = ma.array(ma.getdata(Pattern[1]))
1349        Zero = Parms.get('Zero',[0.,0.])[1]
1350        if PickId:
1351            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
1352            LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
1353            limits = np.array(G2frame.PatternTree.GetItemPyData(LimitId))
1354            excls = limits[2:]
1355            for excl in excls:
1356                xye[0] = ma.masked_inside(xye[0],excl[0],excl[1])
1357        if G2frame.qPlot and 'PWDR' in plottype:
1358            Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, Pattern[2])
1359            X = 2.*np.pi/G2lat.Pos2dsp(Parms,xye[0])
1360        elif G2frame.dPlot and 'PWDR' in plottype:
1361            Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, Pattern[2])
1362            X = G2lat.Pos2dsp(Parms,xye[0])
1363        else:
1364            X = xye[0]
1365        if not lenX:
1366            lenX = len(X)
1367        if 'PWDR' in plottype:
1368            if G2frame.SqrtPlot:
1369                olderr = np.seterr(invalid='ignore') #get around sqrt(-ve) error
1370                Y = np.where(xye[1]>=0.,np.sqrt(xye[1]),-np.sqrt(-xye[1]))
1371                np.seterr(invalid=olderr['invalid'])
1372            else:
1373                Y = xye[1]+offset*N
1374        elif 'SASD' in plottype:
1375            B = xye[5]
1376            if G2frame.sqPlot:
1377                Y = xye[1]*Sample['Scale'][0]*(1.005)**(offset*N)*X**4
1378            else:
1379                Y = xye[1]*Sample['Scale'][0]*(1.005)**(offset*N)
1380        if LimitId and ifpicked:
1381            limits = np.array(G2frame.PatternTree.GetItemPyData(LimitId))
1382            lims = limits[1]
1383            if G2frame.qPlot and 'PWDR' in plottype:
1384                lims = 2.*np.pi/G2lat.Pos2dsp(Parms,lims)
1385            elif G2frame.dPlot and 'PWDR' in plottype:
1386                lims = G2lat.Pos2dsp(Parms,lims)
1387            Lines.append(Plot.axvline(lims[0],color='g',dashes=(5,5),picker=3.))   
1388            Lines.append(Plot.axvline(lims[1],color='r',dashes=(5,5),picker=3.))
1389            for i,item in enumerate(limits[2:]):
1390                Lines.append(Plot.axvline(item[0],color='m',dashes=(5,5),picker=3.))   
1391                Lines.append(Plot.axvline(item[1],color='m',dashes=(5,5),picker=3.))
1392                exclLines += [2*i+2,2*i+3]
1393        if G2frame.Contour:           
1394            if lenX == len(X):
1395                ContourY.append(N)
1396                ContourZ.append(Y)
1397                if 'C' in ParmList[0]['Type'][0]:       
1398                    ContourX = X
1399                else: #'T'OF
1400                    ContourX = range(lenX)
1401                Nseq += 1
1402                Plot.set_ylabel('Data sequence',fontsize=12)
1403        else:
1404            if 'SASD' in plottype and G2frame.logPlot:
1405                X *= (1.01)**(G2frame.Offset[1]*N)
1406            else:
1407                X += G2frame.Offset[1]*.005*N
1408            Xum = ma.getdata(X)
1409            DifLine = ['']
1410            if ifpicked:
1411                if G2frame.SqrtPlot:
1412                    Z = np.where(xye[3]>=0.,np.sqrt(xye[3]),-np.sqrt(-xye[3]))
1413                else:
1414                    Z = xye[3]+offset*N
1415                if 'PWDR' in plottype:
1416                    if G2frame.SqrtPlot:
1417                        W = np.where(xye[4]>=0.,np.sqrt(xye[4]),-np.sqrt(-xye[4]))
1418                        D = np.where(xye[5],(Y-Z),0.)-Ymax*G2frame.delOffset
1419                    else:
1420                        W = xye[4]+offset*N
1421                        D = xye[5]-Ymax*G2frame.delOffset  #powder background
1422                elif 'SASD' in plottype:
1423                    if G2frame.sqPlot:
1424                        W = xye[4]*X**4
1425                        Z = xye[3]*X**4
1426                        B = B*X**4
1427                    else:
1428                        W = xye[4]
1429                    if G2frame.SubBack:
1430                        YB = Y-B
1431                        ZB = Z
1432                    else:
1433                        YB = Y
1434                        ZB = Z+B
1435                    Plot.set_yscale("log",nonposy='mask')
1436                    if np.any(W>0.):
1437                        Plot.set_ylim(bottom=np.min(np.trim_zeros(W))/2.,top=np.max(Y)*2.)
1438                    else:
1439                        Plot.set_ylim(bottom=np.min(np.trim_zeros(YB))/2.,top=np.max(Y)*2.)
1440                if G2frame.logPlot:
1441                    if 'PWDR' in plottype:
1442                        Plot.set_yscale("log",nonposy='mask')
1443                        Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
1444                        Plot.plot(X,Z,colors[(N+1)%6],picker=False)
1445                        Plot.plot(X,W,colors[(N+2)%6],picker=False)     #background
1446                    elif 'SASD' in plottype:
1447                        Plot.set_xscale("log",nonposx='mask')
1448                        Ibeg = np.searchsorted(X,limits[1][0])
1449                        Ifin = np.searchsorted(X,limits[1][1])
1450                        if G2frame.Weight:
1451                            Plot.set_yscale("linear")
1452                            DS = (YB-ZB)*np.sqrt(xye[2])
1453                            Plot.plot(X[Ibeg:Ifin],DS[Ibeg:Ifin],colors[(N+3)%6],picker=False)
1454                            Plot.axhline(0.,color=wx.BLACK)
1455                            Plot.set_ylim(bottom=np.min(DS[Ibeg:Ifin])*1.2,top=np.max(DS[Ibeg:Ifin])*1.2)                                                   
1456                        else:
1457                            Plot.set_yscale("log",nonposy='mask')
1458                            if G2frame.ErrorBars:
1459                                if G2frame.sqPlot:
1460                                    Plot.errorbar(X,YB,yerr=X**4*Sample['Scale'][0]*np.sqrt(1./(Pattern[0]['wtFactor']*xye[2])),
1461                                        ecolor=colors[N%6],picker=3.,clip_on=False)
1462                                else:
1463                                    Plot.errorbar(X,YB,yerr=Sample['Scale'][0]*np.sqrt(1./(Pattern[0]['wtFactor']*xye[2])),
1464                                        ecolor=colors[N%6],picker=3.,clip_on=False)
1465                            else:
1466                                Plot.plot(X,YB,colors[N%6]+'+',picker=3.,clip_on=False)
1467                            Plot.plot(X,W,colors[(N+2)%6],picker=False)     #const. background
1468                            Plot.plot(X,ZB,colors[(N+1)%6],picker=False)
1469                elif G2frame.Weight and 'PWDR' in plottype:
1470                    DY = xye[1]*np.sqrt(xye[2])
1471                    Ymax = max(DY)
1472                    DZ = xye[3]*np.sqrt(xye[2])
1473                    DS = xye[5]*np.sqrt(xye[2])-Ymax*G2frame.delOffset
1474                    ObsLine = Plot.plot(X,DY,colors[N%6]+'+',picker=3.,clip_on=False)         #Io/sig(Io)
1475                    Plot.plot(X,DZ,colors[(N+1)%6],picker=False)                    #Ic/sig(Io)
1476                    DifLine = Plot.plot(X,DS,colors[(N+3)%6],picker=1.)                    #(Io-Ic)/sig(Io)
1477                    Plot.axhline(0.,color=wx.BLACK)
1478                else:
1479                    if G2frame.SubBack:
1480                        if 'PWDR' in plottype:
1481                            Plot.plot(Xum,Y-W,colors[N%6]+'+',picker=False,clip_on=False)  #Io-Ib
1482                            Plot.plot(X,Z-W,colors[(N+1)%6],picker=False)               #Ic-Ib
1483                        else:
1484                            Plot.plot(X,YB,colors[N%6]+'+',picker=3.,clip_on=False)
1485                            Plot.plot(X,ZB,colors[(N+1)%6],picker=False)
1486                    else:
1487                        if 'PWDR' in plottype:
1488                            ObsLine = Plot.plot(Xum,Y,colors[N%6]+'+',picker=3.,clip_on=False)    #Io
1489                            Plot.plot(X,Z,colors[(N+1)%6],picker=False)                 #Ic
1490                        else:
1491                            Plot.plot(X,YB,colors[N%6]+'+',picker=3.,clip_on=False)
1492                            Plot.plot(X,ZB,colors[(N+1)%6],picker=False)
1493                    if 'PWDR' in plottype:
1494                        Plot.plot(X,W,colors[(N+2)%6],picker=False)                 #Ib
1495                        DifLine = Plot.plot(X,D,colors[(N+3)%6],picker=1.)                 #Io-Ic
1496                    Plot.axhline(0.,color=wx.BLACK)
1497                Page.canvas.SetToolTipString('')
1498                if G2frame.PatternTree.GetItemText(PickId) == 'Peak List':
1499                    tip = 'On data point: Pick peak - L or R MB. On line: L-move, R-delete'
1500                    Page.canvas.SetToolTipString(tip)
1501                    data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'))
1502                    for item in data['peaks']:
1503                        if G2frame.qPlot:
1504                            Lines.append(Plot.axvline(2.*np.pi/G2lat.Pos2dsp(Parms,item[0]),color=colors[N%6],picker=2.))
1505                        elif G2frame.dPlot:
1506                            Lines.append(Plot.axvline(G2lat.Pos2dsp(Parms,item[0]),color=colors[N%6],picker=2.))
1507                        else:
1508                            Lines.append(Plot.axvline(item[0],color=colors[N%6],picker=2.))
1509                if G2frame.PatternTree.GetItemText(PickId) == 'Limits':
1510                    tip = 'On data point: Lower limit - L MB; Upper limit - R MB. On limit: MB down to move'
1511                    Page.canvas.SetToolTipString(tip)
1512                    data = G2frame.LimitsTable.GetData()
1513                   
1514            else:   #not picked
1515                if G2frame.logPlot:
1516                    if 'PWDR' in plottype:
1517                        Plot.semilogy(X,Y,colors[N%6],picker=False,nonposy='mask')
1518                    elif 'SASD' in plottype:
1519                        Plot.semilogy(X,Y,colors[N%6],picker=False,nonposy='mask')
1520                else:
1521                    if 'PWDR' in plottype:
1522                        Plot.plot(X,Y,colors[N%6],picker=False)
1523                    elif 'SASD' in plottype:
1524                        Plot.loglog(X,Y,colors[N%6],picker=False,nonposy='mask')
1525                        Plot.set_ylim(bottom=np.min(np.trim_zeros(Y))/2.,top=np.max(Y)*2.)
1526                           
1527                if G2frame.logPlot and 'PWDR' in plottype:
1528                    Plot.set_ylim(bottom=np.min(np.trim_zeros(Y))/2.,top=np.max(Y)*2.)
1529    if PickId and not G2frame.Contour:
1530        Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
1531        if G2frame.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
1532            peaks = np.array((G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))))[0]
1533            for peak in peaks:
1534                if G2frame.qPlot:
1535                    Plot.axvline(2.*np.pi/G2lat.Pos2dsp(Parms,peak[0]),color='b')
1536                if G2frame.dPlot:
1537                    Plot.axvline(G2lat.Pos2dsp(Parms,peak[0]),color='b')
1538                else:
1539                    Plot.axvline(peak[0],color='b')
1540            for hkl in G2frame.HKL:
1541                clr = 'r'
1542                if len(hkl) > 6 and hkl[3]:
1543                    clr = 'g'
1544                if G2frame.qPlot:
1545                    Plot.axvline(2.*np.pi/G2lat.Pos2dsp(Parms,hkl[-1]),color=clr,dashes=(5,5))
1546                if G2frame.dPlot:
1547                    Plot.axvline(G2lat.Pos2dsp(Parms,hkl[-1]),color=clr,dashes=(5,5))
1548                else:
1549                    Plot.axvline(hkl[-1],color=clr,dashes=(5,5))
1550        elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists'] or \
1551            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
1552            refColors=['b','r','c','g','m','k']
1553            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
1554            for pId,phase in enumerate(Phases):
1555                try:    #patch for old style reflection lists
1556                    peaks = Phases[phase]['RefList']
1557                except TypeError:
1558                    peaks = Phases[phase]
1559                if not len(peaks):
1560                    continue
1561                peak = np.array([[peak[4],peak[5]] for peak in peaks])
1562                pos = G2frame.refOffset-pId*Ymax*G2frame.refDelt*np.ones_like(peak)
1563                if G2frame.qPlot:
1564                    Plot.plot(2*np.pi/peak.T[0],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
1565                elif G2frame.dPlot:
1566                    Plot.plot(peak.T[0],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
1567                else:
1568                    Plot.plot(peak.T[1],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
1569            if len(Phases):
1570                handles,legends = Plot.get_legend_handles_labels()  #got double entries in the legends for some reason
1571                if handles:
1572                    Plot.legend(handles[::2],legends[::2],title='Phases',loc='best')    #skip every other one
1573           
1574    if G2frame.Contour:
1575        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
1576        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
1577            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
1578        Page.figure.colorbar(Img)
1579    else:
1580        G2frame.Lines = Lines
1581    if not newPlot:
1582        Page.toolbar.push_current()
1583        Plot.set_xlim(xylim[0])
1584        Plot.set_ylim(xylim[1])
1585#        xylim = []
1586        Page.toolbar.push_current()
1587        Page.toolbar.draw()
1588    else:
1589        Page.canvas.draw()
1590    olderr = np.seterr(invalid='ignore') #ugh - this removes a matplotlib error for mouse clicks in log plots
1591    # and sqrt(-ve) in np.where usage               
1592#    G2frame.Pwdr = True
1593   
1594################################################################################
1595##### PlotDeltSig
1596################################################################################
1597           
1598def PlotDeltSig(G2frame,kind):
1599    'needs a doc string'
1600    try:
1601        plotNum = G2frame.G2plotNB.plotList.index('Error analysis')
1602        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1603        Page.figure.clf()
1604        Plot = Page.figure.gca()          #get a fresh plot after clf()
1605    except ValueError:
1606        newPlot = True
1607        G2frame.Cmax = 1.0
1608        Plot = G2frame.G2plotNB.addMpl('Error analysis').gca()
1609        plotNum = G2frame.G2plotNB.plotList.index('Error analysis')
1610        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1611    Page.Choice = None
1612    PatternId = G2frame.PatternId
1613    Pattern = G2frame.PatternTree.GetItemPyData(PatternId)
1614    Pattern.append(G2frame.PatternTree.GetItemText(PatternId))
1615    wtFactor = Pattern[0]['wtFactor']
1616    if kind == 'PWDR':
1617        limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits'))[1]
1618        xye = np.array(Pattern[1])
1619        xmin = np.searchsorted(xye[0],limits[0])
1620        xmax = np.searchsorted(xye[0],limits[1])
1621        DS = xye[5][xmin:xmax]*np.sqrt(wtFactor*xye[2][xmin:xmax])
1622    elif kind == 'HKLF':
1623        refl = Pattern[1]
1624        DS = []
1625        for ref in refl:
1626            if ref[6] > 0.:
1627                DS.append((ref[5]-ref[7])/ref[6])
1628    Page.SetFocus()
1629    G2frame.G2plotNB.status.DestroyChildren()
1630    DS.sort()
1631    EDS = np.zeros_like(DS)
1632    DX = np.linspace(0.,1.,num=len(DS),endpoint=True)
1633    oldErr = np.seterr(invalid='ignore')    #avoid problem at DX==0
1634    T = np.sqrt(np.log(1.0/DX**2))
1635    top = 2.515517+0.802853*T+0.010328*T**2
1636    bot = 1.0+1.432788*T+0.189269*T**2+0.001308*T**3
1637    EDS = np.where(DX>0,-(T-top/bot),(T-top/bot))
1638    low1 = np.searchsorted(EDS,-1.)
1639    hi1 = np.searchsorted(EDS,1.)
1640    slp,intcp = np.polyfit(EDS[low1:hi1],DS[low1:hi1],deg=1)
1641    frac = 100.*(hi1-low1)/len(DS)
1642    G2frame.G2plotNB.status.SetStatusText(  \
1643        'Over range -1. to 1. :'+' slope = %.3f, intercept = %.3f for %.2f%% of the fitted data'%(slp,intcp,frac),1)
1644    Plot.set_title('Normal probability for '+Pattern[-1])
1645    Plot.set_xlabel(r'expected $\mathsf{\Delta/\sigma}$',fontsize=14)
1646    Plot.set_ylabel(r'observed $\mathsf{\Delta/\sigma}$',fontsize=14)
1647    Plot.plot(EDS,DS,'r+',label='result')
1648    Plot.plot([-2,2],[-2,2],'k',dashes=(5,5),label='ideal')
1649    Plot.legend(loc='upper left')
1650    np.seterr(invalid='warn')
1651    Page.canvas.draw()
1652       
1653################################################################################
1654##### PlotISFG
1655################################################################################
1656           
1657def PlotISFG(G2frame,newPlot=False,type=''):
1658    ''' Plotting package for PDF analysis; displays I(q), S(q), F(q) and G(r) as single
1659    or multiple plots with waterfall and contour plots as options
1660    '''
1661    if not type:
1662        type = G2frame.G2plotNB.plotList[G2frame.G2plotNB.nb.GetSelection()]
1663    if type not in ['I(Q)','S(Q)','F(Q)','G(R)']:
1664        return
1665    superMinusOne = unichr(0xaf)+unichr(0xb9)
1666   
1667    def OnPlotKeyPress(event):
1668        newPlot = False
1669        if event.key == 'u':
1670            if G2frame.Contour:
1671                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
1672            elif G2frame.Offset[0] < 100.:
1673                G2frame.Offset[0] += 1.
1674        elif event.key == 'd':
1675            if G2frame.Contour:
1676                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
1677            elif G2frame.Offset[0] > 0.:
1678                G2frame.Offset[0] -= 1.
1679        elif event.key == 'l':
1680            G2frame.Offset[1] -= 1.
1681        elif event.key == 'r':
1682            G2frame.Offset[1] += 1.
1683        elif event.key == 'o':
1684            G2frame.Offset = [0,0]
1685        elif event.key == 'c':
1686            newPlot = True
1687            G2frame.Contour = not G2frame.Contour
1688            if not G2frame.Contour:
1689                G2frame.SinglePlot = False
1690                G2frame.Offset = [0.,0.]
1691        elif event.key == 's':
1692            if G2frame.Contour:
1693                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
1694                choice.sort()
1695                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
1696                if dlg.ShowModal() == wx.ID_OK:
1697                    sel = dlg.GetSelection()
1698                    G2frame.ContourColor = choice[sel]
1699                else:
1700                    G2frame.ContourColor = 'Paired'
1701                dlg.Destroy()
1702            else:
1703                G2frame.SinglePlot = not G2frame.SinglePlot               
1704        elif event.key == 'i':                  #for smoothing contour plot
1705            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
1706               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
1707               'mitchell','sinc','lanczos']
1708            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
1709            if dlg.ShowModal() == wx.ID_OK:
1710                sel = dlg.GetSelection()
1711                G2frame.Interpolate = choice[sel]
1712            else:
1713                G2frame.Interpolate = 'nearest'
1714            dlg.Destroy()
1715        elif event.key == 't' and not G2frame.Contour:
1716            G2frame.Legend = not G2frame.Legend
1717        PlotISFG(G2frame,newPlot=newPlot,type=type)
1718       
1719    def OnKeyBox(event):
1720        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index(type):
1721            event.key = cb.GetValue()[0]
1722            cb.SetValue(' key press')
1723            wx.CallAfter(OnPlotKeyPress,event)
1724        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
1725                       
1726    def OnMotion(event):
1727        xpos = event.xdata
1728        if xpos:                                        #avoid out of frame mouse position
1729            ypos = event.ydata
1730            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1731            try:
1732                if G2frame.Contour:
1733                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA pattern ID =%5d'%(xpos,int(ypos)),1)
1734                else:
1735                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA %s =%.2f'%(xpos,type,ypos),1)                   
1736            except TypeError:
1737                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
1738   
1739    xylim = []
1740    try:
1741        plotNum = G2frame.G2plotNB.plotList.index(type)
1742        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1743        if not newPlot:
1744            Plot = Page.figure.gca()          #get previous plot & get limits
1745            xylim = Plot.get_xlim(),Plot.get_ylim()
1746        Page.figure.clf()
1747        Plot = Page.figure.gca()
1748    except ValueError:
1749        newPlot = True
1750        G2frame.Cmax = 1.0
1751        Plot = G2frame.G2plotNB.addMpl(type).gca()
1752        plotNum = G2frame.G2plotNB.plotList.index(type)
1753        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1754        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1755        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1756   
1757    Page.SetFocus()
1758    G2frame.G2plotNB.status.DestroyChildren()
1759    if G2frame.Contour:
1760        Page.Choice = (' key press','d: lower contour max','u: raise contour max',
1761            'i: interpolation method','s: color scheme','c: contour off')
1762    else:
1763        Page.Choice = (' key press','l: offset left','r: offset right','d: offset down','u: offset up',
1764            'o: reset offset','t: toggle legend','c: contour on','s: toggle single plot')
1765    Page.keyPress = OnPlotKeyPress
1766    PatternId = G2frame.PatternId
1767    PickId = G2frame.PickId
1768    Plot.set_title(type)
1769    if type == 'G(R)':
1770        Plot.set_xlabel(r'$R,\AA$',fontsize=14)
1771    else:
1772        Plot.set_xlabel(r'$Q,\AA$'+superMinusOne,fontsize=14)
1773    Plot.set_ylabel(r''+type,fontsize=14)
1774    colors=['b','g','r','c','m','k']
1775    name = G2frame.PatternTree.GetItemText(PatternId)[4:]
1776    Pattern = []   
1777    if G2frame.SinglePlot:
1778        name = G2frame.PatternTree.GetItemText(PatternId)
1779        name = type+name[4:]
1780        Id = G2gd.GetPatternTreeItemId(G2frame,PatternId,name)
1781        Pattern = G2frame.PatternTree.GetItemPyData(Id)
1782        if Pattern:
1783            Pattern.append(name)
1784        PlotList = [Pattern,]
1785    else:
1786        PlotList = []
1787        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
1788        while item:
1789            if 'PDF' in G2frame.PatternTree.GetItemText(item):
1790                name = type+G2frame.PatternTree.GetItemText(item)[4:]
1791                Id = G2gd.GetPatternTreeItemId(G2frame,item,name)
1792                Pattern = G2frame.PatternTree.GetItemPyData(Id)
1793                if Pattern:
1794                    Pattern.append(name)
1795                    PlotList.append(Pattern)
1796            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
1797    PDFdata = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'PDF Controls'))
1798    numbDen = G2pwd.GetNumDensity(PDFdata['ElList'],PDFdata['Form Vol'])
1799    Xb = [0.,10.]
1800    Yb = [0.,-40.*np.pi*numbDen]
1801    Ymax = 1.0
1802    lenX = 0
1803    for Pattern in PlotList:
1804        xye = Pattern[1]
1805        Ymax = max(Ymax,max(xye[1]))
1806    offset = G2frame.Offset[0]*Ymax/100.0
1807    if G2frame.Contour:
1808        ContourZ = []
1809        ContourY = []
1810        Nseq = 0
1811    for N,Pattern in enumerate(PlotList):
1812        xye = Pattern[1]
1813        if PickId:
1814            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
1815        X = xye[0]
1816        if not lenX:
1817            lenX = len(X)           
1818        Y = xye[1]+offset*N
1819        if G2frame.Contour:
1820            if lenX == len(X):
1821                ContourY.append(N)
1822                ContourZ.append(Y)
1823                ContourX = X
1824                Nseq += 1
1825                Plot.set_ylabel('Data sequence',fontsize=12)
1826        else:
1827            X = xye[0]+G2frame.Offset[1]*.005*N
1828            if ifpicked:
1829                Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
1830                Page.canvas.SetToolTipString('')
1831            else:
1832                if G2frame.Legend:
1833                    Plot.plot(X,Y,colors[N%6],picker=False,label='Azm:'+Pattern[2].split('=')[1])
1834                else:
1835                    Plot.plot(X,Y,colors[N%6],picker=False)
1836            if type == 'G(R)':
1837                Plot.plot(Xb,Yb,color='k',dashes=(5,5))
1838            elif type == 'F(Q)':
1839                Plot.axhline(0.,color=wx.BLACK)
1840            elif type == 'S(Q)':
1841                Plot.axhline(1.,color=wx.BLACK)
1842    if G2frame.Contour:
1843        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
1844        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
1845            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
1846        Page.figure.colorbar(Img)
1847    elif G2frame.Legend:
1848        Plot.legend(loc='best')
1849    if not newPlot:
1850        Page.toolbar.push_current()
1851        Plot.set_xlim(xylim[0])
1852        Plot.set_ylim(xylim[1])
1853        xylim = []
1854        Page.toolbar.push_current()
1855        Page.toolbar.draw()
1856    else:
1857        Page.canvas.draw()
1858       
1859################################################################################
1860##### PlotCalib
1861################################################################################
1862           
1863def PlotCalib(G2frame,Inst,XY,Sigs,newPlot=False):
1864    '''plot of CW or TOF peak calibration
1865    '''
1866    def OnMotion(event):
1867        xpos = event.xdata
1868        if xpos:                                        #avoid out of frame mouse position
1869            ypos = event.ydata
1870            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1871            try:
1872                G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3g'%(xpos,Title,ypos),1)                   
1873            except TypeError:
1874                G2frame.G2plotNB.status.SetStatusText('Select '+Title+' pattern first',1)
1875            found = []
1876            wid = 1
1877            view = Page.toolbar._views.forward()
1878            if view:
1879                view = view[0][:2]
1880                wid = view[1]-view[0]
1881            found = XY[np.where(np.fabs(XY.T[0]-xpos) < 0.005*wid)]
1882            if len(found):
1883                pos = found[0][1]
1884                if 'C' in Inst['Type'][0]: 
1885                    Page.canvas.SetToolTipString('position=%.4f'%(pos))
1886                else:
1887                    Page.canvas.SetToolTipString('position=%.2f'%(pos))
1888            else:
1889                Page.canvas.SetToolTipString('')
1890
1891    Title = 'Position calibration'
1892    try:
1893        plotNum = G2frame.G2plotNB.plotList.index(Title)
1894        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1895        if not newPlot:
1896            Plot = Page.figure.gca()
1897            xylim = Plot.get_xlim(),Plot.get_ylim()
1898        Page.figure.clf()
1899        Plot = Page.figure.gca()
1900    except ValueError:
1901        newPlot = True
1902        Plot = G2frame.G2plotNB.addMpl(Title).gca()
1903        plotNum = G2frame.G2plotNB.plotList.index(Title)
1904        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1905        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1906   
1907    Page.Choice = None
1908    Page.SetFocus()
1909    G2frame.G2plotNB.status.DestroyChildren()
1910    Plot.set_title(Title)
1911    Plot.set_xlabel(r'd-spacing',fontsize=14)
1912    if 'C' in Inst['Type'][0]:
1913        Plot.set_ylabel(r'$\mathsf{\Delta(2\theta)}$',fontsize=14)
1914    else:
1915        Plot.set_ylabel(r'$\mathsf{\Delta}T/T$',fontsize=14)
1916    for ixy,xyw in enumerate(XY):
1917        if len(xyw) > 2:
1918            X,Y,W = xyw
1919        else:
1920            X,Y = xyw
1921            W = 0.
1922        Yc = G2lat.Dsp2pos(Inst,X)
1923        if 'C' in Inst['Type'][0]:
1924            Y = Y-Yc
1925            E = Sigs[ixy]
1926            bin = W/2.
1927        else:
1928            Y = (Y-Yc)/Yc
1929            E = Sigs[ixy]/Yc
1930            bin = W/(2.*Yc)
1931        if E:
1932            Plot.errorbar(X,Y,ecolor='k',yerr=E)
1933        if ixy:
1934            Plot.plot(X,Y,'kx',picker=3)
1935        else:
1936            Plot.plot(X,Y,'kx',label='peak')
1937        if W:
1938            if ixy:
1939                Plot.plot(X,bin,'b+')
1940            else:
1941                Plot.plot(X,bin,'b+',label='bin width')
1942            Plot.plot(X,-bin,'b+')
1943        Plot.axhline(0.,color='r',linestyle='--')
1944    Plot.legend(loc='best')
1945    if not newPlot:
1946        Page.toolbar.push_current()
1947        Plot.set_xlim(xylim[0])
1948        Plot.set_ylim(xylim[1])
1949#        xylim = []
1950        Page.toolbar.push_current()
1951        Page.toolbar.draw()
1952    else:
1953        Page.canvas.draw()
1954
1955################################################################################
1956##### PlotXY
1957################################################################################
1958           
1959def PlotXY(G2frame,XY,XY2=None,labelX=None,labelY=None,newPlot=False,Title=''):
1960    '''simple plot of xy data, used for diagnostic purposes
1961    '''
1962    def OnMotion(event):
1963        xpos = event.xdata
1964        if xpos:                                        #avoid out of frame mouse position
1965            ypos = event.ydata
1966            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1967            try:
1968                G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,Title,ypos),1)                   
1969            except TypeError:
1970                G2frame.G2plotNB.status.SetStatusText('Select '+Title+' pattern first',1)
1971
1972    try:
1973        plotNum = G2frame.G2plotNB.plotList.index(Title)
1974        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1975        if not newPlot:
1976            Plot = Page.figure.gca()
1977            xylim = Plot.get_xlim(),Plot.get_ylim()
1978        Page.figure.clf()
1979        Plot = Page.figure.gca()
1980    except ValueError:
1981        newPlot = True
1982        Plot = G2frame.G2plotNB.addMpl(Title).gca()
1983        plotNum = G2frame.G2plotNB.plotList.index(Title)
1984        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1985        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1986   
1987    Page.Choice = None
1988    Page.SetFocus()
1989    G2frame.G2plotNB.status.DestroyChildren()
1990    Plot.set_title(Title)
1991    if labelX:
1992        Plot.set_xlabel(r''+labelX,fontsize=14)
1993    else:
1994        Plot.set_xlabel(r'X',fontsize=14)
1995    if labelY:
1996        Plot.set_ylabel(r''+labelY,fontsize=14)
1997    else:
1998        Plot.set_ylabel(r'Y',fontsize=14)
1999    colors=['b','g','r','c','m','k']
2000    for ixy,xy in enumerate(XY):
2001        X,Y = xy
2002        Plot.plot(X,Y,colors[ixy%6]+'+',picker=False)
2003    if len(XY2):
2004        for ixy,xy in enumerate(XY2):
2005            X,Y = xy
2006            Plot.plot(X,Y,colors[ixy%6],picker=False)
2007    if not newPlot:
2008        Page.toolbar.push_current()
2009        Plot.set_xlim(xylim[0])
2010        Plot.set_ylim(xylim[1])
2011        xylim = []
2012        Page.toolbar.push_current()
2013        Page.toolbar.draw()
2014    else:
2015        Page.canvas.draw()
2016
2017################################################################################
2018##### PlotStrain
2019################################################################################
2020           
2021def PlotStrain(G2frame,data,newPlot=False):
2022    '''plot of strain data, used for diagnostic purposes
2023    '''
2024    def OnMotion(event):
2025        xpos = event.xdata
2026        if xpos:                                        #avoid out of frame mouse position
2027            ypos = event.ydata
2028            Page.canvas.SetCursor(wx.CROSS_CURSOR)
2029            try:
2030                G2frame.G2plotNB.status.SetStatusText('d-spacing =%9.5f Azimuth =%9.3f'%(ypos,xpos),1)                   
2031            except TypeError:
2032                G2frame.G2plotNB.status.SetStatusText('Select Strain pattern first',1)
2033
2034    try:
2035        plotNum = G2frame.G2plotNB.plotList.index('Strain')
2036        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2037        if not newPlot:
2038            Plot = Page.figure.gca()
2039            xylim = Plot.get_xlim(),Plot.get_ylim()
2040        Page.figure.clf()
2041        Plot = Page.figure.gca()
2042    except ValueError:
2043        newPlot = True
2044        Plot = G2frame.G2plotNB.addMpl('Strain').gca()
2045        plotNum = G2frame.G2plotNB.plotList.index('Strain')
2046        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2047        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2048   
2049    Page.Choice = None
2050    Page.SetFocus()
2051    G2frame.G2plotNB.status.DestroyChildren()
2052    Plot.set_title('Strain')
2053    Plot.set_ylabel(r'd-spacing',fontsize=14)
2054    Plot.set_xlabel(r'Azimuth',fontsize=14)
2055    colors=['b','g','r','c','m','k']
2056    for N,item in enumerate(data['d-zero']):
2057        Y,X = np.array(item['ImtaObs'])         #plot azimuth as X & d-spacing as Y
2058        Plot.plot(X,Y,colors[N%6]+'+',picker=False)
2059        Y,X = np.array(item['ImtaCalc'])
2060        Plot.plot(X,Y,colors[N%6],picker=False)
2061    if not newPlot:
2062        Page.toolbar.push_current()
2063        Plot.set_xlim(xylim[0])
2064        Plot.set_ylim(xylim[1])
2065        xylim = []
2066        Page.toolbar.push_current()
2067        Page.toolbar.draw()
2068    else:
2069        Page.canvas.draw()
2070       
2071################################################################################
2072##### PlotSASDSizeDist
2073################################################################################
2074           
2075def PlotSASDSizeDist(G2frame):
2076   
2077    def OnPageChanged(event):
2078        PlotText = G2frame.G2plotNB.nb.GetPageText(G2frame.G2plotNB.nb.GetSelection())
2079        if 'Powder' in PlotText:
2080            PlotPatterns(G2frame,plotType='SASD',newPlot=True)
2081        elif 'Size' in PlotText:
2082            PlotSASDSizeDist(G2frame)
2083   
2084    def OnMotion(event):
2085        xpos = event.xdata
2086        if xpos:                                        #avoid out of frame mouse position
2087            ypos = event.ydata
2088            Page.canvas.SetCursor(wx.CROSS_CURSOR)
2089            try:
2090                G2frame.G2plotNB.status.SetStatusText('diameter =%9.3f f(D) =%9.3g'%(xpos,ypos),1)                   
2091            except TypeError:
2092                G2frame.G2plotNB.status.SetStatusText('Select Strain pattern first',1)
2093
2094    try:
2095        plotNum = G2frame.G2plotNB.plotList.index('Size Distribution')
2096        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2097        Page.figure.clf()
2098        Plot = Page.figure.gca()          #get a fresh plot after clf()
2099    except ValueError:
2100        newPlot = True
2101        Plot = G2frame.G2plotNB.addMpl('Size Distribution').gca()
2102        G2frame.G2plotNB.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED,OnPageChanged)
2103        plotNum = G2frame.G2plotNB.plotList.index('Size Distribution')
2104        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2105        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2106    Page.Choice = None
2107    Page.SetFocus()
2108    PatternId = G2frame.PatternId
2109    data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Models'))
2110    Bins,Dbins,BinMag = data['Size']['Distribution']
2111    Plot.set_title('Size Distribution')
2112    Plot.set_xlabel(r'$D, \AA$',fontsize=14)
2113    Plot.set_ylabel(r'$Volume distribution f(D)$',fontsize=14)
2114    if data['Size']['logBins']:
2115        Plot.set_xscale("log",nonposy='mask')
2116        Plot.set_xlim([np.min(2.*Bins)/2.,np.max(2.*Bins)*2.])
2117    Plot.bar(2.*Bins-Dbins,BinMag,2.*Dbins,facecolor='white')       #plot diameters
2118    if 'Size Calc' in data:
2119        Rbins,Dist = data['Size Calc']
2120        for i in range(len(Rbins)):
2121            if len(Rbins[i]):
2122                Plot.plot(2.*Rbins[i],Dist[i])       #plot diameters
2123    Page.canvas.draw()
2124
2125################################################################################
2126##### PlotPowderLines
2127################################################################################
2128           
2129def PlotPowderLines(G2frame):
2130    ''' plotting of powder lines (i.e. no powder pattern) as sticks
2131    '''
2132    global HKL
2133
2134    def OnMotion(event):
2135        xpos = event.xdata
2136        if xpos:                                        #avoid out of frame mouse position
2137            Page.canvas.SetCursor(wx.CROSS_CURSOR)
2138            G2frame.G2plotNB.status.SetFields(['','2-theta =%9.3f '%(xpos,)])
2139            if G2frame.PickId and G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List']:
2140                found = []
2141                if len(HKL):
2142                    view = Page.toolbar._views.forward()[0][:2]
2143                    wid = view[1]-view[0]
2144                    found = HKL[np.where(np.fabs(HKL.T[-1]-xpos) < 0.002*wid)]
2145                if len(found):
2146                    h,k,l = found[0][:3] 
2147                    Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
2148                else:
2149                    Page.canvas.SetToolTipString('')
2150
2151    try:
2152        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
2153        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2154        Page.figure.clf()
2155        Plot = Page.figure.gca()
2156    except ValueError:
2157        Plot = G2frame.G2plotNB.addMpl('Powder Lines').gca()
2158        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
2159        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2160        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2161       
2162    Page.Choice = None
2163    Page.SetFocus()
2164    Plot.set_title('Powder Pattern Lines')
2165    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
2166    PickId = G2frame.PickId
2167    PatternId = G2frame.PatternId
2168    peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))[0]
2169    for peak in peaks:
2170        Plot.axvline(peak[0],color='b')
2171    HKL = np.array(G2frame.HKL)
2172    for hkl in G2frame.HKL:
2173        Plot.axvline(hkl[-1],color='r',dashes=(5,5))
2174    xmin = peaks[0][0]
2175    xmax = peaks[-1][0]
2176    delt = xmax-xmin
2177    xlim = [max(0,xmin-delt/20.),min(180.,xmax+delt/20.)]
2178    Plot.set_xlim(xlim)
2179    Page.canvas.draw()
2180    Page.toolbar.push_current()
2181
2182################################################################################
2183##### PlotPeakWidths
2184################################################################################
2185           
2186def PlotPeakWidths(G2frame):
2187    ''' Plotting of instrument broadening terms as function of 2-theta
2188    Seen when "Instrument Parameters" chosen from powder pattern data tree
2189    '''
2190#    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
2191#    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
2192    gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
2193#    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
2194    PatternId = G2frame.PatternId
2195    limitID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
2196    if limitID:
2197        limits = G2frame.PatternTree.GetItemPyData(limitID)[:2]
2198    else:
2199        return
2200    Parms,Parms2 = G2frame.PatternTree.GetItemPyData( \
2201        G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
2202    if 'C' in Parms['Type'][0]:
2203        lam = G2mth.getWave(Parms)
2204    else:
2205        difC = Parms['difC'][0]
2206    try:  # PATCH: deal with older peak lists, before changed to dict to implement TOF
2207        peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'))['peaks']
2208    except TypeError:
2209        print "Your peak list needs reformatting...",
2210        item = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List')
2211        G2frame.PatternTree.SelectItem(item) 
2212        item = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters')
2213        G2frame.PatternTree.SelectItem(item)
2214        print "done"
2215        return
2216    try:
2217        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
2218        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2219        Page.figure.clf()
2220        Plot = Page.figure.gca()
2221    except ValueError:
2222        Plot = G2frame.G2plotNB.addMpl('Peak Widths').gca()
2223        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
2224        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2225    Page.Choice = None
2226    Page.SetFocus()
2227   
2228    Page.canvas.SetToolTipString('')
2229    colors=['b','g','r','c','m','k']
2230    X = []
2231    Y = []
2232    Z = []
2233    W = []
2234    if 'C' in Parms['Type'][0]:
2235        Plot.set_title('Instrument and sample peak widths')
2236        Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=14)
2237        Plot.set_ylabel(r'$\Delta Q/Q, \Delta d/d$',fontsize=14)
2238        try:
2239            Xmin,Xmax = limits[1]
2240            X = np.linspace(Xmin,Xmax,num=101,endpoint=True)
2241            Q = 4.*np.pi*npsind(X/2.)/lam
2242            Z = np.ones_like(X)
2243            data = G2mth.setPeakparms(Parms,Parms2,X,Z)
2244            s = 1.17741*np.sqrt(data[4])*np.pi/18000.
2245            g = data[6]*np.pi/18000.
2246            G = G2pwd.getgamFW(g,s)
2247            Y = s/nptand(X/2.)
2248            Z = g/nptand(X/2.)
2249            W = G/nptand(X/2.)
2250            Plot.plot(Q,Y,color='r',label='Gaussian')
2251            Plot.plot(Q,Z,color='g',label='Lorentzian')
2252            Plot.plot(Q,W,color='b',label='G+L')
2253           
2254            fit = G2mth.setPeakparms(Parms,Parms2,X,Z,useFit=True)
2255            sf = 1.17741*np.sqrt(fit[4])*np.pi/18000.
2256            gf = fit[6]*np.pi/18000.
2257            Gf = G2pwd.getgamFW(gf,sf)
2258            Yf = sf/nptand(X/2.)
2259            Zf = gf/nptand(X/2.)
2260            Wf = Gf/nptand(X/2.)
2261            Plot.plot(Q,Yf,color='r',dashes=(5,5),label='Gaussian fit')
2262            Plot.plot(Q,Zf,color='g',dashes=(5,5),label='Lorentzian fit')
2263            Plot.plot(Q,Wf,color='b',dashes=(5,5),label='G+L fit')
2264           
2265            X = []
2266            Y = []
2267            Z = []
2268            W = []
2269            V = []
2270            for peak in peaks:
2271                X.append(4.0*math.pi*sind(peak[0]/2.0)/lam)
2272                try:
2273                    s = 1.17741*math.sqrt(peak[4])*math.pi/18000.
2274                except ValueError:
2275                    s = 0.01
2276                g = peak[6]*math.pi/18000.
2277                G = G2pwd.getgamFW(g,s)
2278                Y.append(s/tand(peak[0]/2.))
2279                Z.append(g/tand(peak[0]/2.))
2280                W.append(G/tand(peak[0]/2.))
2281            if len(peaks):
2282                Plot.plot(X,Y,'+',color='r',label='G peak')
2283                Plot.plot(X,Z,'+',color='g',label='L peak')
2284                Plot.plot(X,W,'+',color='b',label='G+L peak')
2285            Plot.legend(loc='best')
2286            Page.canvas.draw()
2287        except ValueError:
2288            print '**** ERROR - default U,V,W profile coefficients yield sqrt of negative value at 2theta =', \
2289                '%.3f'%(2*theta)
2290            G2frame.G2plotNB.Delete('Peak Widths')
2291    else:   #'T'OF
2292        Plot.set_title('Instrument and sample peak coefficients')
2293        Plot.set_xlabel(r'$Q, \AA^{-1}$',fontsize=14)
2294        Plot.set_ylabel(r'$\alpha, \beta, \Delta Q/Q, \Delta d/d$',fontsize=14)
2295        Xmin,Xmax = limits[1]
2296        T = np.linspace(Xmin,Xmax,num=101,endpoint=True)
2297        Z = np.ones_like(T)
2298        data = G2mth.setPeakparms(Parms,Parms2,T,Z)
2299        ds = T/difC
2300        Q = 2.*np.pi/ds
2301        A = data[4]
2302        B = data[6]
2303        S = 1.17741*np.sqrt(data[8])/T
2304        G = data[10]/T
2305        Plot.plot(Q,A,color='r',label='Alpha')
2306        Plot.plot(Q,B,color='g',label='Beta')
2307        Plot.plot(Q,S,color='b',label='Gaussian')
2308        Plot.plot(Q,G,color='m',label='Lorentzian')
2309
2310        fit = G2mth.setPeakparms(Parms,Parms2,T,Z)
2311        ds = T/difC
2312        Q = 2.*np.pi/ds
2313        Af = fit[4]
2314        Bf = fit[6]
2315        Sf = 1.17741*np.sqrt(fit[8])/T
2316        Gf = fit[10]/T
2317        Plot.plot(Q,Af,color='r',dashes=(5,5),label='Alpha fit')
2318        Plot.plot(Q,Bf,color='g',dashes=(5,5),label='Beta fit')
2319        Plot.plot(Q,Sf,color='b',dashes=(5,5),label='Gaussian fit')
2320        Plot.plot(Q,Gf,color='m',dashes=(5,5),label='Lorentzian fit')
2321       
2322        T = []
2323        A = []
2324        B = []
2325        S = []
2326        G = []
2327        W = []
2328        Q = []
2329        V = []
2330        for peak in peaks:
2331            T.append(peak[0])
2332            A.append(peak[4])
2333            B.append(peak[6])
2334            Q.append(2.*np.pi*difC/peak[0])
2335            S.append(1.17741*np.sqrt(peak[8])/peak[0])
2336            G.append(peak[10]/peak[0])
2337           
2338       
2339        Plot.plot(Q,A,'+',color='r',label='Alpha peak')
2340        Plot.plot(Q,B,'+',color='g',label='Beta peak')
2341        Plot.plot(Q,S,'+',color='b',label='Gaussian peak')
2342        Plot.plot(Q,G,'+',color='m',label='Lorentzian peak')
2343        Plot.legend(loc='best')
2344        Page.canvas.draw()
2345
2346   
2347################################################################################
2348##### PlotSizeStrainPO
2349################################################################################
2350           
2351def PlotSizeStrainPO(G2frame,data,Start=False):
2352    '''Plot 3D mustrain/size/preferred orientation figure. In this instance data is for a phase
2353    '''
2354   
2355    PatternId = G2frame.PatternId
2356    generalData = data['General']
2357    SGData = generalData['SGData']
2358    SGLaue = SGData['SGLaue']
2359    if Start:                   #initialize the spherical harmonics qlmn arrays
2360        ptx.pyqlmninit()
2361        Start = False
2362#    MuStrKeys = G2spc.MustrainNames(SGData)
2363    cell = generalData['Cell'][1:]
2364    A,B = G2lat.cell2AB(cell[:6])
2365    Vol = cell[6]
2366    useList = data['Histograms']
2367    phase = generalData['Name']
2368    plotType = generalData['Data plot type']
2369    plotDict = {'Mustrain':'Mustrain','Size':'Size','Preferred orientation':'Pref.Ori.'}
2370    for ptype in plotDict:
2371        G2frame.G2plotNB.Delete(ptype)
2372    if plotType in ['None']:
2373        return       
2374
2375    for item in useList:
2376        if useList[item]['Show']:
2377            break
2378    else:
2379        return            #nothing to show!!
2380   
2381    numPlots = len(useList)
2382
2383    if plotType in ['Mustrain','Size']:
2384        Plot = mp3d.Axes3D(G2frame.G2plotNB.add3D(plotType))
2385    else:
2386        Plot = G2frame.G2plotNB.addMpl(plotType).gca()       
2387    plotNum = G2frame.G2plotNB.plotList.index(plotType)
2388    Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2389    Page.Choice = None
2390    Page.SetFocus()
2391    G2frame.G2plotNB.status.SetStatusText('',1)
2392    if not Page.IsShown():
2393        Page.Show()
2394   
2395    for item in useList:
2396        if useList[item]['Show']:
2397            PHI = np.linspace(0.,360.,30,True)
2398            PSI = np.linspace(0.,180.,30,True)
2399            X = np.outer(npsind(PHI),npsind(PSI))
2400            Y = np.outer(npcosd(PHI),npsind(PSI))
2401            Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
2402            try:        #temp patch instead of 'mustrain' for old files with 'microstrain'
2403                coeff = useList[item][plotDict[plotType]]
2404            except KeyError:
2405                break
2406            if plotType in ['Mustrain','Size']:
2407                if coeff[0] == 'isotropic':
2408                    X *= coeff[1][0]
2409                    Y *= coeff[1][0]
2410                    Z *= coeff[1][0]                               
2411                elif coeff[0] == 'uniaxial':
2412                   
2413                    def uniaxCalc(xyz,iso,aniso,axes):
2414                        Z = np.array(axes)
2415                        cp = abs(np.dot(xyz,Z))
2416                        sp = np.sqrt(1.-cp**2)
2417                        R = iso*aniso/np.sqrt((iso*cp)**2+(aniso*sp)**2)
2418                        return R*xyz
2419                       
2420                    iso,aniso = coeff[1][:2]
2421                    axes = np.inner(A,np.array(coeff[3]))
2422                    axes /= nl.norm(axes)
2423                    Shkl = np.array(coeff[1])
2424                    XYZ = np.dstack((X,Y,Z))
2425                    XYZ = np.nan_to_num(np.apply_along_axis(uniaxCalc,2,XYZ,iso,aniso,axes))
2426                    X,Y,Z = np.dsplit(XYZ,3)
2427                    X = X[:,:,0]
2428                    Y = Y[:,:,0]
2429                    Z = Z[:,:,0]
2430               
2431                elif coeff[0] == 'ellipsoidal':
2432                   
2433                    def ellipseCalc(xyz,E,R):
2434                        XYZ = xyz*E.T
2435                        return np.inner(XYZ.T,R)
2436                       
2437                    S6 = coeff[4]
2438                    Sij = G2lat.U6toUij(S6)
2439                    E,R = nl.eigh(Sij)
2440                    XYZ = np.dstack((X,Y,Z))
2441                    XYZ = np.nan_to_num(np.apply_along_axis(ellipseCalc,2,XYZ,E,R))
2442                    X,Y,Z = np.dsplit(XYZ,3)
2443                    X = X[:,:,0]
2444                    Y = Y[:,:,0]
2445                    Z = Z[:,:,0]
2446                   
2447                elif coeff[0] == 'generalized':
2448                   
2449                    def genMustrain(xyz,SGData,A,Shkl):
2450                        uvw = np.inner(A.T,xyz)
2451                        Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
2452                        Sum = np.sum(np.multiply(Shkl,Strm))
2453                        Sum = np.where(Sum > 0.01,Sum,0.01)
2454                        Sum = np.sqrt(Sum)
2455                        return Sum*xyz
2456                       
2457                    Shkl = np.array(coeff[4])
2458                    if np.any(Shkl):
2459                        XYZ = np.dstack((X,Y,Z))
2460                        XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,SGData,A,Shkl))
2461                        X,Y,Z = np.dsplit(XYZ,3)
2462                        X = X[:,:,0]
2463                        Y = Y[:,:,0]
2464                        Z = Z[:,:,0]
2465                           
2466                if np.any(X) and np.any(Y) and np.any(Z):
2467                    errFlags = np.seterr(all='ignore')
2468                    Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1)
2469                    np.seterr(all='ignore')
2470                    xyzlim = np.array([Plot.get_xlim3d(),Plot.get_ylim3d(),Plot.get_zlim3d()]).T
2471                    XYZlim = [min(xyzlim[0]),max(xyzlim[1])]
2472                    Plot.set_xlim3d(XYZlim)
2473                    Plot.set_ylim3d(XYZlim)
2474                    Plot.set_zlim3d(XYZlim)
2475                    Plot.set_aspect('equal')
2476                if plotType == 'Size':
2477                    Plot.set_title('Crystallite size for '+phase+'\n'+coeff[0]+' model')
2478                    Plot.set_xlabel(r'X, $\mu$m')
2479                    Plot.set_ylabel(r'Y, $\mu$m')
2480                    Plot.set_zlabel(r'Z, $\mu$m')
2481                else:   
2482                    Plot.set_title(r'$\mu$strain for '+phase+'\n'+coeff[0]+' model')
2483                    Plot.set_xlabel(r'X, $\mu$strain')
2484                    Plot.set_ylabel(r'Y, $\mu$strain')
2485                    Plot.set_zlabel(r'Z, $\mu$strain')
2486            else:
2487                h,k,l = generalData['POhkl']
2488                if coeff[0] == 'MD':
2489                    print 'March-Dollase preferred orientation plot'
2490               
2491                else:
2492                    PH = np.array(generalData['POhkl'])
2493                    phi,beta = G2lat.CrsAng(PH,cell[:6],SGData)
2494                    SHCoef = {}
2495                    for item in coeff[5]:
2496                        L,N = eval(item.strip('C'))
2497                        SHCoef['C%d,0,%d'%(L,N)] = coeff[5][item]                       
2498                    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
2499                    X = np.linspace(0,90.0,26)
2500                    Y = G2lat.polfcal(ODFln,'0',X,0.0)
2501                    Plot.plot(X,Y,color='k',label=str(PH))
2502                    Plot.legend(loc='best')
2503                    Plot.set_title('Axial distribution for HKL='+str(PH)+' in '+phase)
2504                    Plot.set_xlabel(r'$\psi$',fontsize=16)
2505                    Plot.set_ylabel('MRD',fontsize=14)
2506    Page.canvas.draw()
2507   
2508################################################################################
2509##### PlotTexture
2510################################################################################
2511           
2512def PlotTexture(G2frame,data,Start=False):
2513    '''Pole figure, inverse pole figure, 3D pole distribution and 3D inverse pole distribution
2514    plotting.
2515    dict generalData contains all phase info needed which is in data
2516    '''
2517
2518    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2519    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2520    PatternId = G2frame.PatternId
2521    generalData = data['General']
2522    SGData = generalData['SGData']
2523    textureData = generalData['SH Texture']
2524    G2frame.G2plotNB.Delete('Texture')
2525    if not textureData['Order']:
2526        return                  #no plot!!
2527    SHData = generalData['SH Texture']
2528    SHCoef = SHData['SH Coeff'][1]
2529    cell = generalData['Cell'][1:7]
2530    Amat,Bmat = G2lat.cell2AB(cell)
2531    sq2 = 1.0/math.sqrt(2.0)
2532   
2533    def rp2xyz(r,p):
2534        z = npcosd(r)
2535        xy = np.sqrt(1.-z**2)
2536        return xy*npsind(p),xy*npcosd(p),z
2537           
2538    def OnMotion(event):
2539        SHData = data['General']['SH Texture']
2540        if event.xdata and event.ydata:                 #avoid out of frame errors
2541            xpos = event.xdata
2542            ypos = event.ydata
2543            if 'Inverse' in SHData['PlotType']:
2544                r = xpos**2+ypos**2
2545                if r <= 1.0:
2546                    if 'equal' in G2frame.Projection: 
2547                        r,p = 2.*npasind(np.sqrt(r)*sq2),npatan2d(ypos,xpos)
2548                    else:
2549                        r,p = 2.*npatand(np.sqrt(r)),npatan2d(ypos,xpos)
2550                    ipf = G2lat.invpolfcal(IODFln,SGData,np.array([r,]),np.array([p,]))
2551                    xyz = np.inner(Amat.T,np.array([rp2xyz(r,p)]))
2552                    y,x,z = list(xyz/np.max(np.abs(xyz)))
2553                   
2554                    G2frame.G2plotNB.status.SetFields(['',
2555                        'psi =%9.3f, beta =%9.3f, MRD =%9.3f hkl=%5.2f,%5.2f,%5.2f'%(r,p,ipf,x,y,z)])
2556                                   
2557            elif 'Axial' in SHData['PlotType']:
2558                pass
2559               
2560            else:                       #ordinary pole figure
2561                z = xpos**2+ypos**2
2562                if z <= 1.0:
2563                    z = np.sqrt(z)
2564                    if 'equal' in G2frame.Projection: 
2565                        r,p = 2.*npasind(z*sq2),npatan2d(ypos,xpos)
2566                    else:
2567                        r,p = 2.*npatand(z),npatan2d(ypos,xpos)
2568                    pf = G2lat.polfcal(ODFln,SamSym[textureData['Model']],np.array([r,]),np.array([p,]))
2569                    G2frame.G2plotNB.status.SetFields(['','phi =%9.3f, gam =%9.3f, MRD =%9.3f'%(r,p,pf)])
2570   
2571    try:
2572        plotNum = G2frame.G2plotNB.plotList.index('Texture')
2573        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2574        Page.figure.clf()
2575        Plot = Page.figure.gca()
2576        if not Page.IsShown():
2577            Page.Show()
2578    except ValueError:
2579        Plot = G2frame.G2plotNB.addMpl('Texture').gca()
2580        plotNum = G2frame.G2plotNB.plotList.index('Texture')
2581        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2582        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2583
2584    Page.Choice = None
2585    Page.SetFocus()
2586    G2frame.G2plotNB.status.SetFields(['',''])   
2587    PH = np.array(SHData['PFhkl'])
2588    phi,beta = G2lat.CrsAng(PH,cell,SGData)
2589    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
2590    PX = np.array(SHData['PFxyz'])
2591    gam = atan2d(PX[0],PX[1])
2592    xy = math.sqrt(PX[0]**2+PX[1]**2)
2593    xyz = math.sqrt(PX[0]**2+PX[1]**2+PX[2]**2)
2594    psi = asind(xy/xyz)
2595    IODFln = G2lat.Glnh(Start,SHCoef,psi,gam,SamSym[textureData['Model']])
2596    if 'Axial' in SHData['PlotType']:
2597        X = np.linspace(0,90.0,26)
2598        Y = G2lat.polfcal(ODFln,SamSym[textureData['Model']],X,0.0)
2599        Plot.plot(X,Y,color='k',label=str(SHData['PFhkl']))
2600        Plot.legend(loc='best')
2601        Plot.set_title('Axial distribution for HKL='+str(SHData['PFhkl']))
2602        Plot.set_xlabel(r'$\psi$',fontsize=16)
2603        Plot.set_ylabel('MRD',fontsize=14)
2604       
2605    else:       
2606        npts = 201
2607        if 'Inverse' in SHData['PlotType']:
2608            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
2609            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
2610            if 'equal' in G2frame.Projection:
2611                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
2612            else:
2613                R = np.where(R <= 1.,2.*npatand(R),0.0)
2614            Z = np.zeros_like(R)
2615            Z = G2lat.invpolfcal(IODFln,SGData,R,P)
2616            Z = np.reshape(Z,(npts,npts))
2617            CS = Plot.contour(Y,X,Z,aspect='equal')
2618            Plot.clabel(CS,fontsize=9,inline=1)
2619            try:
2620                Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
2621            except ValueError:
2622                pass
2623            Page.figure.colorbar(Img)
2624            Plot.set_title('Inverse pole figure for XYZ='+str(SHData['PFxyz']))
2625            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
2626                       
2627        else:
2628            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
2629            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
2630            if 'equal' in G2frame.Projection:
2631                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
2632            else:
2633                R = np.where(R <= 1.,2.*npatand(R),0.0)
2634            Z = np.zeros_like(R)
2635            Z = G2lat.polfcal(ODFln,SamSym[textureData['Model']],R,P)
2636            Z = np.reshape(Z,(npts,npts))
2637            try:
2638                CS = Plot.contour(Y,X,Z,aspect='equal')
2639                Plot.clabel(CS,fontsize=9,inline=1)
2640            except ValueError:
2641                pass
2642            Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
2643            Page.figure.colorbar(Img)
2644            Plot.set_title('Pole figure for HKL='+str(SHData['PFhkl']))
2645            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
2646    Page.canvas.draw()
2647
2648################################################################################
2649##### PlotCovariance
2650################################################################################
2651           
2652def PlotCovariance(G2frame,Data):
2653    'needs a doc string'
2654    if not Data:
2655        print 'No covariance matrix available'
2656        return
2657    varyList = Data['varyList']
2658    values = Data['variables']
2659    Xmax = len(varyList)
2660    covMatrix = Data['covMatrix']
2661    sig = np.sqrt(np.diag(covMatrix))
2662    xvar = np.outer(sig,np.ones_like(sig))
2663    covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
2664    title = ' for\n'+Data['title']
2665    newAtomDict = Data.get('newAtomDict',{})
2666   
2667
2668    def OnPlotKeyPress(event):
2669        newPlot = False
2670        if event.key == 's':
2671            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
2672            choice.sort()
2673            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
2674            if dlg.ShowModal() == wx.ID_OK:
2675                sel = dlg.GetSelection()
2676                G2frame.VcovColor = choice[sel]
2677            else:
2678                G2frame.VcovColor = 'RdYlGn'
2679            dlg.Destroy()
2680        PlotCovariance(G2frame,Data)
2681
2682    def OnMotion(event):
2683        #there is a problem here - reports wrong values
2684        if event.button:
2685            ytics = imgAx.get_yticks()
2686            ytics = np.where(ytics<len(varyList),ytics,-1)
2687            ylabs = [np.where(0<=i ,varyList[int(i)],' ') for i in ytics]
2688            imgAx.set_yticklabels(ylabs)           
2689        if event.xdata and event.ydata:                 #avoid out of frame errors
2690            xpos = int(event.xdata+.5)
2691            ypos = int(event.ydata+.5)
2692            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
2693                if xpos == ypos:
2694                    value = values[xpos]
2695                    name = varyList[xpos]
2696                    if varyList[xpos] in newAtomDict:
2697                        name,value = newAtomDict[name]                       
2698                    msg = '%s value = %.4g, esd = %.4g'%(name,value,sig[xpos])
2699                else:
2700                    msg = '%s - %s: %5.3f'%(varyList[xpos],varyList[ypos],covArray[xpos][ypos])
2701                Page.canvas.SetToolTipString(msg)
2702                G2frame.G2plotNB.status.SetFields(['',msg])
2703               
2704    try:
2705        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
2706        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2707        Page.figure.clf()
2708        Plot = Page.figure.gca()
2709        if not Page.IsShown():
2710            Page.Show()
2711    except ValueError:
2712        Plot = G2frame.G2plotNB.addMpl('Covariance').gca()
2713        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
2714        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2715        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2716        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
2717    Page.Choice = ['s: to change colors']
2718    Page.keyPress = OnPlotKeyPress
2719    Page.SetFocus()
2720    G2frame.G2plotNB.status.SetFields(['',''])   
2721    acolor = mpl.cm.get_cmap(G2frame.VcovColor)
2722    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower',
2723        vmin=-1.,vmax=1.)
2724    imgAx = Img.get_axes()
2725    ytics = imgAx.get_yticks()
2726    ylabs = [varyList[int(i)] for i in ytics[:-1]]
2727    imgAx.set_yticklabels(ylabs)
2728    colorBar = Page.figure.colorbar(Img)
2729    Plot.set_title('V-Cov matrix'+title)
2730    Plot.set_xlabel('Variable number')
2731    Plot.set_ylabel('Variable name')
2732    Page.canvas.draw()
2733   
2734################################################################################
2735##### PlotTorsion
2736################################################################################
2737
2738def PlotTorsion(G2frame,phaseName,Torsion,TorName,Names=[],Angles=[],Coeff=[]):
2739    'needs a doc string'
2740   
2741    global names
2742    names = Names
2743    sum = np.sum(Torsion)
2744    torsion = np.log(2*Torsion+1.)/sum
2745    tMin = np.min(torsion)
2746    tMax = np.max(torsion)
2747    torsion = 3.*(torsion-tMin)/(tMax-tMin)
2748    X = np.linspace(0.,360.,num=45)
2749   
2750    def OnPick(event):
2751        ind = event.ind[0]
2752        msg = 'atoms:'+names[ind]
2753        Page.canvas.SetToolTipString(msg)
2754        try:
2755            page = G2frame.dataDisplay.GetSelection()
2756        except:
2757            return
2758        if G2frame.dataDisplay.GetPageText(page) == 'Torsion restraints':
2759            torGrid = G2frame.dataDisplay.GetPage(page).Torsions
2760            torGrid.ClearSelection()
2761            for row in range(torGrid.GetNumberRows()):
2762                if names[ind] in torGrid.GetCellValue(row,0):
2763                    torGrid.SelectRow(row)
2764            torGrid.ForceRefresh()
2765               
2766    def OnMotion(event):
2767        if event.xdata and event.ydata:                 #avoid out of frame errors
2768            xpos = event.xdata
2769            ypos = event.ydata
2770            msg = 'torsion,energy: %5.3f %5.3f'%(xpos,ypos)
2771            Page.canvas.SetToolTipString(msg)
2772
2773    try:
2774        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
2775        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2776        Page.figure.clf()
2777        Plot = Page.figure.gca()
2778        if not Page.IsShown():
2779            Page.Show()
2780    except ValueError:
2781        Plot = G2frame.G2plotNB.addMpl('Torsion').gca()
2782        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
2783        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2784        Page.canvas.mpl_connect('pick_event', OnPick)
2785        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2786   
2787    Page.SetFocus()
2788    G2frame.G2plotNB.status.SetFields(['','Use mouse LB to identify torsion atoms'])
2789    Plot.plot(X,torsion,'b+')
2790    if len(Coeff):
2791        X2 = np.linspace(0.,360.,45)
2792        Y2 = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in X2])
2793        Plot.plot(X2,Y2,'r')
2794    if len(Angles):
2795        Eval = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in Angles])
2796        Plot.plot(Angles,Eval,'ro',picker=5)
2797    Plot.set_xlim((0.,360.))
2798    Plot.set_title('Torsion angles for '+TorName+' in '+phaseName)
2799    Plot.set_xlabel('angle',fontsize=16)
2800    Plot.set_ylabel('Energy',fontsize=16)
2801    Page.canvas.draw()
2802   
2803################################################################################
2804##### PlotRama
2805################################################################################
2806
2807def PlotRama(G2frame,phaseName,Rama,RamaName,Names=[],PhiPsi=[],Coeff=[]):
2808    'needs a doc string'
2809
2810    global names
2811    names = Names
2812    rama = np.log(2*Rama+1.)
2813    ramaMax = np.max(rama)
2814    rama = np.reshape(rama,(45,45))
2815    global Phi,Psi
2816    Phi = []
2817    Psi = []
2818
2819    def OnPlotKeyPress(event):
2820        newPlot = False
2821        if event.key == 's':
2822            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
2823            choice.sort()
2824            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
2825            if dlg.ShowModal() == wx.ID_OK:
2826                sel = dlg.GetSelection()
2827                G2frame.RamaColor = choice[sel]
2828            else:
2829                G2frame.RamaColor = 'RdYlGn'
2830            dlg.Destroy()
2831        PlotRama(G2frame,phaseName,Rama)
2832
2833    def OnPick(event):
2834        ind = event.ind[0]
2835        msg = 'atoms:'+names[ind]
2836        Page.canvas.SetToolTipString(msg)
2837        try:
2838            page = G2frame.dataDisplay.GetSelection()
2839        except:
2840            return
2841        if G2frame.dataDisplay.GetPageText(page) == 'Ramachandran restraints':
2842            ramaGrid = G2frame.dataDisplay.GetPage(page).Ramas
2843            ramaGrid.ClearSelection()
2844            for row in range(ramaGrid.GetNumberRows()):
2845                if names[ind] in ramaGrid.GetCellValue(row,0):
2846                    ramaGrid.SelectRow(row)
2847            ramaGrid.ForceRefresh()
2848
2849    def OnMotion(event):
2850        if event.xdata and event.ydata:                 #avoid out of frame errors
2851            xpos = event.xdata
2852            ypos = event.ydata
2853            msg = 'phi/psi: %5.3f %5.3f'%(xpos,ypos)
2854            Page.canvas.SetToolTipString(msg)
2855           
2856    try:
2857        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
2858        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2859        Page.figure.clf()
2860        Plot = Page.figure.gca()
2861        if not Page.IsShown():
2862            Page.Show()
2863    except ValueError:
2864        Plot = G2frame.G2plotNB.addMpl('Ramachandran').gca()
2865        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
2866        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2867        Page.canvas.mpl_connect('pick_event', OnPick)
2868        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2869        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
2870
2871    Page.Choice = ['s: to change colors']
2872    Page.keyPress = OnPlotKeyPress
2873    Page.SetFocus()
2874    G2frame.G2plotNB.status.SetFields(['','Use mouse LB to identify phi/psi atoms'])
2875    acolor = mpl.cm.get_cmap(G2frame.RamaColor)
2876    if RamaName == 'All' or '-1' in RamaName:
2877        if len(Coeff): 
2878            X,Y = np.meshgrid(np.linspace(-180.,180.,45),np.linspace(-180.,180.,45))
2879            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
2880            Plot.contour(X,Y,np.reshape(Z,(45,45)))
2881        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
2882            extent=[-180,180,-180,180],origin='lower')
2883        if len(PhiPsi):
2884            Phi,Psi = PhiPsi.T
2885            Phi = np.where(Phi>180.,Phi-360.,Phi)
2886            Psi = np.where(Psi>180.,Psi-360.,Psi)
2887            Plot.plot(Phi,Psi,'ro',picker=5)
2888        Plot.set_xlim((-180.,180.))
2889        Plot.set_ylim((-180.,180.))
2890    else:
2891        if len(Coeff): 
2892            X,Y = np.meshgrid(np.linspace(0.,360.,45),np.linspace(0.,360.,45))
2893            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
2894            Plot.contour(X,Y,np.reshape(Z,(45,45)))
2895        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
2896            extent=[0,360,0,360],origin='lower')
2897        if len(PhiPsi):
2898            Phi,Psi = PhiPsi.T
2899            Plot.plot(Phi,Psi,'ro',picker=5)
2900        Plot.set_xlim((0.,360.))
2901        Plot.set_ylim((0.,360.))
2902    Plot.set_title('Ramachandran for '+RamaName+' in '+phaseName)
2903    Plot.set_xlabel(r'$\phi$',fontsize=16)
2904    Plot.set_ylabel(r'$\psi$',fontsize=16)
2905    colorBar = Page.figure.colorbar(Img)
2906    Page.canvas.draw()
2907
2908
2909################################################################################
2910##### PlotSeq
2911################################################################################
2912def PlotSelectedSequence(G2frame,ColumnList,TableGet,SelectX,fitnum=None,fitvals=None):
2913    '''Plot a result from a sequential refinement
2914
2915    :param wx.Frame G2frame: The main GSAS-II tree "window"
2916    :param list ColumnList: list of int values corresponding to columns
2917      selected as y values
2918    :param function TableGet: a function that takes a column number
2919      as argument and returns the column label, the values and there ESDs (or None)
2920    :param function SelectX: a function that returns a selected column
2921      number (or None) as the X-axis selection
2922    '''
2923    global Title,xLabel,yLabel
2924    xLabel = yLabel = Title = ''
2925    def OnMotion(event):
2926        if event.xdata and event.ydata:                 #avoid out of frame errors
2927            xpos = event.xdata
2928            ypos = event.ydata
2929            msg = '%5.3f %.6g'%(xpos,ypos)
2930            Page.canvas.SetToolTipString(msg)
2931
2932    def OnKeyPress(event):
2933        global Title,xLabel,yLabel
2934        if event.key == 's':
2935            G2frame.seqXaxis = G2frame.seqXselect()
2936            Draw()
2937        elif event.key == 't':
2938            dlg = G2gd.MultiStringDialog(G2frame,'Set titles & labels',[' Title ',' x-Label ',' y-Label '],
2939                [Title,xLabel,yLabel])
2940            if dlg.Show():
2941                Title,xLabel,yLabel = dlg.GetValues()
2942            dlg.Destroy()
2943            Draw()
2944           
2945    def Draw():
2946        global Title,xLabel,yLabel
2947        Page.SetFocus()
2948        G2frame.G2plotNB.status.SetStatusText('press s to select X axis, t to change titles',1)
2949        Plot.clear()
2950        if G2frame.seqXaxis is not None:   
2951            xName,X,Xsig = Page.seqTableGet(G2frame.seqXaxis)
2952        else:
2953            X = np.arange(0,G2frame.SeqTable.GetNumberRows(),1)
2954            xName = 'Data sequence number'
2955        for col in Page.seqYaxisList:
2956            name,Y,sig = Page.seqTableGet(col)
2957            if sig:
2958                if G2frame.seqReverse and not G2frame.seqXaxis:
2959                    Y = Y[::-1]
2960                    sig = sig[::-1]
2961                Plot.errorbar(X,Y,yerr=sig,label=name)
2962            else:
2963                if G2frame.seqReverse and not G2frame.seqXaxis:
2964                    Y = Y[::-1]
2965                Plot.plot(X,Y)
2966                Plot.plot(X,Y,'o',label=name)
2967        if Page.fitvals:
2968            if G2frame.seqReverse and not G2frame.seqXaxis:
2969                Page.fitvals = Page.fitvals[::-1]
2970            Plot.plot(X,Page.fitvals,label='Fit')
2971           
2972        Plot.legend(loc='best')
2973        print Title,xLabel,yLabel
2974        if Title:
2975            Plot.set_title(Title)
2976        else:
2977            Plot.set_title('')
2978        if xLabel:
2979            Plot.set_xlabel(xLabel)
2980        else:
2981            Plot.set_xlabel(xName)
2982        if yLabel:
2983            Plot.set_ylabel(yLabel)
2984        else:
2985            Plot.set_ylabel('Parameter values')
2986        Page.canvas.draw()           
2987           
2988    G2frame.seqXselect = SelectX
2989    try:
2990        G2frame.seqXaxis
2991    except:
2992        G2frame.seqXaxis = None
2993
2994    if fitnum is None:
2995        label = 'Sequential refinement'
2996    else:
2997        label = 'Parametric fit #'+str(fitnum+1)
2998    try:
2999        plotNum = G2frame.G2plotNB.plotList.index(label)
3000        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3001        Page.figure.clf()
3002        Plot = Page.figure.gca()
3003        if not Page.IsShown():
3004            Page.Show()
3005    except ValueError:
3006        Plot = G2frame.G2plotNB.addMpl(label).gca()
3007        plotNum = G2frame.G2plotNB.plotList.index(label)
3008        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3009        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
3010        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3011    Page.Choice = ['s - select x-axis','t - change titles',]
3012    Page.keyPress = OnKeyPress
3013    Page.seqYaxisList = ColumnList
3014    Page.seqTableGet = TableGet
3015    Page.fitvals = fitvals
3016       
3017    Draw()
3018    G2frame.G2plotNB.nb.SetSelection(plotNum) # raises plot tab
3019               
3020################################################################################
3021##### PlotExposedImage & PlotImage
3022################################################################################
3023           
3024def PlotExposedImage(G2frame,newPlot=False,event=None):
3025    '''General access module for 2D image plotting
3026    '''
3027    plotNo = G2frame.G2plotNB.nb.GetSelection()
3028    if G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
3029        PlotImage(G2frame,newPlot,event,newImage=True)
3030    elif G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
3031        PlotIntegration(G2frame,newPlot,event)
3032
3033def OnStartMask(G2frame):
3034    '''Initiate the start of a Frame or Polygon map
3035
3036    :param wx.Frame G2frame: The main GSAS-II tree "window"
3037    :param str eventkey: a single letter ('f' or 'p') that
3038      determines what type of mask is created.   
3039    '''
3040    Masks = G2frame.PatternTree.GetItemPyData(
3041        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
3042    if G2frame.MaskKey == 'f':
3043        Masks['Frames'] = []
3044    elif G2frame.MaskKey == 'p':
3045        Masks['Polygons'].append([])
3046    elif G2frame.MaskKey == 's':
3047        Masks['Points'].append([])
3048    elif G2frame.MaskKey == 'a':
3049        Masks['Arcs'].append([])
3050    elif G2frame.MaskKey == 'r':
3051        Masks['Rings'].append([])
3052    G2imG.UpdateMasks(G2frame,Masks)
3053    PlotImage(G2frame,newImage=True)
3054   
3055def OnStartNewDzero(G2frame):
3056    '''Initiate the start of adding a new d-zero to a strain data set
3057
3058    :param wx.Frame G2frame: The main GSAS-II tree "window"
3059    :param str eventkey: a single letter ('a') that
3060      triggers the addition of a d-zero.   
3061    '''
3062    G2frame.dataFrame.GetStatusBar().SetStatusText('Add strain ring active - LB pick d-zero value',0)
3063    G2frame.PickId = G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain')
3064    data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
3065    return data
3066
3067def PlotImage(G2frame,newPlot=False,event=None,newImage=True):
3068    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
3069    masks, etc.
3070    '''
3071    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
3072    import numpy.ma as ma
3073    Dsp = lambda tth,wave: wave/(2.*npsind(tth/2.))
3074    global Data,Masks,StrSta
3075    colors=['b','g','r','c','m','k']
3076    Data = G2frame.PatternTree.GetItemPyData(
3077        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3078# patch
3079    if 'invert_x' not in Data:
3080        Data['invert_x'] = False
3081        Data['invert_y'] = True
3082# end patch
3083    Masks = G2frame.PatternTree.GetItemPyData(
3084        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
3085    try:    #may be absent
3086        StrSta = G2frame.PatternTree.GetItemPyData(
3087            G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain'))
3088    except TypeError:   #is missing
3089        StrSta = {}
3090
3091    def OnImMotion(event):
3092        Page.canvas.SetToolTipString('')
3093        sizexy = Data['size']
3094        if event.xdata and event.ydata and len(G2frame.ImageZ):                 #avoid out of frame errors
3095            Page.canvas.SetToolTipString('%8.2f %8.2fmm'%(event.xdata,event.ydata))
3096            Page.canvas.SetCursor(wx.CROSS_CURSOR)
3097            item = G2frame.itemPicked
3098            pixelSize = Data['pixelSize']
3099            scalex = 1000./pixelSize[0]
3100            scaley = 1000./pixelSize[1]
3101            if item and G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Image Controls':
3102                if 'Text' in str(item):
3103                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
3104                else:
3105                    xcent,ycent = Data['center']
3106                    xpos = event.xdata-xcent
3107                    ypos = event.ydata-ycent
3108                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
3109                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
3110                        Page.canvas.SetToolTipString('%6d deg'%(azm))
3111                    elif 'line1' in  str(item) or 'line2' in str(item):
3112                        Page.canvas.SetToolTipString('%8.3f deg'%(tth))                           
3113            else:
3114                xpos = event.xdata
3115                ypos = event.ydata
3116                xpix = xpos*scalex
3117                ypix = ypos*scaley
3118                Int = 0
3119                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
3120                    Int = G2frame.ImageZ[ypix][xpix]
3121                tth,azm,D,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
3122                Q = 2.*math.pi/dsp
3123                fields = ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)]
3124                if G2frame.MaskKey in ['p','f']:
3125                    fields[1] = 'Polygon/frame mask pick - LB next point, RB close polygon'
3126                elif G2frame.StrainKey:
3127                    fields[0] = 'd-zero pick active'
3128                G2frame.G2plotNB.status.SetFields(fields)
3129
3130    def OnImPlotKeyPress(event):
3131        try:
3132            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
3133        except TypeError:
3134            return
3135        if PickName == 'Masks':
3136            if event.key in ['l','p','f','s','a','r']:
3137                G2frame.MaskKey = event.key
3138                OnStartMask(G2frame)
3139                PlotImage(G2frame,newPlot=False)
3140               
3141        elif PickName == 'Stress/Strain':
3142            if event.key in ['a',]:
3143                G2frame.StrainKey = event.key
3144                StrSta = OnStartNewDzero(G2frame)
3145                PlotImage(G2frame,newPlot=False)
3146               
3147        elif PickName == 'Image Controls':
3148            if event.key in ['c',]:
3149                Xpos = event.xdata
3150                if not Xpos:            #got point out of frame
3151                    return
3152                Ypos = event.ydata
3153                dlg = wx.MessageDialog(G2frame,'Are you sure you want to change the center?',
3154                    'Center change',style=wx.OK|wx.CANCEL)
3155                try:
3156                    if dlg.ShowModal() == wx.ID_OK:
3157                        print 'move center to: ',Xpos,Ypos
3158                        Data['center'] = [Xpos,Ypos]
3159                        G2imG.UpdateImageControls(G2frame,Data,Masks)
3160                        PlotImage(G2frame,newPlot=False)
3161                finally:
3162                    dlg.Destroy()
3163                return
3164            elif event.key == 'l':
3165                G2frame.logPlot = not G2frame.logPlot
3166            elif event.key in ['x',]:
3167                Data['invert_x'] = not Data['invert_x']
3168            elif event.key in ['y',]:
3169                Data['invert_y'] = not Data['invert_y']
3170            PlotImage(G2frame,newPlot=True)
3171           
3172    def OnKeyBox(event):
3173        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):
3174            event.key = cb.GetValue()[0]
3175            cb.SetValue(' key press')
3176            if event.key in ['l','s','a','r','p','x','y']:
3177                wx.CallAfter(OnImPlotKeyPress,event)
3178        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
3179                       
3180    def OnImPick(event):
3181        if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']:
3182            return
3183        if G2frame.itemPicked is not None: return
3184        G2frame.itemPicked = event.artist
3185        G2frame.mousePicked = event.mouseevent
3186       
3187    def OnImRelease(event):
3188        try:
3189            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
3190        except TypeError:
3191            return
3192        if PickName not in ['Image Controls','Masks','Stress/Strain']:
3193            return
3194        pixelSize = Data['pixelSize']
3195        scalex = 1000./pixelSize[0]
3196        scaley = 1000./pixelSize[1]
3197#        pixLimit = Data['pixLimit']    #can be too tight
3198        pixLimit = 20       #this makes the search box 40x40 pixels
3199        if G2frame.itemPicked is None and PickName == 'Image Controls' and len(G2frame.ImageZ):
3200            Xpos = event.xdata
3201            if not (Xpos and G2frame.ifGetRing):                   #got point out of frame
3202                return
3203            Ypos = event.ydata
3204            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
3205                if event.button == 1:
3206                    Xpix = Xpos*scalex
3207                    Ypix = Ypos*scaley
3208                    xpos,ypos,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
3209                    if I and J:
3210                        xpos += .5                              #shift to pixel center
3211                        ypos += .5
3212                        xpos /= scalex                          #convert to mm
3213                        ypos /= scaley
3214                        Data['ring'].append([xpos,ypos])
3215                elif event.button == 3:
3216                    G2frame.dataFrame.GetStatusBar().SetStatusText('Calibrating...',0)
3217                    if G2img.ImageCalibrate(G2frame,Data):
3218                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check',0)
3219                        print 'Calibration successful'
3220                    else:
3221                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose',0)
3222                        print 'Calibration failed'
3223                    G2frame.ifGetRing = False
3224                    G2imG.UpdateImageControls(G2frame,Data,Masks)
3225                    return
3226                PlotImage(G2frame,newImage=False)
3227            return
3228        elif G2frame.MaskKey and PickName == 'Masks':
3229            Xpos,Ypos = [event.xdata,event.ydata]
3230            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3231                return
3232            if G2frame.MaskKey == 's' and event.button == 1:
3233                Masks['Points'][-1] = [Xpos,Ypos,1.]
3234                G2frame.MaskKey = ''               
3235            elif G2frame.MaskKey == 'r' and event.button == 1:
3236                tth = G2img.GetTth(Xpos,Ypos,Data)
3237                Masks['Rings'][-1] = [tth,0.1]
3238                G2frame.MaskKey = ''               
3239            elif G2frame.MaskKey == 'a' and event.button == 1:
3240                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
3241                azm = int(azm)               
3242                Masks['Arcs'][-1] = [tth,[azm-5,azm+5],0.1]
3243                G2frame.MaskKey = ''               
3244            elif G2frame.MaskKey =='p':
3245                polygon = Masks['Polygons'][-1]
3246                if len(polygon) > 2 and event.button == 3:
3247                    x0,y0 = polygon[0]
3248                    polygon.append([x0,y0])
3249                    G2frame.MaskKey = ''
3250                    G2frame.G2plotNB.status.SetFields(['','Polygon closed'])
3251                else:
3252                    G2frame.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(Xpos,Ypos)])
3253                    polygon.append([Xpos,Ypos])
3254            elif G2frame.MaskKey =='f':
3255                frame = Masks['Frames']
3256                if len(frame) > 2 and event.button == 3:
3257                    x0,y0 = frame[0]
3258                    frame.append([x0,y0])
3259                    G2frame.MaskKey = ''
3260                    G2frame.G2plotNB.status.SetFields(['','Frame closed'])
3261                else:
3262                    G2frame.G2plotNB.status.SetFields(['','New frame point: %.1f,%.1f'%(Xpos,Ypos)])
3263                    frame.append([Xpos,Ypos])
3264            G2imG.UpdateMasks(G2frame,Masks)
3265            PlotImage(G2frame,newImage=False)
3266        elif PickName == 'Stress/Strain' and G2frame.StrainKey:
3267            Xpos,Ypos = [event.xdata,event.ydata]
3268            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3269                return
3270            dsp = G2img.GetDsp(Xpos,Ypos,Data)
3271            StrSta['d-zero'].append({'Dset':dsp,'Dcalc':0.0,'pixLimit':10,'cutoff':1.0,
3272                'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]],'Emat':[1.0,1.0,1.0]})
3273            R,r = G2img.MakeStrStaRing(StrSta['d-zero'][-1],G2frame.ImageZ,Data)
3274            if not len(R):
3275                del StrSta['d-zero'][-1]
3276                G2frame.ErrorDialog('Strain peak selection','WARNING - No points found for this ring selection')
3277            StrSta['d-zero'] = G2mth.sortArray(StrSta['d-zero'],'Dset',reverse=True)
3278            G2frame.StrainKey = ''
3279            G2imG.UpdateStressStrain(G2frame,StrSta)
3280            PlotStrain(G2frame,StrSta)
3281            PlotImage(G2frame,newPlot=False)           
3282        else:
3283            Xpos,Ypos = [event.xdata,event.ydata]
3284            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3285                return
3286            if G2frame.ifGetRing:                          #delete a calibration ring pick
3287                xypos = [Xpos,Ypos]
3288                rings = Data['ring']
3289                for ring in rings:
3290                    if np.allclose(ring,xypos,.01,0):
3291                        rings.remove(ring)
3292            else:
3293                tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data)[:3]
3294                itemPicked = str(G2frame.itemPicked)
3295                if 'Line2D' in itemPicked and PickName == 'Image Controls':
3296                    if 'line1' in itemPicked:
3297                        Data['IOtth'][0] = max(tth,0.001)
3298                    elif 'line2' in itemPicked:
3299                        Data['IOtth'][1] = tth
3300                    elif 'line3' in itemPicked:
3301                        Data['LRazimuth'][0] = int(azm)
3302                    elif 'line4' in itemPicked and not Data['fullIntegrate']:
3303                        Data['LRazimuth'][1] = int(azm)
3304                   
3305                    Data['LRazimuth'][0] %= 360
3306                    Data['LRazimuth'][1] %= 360
3307                    if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
3308                        Data['LRazimuth'][1] += 360                       
3309                    if Data['fullIntegrate']:
3310                        Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
3311                       
3312                    if  Data['IOtth'][0] > Data['IOtth'][1]:
3313                        Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
3314                       
3315                    G2frame.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
3316                    G2frame.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
3317                    G2frame.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
3318                    G2frame.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
3319                elif 'Circle' in itemPicked and PickName == 'Masks':
3320                    spots = Masks['Points']
3321                    newPos = itemPicked.split(')')[0].split('(')[2].split(',')
3322                    newPos = np.array([float(newPos[0]),float(newPos[1])])
3323                    for spot in spots:
3324                        if np.allclose(np.array([spot[:2]]),newPos):
3325                            spot[:2] = Xpos,Ypos
3326                    G2imG.UpdateMasks(G2frame,Masks)
3327                elif 'Line2D' in itemPicked and PickName == 'Masks':
3328                    Obj = G2frame.itemPicked.findobj()
3329                    rings = Masks['Rings']
3330                    arcs = Masks['Arcs']
3331                    polygons = Masks['Polygons']
3332                    frame = Masks['Frames']
3333                    for ring in G2frame.ringList:
3334                        if Obj == ring[0]:
3335                            rN = ring[1]
3336                            if ring[2] == 'o':
3337                                rings[rN][0] = G2img.GetTth(Xpos,Ypos,Data)-rings[rN][1]/2.
3338                            else:
3339                                rings[rN][0] = G2img.GetTth(Xpos,Ypos,Data)+rings[rN][1]/2.
3340                    for arc in G2frame.arcList:
3341                        if Obj == arc[0]:
3342                            aN = arc[1]
3343                            if arc[2] == 'o':
3344                                arcs[aN][0] = G2img.GetTth(Xpos,Ypos,Data)-arcs[aN][2]/2
3345                            elif arc[2] == 'i':
3346                                arcs[aN][0] = G2img.GetTth(Xpos,Ypos,Data)+arcs[aN][2]/2
3347                            elif arc[2] == 'l':
3348                                arcs[aN][1][0] = int(G2img.GetAzm(Xpos,Ypos,Data))
3349                            else:
3350                                arcs[aN][1][1] = int(G2img.GetAzm(Xpos,Ypos,Data))
3351                    for poly in G2frame.polyList:   #merging points problem here?
3352                        if Obj == poly[0]:
3353                            ind = G2frame.itemPicked.contains(G2frame.mousePicked)[1]['ind'][0]
3354                            oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
3355                            pN = poly[1]
3356                            for i,xy in enumerate(polygons[pN]):
3357                                if np.allclose(np.array([xy]),oldPos,atol=1.0):
3358                                    if event.button == 1:
3359                                        polygons[pN][i] = Xpos,Ypos
3360                                    elif event.button == 3:
3361                                        polygons[pN].insert(i,[Xpos,Ypos])
3362                                        break
3363                    if frame:
3364                        oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
3365                        for i,xy in enumerate(frame):
3366                            if np.allclose(np.array([xy]),oldPos,atol=1.0):
3367                                if event.button == 1:
3368                                    frame[i] = Xpos,Ypos
3369                                elif event.button == 3:
3370                                    frame.insert(i,[Xpos,Ypos])
3371                                    break
3372                    G2imG.UpdateMasks(G2frame,Masks)
3373#                else:                  #keep for future debugging
3374#                    print str(G2frame.itemPicked),event.xdata,event.ydata,event.button
3375            PlotImage(G2frame,newImage=True)
3376            G2frame.itemPicked = None
3377           
3378    try:
3379        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
3380        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3381        if not newPlot:
3382            Plot = Page.figure.gca()          #get previous powder plot & get limits
3383            xylim = Plot.get_xlim(),Plot.get_ylim()
3384        if newImage:
3385            Page.figure.clf()
3386            Plot = Page.figure.gca()          #get a fresh plot after clf()
3387    except ValueError:
3388        Plot = G2frame.G2plotNB.addMpl('2D Powder Image').gca()
3389        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
3390        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3391        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
3392        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
3393        Page.canvas.mpl_connect('pick_event', OnImPick)
3394        Page.canvas.mpl_connect('button_release_event', OnImRelease)
3395        xylim = []
3396    Page.Choice = None
3397    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
3398        Page.SetFocus()
3399    Title = G2frame.PatternTree.GetItemText(G2frame.Image)[4:]
3400    G2frame.G2plotNB.status.DestroyChildren()
3401    if G2frame.logPlot:
3402        Title = 'log('+Title+')'
3403    Plot.set_title(Title)
3404    try:
3405        if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Image Controls',]:
3406            Page.Choice = (' key press','l: log(I) on','x: flip x','y: flip y',)
3407            if G2frame.logPlot:
3408                Page.Choice[1] = 'l: log(I) off'
3409            Page.keyPress = OnImPlotKeyPress
3410        elif G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Masks',]:
3411            Page.Choice = (' key press','l: log(I) on','s: spot mask','a: arc mask','r: ring mask',
3412                'p: polygon mask','f: frame mask',)
3413            if G2frame.logPlot:
3414                Page.Choice[1] = 'l: log(I) off'
3415            Page.keyPress = OnImPlotKeyPress
3416        elif G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Stress/Strain',]:
3417            Page.Choice = (' key press','a: add new ring',)
3418            Page.keyPress = OnImPlotKeyPress
3419    except TypeError:
3420        pass
3421    size,imagefile = G2frame.PatternTree.GetItemPyData(G2frame.Image)
3422    dark = Data.get('dark image',[0,''])
3423    if dark[0]:
3424        darkfile = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, 
3425            G2frame.root,dark[0]))[1]
3426    if imagefile != G2frame.oldImagefile:
3427        imagefile = G2IO.CheckImageFile(G2frame,imagefile)
3428        if not imagefile:
3429            G2frame.G2plotNB.Delete('2D Powder Image')
3430            return
3431        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
3432        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
3433        if dark[0]:
3434            darkImg = G2IO.GetImageData(G2frame,darkfile,imageOnly=True)
3435            G2frame.ImageZ += dark[1]*darkImg
3436        G2frame.oldImagefile = imagefile
3437
3438    imScale = 1
3439    if len(G2frame.ImageZ) > 1024:
3440        imScale = len(G2frame.ImageZ)/1024
3441    sizexy = Data['size']
3442    pixelSize = Data['pixelSize']
3443    scalex = 1000./pixelSize[0]
3444    scaley = 1000./pixelSize[1]
3445    Xmax = sizexy[0]*pixelSize[0]/1000.
3446    Ymax = sizexy[1]*pixelSize[1]/1000.
3447    xlim = (0,Xmax)
3448    ylim = (Ymax,0)
3449    Imin,Imax = Data['range'][1]
3450    acolor = mpl.cm.get_cmap(Data['color'])
3451    xcent,ycent = Data['center']
3452    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
3453    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
3454    #do threshold mask - "real" mask - others are just bondaries
3455    Zlim = Masks['Thresholds'][1]
3456    wx.BeginBusyCursor()
3457    try:
3458           
3459        if newImage:                   
3460            Imin,Imax = Data['range'][1]
3461            MA = ma.masked_greater(ma.masked_less(G2frame.ImageZ,Zlim[0]),Zlim[1])
3462            MaskA = ma.getmaskarray(MA)
3463            A = G2img.ImageCompress(MA,imScale)
3464            AM = G2img.ImageCompress(MaskA,imScale)
3465            if G2frame.logPlot:
3466                A = np.where(A>Imin,np.where(A<Imax,A,0),0)
3467                A = np.where(A>0,np.log(A),0)
3468                AM = np.where(AM>0,np.log(AM),0)
3469                Imin,Imax = [np.amin(A),np.amax(A)]
3470            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
3471                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Ymax,0])
3472            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
3473                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
3474   
3475        Plot.plot(xcent,ycent,'x')
3476        #G2frame.PatternTree.GetItemText(item)
3477        if Data['showLines']:
3478            LRAzim = Data['LRazimuth']                  #NB: integers
3479            Nazm = Data['outAzimuths']
3480            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
3481            AzmthOff = Data['azmthOff']
3482            IOtth = Data['IOtth']
3483            wave = Data['wavelength']
3484            dspI = wave/(2.0*sind(IOtth[0]/2.0))
3485            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
3486            dspO = wave/(2.0*sind(IOtth[1]/2.0))
3487            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
3488            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
3489            if ellI:
3490                xyI = []
3491                for azm in Azm:
3492                    xy = G2img.GetDetectorXY(dspI,azm,Data)
3493                    if np.any(xy):
3494                        xyI.append(xy)
3495                if len(xyI):
3496                    xyI = np.array(xyI)
3497                    arcxI,arcyI = xyI.T
3498                    Plot.plot(arcxI,arcyI,picker=3)
3499            if ellO:
3500                xyO = []
3501                for azm in Azm:
3502                    xy = G2img.GetDetectorXY(dspO,azm,Data)
3503                    if np.any(xy):
3504                        xyO.append(xy)
3505                if len(xyO):
3506                    xyO = np.array(xyO)
3507                    arcxO,arcyO = xyO.T               
3508                    Plot.plot(arcxO,arcyO,picker=3)
3509            if ellO and ellI:
3510                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
3511                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
3512            for i in range(Nazm):
3513                cake = LRAzim[0]+i*delAzm-AzmthOff
3514                if Data['centerAzm']:
3515                    cake += delAzm/2.
3516                ind = np.searchsorted(Azm,cake)
3517                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
3518                   
3519        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Image Controls':
3520            for xring,yring in Data['ring']:
3521                Plot.plot(xring,yring,'r+',picker=3)
3522            if Data['setRings']:
3523                N = 0
3524                for ring in Data['rings']:
3525                    xring,yring = np.array(ring).T[:2]
3526                    Plot.plot(xring,yring,'.',color=colors[N%6])
3527                    N += 1           
3528            for ellipse in Data['ellipses']:      #what about hyperbola?
3529                cent,phi,[width,height],col = ellipse
3530                if width > 0:       #ellipses
3531                    Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
3532                    Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
3533        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Stress/Strain':
3534            for N,ring in enumerate(StrSta['d-zero']):
3535                xring,yring = ring['ImxyObs']
3536                Plot.plot(xring,yring,colors[N%6]+'.')
3537        #masks - mask lines numbered after integration limit lines
3538        spots = Masks['Points']
3539        rings = Masks['Rings']
3540        arcs = Masks['Arcs']
3541        polygons = Masks['Polygons']
3542        if 'Frames' not in Masks:
3543            Masks['Frames'] = []
3544        frame = Masks['Frames']
3545        for spot in spots:
3546            if spot:
3547                x,y,d = spot
3548                Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
3549        G2frame.ringList = []
3550        for iring,ring in enumerate(rings):
3551            if ring:
3552                tth,thick = ring
3553                wave = Data['wavelength']
3554                xy1 = []
3555                xy2 = []
3556                Azm = np.linspace(0,362,181)
3557                for azm in Azm:
3558                    xy1.append(G2img.GetDetectorXY(Dsp(tth+thick/2.,wave),azm,Data))      #what about hyperbola
3559                    xy2.append(G2img.GetDetectorXY(Dsp(tth-thick/2.,wave),azm,Data))      #what about hyperbola
3560                x1,y1 = np.array(xy1).T
3561                x2,y2 = np.array(xy2).T
3562                G2frame.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
3563                G2frame.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
3564        G2frame.arcList = []
3565        for iarc,arc in enumerate(arcs):
3566            if arc:
3567                tth,azm,thick = arc           
3568                wave = Data['wavelength']
3569                xy1 = []
3570                xy2 = []
3571                aR = azm[0],azm[1],azm[1]-azm[0]
3572                if azm[1]-azm[0] > 180:
3573                    aR[2] /= 2
3574                Azm = np.linspace(aR[0],aR[1],aR[2])
3575                for azm in Azm:
3576                    xy1.append(G2img.GetDetectorXY(Dsp(tth+thick/2.,wave),azm,Data))      #what about hyperbola
3577                    xy2.append(G2img.GetDetectorXY(Dsp(tth-thick/2.,wave),azm,Data))      #what about hyperbola
3578                x1,y1 = np.array(xy1).T
3579                x2,y2 = np.array(xy2).T
3580                G2frame.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
3581                G2frame.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
3582                G2frame.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
3583                G2frame.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
3584        G2frame.polyList = []
3585        for ipoly,polygon in enumerate(polygons):
3586            if polygon:
3587                x,y = np.hsplit(np.array(polygon),2)
3588                G2frame.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
3589                Plot.plot(x,y,'r')           
3590        G2frame.frameList = []
3591        if frame:
3592            x,y = np.hsplit(np.array(frame),2)
3593            G2frame.frameList.append([Plot.plot(x,y,'g+',picker=10),0])
3594            Plot.plot(x,y,'g')           
3595        if newImage:
3596            colorBar = Page.figure.colorbar(Img)
3597        Plot.set_xlim(xlim)
3598        Plot.set_ylim(ylim)
3599        if Data['invert_x']:
3600            Plot.invert_xaxis()
3601        if Data['invert_y']:
3602            Plot.invert_yaxis()
3603        if not newPlot and xylim:
3604            Page.toolbar.push_current()
3605            Plot.set_xlim(xylim[0])
3606            Plot.set_ylim(xylim[1])
3607            xylim = []
3608            Page.toolbar.push_current()
3609            Page.toolbar.draw()
3610            # patch for wx 2.9 on Mac, to force a redraw
3611            i,j= wx.__version__.split('.')[0:2]
3612            if int(i)+int(j)/10. > 2.8 and 'wxOSX' in wx.PlatformInfo:
3613                Page.canvas.draw()
3614        else:
3615            Page.canvas.draw()
3616    finally:
3617        wx.EndBusyCursor()
3618       
3619################################################################################
3620##### PlotIntegration
3621################################################################################
3622           
3623def PlotIntegration(G2frame,newPlot=False,event=None):
3624    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
3625    '''
3626           
3627    def OnMotion(event):
3628        Page.canvas.SetToolTipString('')
3629        Page.canvas.SetCursor(wx.CROSS_CURSOR)
3630        azm = event.ydata
3631        tth = event.xdata
3632        if azm and tth:
3633            G2frame.G2plotNB.status.SetFields(\
3634                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
3635                               
3636    try:
3637        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
3638        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3639        if not newPlot:
3640            Plot = Page.figure.gca()          #get previous plot & get limits
3641            xylim = Plot.get_xlim(),Plot.get_ylim()
3642        Page.figure.clf()
3643        Plot = Page.figure.gca()          #get a fresh plot after clf()
3644       
3645    except ValueError:
3646        Plot = G2frame.G2plotNB.addMpl('2D Integration').gca()
3647        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
3648        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3649        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3650        Page.views = False
3651        view = False
3652    Page.Choice = None
3653    if not event:
3654        Page.SetFocus()
3655       
3656    Data = G2frame.PatternTree.GetItemPyData(
3657        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3658    image = G2frame.Integrate[0]
3659    xsc = G2frame.Integrate[1]
3660    ysc = G2frame.Integrate[2]
3661    Imin,Imax = Data['range'][1]
3662    acolor = mpl.cm.get_cmap(Data['color'])
3663    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
3664    Plot.set_ylabel('azimuth',fontsize=12)
3665    Plot.set_xlabel('2-theta',fontsize=12)
3666    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
3667        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
3668    colorBar = Page.figure.colorbar(Img)
3669#    if Data['ellipses']:           
3670#        for ellipse in Data['ellipses']:
3671#            x,y = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
3672#            tth,azm = G2img.GetTthAzm(x,y,Data)
3673##            azm = np.where(azm < 0.,azm+360,azm)
3674#            Plot.plot(tth,azm,'b,')
3675    if not newPlot:
3676        Page.toolbar.push_current()
3677        Plot.set_xlim(xylim[0])
3678        Plot.set_ylim(xylim[1])
3679        xylim = []
3680        Page.toolbar.push_current()
3681        Page.toolbar.draw()
3682    else:
3683        Page.canvas.draw()
3684               
3685################################################################################
3686##### PlotTRImage
3687################################################################################
3688           
3689def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):
3690    '''a test plot routine - not normally used
3691    ''' 
3692           
3693    def OnMotion(event):
3694        Page.canvas.SetToolTipString('')
3695        Page.canvas.SetCursor(wx.CROSS_CURSOR)
3696        azm = event.xdata
3697        tth = event.ydata
3698        if azm and tth:
3699            G2frame.G2plotNB.status.SetFields(\
3700                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
3701                               
3702    try:
3703        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
3704        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3705        if not newPlot:
3706            Plot = Page.figure.gca()          #get previous plot & get limits
3707            xylim = Plot.get_xlim(),Plot.get_ylim()
3708        Page.figure.clf()
3709        Plot = Page.figure.gca()          #get a fresh plot after clf()
3710       
3711    except ValueError:
3712        Plot = G2frame.G2plotNB.addMpl('2D Transformed Powder Image').gca()
3713        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
3714        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3715        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3716        Page.views = False
3717        view = False
3718    Page.Choice = None
3719    Page.SetFocus()
3720       
3721    Data = G2frame.PatternTree.GetItemPyData(
3722        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3723    Imin,Imax = Data['range'][1]
3724    step = (Imax-Imin)/5.
3725    V = np.arange(Imin,Imax,step)
3726    acolor = mpl.cm.get_cmap(Data['color'])
3727    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
3728    Plot.set_xlabel('azimuth',fontsize=12)
3729    Plot.set_ylabel('2-theta',fontsize=12)
3730    Plot.contour(tax,tay,taz,V,cmap=acolor)
3731    if Data['showLines']:
3732        IOtth = Data['IOtth']
3733        if Data['fullIntegrate']:
3734            LRAzim = [-180,180]
3735        else:
3736            LRAzim = Data['LRazimuth']                  #NB: integers
3737        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
3738        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
3739        if not Data['fullIntegrate']:
3740            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
3741            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
3742    if Data['setRings']:
3743        rings = np.concatenate((Data['rings']),axis=0)
3744        for xring,yring,dsp in rings:
3745            x,y = G2img.GetTthAzm(xring,yring,Data)
3746            Plot.plot(y,x,'r+')           
3747    if Data['ellipses']:           
3748        for ellipse in Data['ellipses']:
3749            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
3750            x,y = np.hsplit(ring,2)
3751            tth,azm = G2img.GetTthAzm(x,y,Data)
3752            Plot.plot(azm,tth,'b,')
3753    if not newPlot:
3754        Page.toolbar.push_current()
3755        Plot.set_xlim(xylim[0])
3756        Plot.set_ylim(xylim[1])
3757        xylim = []
3758        Page.toolbar.push_current()
3759        Page.toolbar.draw()
3760    else:
3761        Page.canvas.draw()
3762       
3763################################################################################
3764##### PlotStructure
3765################################################################################
3766           
3767def PlotStructure(G2frame,data,firstCall=False):
3768    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
3769    thermal motion ellipsoids and polyhedra
3770    '''
3771
3772    def FindPeaksBonds(XYZ):
3773        rFact = data['Drawing'].get('radiusFactor',0.85)    #data['Drawing'] could be empty!
3774        Bonds = [[] for x in XYZ]
3775        for i,xyz in enumerate(XYZ):
3776            Dx = XYZ-xyz
3777            dist = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=1))
3778            IndB = ma.nonzero(ma.masked_greater(dist,rFact*2.2))
3779            for j in IndB[0]:
3780                Bonds[i].append(Dx[j]/2.)
3781                Bonds[j].append(-Dx[j]/2.)
3782        return Bonds
3783
3784    # PlotStructure initialization here
3785    ForthirdPI = 4.0*math.pi/3.0
3786    generalData = data['General']
3787    cell = generalData['Cell'][1:7]
3788    Vol = generalData['Cell'][7:8][0]
3789    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
3790    Gmat,gmat = G2lat.cell2Gmat(cell)
3791    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3792    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3793    SGData = generalData['SGData']
3794    Mydir = generalData['Mydir']
3795    atomData = data['Atoms']
3796    mapPeaks = []
3797    drawingData = data['Drawing']
3798    if not drawingData:
3799        return          #nothing setup, nothing to draw   
3800    if 'Map Peaks' in data:
3801        mapPeaks = np.array(data['Map Peaks'])
3802        peakMax = 100.
3803        if len(mapPeaks):
3804            peakMax = np.max(mapPeaks.T[0])
3805    resRBData = data['RBModels'].get('Residue',[])
3806    vecRBData = data['RBModels'].get('Vector',[])
3807    rbAtmDict = {}
3808    for rbObj in resRBData+vecRBData:
3809        exclList = ['X' for i in range(len(rbObj['Ids']))]
3810        rbAtmDict.update(dict(zip(rbObj['Ids'],exclList)))
3811    testRBObj = data.get('testRBObj',{})
3812    rbObj = testRBObj.get('rbObj',{})
3813    MCSA = data.get('MCSA',{})
3814    mcsaModels = MCSA.get('Models',[])
3815    if len(mcsaModels) > 1:
3816        XYZs,Types = G2mth.UpdateMCSAxyz(Bmat,MCSA)
3817        mcsaXYZ = []
3818        mcsaTypes = []
3819        neqv = 0
3820        for xyz,atyp in zip(XYZs,Types):
3821            equiv = G2spc.GenAtom(xyz,SGData,All=True,Move=False)
3822            neqv = max(neqv,len(equiv))
3823            for item in equiv:
3824                mcsaXYZ.append(item[0]) 
3825                mcsaTypes.append(atyp)
3826        mcsaXYZ = np.array(mcsaXYZ)
3827        mcsaTypes = np.array(mcsaTypes)
3828        nuniq = mcsaXYZ.shape[0]/neqv
3829        mcsaXYZ = np.reshape(mcsaXYZ,(nuniq,neqv,3))
3830        mcsaTypes = np.reshape(mcsaTypes,(nuniq,neqv))
3831        cent = np.fix(np.sum(mcsaXYZ+2.,axis=0)/nuniq)-2
3832        cent[0] = [0,0,0]   #make sure 1st one isn't moved
3833        mcsaXYZ = np.swapaxes(mcsaXYZ,0,1)-cent[:,np.newaxis,:]
3834        mcsaTypes = np.swapaxes(mcsaTypes,0,1)
3835        mcsaXYZ = np.reshape(mcsaXYZ,(nuniq*neqv,3))
3836        mcsaTypes = np.reshape(mcsaTypes,(nuniq*neqv))                       
3837        mcsaBonds = FindPeaksBonds(mcsaXYZ)       
3838    drawAtoms = drawingData.get('Atoms',[])
3839    mapData = {}
3840    flipData = {}
3841    rhoXYZ = []
3842    showBonds = False
3843    if 'Map' in generalData:
3844        mapData = generalData['Map']
3845        showBonds = mapData.get('Show bonds',False)
3846    if 'Flip' in generalData:
3847        flipData = generalData['Flip']                       
3848        flipData['mapRoll'] = [0,0,0]
3849    Wt = np.array([255,255,255])
3850    Rd = np.array([255,0,0])
3851    Gr = np.array([0,255,0])
3852    wxGreen = wx.Colour(0,255,0)
3853    Bl = np.array([0,0,255])
3854    Or = np.array([255,128,0])
3855    wxOrange = wx.Colour(255,128,0)
3856    uBox = np.array([[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]])
3857    uEdges = np.array([
3858        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
3859        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
3860        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
3861    mD = 0.1
3862    mV = np.array([[[-mD,0,0],[mD,0,0]],[[0,-mD,0],[0,mD,0]],[[0,0,-mD],[0,0,mD]]])
3863    mapPeakVecs = np.inner(mV,Bmat)
3864
3865    backColor = np.array(list(drawingData['backColor'])+[0,])
3866    Bc = np.array(list(drawingData['backColor']))
3867    uColors = [Rd,Gr,Bl,Wt-Bc, Wt-Bc,Wt-Bc,Wt-Bc,Wt-Bc, Wt-Bc,Wt-Bc,Wt-Bc,Wt-Bc]
3868    altDown = False
3869    shiftDown = False
3870    ctrlDown = False
3871   
3872    def OnKeyBox(event):
3873#        Draw()                          #make sure plot is fresh!!
3874        mode = cb.GetValue()
3875        if mode in ['jpeg','bmp','tiff',]:
3876            try:
3877                import Image as Im
3878            except ImportError:
3879                try:
3880                    from PIL import Image as Im
3881                except ImportError:
3882                    print "PIL/pillow Image module not present. Cannot save images without this"
3883                    raise Exception("PIL/pillow Image module not found")
3884            Fname = os.path.join(Mydir,generalData['Name']+'.'+mode)
3885            size = Page.canvas.GetSize()
3886            glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
3887            if mode in ['jpeg',]:
3888                Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
3889                im = Image.new("RGBA", (size[0],size[1]))
3890            else:
3891                Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
3892                im = Image.new("RGB", (size[0],size[1]))
3893            im.fromstring(Pix)
3894            im.save(Fname,mode)
3895            cb.SetValue(' save as/key:')
3896            G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
3897        else:
3898            event.key = cb.GetValue()[0]
3899            cb.SetValue(' save as/key:')
3900            wx.CallAfter(OnKey,event)
3901        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
3902
3903    def OnKey(event):           #on key UP!!
3904#        Draw()                          #make sure plot is fresh!!
3905        try:
3906            keyCode = event.GetKeyCode()
3907            if keyCode > 255:
3908                keyCode = 0
3909            key = chr(keyCode)
3910        except AttributeError:       #if from OnKeyBox above
3911            key = str(event.key).upper()
3912        indx = drawingData['selectedAtoms']
3913        cx,ct = drawingData['atomPtrs'][:2]
3914        if key in ['C']:
3915            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
3916            drawingData['viewDir'] = [0,0,1]
3917            drawingData['oldxy'] = []
3918            V0 = np.array([0,0,1])
3919            V = np.inner(Amat,V0)
3920            V /= np.sqrt(np.sum(V**2))
3921            A = np.arccos(np.sum(V*V0))
3922            Q = G2mth.AV2Q(A,[0,1,0])
3923            drawingData['Quaternion'] = Q
3924            SetViewPointText(drawingData['viewPoint'][0])
3925            SetViewDirText(drawingData['viewDir'])
3926            Q = drawingData['Quaternion']
3927            G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
3928        elif key in ['N']:
3929            drawAtoms = drawingData['Atoms']
3930            if not len(drawAtoms):      #no atoms
3931                return
3932            pI = drawingData['viewPoint'][1]
3933            if not len(pI):
3934                pI = [0,0]
3935            if indx:
3936                pI[0] = indx[pI[1]]
3937                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
3938                pI[1] += 1
3939                if pI[1] >= len(indx):
3940                    pI[1] = 0
3941            else:
3942                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
3943                pI[0] += 1
3944                if pI[0] >= len(drawAtoms):
3945                    pI[0] = 0
3946            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
3947            SetViewPointText(drawingData['viewPoint'][0])
3948            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
3949               
3950        elif key in ['P']:
3951            drawAtoms = drawingData['Atoms']
3952            if not len(drawAtoms):      #no atoms
3953                return
3954            pI = drawingData['viewPoint'][1]
3955            if not len(pI):
3956                pI = [0,0]
3957            if indx:
3958                pI[0] = indx[pI[1]]
3959                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
3960                pI[1] -= 1
3961                if pI[1] < 0:
3962                    pI[1] = len(indx)-1
3963            else:
3964                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
3965                pI[0] -= 1
3966                if pI[0] < 0:
3967                    pI[0] = len(drawAtoms)-1
3968            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
3969            SetViewPointText(drawingData['viewPoint'][0])           
3970            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
3971        elif key in ['U','D','L','R'] and mapData['Flip'] == True:
3972            dirDict = {'U':[0,1],'D':[0,-1],'L':[-1,0],'R':[1,0]}
3973            SetMapRoll(dirDict[key])
3974            SetPeakRoll(dirDict[key])
3975            SetMapPeaksText(mapPeaks)
3976        Draw('key')
3977           
3978    def GetTruePosition(xy,Add=False):
3979        View = glGetIntegerv(GL_VIEWPORT)
3980        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
3981        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
3982        Zmax = 1.
3983        if Add:
3984            Indx = GetSelectedAtoms()
3985        if G2frame.dataDisplay.GetPageText(getSelection()) == 'Map peaks':
3986            for i,peak in enumerate(mapPeaks):
3987                x,y,z = peak[1:4]
3988                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
3989                XY = [int(X),int(View[3]-Y)]
3990                if np.allclose(xy,XY,atol=10) and Z < Zmax:
3991                    Zmax = Z
3992                    try:
3993                        Indx.remove(i)
3994                        ClearSelectedAtoms()
3995                        for id in Indx:
3996                            SetSelectedAtoms(id,Add)
3997                    except:
3998                        SetSelectedAtoms(i,Add)
3999        else:
4000            cx = drawingData['atomPtrs'][0]
4001            for i,atom in enumerate(drawAtoms):
4002                x,y,z = atom[cx:cx+3]
4003                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
4004                XY = [int(X),int(View[3]-Y)]
4005                if np.allclose(xy,XY,atol=10) and Z < Zmax:
4006                    Zmax = Z
4007                    try:
4008                        Indx.remove(i)
4009                        ClearSelectedAtoms()
4010                        for id in Indx:
4011                            SetSelectedAtoms(id,Add)
4012                    except:
4013                        SetSelectedAtoms(i,Add)
4014                                       
4015    def OnMouseDown(event):
4016        xy = event.GetPosition()
4017        if event.ShiftDown():
4018            if event.LeftIsDown():
4019                GetTruePosition(xy)
4020            elif event.RightIsDown():
4021                GetTruePosition(xy,True)
4022        else:
4023            drawingData['oldxy'] = list(xy)
4024       
4025    def OnMouseMove(event):
4026        if event.ShiftDown():           #don't want any inadvertant moves when picking
4027            return
4028        newxy = event.GetPosition()
4029                               
4030        if event.Dragging():
4031            if event.AltDown() and rbObj:
4032                if event.LeftIsDown():
4033                    SetRBRotation(newxy)
4034                    Q = rbObj['Orient'][0]
4035                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4036                elif event.RightIsDown():
4037                    SetRBTranslation(newxy)
4038                    Tx,Ty,Tz = rbObj['Orig'][0]
4039                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
4040                elif event.MiddleIsDown():
4041                    SetRBRotationZ(newxy)
4042                    Q = rbObj['Orient'][0]
4043                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4044                Draw('move')
4045            elif not event.ControlDown():
4046                if event.LeftIsDown():
4047                    SetRotation(newxy)
4048                    Q = drawingData['Quaternion']
4049                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4050                elif event.RightIsDown():
4051                    SetTranslation(newxy)
4052                    Tx,Ty,Tz = drawingData['viewPoint'][0]
4053                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
4054                elif event.MiddleIsDown():
4055                    SetRotationZ(newxy)
4056                    Q = drawingData['Quaternion']
4057                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4058                Draw('move')
4059           
4060       
4061    def OnMouseWheel(event):
4062        if event.ShiftDown():
4063            return
4064        drawingData['cameraPos'] += event.GetWheelRotation()/24.
4065        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
4066        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(drawingData['cameraPos']),1)
4067        page = getSelection()
4068        if page:
4069            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
4070                G2frame.dataDisplay.cameraPosTxt.SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
4071                G2frame.dataDisplay.cameraSlider.SetValue(drawingData['cameraPos'])
4072            Draw('wheel')
4073       
4074    def getSelection():
4075        try:
4076            return G2frame.dataDisplay.GetSelection()
4077        except AttributeError:
4078            G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!',1)
4079            return 0
4080           
4081    def SetViewPointText(VP):
4082        page = getSelection()
4083        if page:
4084            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
4085                G2frame.dataDisplay.viewPoint.SetValue('%.3f %.3f %.3f'%(VP[0],VP[1],VP[2]))
4086               
4087    def SetRBOrigText():
4088        page = getSelection()
4089        if page:
4090            if G2frame.dataDisplay.GetPageText(page) == 'RB Models':
4091                for i,sizer in enumerate(testRBObj['Sizers']['Xsizers']):
4092                    sizer.SetValue('%8.5f'%(testRBObj['rbObj']['Orig'][0][i]))
4093                   
4094    def SetRBOrienText():
4095        page = getSelection()
4096        if page:
4097            if G2frame.dataDisplay.GetPageText(page) == 'RB Models':
4098                for i,sizer in enumerate(testRBObj['Sizers']['Osizers']):
4099                    sizer.SetValue('%8.5f'%(testRBObj['rbObj']['Orient'][0][i]))
4100               
4101    def SetViewDirText(VD):
4102        page = getSelection()
4103        if page:
4104            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
4105                G2frame.dataDisplay.viewDir.SetValue('%.3f %.3f %.3f'%(VD[0],VD[1],VD[2]))
4106               
4107    def SetMapPeaksText(mapPeaks):
4108        page = getSelection()
4109        if page:
4110            if G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
4111                G2frame.MapPeaksTable.SetData(mapPeaks)
4112                panel = G2frame.dataDisplay.GetPage(page).GetChildren()
4113                names = [child.GetName() for child in panel]
4114                panel[names.index('grid window')].Refresh()
4115           
4116    def ClearSelectedAtoms():
4117        page = getSelection()
4118        if page:
4119            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
4120                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
4121            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
4122                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
4123            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
4124                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
4125               
4126                   
4127    def SetSelectedAtoms(ind,Add=False):
4128        page = getSelection()
4129        if page:
4130            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
4131                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)      #this is the Atoms grid in Draw Atoms
4132            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
4133                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)                 
4134            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
4135                Id = drawAtoms[ind][-3]
4136                for i,atom in enumerate(atomData):
4137                    if atom[-1] == Id:
4138                        G2frame.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
4139                 
4140    def GetSelectedAtoms():
4141        page = getSelection()
4142        Ind = []
4143        if page:
4144            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
4145                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
4146            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
4147                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()
4148            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
4149                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
4150        return Ind
4151                                       
4152    def SetBackground():
4153        R,G,B,A = Page.camera['backColor']
4154        glClearColor(R,G,B,A)
4155        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
4156       
4157    def SetLights():
4158        glEnable(GL_DEPTH_TEST)
4159        glShadeModel(GL_SMOOTH)
4160        glEnable(GL_LIGHTING)
4161        glEnable(GL_LIGHT0)
4162        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
4163        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
4164        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
4165       
4166    def GetRoll(newxy,rho):
4167        Q = drawingData['Quaternion']
4168        dxy = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,newxy+[0,]))
4169        dxy = np.array(dxy*rho.shape)       
4170        roll = np.where(dxy>0.5,1,np.where(dxy<-.5,-1,0))
4171        return roll
4172               
4173    def SetMapRoll(newxy):
4174        rho = mapData['rho']
4175        roll = GetRoll(newxy,rho)
4176        mapData['rho'] = np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
4177        drawingData['oldxy'] = list(newxy)
4178       
4179    def SetPeakRoll(newxy):
4180        rho = mapData['rho']
4181        roll = GetRoll(newxy,rho)
4182        steps = 1./np.array(rho.shape)
4183        dxy = roll*steps
4184        for peak in mapPeaks:
4185            peak[1:4] += dxy
4186            peak[1:4] %= 1.
4187            peak[4] = np.sqrt(np.sum(np.inner(Amat,peak[1:4])**2))
4188               
4189    def SetTranslation(newxy):
4190#first get translation vector in screen coords.       
4191        oldxy = drawingData['oldxy']
4192        if not len(oldxy): oldxy = list(newxy)
4193        dxy = newxy-oldxy
4194        drawingData['oldxy'] = list(newxy)
4195        V = np.array([-dxy[0],dxy[1],0.])
4196#then transform to rotated crystal coordinates & apply to view point       
4197        Q = drawingData['Quaternion']
4198        V = np.inner(Bmat,G2mth.prodQVQ(G2mth.invQ(Q),V))
4199        Tx,Ty,Tz = drawingData['viewPoint'][0]
4200        Tx += V[0]*0.01
4201        Ty += V[1]*0.01
4202        Tz += V[2]*0.01
4203        drawingData['viewPoint'][0] =  Tx,Ty,Tz
4204        SetViewPointText([Tx,Ty,Tz])
4205       
4206    def SetRBTranslation(newxy):
4207#first get translation vector in screen coords.       
4208        oldxy = drawingData['oldxy']
4209        if not len(oldxy): oldxy = list(newxy)
4210        dxy = newxy-oldxy
4211        drawingData['oldxy'] = list(newxy)
4212        V = np.array([-dxy[0],dxy[1],0.])
4213#then transform to rotated crystal coordinates & apply to RB origin       
4214        Q = drawingData['Quaternion']
4215        V = np.inner(Bmat,G2mth.prodQVQ(G2mth.invQ(Q),V))
4216        Tx,Ty,Tz = rbObj['Orig'][0]
4217        Tx -= V[0]*0.01
4218        Ty -= V[1]*0.01
4219        Tz -= V[2]*0.01
4220        rbObj['Orig'][0] =  Tx,Ty,Tz
4221        SetRBOrigText()
4222       
4223    def SetRotation(newxy):
4224        'Perform a rotation in x-y space due to a left-mouse drag'
4225    #first get rotation vector in screen coords. & angle increment       
4226        oldxy = drawingData['oldxy']
4227        if not len(oldxy): oldxy = list(newxy)
4228        dxy = newxy-oldxy
4229        drawingData['oldxy'] = list(newxy)
4230        V = np.array([dxy[1],dxy[0],0.])
4231        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
4232        if not A: return # nothing changed, nothing to do
4233    # next transform vector back to xtal coordinates via inverse quaternion
4234    # & make new quaternion
4235        Q = drawingData['Quaternion']
4236        V = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,V))
4237        DQ = G2mth.AVdeg2Q(A,V)
4238        Q = G2mth.prodQQ(Q,DQ)
4239        drawingData['Quaternion'] = Q
4240    # finally get new view vector - last row of rotation matrix
4241        VD = np.inner(Bmat,G2mth.Q2Mat(Q)[2])
4242        VD /= np.sqrt(np.sum(VD**2))
4243        drawingData['viewDir'] = VD
4244        SetViewDirText(VD)
4245       
4246    def SetRotationZ(newxy):                       
4247#first get rotation vector (= view vector) in screen coords. & angle increment       
4248        View = glGetIntegerv(GL_VIEWPORT)
4249        cent = [View[2]/2,View[3]/2]
4250        oldxy = drawingData['oldxy']
4251        if not len(oldxy): oldxy = list(newxy)
4252        dxy = newxy-oldxy
4253        drawingData['oldxy'] = list(newxy)
4254        V = drawingData['viewDir']
4255        A = [0,0]
4256        A[0] = dxy[1]*.25
4257        A[1] = dxy[0]*.25
4258        if newxy[0] > cent[0]:
4259            A[0] *= -1
4260        if newxy[1] < cent[1]:
4261            A[1] *= -1       
4262# next transform vector back to xtal coordinates & make new quaternion
4263        Q = drawingData['Quaternion']
4264        V = np.inner(Amat,V)
4265        Qx = G2mth.AVdeg2Q(A[0],V)
4266        Qy = G2mth.AVdeg2Q(A[1],V)
4267        Q = G2mth.prodQQ(Q,Qx)
4268        Q = G2mth.prodQQ(Q,Qy)
4269        drawingData['Quaternion'] = Q
4270
4271    def SetRBRotation(newxy):
4272#first get rotation vector in screen coords. & angle increment       
4273        oldxy = drawingData['oldxy']
4274        if not len(oldxy): oldxy = list(newxy)
4275        dxy = newxy-oldxy
4276        drawingData['oldxy'] = list(newxy)
4277        V = np.array([dxy[1],dxy[0],0.])
4278        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
4279# next transform vector back to xtal coordinates via inverse quaternion
4280# & make new quaternion
4281        Q = rbObj['Orient'][0]              #rotate RB to Cart
4282        QC = drawingData['Quaternion']      #rotate Cart to drawing
4283        V = G2mth.prodQVQ(G2mth.invQ(QC),V)
4284        V = G2mth.prodQVQ(G2mth.invQ(Q),V)
4285        DQ = G2mth.AVdeg2Q(A,V)
4286        Q = G2mth.prodQQ(Q,DQ)
4287        rbObj['Orient'][0] = Q
4288        SetRBOrienText()
4289       
4290    def SetRBRotationZ(newxy):                       
4291#first get rotation vector (= view vector) in screen coords. & angle increment       
4292        View = glGetIntegerv(GL_VIEWPORT)
4293        cent = [View[2]/2,View[3]/2]
4294        oldxy = drawingData['oldxy']
4295        if not len(oldxy): oldxy = list(newxy)
4296        dxy = newxy-oldxy
4297        drawingData['oldxy'] = list(newxy)
4298        V = drawingData['viewDir']
4299        A = [0,0]
4300        A[0] = dxy[1]*.25
4301        A[1] = dxy[0]*.25
4302        if newxy[0] < cent[0]:
4303            A[0] *= -1
4304        if newxy[1] > cent[1]:
4305            A[1] *= -1       
4306# next transform vector back to RB coordinates & make new quaternion
4307        Q = rbObj['Orient'][0]              #rotate RB to cart
4308        V = np.inner(Amat,V)
4309        V = -G2mth.prodQVQ(G2mth.invQ(Q),V)
4310        Qx = G2mth.AVdeg2Q(A[0],V)
4311        Qy = G2mth.AVdeg2Q(A[1],V)
4312        Q = G2mth.prodQQ(Q,Qx)
4313        Q = G2mth.prodQQ(Q,Qy)
4314        rbObj['Orient'][0] = Q
4315        SetRBOrienText()
4316
4317    def RenderBox():
4318        glEnable(GL_COLOR_MATERIAL)
4319        glLineWidth(2)
4320        glEnable(GL_BLEND)
4321        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
4322        glEnable(GL_LINE_SMOOTH)
4323        glBegin(GL_LINES)
4324        for line,color in zip(uEdges,uColors):
4325            glColor3ubv(color<