source: trunk/GSASIIplot.py @ 1583

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

cleanup masks before integration; add zoom buttons to mpl window; wx2.9+ fix for seq ref window

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