source: trunk/GSASIIplot.py @ 1628

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

fix problem with U/D offset of SASD plots; L/R worked OK.

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