source: trunk/GSASIIplot.py @ 1635

Last change on this file since 1635 was 1635, checked in by vondreele, 8 years ago

modulation functions & constraints

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