source: trunk/GSASIIplot.py @ 1448

Last change on this file since 1448 was 1448, checked in by toby, 9 years ago

quick wx2.9+ fixes

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