source: trunk/GSASIIplot.py @ 1786

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

fix anomalous dispersion calcs.
fix SC weights on derivatives

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