source: trunk/GSASIIplot.py @ 1619

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

add columnsorter for hetrogeneous seq ref; reorg to start moving widgets out of g2grid

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 218.9 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIplot: plotting routines*
4===============================
5
6'''
7########### SVN repository information ###################
8# $Date: 2014-12-27 03:58:22 +0000 (Sat, 27 Dec 2014) $
9# $Author: toby $
10# $Revision: 1619 $
11# $URL: trunk/GSASIIplot.py $
12# $Id: GSASIIplot.py 1619 2014-12-27 03:58:22Z toby $
13########### SVN repository information ###################
14import math
15import time
16import copy
17import sys
18import os.path
19import numpy as np
20import numpy.ma as ma
21import numpy.linalg as nl
22import wx
23import wx.aui
24import wx.glcanvas
25import matplotlib as mpl
26import mpl_toolkits.mplot3d.axes3d as mp3d
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 1619 $")
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:
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:  #will be trouble for SS reflection lists - will need peak[:7]
1389                try:
1390                    for peak in Phases[G2frame.RefList]['RefList']:
1391                        if len(peak) > 15:
1392                            HKL.append(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:
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    print 'modulation plot for',atom[0]
2779    try:
2780        plotNum = G2frame.G2plotNB.plotList.index('Modulation')
2781        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2782        Page.figure.clf()
2783        Plot = Page.figure.gca()
2784        if not Page.IsShown():
2785            Page.Show()
2786    except ValueError:
2787        Plot = G2frame.G2plotNB.addMpl('Modulation').gca()
2788        plotNum = G2frame.G2plotNB.plotList.index('Modulation')
2789        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2790#        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2791
2792    General = data['General']
2793    cx,ct,cs,cia = General['AtomPtrs']
2794    Map = General['4DmapData']
2795    MapType = Map['MapType']
2796    rhoSize = np.array(Map['rho'].shape)
2797    atxyz = np.array(atom[cx:cx+3])
2798    Title = MapType+' modulation map for atom '+atom[0]+    \
2799        ' at %.4f %.4f %.4f'%(atxyz[0],atxyz[1],atxyz[2])
2800    ix = -np.array(np.rint(rhoSize[:3]*atxyz),dtype='i')
2801    ix += (rhoSize[:3]/2)
2802    ix = ix%rhoSize[:3]
2803    print 'roll',ix,atxyz
2804    rho = np.roll(np.roll(np.roll(Map['rho'],ix[0],axis=0),ix[1],axis=1),ix[2],axis=2)
2805    ix = rhoSize[:3]/2
2806    print 'slab',ix,rhoSize
2807    ib = 4
2808    if Ax == 'x':
2809        slab = np.sum(np.sum(rho[:,ix[1]-ib:ix[1]+ib,ix[2]-ib:ix[2]+ib,:],axis=2),axis=1)
2810    elif Ax == 'y':
2811        slab = np.sum(np.sum(rho[ix[0]-ib:ix[0]+ib,:,ix[2]-ib:ix[2]+ib,:],axis=2),axis=0)
2812    elif Ax == 'z':
2813        slab = np.sum(np.sum(rho[ix[0]-ib:ix[0]+ib,ix[1]-ib:ix[1]+ib,:,:],axis=1),axis=0)
2814    Plot.set_title(Title)
2815    Plot.set_xlabel('t')
2816    Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax))
2817    Plot.contour(slab,20,extent=(0.,1.,-.5,.5))
2818    Page.canvas.draw()
2819   
2820################################################################################
2821##### PlotCovariance
2822################################################################################
2823           
2824def PlotCovariance(G2frame,Data):
2825    'needs a doc string'
2826    if not Data:
2827        print 'No covariance matrix available'
2828        return
2829    varyList = Data['varyList']
2830    values = Data['variables']
2831    Xmax = len(varyList)
2832    covMatrix = Data['covMatrix']
2833    sig = np.sqrt(np.diag(covMatrix))
2834    xvar = np.outer(sig,np.ones_like(sig))
2835    covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
2836    title = ' for\n'+Data['title']
2837    newAtomDict = Data.get('newAtomDict',{})
2838   
2839
2840    def OnPlotKeyPress(event):
2841        newPlot = False
2842        if event.key == 's':
2843            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
2844            choice.sort()
2845            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
2846            if dlg.ShowModal() == wx.ID_OK:
2847                sel = dlg.GetSelection()
2848                G2frame.VcovColor = choice[sel]
2849            else:
2850                G2frame.VcovColor = 'RdYlGn'
2851            dlg.Destroy()
2852        PlotCovariance(G2frame,Data)
2853
2854    def OnMotion(event):
2855        #there is a problem here - reports wrong values
2856        if event.button:
2857            ytics = imgAx.get_yticks()
2858            ytics = np.where(ytics<len(varyList),ytics,-1)
2859            ylabs = [np.where(0<=i ,varyList[int(i)],' ') for i in ytics]
2860            imgAx.set_yticklabels(ylabs)           
2861        if event.xdata and event.ydata:                 #avoid out of frame errors
2862            xpos = int(event.xdata+.5)
2863            ypos = int(event.ydata+.5)
2864            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
2865                if xpos == ypos:
2866                    value = values[xpos]
2867                    name = varyList[xpos]
2868                    if varyList[xpos] in newAtomDict:
2869                        name,value = newAtomDict[name]                       
2870                    msg = '%s value = %.4g, esd = %.4g'%(name,value,sig[xpos])
2871                else:
2872                    msg = '%s - %s: %5.3f'%(varyList[xpos],varyList[ypos],covArray[xpos][ypos])
2873                Page.canvas.SetToolTipString(msg)
2874                G2frame.G2plotNB.status.SetFields(['',msg])
2875               
2876    try:
2877        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
2878        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2879        Page.figure.clf()
2880        Plot = Page.figure.gca()
2881        if not Page.IsShown():
2882            Page.Show()
2883    except ValueError:
2884        Plot = G2frame.G2plotNB.addMpl('Covariance').gca()
2885        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
2886        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2887        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2888        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
2889    Page.Choice = ['s: to change colors']
2890    Page.keyPress = OnPlotKeyPress
2891    Page.SetFocus()
2892    G2frame.G2plotNB.status.SetFields(['',''])   
2893    acolor = mpl.cm.get_cmap(G2frame.VcovColor)
2894    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower',
2895        vmin=-1.,vmax=1.)
2896    imgAx = Img.get_axes()
2897    ytics = imgAx.get_yticks()
2898    ylabs = [varyList[int(i)] for i in ytics[:-1]]
2899    imgAx.set_yticklabels(ylabs)
2900    colorBar = Page.figure.colorbar(Img)
2901    Plot.set_title('V-Cov matrix'+title)
2902    Plot.set_xlabel('Variable number')
2903    Plot.set_ylabel('Variable name')
2904    Page.canvas.draw()
2905   
2906################################################################################
2907##### PlotTorsion
2908################################################################################
2909
2910def PlotTorsion(G2frame,phaseName,Torsion,TorName,Names=[],Angles=[],Coeff=[]):
2911    'needs a doc string'
2912   
2913    global names
2914    names = Names
2915    sum = np.sum(Torsion)
2916    torsion = np.log(2*Torsion+1.)/sum
2917    tMin = np.min(torsion)
2918    tMax = np.max(torsion)
2919    torsion = 3.*(torsion-tMin)/(tMax-tMin)
2920    X = np.linspace(0.,360.,num=45)
2921   
2922    def OnPick(event):
2923        ind = event.ind[0]
2924        msg = 'atoms:'+names[ind]
2925        Page.canvas.SetToolTipString(msg)
2926        try:
2927            page = G2frame.dataDisplay.GetSelection()
2928        except:
2929            return
2930        if G2frame.dataDisplay.GetPageText(page) == 'Torsion restraints':
2931            torGrid = G2frame.dataDisplay.GetPage(page).Torsions
2932            torGrid.ClearSelection()
2933            for row in range(torGrid.GetNumberRows()):
2934                if names[ind] in torGrid.GetCellValue(row,0):
2935                    torGrid.SelectRow(row)
2936            torGrid.ForceRefresh()
2937               
2938    def OnMotion(event):
2939        if event.xdata and event.ydata:                 #avoid out of frame errors
2940            xpos = event.xdata
2941            ypos = event.ydata
2942            msg = 'torsion,energy: %5.3f %5.3f'%(xpos,ypos)
2943            Page.canvas.SetToolTipString(msg)
2944
2945    try:
2946        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
2947        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2948        Page.figure.clf()
2949        Plot = Page.figure.gca()
2950        if not Page.IsShown():
2951            Page.Show()
2952    except ValueError:
2953        Plot = G2frame.G2plotNB.addMpl('Torsion').gca()
2954        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
2955        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2956        Page.canvas.mpl_connect('pick_event', OnPick)
2957        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2958   
2959    Page.SetFocus()
2960    G2frame.G2plotNB.status.SetFields(['','Use mouse LB to identify torsion atoms'])
2961    Plot.plot(X,torsion,'b+')
2962    if len(Coeff):
2963        X2 = np.linspace(0.,360.,45)
2964        Y2 = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in X2])
2965        Plot.plot(X2,Y2,'r')
2966    if len(Angles):
2967        Eval = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in Angles])
2968        Plot.plot(Angles,Eval,'ro',picker=5)
2969    Plot.set_xlim((0.,360.))
2970    Plot.set_title('Torsion angles for '+TorName+' in '+phaseName)
2971    Plot.set_xlabel('angle',fontsize=16)
2972    Plot.set_ylabel('Energy',fontsize=16)
2973    Page.canvas.draw()
2974   
2975################################################################################
2976##### PlotRama
2977################################################################################
2978
2979def PlotRama(G2frame,phaseName,Rama,RamaName,Names=[],PhiPsi=[],Coeff=[]):
2980    'needs a doc string'
2981
2982    global names
2983    names = Names
2984    rama = np.log(2*Rama+1.)
2985    ramaMax = np.max(rama)
2986    rama = np.reshape(rama,(45,45))
2987    global Phi,Psi
2988    Phi = []
2989    Psi = []
2990
2991    def OnPlotKeyPress(event):
2992        newPlot = False
2993        if event.key == 's':
2994            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
2995            choice.sort()
2996            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
2997            if dlg.ShowModal() == wx.ID_OK:
2998                sel = dlg.GetSelection()
2999                G2frame.RamaColor = choice[sel]
3000            else:
3001                G2frame.RamaColor = 'RdYlGn'
3002            dlg.Destroy()
3003        PlotRama(G2frame,phaseName,Rama)
3004
3005    def OnPick(event):
3006        ind = event.ind[0]
3007        msg = 'atoms:'+names[ind]
3008        Page.canvas.SetToolTipString(msg)
3009        try:
3010            page = G2frame.dataDisplay.GetSelection()
3011        except:
3012            return
3013        if G2frame.dataDisplay.GetPageText(page) == 'Ramachandran restraints':
3014            ramaGrid = G2frame.dataDisplay.GetPage(page).Ramas
3015            ramaGrid.ClearSelection()
3016            for row in range(ramaGrid.GetNumberRows()):
3017                if names[ind] in ramaGrid.GetCellValue(row,0):
3018                    ramaGrid.SelectRow(row)
3019            ramaGrid.ForceRefresh()
3020
3021    def OnMotion(event):
3022        if event.xdata and event.ydata:                 #avoid out of frame errors
3023            xpos = event.xdata
3024            ypos = event.ydata
3025            msg = 'phi/psi: %5.3f %5.3f'%(xpos,ypos)
3026            Page.canvas.SetToolTipString(msg)
3027           
3028    try:
3029        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
3030        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3031        Page.figure.clf()
3032        Plot = Page.figure.gca()
3033        if not Page.IsShown():
3034            Page.Show()
3035    except ValueError:
3036        Plot = G2frame.G2plotNB.addMpl('Ramachandran').gca()
3037        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
3038        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3039        Page.canvas.mpl_connect('pick_event', OnPick)
3040        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3041        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
3042
3043    Page.Choice = ['s: to change colors']
3044    Page.keyPress = OnPlotKeyPress
3045    Page.SetFocus()
3046    G2frame.G2plotNB.status.SetFields(['','Use mouse LB to identify phi/psi atoms'])
3047    acolor = mpl.cm.get_cmap(G2frame.RamaColor)
3048    if RamaName == 'All' or '-1' in RamaName:
3049        if len(Coeff): 
3050            X,Y = np.meshgrid(np.linspace(-180.,180.,45),np.linspace(-180.,180.,45))
3051            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
3052            Plot.contour(X,Y,np.reshape(Z,(45,45)))
3053        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
3054            extent=[-180,180,-180,180],origin='lower')
3055        if len(PhiPsi):
3056            Phi,Psi = PhiPsi.T
3057            Phi = np.where(Phi>180.,Phi-360.,Phi)
3058            Psi = np.where(Psi>180.,Psi-360.,Psi)
3059            Plot.plot(Phi,Psi,'ro',picker=5)
3060        Plot.set_xlim((-180.,180.))
3061        Plot.set_ylim((-180.,180.))
3062    else:
3063        if len(Coeff): 
3064            X,Y = np.meshgrid(np.linspace(0.,360.,45),np.linspace(0.,360.,45))
3065            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
3066            Plot.contour(X,Y,np.reshape(Z,(45,45)))
3067        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
3068            extent=[0,360,0,360],origin='lower')
3069        if len(PhiPsi):
3070            Phi,Psi = PhiPsi.T
3071            Plot.plot(Phi,Psi,'ro',picker=5)
3072        Plot.set_xlim((0.,360.))
3073        Plot.set_ylim((0.,360.))
3074    Plot.set_title('Ramachandran for '+RamaName+' in '+phaseName)
3075    Plot.set_xlabel(r'$\phi$',fontsize=16)
3076    Plot.set_ylabel(r'$\psi$',fontsize=16)
3077    colorBar = Page.figure.colorbar(Img)
3078    Page.canvas.draw()
3079
3080
3081################################################################################
3082##### PlotSeq
3083################################################################################
3084def PlotSelectedSequence(G2frame,ColumnList,TableGet,SelectX,fitnum=None,fitvals=None):
3085    '''Plot a result from a sequential refinement
3086
3087    :param wx.Frame G2frame: The main GSAS-II tree "window"
3088    :param list ColumnList: list of int values corresponding to columns
3089      selected as y values
3090    :param function TableGet: a function that takes a column number
3091      as argument and returns the column label, the values and there ESDs (or None)
3092    :param function SelectX: a function that returns a selected column
3093      number (or None) as the X-axis selection
3094    '''
3095    global Title,xLabel,yLabel
3096    xLabel = yLabel = Title = ''
3097    def OnMotion(event):
3098        if event.xdata and event.ydata:                 #avoid out of frame errors
3099            xpos = event.xdata
3100            ypos = event.ydata
3101            msg = '%5.3f %.6g'%(xpos,ypos)
3102            Page.canvas.SetToolTipString(msg)
3103
3104    def OnKeyPress(event):
3105        global Title,xLabel,yLabel
3106        if event.key == 's':
3107            G2frame.seqXaxis = G2frame.seqXselect()
3108            Draw()
3109        elif event.key == 't':
3110            dlg = G2gd.MultiStringDialog(G2frame,'Set titles & labels',[' Title ',' x-Label ',' y-Label '],
3111                [Title,xLabel,yLabel])
3112            if dlg.Show():
3113                Title,xLabel,yLabel = dlg.GetValues()
3114            dlg.Destroy()
3115            Draw()
3116           
3117    def Draw():
3118        global Title,xLabel,yLabel
3119        Page.SetFocus()
3120        G2frame.G2plotNB.status.SetStatusText('press s to select X axis, t to change titles',1)
3121        Plot.clear()
3122        if G2frame.seqXaxis is not None:   
3123            xName,X,Xsig = Page.seqTableGet(G2frame.seqXaxis)
3124        else:
3125            X = np.arange(0,G2frame.SeqTable.GetNumberRows(),1)
3126            xName = 'Data sequence number'
3127        for col in Page.seqYaxisList:
3128            name,Y,sig = Page.seqTableGet(col)
3129            # deal with missing (None) values
3130            Xnew = []
3131            Ynew = []
3132            Ysnew = []
3133            for i in range(len(X)):
3134                if X[i] is None or Y[i] is None: continue
3135                Xnew.append(X[i])
3136                Ynew.append(Y[i])
3137                if sig: Ysnew.append(sig[i])
3138            if Ysnew:
3139                if G2frame.seqReverse and not G2frame.seqXaxis:
3140                    Ynew = Ynew[::-1]
3141                    Ysnew = Ysnew[::-1]
3142                Plot.errorbar(Xnew,Ynew,yerr=Ysnew,label=name)
3143            else:
3144                if G2frame.seqReverse and not G2frame.seqXaxis:
3145                    Ynew = Ynew[::-1]
3146                Plot.plot(Xnew,Ynew)
3147                Plot.plot(Xnew,Ynew,'o',label=name)
3148        if Page.fitvals: # TODO: deal with fitting of None values
3149            if G2frame.seqReverse and not G2frame.seqXaxis:
3150                Page.fitvals = Page.fitvals[::-1]
3151            Plot.plot(X,Page.fitvals,label='Fit')
3152           
3153        Plot.legend(loc='best')
3154        print Title,xLabel,yLabel
3155        if Title:
3156            Plot.set_title(Title)
3157        else:
3158            Plot.set_title('')
3159        if xLabel:
3160            Plot.set_xlabel(xLabel)
3161        else:
3162            Plot.set_xlabel(xName)
3163        if yLabel:
3164            Plot.set_ylabel(yLabel)
3165        else:
3166            Plot.set_ylabel('Parameter values')
3167        Page.canvas.draw()           
3168           
3169    G2frame.seqXselect = SelectX
3170    try:
3171        G2frame.seqXaxis
3172    except:
3173        G2frame.seqXaxis = None
3174
3175    if fitnum is None:
3176        label = 'Sequential refinement'
3177    else:
3178        label = 'Parametric fit #'+str(fitnum+1)
3179    try:
3180        plotNum = G2frame.G2plotNB.plotList.index(label)
3181        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3182        Page.figure.clf()
3183        Plot = Page.figure.gca()
3184        if not Page.IsShown():
3185            Page.Show()
3186    except ValueError:
3187        Plot = G2frame.G2plotNB.addMpl(label).gca()
3188        plotNum = G2frame.G2plotNB.plotList.index(label)
3189        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3190        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
3191        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3192    Page.Choice = ['s - select x-axis','t - change titles',]
3193    Page.keyPress = OnKeyPress
3194    Page.seqYaxisList = ColumnList
3195    Page.seqTableGet = TableGet
3196    Page.fitvals = fitvals
3197       
3198    Draw()
3199    G2frame.G2plotNB.nb.SetSelection(plotNum) # raises plot tab
3200               
3201################################################################################
3202##### PlotExposedImage & PlotImage
3203################################################################################
3204           
3205def PlotExposedImage(G2frame,newPlot=False,event=None):
3206    '''General access module for 2D image plotting
3207    '''
3208    plotNo = G2frame.G2plotNB.nb.GetSelection()
3209    if G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
3210        PlotImage(G2frame,newPlot,event,newImage=True)
3211    elif G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
3212        PlotIntegration(G2frame,newPlot,event)
3213
3214def OnStartMask(G2frame):
3215    '''Initiate the start of a Frame or Polygon map
3216
3217    :param wx.Frame G2frame: The main GSAS-II tree "window"
3218    :param str eventkey: a single letter ('f' or 'p') that
3219      determines what type of mask is created.   
3220    '''
3221    Masks = G2frame.PatternTree.GetItemPyData(
3222        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
3223    if G2frame.MaskKey == 'f':
3224        Masks['Frames'] = []
3225    elif G2frame.MaskKey == 'p':
3226        Masks['Polygons'].append([])
3227    elif G2frame.MaskKey == 's':
3228        Masks['Points'].append([])
3229    elif G2frame.MaskKey == 'a':
3230        Masks['Arcs'].append([])
3231    elif G2frame.MaskKey == 'r':
3232        Masks['Rings'].append([])
3233    G2imG.UpdateMasks(G2frame,Masks)
3234    PlotImage(G2frame,newImage=True)
3235   
3236def OnStartNewDzero(G2frame):
3237    '''Initiate the start of adding a new d-zero to a strain data set
3238
3239    :param wx.Frame G2frame: The main GSAS-II tree "window"
3240    :param str eventkey: a single letter ('a') that
3241      triggers the addition of a d-zero.   
3242    '''
3243    G2frame.dataFrame.GetStatusBar().SetStatusText('Add strain ring active - LB pick d-zero value',0)
3244    G2frame.PickId = G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain')
3245    data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
3246    return data
3247
3248def PlotImage(G2frame,newPlot=False,event=None,newImage=True):
3249    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
3250    masks, etc.
3251    '''
3252    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
3253    import numpy.ma as ma
3254    Dsp = lambda tth,wave: wave/(2.*npsind(tth/2.))
3255    global Data,Masks,StrSta
3256    colors=['b','g','r','c','m','k']
3257    Data = G2frame.PatternTree.GetItemPyData(
3258        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3259# patch
3260    if 'invert_x' not in Data:
3261        Data['invert_x'] = False
3262        Data['invert_y'] = True
3263# end patch
3264    Masks = G2frame.PatternTree.GetItemPyData(
3265        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
3266    try:    #may be absent
3267        StrSta = G2frame.PatternTree.GetItemPyData(
3268            G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain'))
3269    except TypeError:   #is missing
3270        StrSta = {}
3271
3272    def OnImMotion(event):
3273        Page.canvas.SetToolTipString('')
3274        sizexy = Data['size']
3275        if event.xdata and event.ydata and len(G2frame.ImageZ):                 #avoid out of frame errors
3276            Page.canvas.SetToolTipString('%8.2f %8.2fmm'%(event.xdata,event.ydata))
3277            Page.canvas.SetCursor(wx.CROSS_CURSOR)
3278            item = G2frame.itemPicked
3279            pixelSize = Data['pixelSize']
3280            scalex = 1000./pixelSize[0]
3281            scaley = 1000./pixelSize[1]
3282            if item and G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Image Controls':
3283                if 'Text' in str(item):
3284                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
3285                else:
3286                    xcent,ycent = Data['center']
3287                    xpos = event.xdata-xcent
3288                    ypos = event.ydata-ycent
3289                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
3290                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
3291                        Page.canvas.SetToolTipString('%6d deg'%(azm))
3292                    elif 'line1' in  str(item) or 'line2' in str(item):
3293                        Page.canvas.SetToolTipString('%8.3f deg'%(tth))                           
3294            else:
3295                xpos = event.xdata
3296                ypos = event.ydata
3297                xpix = xpos*scalex
3298                ypix = ypos*scaley
3299                Int = 0
3300                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
3301                    Int = G2frame.ImageZ[ypix][xpix]
3302                tth,azm,D,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
3303                Q = 2.*math.pi/dsp
3304                fields = ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)]
3305                if G2frame.MaskKey in ['p','f']:
3306                    fields[1] = 'Polygon/frame mask pick - LB next point, RB close polygon'
3307                elif G2frame.StrainKey:
3308                    fields[0] = 'd-zero pick active'
3309                G2frame.G2plotNB.status.SetFields(fields)
3310
3311    def OnImPlotKeyPress(event):
3312        try:
3313            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
3314        except TypeError:
3315            return
3316        if PickName == 'Masks':
3317            if event.key in ['l','p','f','s','a','r']:
3318                G2frame.MaskKey = event.key
3319                OnStartMask(G2frame)
3320                PlotImage(G2frame,newPlot=False)
3321               
3322        elif PickName == 'Stress/Strain':
3323            if event.key in ['a',]:
3324                G2frame.StrainKey = event.key
3325                StrSta = OnStartNewDzero(G2frame)
3326                PlotImage(G2frame,newPlot=False)
3327               
3328        elif PickName == 'Image Controls':
3329            if event.key in ['c',]:
3330                Xpos = event.xdata
3331                if not Xpos:            #got point out of frame
3332                    return
3333                Ypos = event.ydata
3334                dlg = wx.MessageDialog(G2frame,'Are you sure you want to change the center?',
3335                    'Center change',style=wx.OK|wx.CANCEL)
3336                try:
3337                    if dlg.ShowModal() == wx.ID_OK:
3338                        print 'move center to: ',Xpos,Ypos
3339                        Data['center'] = [Xpos,Ypos]
3340                        G2imG.UpdateImageControls(G2frame,Data,Masks)
3341                        PlotImage(G2frame,newPlot=False)
3342                finally:
3343                    dlg.Destroy()
3344                return
3345            elif event.key == 'l':
3346                G2frame.logPlot = not G2frame.logPlot
3347            elif event.key in ['x',]:
3348                Data['invert_x'] = not Data['invert_x']
3349            elif event.key in ['y',]:
3350                Data['invert_y'] = not Data['invert_y']
3351            PlotImage(G2frame,newPlot=True)
3352           
3353    def OnKeyBox(event):
3354        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):
3355            event.key = cb.GetValue()[0]
3356            cb.SetValue(' key press')
3357            if event.key in ['l','s','a','r','p','x','y']:
3358                wx.CallAfter(OnImPlotKeyPress,event)
3359        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
3360                       
3361    def OnImPick(event):
3362        if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']:
3363            return
3364        if G2frame.itemPicked is not None: return
3365        G2frame.itemPicked = event.artist
3366        G2frame.mousePicked = event.mouseevent
3367       
3368    def OnImRelease(event):
3369        try:
3370            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
3371        except TypeError:
3372            return
3373        if PickName not in ['Image Controls','Masks','Stress/Strain']:
3374            return
3375        pixelSize = Data['pixelSize']
3376        scalex = 1000./pixelSize[0]
3377        scaley = 1000./pixelSize[1]
3378#        pixLimit = Data['pixLimit']    #can be too tight
3379        pixLimit = 20       #this makes the search box 40x40 pixels
3380        if G2frame.itemPicked is None and PickName == 'Image Controls' and len(G2frame.ImageZ):
3381            Xpos = event.xdata
3382            if not (Xpos and G2frame.ifGetRing):                   #got point out of frame
3383                return
3384            Ypos = event.ydata
3385            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
3386                if event.button == 1:
3387                    Xpix = Xpos*scalex
3388                    Ypix = Ypos*scaley
3389                    xpos,ypos,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
3390                    if I and J:
3391                        xpos += .5                              #shift to pixel center
3392                        ypos += .5
3393                        xpos /= scalex                          #convert to mm
3394                        ypos /= scaley
3395                        Data['ring'].append([xpos,ypos])
3396                elif event.button == 3:
3397                    G2frame.dataFrame.GetStatusBar().SetStatusText('Calibrating...',0)
3398                    if G2img.ImageCalibrate(G2frame,Data):
3399                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check',0)
3400                        print 'Calibration successful'
3401                    else:
3402                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose',0)
3403                        print 'Calibration failed'
3404                    G2frame.ifGetRing = False
3405                    G2imG.UpdateImageControls(G2frame,Data,Masks)
3406                    return
3407                PlotImage(G2frame,newImage=False)
3408            return
3409        elif G2frame.MaskKey and PickName == 'Masks':
3410            Xpos,Ypos = [event.xdata,event.ydata]
3411            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3412                return
3413            if G2frame.MaskKey == 's' and event.button == 1:
3414                Masks['Points'][-1] = [Xpos,Ypos,1.]
3415                G2frame.MaskKey = ''               
3416            elif G2frame.MaskKey == 'r' and event.button == 1:
3417                tth = G2img.GetTth(Xpos,Ypos,Data)
3418                Masks['Rings'][-1] = [tth,0.1]
3419                G2frame.MaskKey = ''               
3420            elif G2frame.MaskKey == 'a' and event.button == 1:
3421                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
3422                azm = int(azm)               
3423                Masks['Arcs'][-1] = [tth,[azm-5,azm+5],0.1]
3424                G2frame.MaskKey = ''               
3425            elif G2frame.MaskKey =='p':
3426                polygon = Masks['Polygons'][-1]
3427                if len(polygon) > 2 and event.button == 3:
3428                    x0,y0 = polygon[0]
3429                    polygon.append([x0,y0])
3430                    G2frame.MaskKey = ''
3431                    G2frame.G2plotNB.status.SetFields(['','Polygon closed'])
3432                else:
3433                    G2frame.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(Xpos,Ypos)])
3434                    polygon.append([Xpos,Ypos])
3435            elif G2frame.MaskKey =='f':
3436                frame = Masks['Frames']
3437                if len(frame) > 2 and event.button == 3:
3438                    x0,y0 = frame[0]
3439                    frame.append([x0,y0])
3440                    G2frame.MaskKey = ''
3441                    G2frame.G2plotNB.status.SetFields(['','Frame closed'])
3442                else:
3443                    G2frame.G2plotNB.status.SetFields(['','New frame point: %.1f,%.1f'%(Xpos,Ypos)])
3444                    frame.append([Xpos,Ypos])
3445            G2imG.UpdateMasks(G2frame,Masks)
3446            PlotImage(G2frame,newImage=False)
3447        elif PickName == 'Stress/Strain' and G2frame.StrainKey:
3448            Xpos,Ypos = [event.xdata,event.ydata]
3449            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3450                return
3451            dsp = G2img.GetDsp(Xpos,Ypos,Data)
3452            StrSta['d-zero'].append({'Dset':dsp,'Dcalc':0.0,'pixLimit':10,'cutoff':1.0,
3453                'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]],'Emat':[1.0,1.0,1.0]})
3454            R,r = G2img.MakeStrStaRing(StrSta['d-zero'][-1],G2frame.ImageZ,Data)
3455            if not len(R):
3456                del StrSta['d-zero'][-1]
3457                G2frame.ErrorDialog('Strain peak selection','WARNING - No points found for this ring selection')
3458            StrSta['d-zero'] = G2mth.sortArray(StrSta['d-zero'],'Dset',reverse=True)
3459            G2frame.StrainKey = ''
3460            G2imG.UpdateStressStrain(G2frame,StrSta)
3461            PlotStrain(G2frame,StrSta)
3462            PlotImage(G2frame,newPlot=False)           
3463        else:
3464            Xpos,Ypos = [event.xdata,event.ydata]
3465            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3466                return
3467            if G2frame.ifGetRing:                          #delete a calibration ring pick
3468                xypos = [Xpos,Ypos]
3469                rings = Data['ring']
3470                for ring in rings:
3471                    if np.allclose(ring,xypos,.01,0):
3472                        rings.remove(ring)
3473            else:
3474                tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data)[:3]
3475                itemPicked = str(G2frame.itemPicked)
3476                if 'Line2D' in itemPicked and PickName == 'Image Controls':
3477                    if 'line1' in itemPicked:
3478                        Data['IOtth'][0] = max(tth,0.001)
3479                    elif 'line2' in itemPicked:
3480                        Data['IOtth'][1] = tth
3481                    elif 'line3' in itemPicked:
3482                        Data['LRazimuth'][0] = int(azm)
3483                    elif 'line4' in itemPicked and not Data['fullIntegrate']:
3484                        Data['LRazimuth'][1] = int(azm)
3485                   
3486                    Data['LRazimuth'][0] %= 360
3487                    Data['LRazimuth'][1] %= 360
3488                    if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
3489                        Data['LRazimuth'][1] += 360                       
3490                    if Data['fullIntegrate']:
3491                        Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
3492                       
3493                    if  Data['IOtth'][0] > Data['IOtth'][1]:
3494                        Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
3495                       
3496                    G2frame.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
3497                    G2frame.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
3498                    G2frame.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
3499                    G2frame.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
3500                elif 'Circle' in itemPicked and PickName == 'Masks':
3501                    spots = Masks['Points']
3502                    newPos = itemPicked.split(')')[0].split('(')[2].split(',')
3503                    newPos = np.array([float(newPos[0]),float(newPos[1])])
3504                    for spot in spots:
3505                        if spot and np.allclose(np.array([spot[:2]]),newPos):
3506                            spot[:2] = Xpos,Ypos
3507                    G2imG.UpdateMasks(G2frame,Masks)
3508                elif 'Line2D' in itemPicked and PickName == 'Masks':
3509                    Obj = G2frame.itemPicked.findobj()
3510                    rings = Masks['Rings']
3511                    arcs = Masks['Arcs']
3512                    polygons = Masks['Polygons']
3513                    frame = Masks['Frames']
3514                    for ring in G2frame.ringList:
3515                        if Obj == ring[0]:
3516                            rN = ring[1]
3517                            if ring[2] == 'o':
3518                                rings[rN][0] = G2img.GetTth(Xpos,Ypos,Data)-rings[rN][1]/2.
3519                            else:
3520                                rings[rN][0] = G2img.GetTth(Xpos,Ypos,Data)+rings[rN][1]/2.
3521                    for arc in G2frame.arcList:
3522                        if Obj == arc[0]:
3523                            aN = arc[1]
3524                            if arc[2] == 'o':
3525                                arcs[aN][0] = G2img.GetTth(Xpos,Ypos,Data)-arcs[aN][2]/2
3526                            elif arc[2] == 'i':
3527                                arcs[aN][0] = G2img.GetTth(Xpos,Ypos,Data)+arcs[aN][2]/2
3528                            elif arc[2] == 'l':
3529                                arcs[aN][1][0] = int(G2img.GetAzm(Xpos,Ypos,Data))
3530                            else:
3531                                arcs[aN][1][1] = int(G2img.GetAzm(Xpos,Ypos,Data))
3532                    for poly in G2frame.polyList:   #merging points problem here?
3533                        if Obj == poly[0]:
3534                            ind = G2frame.itemPicked.contains(G2frame.mousePicked)[1]['ind'][0]
3535                            oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
3536                            pN = poly[1]
3537                            for i,xy in enumerate(polygons[pN]):
3538                                if np.allclose(np.array([xy]),oldPos,atol=1.0):
3539                                    if event.button == 1:
3540                                        polygons[pN][i] = Xpos,Ypos
3541                                    elif event.button == 3:
3542                                        polygons[pN].insert(i,[Xpos,Ypos])
3543                                        break
3544                    if frame:
3545                        oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
3546                        for i,xy in enumerate(frame):
3547                            if np.allclose(np.array([xy]),oldPos,atol=1.0):
3548                                if event.button == 1:
3549                                    frame[i] = Xpos,Ypos
3550                                elif event.button == 3:
3551                                    frame.insert(i,[Xpos,Ypos])
3552                                    break
3553                    G2imG.UpdateMasks(G2frame,Masks)
3554#                else:                  #keep for future debugging
3555#                    print str(G2frame.itemPicked),event.xdata,event.ydata,event.button
3556            PlotImage(G2frame,newImage=True)
3557            G2frame.itemPicked = None
3558           
3559    try:
3560        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
3561        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3562        if not newPlot:
3563            Plot = Page.figure.gca()          #get previous powder plot & get limits
3564            xylim = Plot.get_xlim(),Plot.get_ylim()
3565        if newImage:
3566            Page.figure.clf()
3567            Plot = Page.figure.gca()          #get a fresh plot after clf()
3568    except ValueError:
3569        Plot = G2frame.G2plotNB.addMpl('2D Powder Image').gca()
3570        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
3571        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3572        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
3573        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
3574        Page.canvas.mpl_connect('pick_event', OnImPick)
3575        Page.canvas.mpl_connect('button_release_event', OnImRelease)
3576        xylim = []
3577    Page.Choice = None
3578    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
3579        Page.SetFocus()
3580    Title = G2frame.PatternTree.GetItemText(G2frame.Image)[4:]
3581    G2frame.G2plotNB.status.DestroyChildren()
3582    if G2frame.logPlot:
3583        Title = 'log('+Title+')'
3584    Plot.set_title(Title)
3585    try:
3586        if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Image Controls',]:
3587            Page.Choice = (' key press','l: log(I) on','x: flip x','y: flip y',)
3588            if G2frame.logPlot:
3589                Page.Choice[1] = 'l: log(I) off'
3590            Page.keyPress = OnImPlotKeyPress
3591        elif G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Masks',]:
3592            Page.Choice = (' key press','l: log(I) on','s: spot mask','a: arc mask','r: ring mask',
3593                'p: polygon mask','f: frame mask',)
3594            if G2frame.logPlot:
3595                Page.Choice[1] = 'l: log(I) off'
3596            Page.keyPress = OnImPlotKeyPress
3597        elif G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Stress/Strain',]:
3598            Page.Choice = (' key press','a: add new ring',)
3599            Page.keyPress = OnImPlotKeyPress
3600    except TypeError:
3601        pass
3602    size,imagefile = G2frame.PatternTree.GetItemPyData(G2frame.Image)
3603    dark = Data.get('dark image',[0,''])
3604    if dark[0]:
3605        darkfile = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, 
3606            G2frame.root,dark[0]))[1]
3607    if imagefile != G2frame.oldImagefile:
3608        imagefile = G2IO.CheckImageFile(G2frame,imagefile)
3609        if not imagefile:
3610            G2frame.G2plotNB.Delete('2D Powder Image')
3611            return
3612        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
3613        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
3614        if dark[0]:
3615            darkImg = G2IO.GetImageData(G2frame,darkfile,imageOnly=True)
3616            G2frame.ImageZ += dark[1]*darkImg
3617        G2frame.oldImagefile = imagefile
3618
3619    imScale = 1
3620    if len(G2frame.ImageZ) > 1024:
3621        imScale = len(G2frame.ImageZ)/1024
3622    sizexy = Data['size']
3623    pixelSize = Data['pixelSize']
3624    scalex = 1000./pixelSize[0]
3625    scaley = 1000./pixelSize[1]
3626    Xmax = sizexy[0]*pixelSize[0]/1000.
3627    Ymax = sizexy[1]*pixelSize[1]/1000.
3628    xlim = (0,Xmax)
3629    ylim = (Ymax,0)
3630    Imin,Imax = Data['range'][1]
3631    acolor = mpl.cm.get_cmap(Data['color'])
3632    xcent,ycent = Data['center']
3633    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
3634    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
3635    #do threshold mask - "real" mask - others are just bondaries
3636    Zlim = Masks['Thresholds'][1]
3637    wx.BeginBusyCursor()
3638    try:
3639           
3640        if newImage:                   
3641            Imin,Imax = Data['range'][1]
3642            MA = ma.masked_greater(ma.masked_less(G2frame.ImageZ,Zlim[0]),Zlim[1])
3643            MaskA = ma.getmaskarray(MA)
3644            A = G2img.ImageCompress(MA,imScale)
3645            AM = G2img.ImageCompress(MaskA,imScale)
3646            if G2frame.logPlot:
3647                A = np.where(A>Imin,np.where(A<Imax,A,0),0)
3648                A = np.where(A>0,np.log(A),0)
3649                AM = np.where(AM>0,np.log(AM),0)
3650                Imin,Imax = [np.amin(A),np.amax(A)]
3651            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
3652                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Ymax,0])
3653            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
3654                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
3655   
3656        Plot.plot(xcent,ycent,'x')
3657        #G2frame.PatternTree.GetItemText(item)
3658        if Data['showLines']:
3659            LRAzim = Data['LRazimuth']                  #NB: integers
3660            Nazm = Data['outAzimuths']
3661            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
3662            AzmthOff = Data['azmthOff']
3663            IOtth = Data['IOtth']
3664            wave = Data['wavelength']
3665            dspI = wave/(2.0*sind(IOtth[0]/2.0))
3666            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
3667            dspO = wave/(2.0*sind(IOtth[1]/2.0))
3668            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
3669            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
3670            if ellI:
3671                xyI = []
3672                for azm in Azm:
3673                    xy = G2img.GetDetectorXY(dspI,azm,Data)
3674                    if np.any(xy):
3675                        xyI.append(xy)
3676                if len(xyI):
3677                    xyI = np.array(xyI)
3678                    arcxI,arcyI = xyI.T
3679                    Plot.plot(arcxI,arcyI,picker=3)
3680            if ellO:
3681                xyO = []
3682                for azm in Azm:
3683                    xy = G2img.GetDetectorXY(dspO,azm,Data)
3684                    if np.any(xy):
3685                        xyO.append(xy)
3686                if len(xyO):
3687                    xyO = np.array(xyO)
3688                    arcxO,arcyO = xyO.T               
3689                    Plot.plot(arcxO,arcyO,picker=3)
3690            if ellO and ellI:
3691                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
3692                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
3693            for i in range(Nazm):
3694                cake = LRAzim[0]+i*delAzm-AzmthOff
3695                if Data['centerAzm']:
3696                    cake += delAzm/2.
3697                ind = np.searchsorted(Azm,cake)
3698                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
3699                   
3700        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Image Controls':
3701            for xring,yring in Data['ring']:
3702                Plot.plot(xring,yring,'r+',picker=3)
3703            if Data['setRings']:
3704                N = 0
3705                for ring in Data['rings']:
3706                    xring,yring = np.array(ring).T[:2]
3707                    Plot.plot(xring,yring,'.',color=colors[N%6])
3708                    N += 1           
3709            for ellipse in Data['ellipses']:      #what about hyperbola?
3710                cent,phi,[width,height],col = ellipse
3711                if width > 0:       #ellipses
3712                    Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
3713                    Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
3714        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Stress/Strain':
3715            for N,ring in enumerate(StrSta['d-zero']):
3716                xring,yring = ring['ImxyObs']
3717                Plot.plot(xring,yring,colors[N%6]+'.')
3718        #masks - mask lines numbered after integration limit lines
3719        spots = Masks['Points']
3720        rings = Masks['Rings']
3721        arcs = Masks['Arcs']
3722        polygons = Masks['Polygons']
3723        if 'Frames' not in Masks:
3724            Masks['Frames'] = []
3725        frame = Masks['Frames']
3726        for spot in spots:
3727            if spot:
3728                x,y,d = spot
3729                Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
3730        G2frame.ringList = []
3731        for iring,ring in enumerate(rings):
3732            if ring:
3733                tth,thick = ring
3734                wave = Data['wavelength']
3735                xy1 = []
3736                xy2 = []
3737                Azm = np.linspace(0,362,181)
3738                for azm in Azm:
3739                    xy1.append(G2img.GetDetectorXY(Dsp(tth+thick/2.,wave),azm,Data))      #what about hyperbola
3740                    xy2.append(G2img.GetDetectorXY(Dsp(tth-thick/2.,wave),azm,Data))      #what about hyperbola
3741                x1,y1 = np.array(xy1).T
3742                x2,y2 = np.array(xy2).T
3743                G2frame.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
3744                G2frame.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
3745        G2frame.arcList = []
3746        for iarc,arc in enumerate(arcs):
3747            if arc:
3748                tth,azm,thick = arc           
3749                wave = Data['wavelength']
3750                xy1 = []
3751                xy2 = []
3752                aR = azm[0],azm[1],azm[1]-azm[0]
3753                if azm[1]-azm[0] > 180:
3754                    aR[2] /= 2
3755                Azm = np.linspace(aR[0],aR[1],aR[2])
3756                for azm in Azm:
3757                    xy1.append(G2img.GetDetectorXY(Dsp(tth+thick/2.,wave),azm,Data))      #what about hyperbola
3758                    xy2.append(G2img.GetDetectorXY(Dsp(tth-thick/2.,wave),azm,Data))      #what about hyperbola
3759                x1,y1 = np.array(xy1).T
3760                x2,y2 = np.array(xy2).T
3761                G2frame.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
3762                G2frame.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
3763                G2frame.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
3764                G2frame.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
3765        G2frame.polyList = []
3766        for ipoly,polygon in enumerate(polygons):
3767            if polygon:
3768                x,y = np.hsplit(np.array(polygon),2)
3769                G2frame.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
3770                Plot.plot(x,y,'r')           
3771        G2frame.frameList = []
3772        if frame:
3773            x,y = np.hsplit(np.array(frame),2)
3774            G2frame.frameList.append([Plot.plot(x,y,'g+',picker=10),0])
3775            Plot.plot(x,y,'g')           
3776        if newImage:
3777            colorBar = Page.figure.colorbar(Img)
3778        Plot.set_xlim(xlim)
3779        Plot.set_ylim(ylim)
3780        if Data['invert_x']:
3781            Plot.invert_xaxis()
3782        if Data['invert_y']:
3783            Plot.invert_yaxis()
3784        if not newPlot and xylim:
3785            Page.toolbar.push_current()
3786            Plot.set_xlim(xylim[0])
3787            Plot.set_ylim(xylim[1])
3788            xylim = []
3789            Page.toolbar.push_current()
3790            Page.toolbar.draw()
3791            # patch for wx 2.9 on Mac, to force a redraw
3792            i,j= wx.__version__.split('.')[0:2]
3793            if int(i)+int(j)/10. > 2.8 and 'wxOSX' in wx.PlatformInfo:
3794                Page.canvas.draw()
3795        else:
3796            Page.canvas.draw()
3797    finally:
3798        wx.EndBusyCursor()
3799       
3800################################################################################
3801##### PlotIntegration
3802################################################################################
3803           
3804def PlotIntegration(G2frame,newPlot=False,event=None):
3805    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
3806    '''
3807           
3808    def OnMotion(event):
3809        Page.canvas.SetToolTipString('')
3810        Page.canvas.SetCursor(wx.CROSS_CURSOR)
3811        azm = event.ydata
3812        tth = event.xdata
3813        if azm and tth:
3814            G2frame.G2plotNB.status.SetFields(\
3815                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
3816                               
3817    try:
3818        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
3819        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3820        if not newPlot:
3821            Plot = Page.figure.gca()          #get previous plot & get limits
3822            xylim = Plot.get_xlim(),Plot.get_ylim()
3823        Page.figure.clf()
3824        Plot = Page.figure.gca()          #get a fresh plot after clf()
3825       
3826    except ValueError:
3827        Plot = G2frame.G2plotNB.addMpl('2D Integration').gca()
3828        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
3829        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3830        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3831        Page.views = False
3832        view = False
3833    Page.Choice = None
3834    if not event:
3835        Page.SetFocus()
3836       
3837    Data = G2frame.PatternTree.GetItemPyData(
3838        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3839    image = G2frame.Integrate[0]
3840    xsc = G2frame.Integrate[1]
3841    ysc = G2frame.Integrate[2]
3842    Imin,Imax = Data['range'][1]
3843    acolor = mpl.cm.get_cmap(Data['color'])
3844    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
3845    Plot.set_ylabel('azimuth',fontsize=12)
3846    Plot.set_xlabel('2-theta',fontsize=12)
3847    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
3848        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
3849    colorBar = Page.figure.colorbar(Img)
3850#    if Data['ellipses']:           
3851#        for ellipse in Data['ellipses']:
3852#            x,y = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
3853#            tth,azm = G2img.GetTthAzm(x,y,Data)
3854##            azm = np.where(azm < 0.,azm+360,azm)
3855#            Plot.plot(tth,azm,'b,')
3856    if not newPlot:
3857        Page.toolbar.push_current()
3858        Plot.set_xlim(xylim[0])
3859        Plot.set_ylim(xylim[1])
3860        xylim = []
3861        Page.toolbar.push_current()
3862        Page.toolbar.draw()
3863    else:
3864        Page.canvas.draw()
3865               
3866################################################################################
3867##### PlotTRImage
3868################################################################################
3869           
3870def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):
3871    '''a test plot routine - not normally used
3872    ''' 
3873           
3874    def OnMotion(event):
3875        Page.canvas.SetToolTipString('')
3876        Page.canvas.SetCursor(wx.CROSS_CURSOR)
3877        azm = event.xdata
3878        tth = event.ydata
3879        if azm and tth:
3880            G2frame.G2plotNB.status.SetFields(\
3881                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
3882                               
3883    try:
3884        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
3885        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3886        if not newPlot:
3887            Plot = Page.figure.gca()          #get previous plot & get limits
3888            xylim = Plot.get_xlim(),Plot.get_ylim()
3889        Page.figure.clf()
3890        Plot = Page.figure.gca()          #get a fresh plot after clf()
3891       
3892    except ValueError:
3893        Plot = G2frame.G2plotNB.addMpl('2D Transformed Powder Image').gca()
3894        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
3895        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3896        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3897        Page.views = False
3898        view = False
3899    Page.Choice = None
3900    Page.SetFocus()
3901       
3902    Data = G2frame.PatternTree.GetItemPyData(
3903        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3904    Imin,Imax = Data['range'][1]
3905    step = (Imax-Imin)/5.
3906    V = np.arange(Imin,Imax,step)
3907    acolor = mpl.cm.get_cmap(Data['color'])
3908    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
3909    Plot.set_xlabel('azimuth',fontsize=12)
3910    Plot.set_ylabel('2-theta',fontsize=12)
3911    Plot.contour(tax,tay,taz,V,cmap=acolor)
3912    if Data['showLines']:
3913        IOtth = Data['IOtth']
3914        if Data['fullIntegrate']:
3915            LRAzim = [-180,180]
3916        else:
3917            LRAzim = Data['LRazimuth']                  #NB: integers
3918        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
3919        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
3920        if not Data['fullIntegrate']:
3921            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
3922            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
3923    if Data['setRings']:
3924        rings = np.concatenate((Data['rings']),axis=0)
3925        for xring,yring,dsp in rings:
3926            x,y = G2img.GetTthAzm(xring,yring,Data)
3927            Plot.plot(y,x,'r+')           
3928    if Data['ellipses']:           
3929        for ellipse in Data['ellipses']:
3930            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
3931            x,y = np.hsplit(ring,2)
3932            tth,azm = G2img.GetTthAzm(x,y,Data)
3933            Plot.plot(azm,tth,'b,')
3934    if not newPlot:
3935        Page.toolbar.push_current()
3936        Plot.set_xlim(xylim[0])
3937        Plot.set_ylim(xylim[1])
3938        xylim = []
3939        Page.toolbar.push_current()
3940        Page.toolbar.draw()
3941    else:
3942        Page.canvas.draw()
3943       
3944################################################################################
3945##### PlotStructure
3946################################################################################
3947           
3948def PlotStructure(G2frame,data,firstCall=False):
3949    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
3950    thermal motion ellipsoids and polyhedra
3951    '''
3952
3953    def FindPeaksBonds(XYZ):
3954        rFact = data['Drawing'].get('radiusFactor',0.85)    #data['Drawing'] could be empty!
3955        Bonds = [[] for x in XYZ]
3956        for i,xyz in enumerate(XYZ):
3957            Dx = XYZ-xyz
3958            dist = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=1))
3959            IndB = ma.nonzero(ma.masked_greater(dist,rFact*2.2))
3960            for j in IndB[0]:
3961                Bonds[i].append(Dx[j]/2.)
3962                Bonds[j].append(-Dx[j]/2.)
3963        return Bonds
3964
3965    # PlotStructure initialization here
3966    ForthirdPI = 4.0*math.pi/3.0
3967    generalData = data['General']
3968    cell = generalData['Cell'][1:7]
3969    Vol = generalData['Cell'][7:8][0]
3970    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
3971    Gmat,gmat = G2lat.cell2Gmat(cell)
3972    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3973    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3974    SGData = generalData['SGData']
3975    Mydir = generalData['Mydir']
3976    atomData = data['Atoms']
3977    mapPeaks = []
3978    drawingData = data['Drawing']
3979    if not drawingData:
3980        return          #nothing setup, nothing to draw   
3981    if 'Map Peaks' in data:
3982        mapPeaks = np.array(data['Map Peaks'])
3983        peakMax = 100.
3984        if len(mapPeaks):
3985            peakMax = np.max(mapPeaks.T[0])
3986    resRBData = data['RBModels'].get('Residue',[])
3987    vecRBData = data['RBModels'].get('Vector',[])
3988    rbAtmDict = {}
3989    for rbObj in resRBData+vecRBData:
3990        exclList = ['X' for i in range(len(rbObj['Ids']))]
3991        rbAtmDict.update(dict(zip(rbObj['Ids'],exclList)))
3992    testRBObj = data.get('testRBObj',{})
3993    rbObj = testRBObj.get('rbObj',{})
3994    MCSA = data.get('MCSA',{})
3995    mcsaModels = MCSA.get('Models',[])
3996    if len(mcsaModels) > 1:
3997        XYZs,Types = G2mth.UpdateMCSAxyz(Bmat,MCSA)
3998        mcsaXYZ = []
3999        mcsaTypes = []
4000        neqv = 0
4001        for xyz,atyp in zip(XYZs,Types):
4002            equiv = G2spc.GenAtom(xyz,SGData,All=True,Move=False)
4003            neqv = max(neqv,len(equiv))
4004            for item in equiv:
4005                mcsaXYZ.append(item[0]) 
4006                mcsaTypes.append(atyp)
4007        mcsaXYZ = np.array(mcsaXYZ)
4008        mcsaTypes = np.array(mcsaTypes)
4009        nuniq = mcsaXYZ.shape[0]/neqv
4010        mcsaXYZ = np.reshape(mcsaXYZ,(nuniq,neqv,3))
4011        mcsaTypes = np.reshape(mcsaTypes,(nuniq,neqv))
4012        cent = np.fix(np.sum(mcsaXYZ+2.,axis=0)/nuniq)-2
4013        cent[0] = [0,0,0]   #make sure 1st one isn't moved
4014        mcsaXYZ = np.swapaxes(mcsaXYZ,0,1)-cent[:,np.newaxis,:]
4015        mcsaTypes = np.swapaxes(mcsaTypes,0,1)
4016        mcsaXYZ = np.reshape(mcsaXYZ,(nuniq*neqv,3))
4017        mcsaTypes = np.reshape(mcsaTypes,(nuniq*neqv))                       
4018        mcsaBonds = FindPeaksBonds(mcsaXYZ)       
4019    drawAtoms = drawingData.get('Atoms',[])
4020    mapData = {}
4021    flipData = {}
4022    rhoXYZ = []
4023    showBonds = False
4024    if 'Map' in generalData:
4025        mapData = generalData['Map']
4026        showBonds = mapData.get('Show bonds',False)
4027    if 'Flip' in generalData:
4028        flipData = generalData['Flip']                       
4029        flipData['mapRoll'] = [0,0,0]
4030    Wt = np.array([255,255,255])
4031    Rd = np.array([255,0,0])
4032    Gr = np.array([0,255,0])
4033    wxGreen = wx.Colour(0,255,0)
4034    Bl = np.array([0,0,255])
4035    Or = np.array([255,128,0])
4036    wxOrange = wx.Colour(255,128,0)
4037    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]])
4038    uEdges = np.array([
4039        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
4040        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
4041        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
4042    mD = 0.1
4043    mV = np.array([[[-mD,0,0],[mD,0,0]],[[0,-mD,0],[0,mD,0]],[[0,0,-mD],[0,0,mD]]])
4044    mapPeakVecs = np.inner(mV,Bmat)
4045
4046    backColor = np.array(list(drawingData['backColor'])+[0,])
4047    Bc = np.array(list(drawingData['backColor']))
4048    uColors = [Rd,Gr,Bl,Wt-Bc, Wt-Bc,Wt-Bc,Wt-Bc,Wt-Bc, Wt-Bc,Wt-Bc,Wt-Bc,Wt-Bc]
4049    altDown = False
4050    shiftDown = False
4051    ctrlDown = False
4052   
4053    def OnKeyBox(event):
4054#        Draw()                          #make sure plot is fresh!!
4055        mode = cb.GetValue()
4056        if mode in ['jpeg','bmp','tiff',]:
4057            try:
4058                import Image as Im
4059            except ImportError:
4060                try:
4061                    from PIL import Image as Im
4062                except ImportError:
4063                    print "PIL/pillow Image module not present. Cannot save images without this"
4064                    raise Exception("PIL/pillow Image module not found")
4065            Fname = os.path.join(Mydir,generalData['Name']+'.'+mode)
4066            size = Page.canvas.GetSize()
4067            glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
4068            if mode in ['jpeg',]:
4069                Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
4070                im = Im.new("RGBA", (size[0],size[1]))
4071            else:
4072                Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
4073                im = Im.new("RGB", (size[0],size[1]))
4074            im.fromstring(Pix)
4075            im.save(Fname,mode)
4076            cb.SetValue(' save as/key:')
4077            G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
4078        else:
4079            event.key = cb.GetValue()[0]
4080            cb.SetValue(' save as/key:')
4081            wx.CallAfter(OnKey,event)
4082        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
4083
4084    def OnKey(event):           #on key UP!!
4085#        Draw()                          #make sure plot is fresh!!
4086        try:
4087            keyCode = event.GetKeyCode()
4088            if keyCode > 255:
4089                keyCode = 0
4090            key = chr(keyCode)
4091        except AttributeError:       #if from OnKeyBox above
4092            key = str(event.key).upper()
4093        indx = drawingData['selectedAtoms']
4094        cx,ct = drawingData['atomPtrs'][:2]
4095        if key in ['C']:
4096            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
4097            drawingData['viewDir'] = [0,0,1]
4098            drawingData['oldxy'] = []
4099            V0 = np.array([0,0,1])
4100            V = np.inner(Amat,V0)
4101            V /= np.sqrt(np.sum(V**2))
4102            A = np.arccos(np.sum(V*V0))
4103            Q = G2mth.AV2Q(A,[0,1,0])
4104            drawingData['Quaternion'] = Q
4105            SetViewPointText(drawingData['viewPoint'][0])
4106            SetViewDirText(drawingData['viewDir'])
4107            Q = drawingData['Quaternion']
4108            G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4109        elif key in ['N']:
4110            drawAtoms = drawingData['Atoms']
4111            if not len(drawAtoms):      #no atoms
4112                return
4113            pI = drawingData['viewPoint'][1]
4114            if not len(pI):
4115                pI = [0,0]
4116            if indx:
4117                pI[0] = indx[pI[1]]
4118                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
4119                pI[1] += 1
4120                if pI[1] >= len(indx):
4121                    pI[1] = 0
4122            else:
4123                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
4124                pI[0] += 1
4125                if pI[0] >= len(drawAtoms):
4126                    pI[0] = 0
4127            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
4128            SetViewPointText(drawingData['viewPoint'][0])
4129            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
4130               
4131        elif key in ['P']:
4132            drawAtoms = drawingData['Atoms']
4133            if not len(drawAtoms):      #no atoms
4134                return
4135            pI = drawingData['viewPoint'][1]
4136            if not len(pI):
4137                pI = [0,0]
4138            if indx:
4139                pI[0] = indx[pI[1]]
4140                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
4141                pI[1] -= 1
4142                if pI[1] < 0:
4143                    pI[1] = len(indx)-1
4144            else:
4145                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
4146                pI[0] -= 1
4147                if pI[0] < 0:
4148                    pI[0] = len(drawAtoms)-1
4149            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
4150            SetViewPointText(drawingData['viewPoint'][0])           
4151            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
4152        elif key in ['U','D','L','R'] and mapData['Flip'] == True:
4153            dirDict = {'U':[0,1],'D':[0,-1],'L':[-1,0],'R':[1,0]}
4154            SetMapRoll(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 = []
4324        if page:
4325            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
4326