source: trunk/GSASIIplot.py @ 1623

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

remove stray prints
double up modulation plot
set roll of 4D map to match that for 3D map
remove PlotPattern? call after reflection list selection - not needed & crashes for single crystal data
4D charge flipping reproduces J2K modulation maps!

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