source: trunk/GSASIIplot.py @ 1579

Last change on this file since 1579 was 1579, checked in by toby, 8 years ago

revamp masks GUI

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