source: trunk/GSASIIplot.py @ 1443

Last change on this file since 1443 was 1443, checked in by vondreele, 7 years ago

add calibration of lam, difC, etc. from index peak positions
new plotCalib routine

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