source: trunk/GSASIIplot.py @ 1779

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

Add FWHM to exported peak lists
change default F2 refinement to False
Add histogram name to axial dist plots from DData

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