source: trunk/GSASIIplot.py @ 1604

Last change on this file since 1604 was 1604, checked in by vondreele, 9 years ago

make modulation wave maps
fix atom index bugs
begin modulated structure imput to least squares

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