source: trunk/GSASIIplot.py @ 1625

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

add a parameter to result from G2stIO.GetPhaseData?
add modulation functions to G2Math
add modulation names to G2obj
implement various wave types for modulations
plot position modulation wave function on contour plots
implement shift of modulation plot by +/-/0 keys
temporarily default G2spc.GetSSfxuinel to 1,-1 site symm.
move maxSSwave dict out of parmDict - now in controlDict
implement import of Sawtooth parms from J2K files.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 221.1 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIplot: plotting routines*
4===============================
5
6'''
7########### SVN repository information ###################
8# $Date: 2015-01-08 22:24:32 +0000 (Thu, 08 Jan 2015) $
9# $Author: vondreele $
10# $Revision: 1625 $
11# $URL: trunk/GSASIIplot.py $
12# $Id: GSASIIplot.py 1625 2015-01-08 22:24:32Z 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: 1625 $")
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
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.005)**(offset*N)*X**4
1491            else:
1492                Y = xye[1]*Sample['Scale'][0]*(1.005)**(offset*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        wx.CallAfter(ModulationPlot,G2frame,data,Atom,Ax,Off)
2802
2803    try:
2804        plotNum = G2frame.G2plotNB.plotList.index('Modulation')
2805        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2806        Page.figure.clf()
2807        Plot = Page.figure.gca()
2808        if not Page.IsShown():
2809            Page.Show()
2810    except ValueError:
2811        Plot = G2frame.G2plotNB.addMpl('Modulation').gca()
2812        plotNum = G2frame.G2plotNB.plotList.index('Modulation')
2813        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2814        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2815        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
2816       
2817    Page.Choice = ['+: shift up','-: shift down','0: reset']
2818    Page.keyPress = OnPlotKeyPress
2819    Page.SetFocus()
2820    General = data['General']
2821    cx,ct,cs,cia = General['AtomPtrs']
2822    Map = General['4DmapData']
2823    MapType = Map['MapType']
2824    rhoSize = np.array(Map['rho'].shape)
2825    atxyz = np.array(atom[cx:cx+3])
2826    waveType = atom[-1]['SS1']['waveType']
2827    Spos = atom[-1]['SS1']['Spos']
2828    tau = np.linspace(0.,2.,101)
2829    wave = np.zeros((3,101))
2830    if len(Spos):
2831        scof = []
2832        ccof = []
2833        for i,spos in enumerate(Spos):
2834            if waveType in ['Sawtooth','ZigZag'] and not i:
2835                Toff = spos[0][0]
2836                slopes = np.array(spos[0][1:])
2837                if waveType == 'Sawtooth':
2838                    wave = G2mth.posSawtooth(tau,Toff,slopes)
2839                elif waveType == 'ZigZag':
2840                    wave = G2mth.posZigZag(tau,Toff,slopes)
2841            else:
2842                scof.append(spos[0][:3])
2843                ccof.append(spos[0][3:])
2844        wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof))
2845    Title = MapType+' modulation map for atom '+atom[0]+    \
2846        ' at %.4f %.4f %.4f'%(atxyz[0],atxyz[1],atxyz[2])
2847    ix = -np.array(np.rint(rhoSize[:3]*atxyz),dtype='i')
2848    ix += (rhoSize[:3]/2)
2849    ix = ix%rhoSize[:3]
2850    rho = np.roll(np.roll(np.roll(Map['rho'],ix[0],axis=0),ix[1],axis=1),ix[2],axis=2)
2851    ix = rhoSize[:3]/2
2852    ib = 4
2853    if Ax == 'x':
2854        slab = np.sum(np.sum(rho[:,ix[1]-ib:ix[1]+ib,ix[2]-ib:ix[2]+ib,:],axis=2),axis=1)
2855        Plot.plot(wave[0],tau)
2856    elif Ax == 'y':
2857        slab = np.sum(np.sum(rho[ix[0]-ib:ix[0]+ib,:,ix[2]-ib:ix[2]+ib,:],axis=2),axis=0)
2858        Plot.plot(wave[1],tau)
2859    elif Ax == 'z':
2860        slab = np.sum(np.sum(rho[ix[0]-ib:ix[0]+ib,ix[1]-ib:ix[1]+ib,:,:],axis=1),axis=0)
2861        Plot.plot(wave[2],tau)
2862    Plot.set_title(Title)
2863    Plot.set_xlabel('t')
2864    Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax))
2865    Slab = np.concatenate((slab,slab),axis=1).T
2866    Plot.contour(Slab,20,extent=(-.5+Off*.005,.5+Off*.005,0.,2.))
2867   
2868    Page.canvas.draw()
2869   
2870################################################################################
2871##### PlotCovariance
2872################################################################################
2873           
2874def PlotCovariance(G2frame,Data):
2875    'needs a doc string'
2876    if not Data:
2877        print 'No covariance matrix available'
2878        return
2879    varyList = Data['varyList']
2880    values = Data['variables']
2881    Xmax = len(varyList)
2882    covMatrix = Data['covMatrix']
2883    sig = np.sqrt(np.diag(covMatrix))
2884    xvar = np.outer(sig,np.ones_like(sig))
2885    covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
2886    title = ' for\n'+Data['title']
2887    newAtomDict = Data.get('newAtomDict',{})
2888   
2889
2890    def OnPlotKeyPress(event):
2891        newPlot = False
2892        if event.key == 's':
2893            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
2894            choice.sort()
2895            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
2896            if dlg.ShowModal() == wx.ID_OK:
2897                sel = dlg.GetSelection()
2898                G2frame.VcovColor = choice[sel]
2899            else:
2900                G2frame.VcovColor = 'RdYlGn'
2901            dlg.Destroy()
2902        PlotCovariance(G2frame,Data)
2903
2904    def OnMotion(event):
2905        #there is a problem here - reports wrong values
2906        if event.button:
2907            ytics = imgAx.get_yticks()
2908            ytics = np.where(ytics<len(varyList),ytics,-1)
2909            ylabs = [np.where(0<=i ,varyList[int(i)],' ') for i in ytics]
2910            imgAx.set_yticklabels(ylabs)           
2911        if event.xdata and event.ydata:                 #avoid out of frame errors
2912            xpos = int(event.xdata+.5)
2913            ypos = int(event.ydata+.5)
2914            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
2915                if xpos == ypos:
2916                    value = values[xpos]
2917                    name = varyList[xpos]
2918                    if varyList[xpos] in newAtomDict:
2919                        name,value = newAtomDict[name]                       
2920                    msg = '%s value = %.4g, esd = %.4g'%(name,value,sig[xpos])
2921                else:
2922                    msg = '%s - %s: %5.3f'%(varyList[xpos],varyList[ypos],covArray[xpos][ypos])
2923                Page.canvas.SetToolTipString(msg)
2924                G2frame.G2plotNB.status.SetFields(['',msg])
2925               
2926    try:
2927        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
2928        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2929        Page.figure.clf()
2930        Plot = Page.figure.gca()
2931        if not Page.IsShown():
2932            Page.Show()
2933    except ValueError:
2934        Plot = G2frame.G2plotNB.addMpl('Covariance').gca()
2935        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
2936        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2937        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2938        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
2939    Page.Choice = ['s: to change colors']
2940    Page.keyPress = OnPlotKeyPress
2941    Page.SetFocus()
2942    G2frame.G2plotNB.status.SetFields(['',''])   
2943    acolor = mpl.cm.get_cmap(G2frame.VcovColor)
2944    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower',
2945        vmin=-1.,vmax=1.)
2946    imgAx = Img.get_axes()
2947    ytics = imgAx.get_yticks()
2948    ylabs = [varyList[int(i)] for i in ytics[:-1]]
2949    imgAx.set_yticklabels(ylabs)
2950    colorBar = Page.figure.colorbar(Img)
2951    Plot.set_title('V-Cov matrix'+title)
2952    Plot.set_xlabel('Variable number')
2953    Plot.set_ylabel('Variable name')
2954    Page.canvas.draw()
2955   
2956################################################################################
2957##### PlotTorsion
2958################################################################################
2959
2960def PlotTorsion(G2frame,phaseName,Torsion,TorName,Names=[],Angles=[],Coeff=[]):
2961    'needs a doc string'
2962   
2963    global names
2964    names = Names
2965    sum = np.sum(Torsion)
2966    torsion = np.log(2*Torsion+1.)/sum
2967    tMin = np.min(torsion)
2968    tMax = np.max(torsion)
2969    torsion = 3.*(torsion-tMin)/(tMax-tMin)
2970    X = np.linspace(0.,360.,num=45)
2971   
2972    def OnPick(event):
2973        ind = event.ind[0]
2974        msg = 'atoms:'+names[ind]
2975        Page.canvas.SetToolTipString(msg)
2976        try:
2977            page = G2frame.dataDisplay.GetSelection()
2978        except:
2979            return
2980        if G2frame.dataDisplay.GetPageText(page) == 'Torsion restraints':
2981            torGrid = G2frame.dataDisplay.GetPage(page).Torsions
2982            torGrid.ClearSelection()
2983            for row in range(torGrid.GetNumberRows()):
2984                if names[ind] in torGrid.GetCellValue(row,0):
2985                    torGrid.SelectRow(row)
2986            torGrid.ForceRefresh()
2987               
2988    def OnMotion(event):
2989        if event.xdata and event.ydata:                 #avoid out of frame errors
2990            xpos = event.xdata
2991            ypos = event.ydata
2992            msg = 'torsion,energy: %5.3f %5.3f'%(xpos,ypos)
2993            Page.canvas.SetToolTipString(msg)
2994
2995    try:
2996        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
2997        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2998        Page.figure.clf()
2999        Plot = Page.figure.gca()
3000        if not Page.IsShown():
3001            Page.Show()
3002    except ValueError:
3003        Plot = G2frame.G2plotNB.addMpl('Torsion').gca()
3004        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
3005        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3006        Page.canvas.mpl_connect('pick_event', OnPick)
3007        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3008   
3009    Page.SetFocus()
3010    G2frame.G2plotNB.status.SetFields(['','Use mouse LB to identify torsion atoms'])
3011    Plot.plot(X,torsion,'b+')
3012    if len(Coeff):
3013        X2 = np.linspace(0.,360.,45)
3014        Y2 = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in X2])
3015        Plot.plot(X2,Y2,'r')
3016    if len(Angles):
3017        Eval = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in Angles])
3018        Plot.plot(Angles,Eval,'ro',picker=5)
3019    Plot.set_xlim((0.,360.))
3020    Plot.set_title('Torsion angles for '+TorName+' in '+phaseName)
3021    Plot.set_xlabel('angle',fontsize=16)
3022    Plot.set_ylabel('Energy',fontsize=16)
3023    Page.canvas.draw()
3024   
3025################################################################################
3026##### PlotRama
3027################################################################################
3028
3029def PlotRama(G2frame,phaseName,Rama,RamaName,Names=[],PhiPsi=[],Coeff=[]):
3030    'needs a doc string'
3031
3032    global names
3033    names = Names
3034    rama = np.log(2*Rama+1.)
3035    ramaMax = np.max(rama)
3036    rama = np.reshape(rama,(45,45))
3037    global Phi,Psi
3038    Phi = []
3039    Psi = []
3040
3041    def OnPlotKeyPress(event):
3042        newPlot = False
3043        if event.key == 's':
3044            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
3045            choice.sort()
3046            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
3047            if dlg.ShowModal() == wx.ID_OK:
3048                sel = dlg.GetSelection()
3049                G2frame.RamaColor = choice[sel]
3050            else:
3051                G2frame.RamaColor = 'RdYlGn'
3052            dlg.Destroy()
3053        PlotRama(G2frame,phaseName,Rama)
3054
3055    def OnPick(event):
3056        ind = event.ind[0]
3057        msg = 'atoms:'+names[ind]
3058        Page.canvas.SetToolTipString(msg)
3059        try:
3060            page = G2frame.dataDisplay.GetSelection()
3061        except:
3062            return
3063        if G2frame.dataDisplay.GetPageText(page) == 'Ramachandran restraints':
3064            ramaGrid = G2frame.dataDisplay.GetPage(page).Ramas
3065            ramaGrid.ClearSelection()
3066            for row in range(ramaGrid.GetNumberRows()):
3067                if names[ind] in ramaGrid.GetCellValue(row,0):
3068                    ramaGrid.SelectRow(row)
3069            ramaGrid.ForceRefresh()
3070
3071    def OnMotion(event):
3072        if event.xdata and event.ydata:                 #avoid out of frame errors
3073            xpos = event.xdata
3074            ypos = event.ydata
3075            msg = 'phi/psi: %5.3f %5.3f'%(xpos,ypos)
3076            Page.canvas.SetToolTipString(msg)
3077           
3078    try:
3079        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
3080        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3081        Page.figure.clf()
3082        Plot = Page.figure.gca()
3083        if not Page.IsShown():
3084            Page.Show()
3085    except ValueError:
3086        Plot = G2frame.G2plotNB.addMpl('Ramachandran').gca()
3087        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
3088        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3089        Page.canvas.mpl_connect('pick_event', OnPick)
3090        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3091        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
3092
3093    Page.Choice = ['s: to change colors']
3094    Page.keyPress = OnPlotKeyPress
3095    Page.SetFocus()
3096    G2frame.G2plotNB.status.SetFields(['','Use mouse LB to identify phi/psi atoms'])
3097    acolor = mpl.cm.get_cmap(G2frame.RamaColor)
3098    if RamaName == 'All' or '-1' in RamaName:
3099        if len(Coeff): 
3100            X,Y = np.meshgrid(np.linspace(-180.,180.,45),np.linspace(-180.,180.,45))
3101            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
3102            Plot.contour(X,Y,np.reshape(Z,(45,45)))
3103        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
3104            extent=[-180,180,-180,180],origin='lower')
3105        if len(PhiPsi):
3106            Phi,Psi = PhiPsi.T
3107            Phi = np.where(Phi>180.,Phi-360.,Phi)
3108            Psi = np.where(Psi>180.,Psi-360.,Psi)
3109            Plot.plot(Phi,Psi,'ro',picker=5)
3110        Plot.set_xlim((-180.,180.))
3111        Plot.set_ylim((-180.,180.))
3112    else:
3113        if len(Coeff): 
3114            X,Y = np.meshgrid(np.linspace(0.,360.,45),np.linspace(0.,360.,45))
3115            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
3116            Plot.contour(X,Y,np.reshape(Z,(45,45)))
3117        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
3118            extent=[0,360,0,360],origin='lower')
3119        if len(PhiPsi):
3120            Phi,Psi = PhiPsi.T
3121            Plot.plot(Phi,Psi,'ro',picker=5)
3122        Plot.set_xlim((0.,360.))
3123        Plot.set_ylim((0.,360.))
3124    Plot.set_title('Ramachandran for '+RamaName+' in '+phaseName)
3125    Plot.set_xlabel(r'$\phi$',fontsize=16)
3126    Plot.set_ylabel(r'$\psi$',fontsize=16)
3127    colorBar = Page.figure.colorbar(Img)
3128    Page.canvas.draw()
3129
3130
3131################################################################################
3132##### PlotSeq
3133################################################################################
3134def PlotSelectedSequence(G2frame,ColumnList,TableGet,SelectX,fitnum=None,fitvals=None):
3135    '''Plot a result from a sequential refinement
3136
3137    :param wx.Frame G2frame: The main GSAS-II tree "window"
3138    :param list ColumnList: list of int values corresponding to columns
3139      selected as y values
3140    :param function TableGet: a function that takes a column number
3141      as argument and returns the column label, the values and there ESDs (or None)
3142    :param function SelectX: a function that returns a selected column
3143      number (or None) as the X-axis selection
3144    '''
3145    global Title,xLabel,yLabel
3146    xLabel = yLabel = Title = ''
3147    def OnMotion(event):
3148        if event.xdata and event.ydata:                 #avoid out of frame errors
3149            xpos = event.xdata
3150            ypos = event.ydata
3151            msg = '%5.3f %.6g'%(xpos,ypos)
3152            Page.canvas.SetToolTipString(msg)
3153
3154    def OnKeyPress(event):
3155        global Title,xLabel,yLabel
3156        if event.key == 's':
3157            G2frame.seqXaxis = G2frame.seqXselect()
3158            Draw()
3159        elif event.key == 't':
3160            dlg = G2gd.MultiStringDialog(G2frame,'Set titles & labels',[' Title ',' x-Label ',' y-Label '],
3161                [Title,xLabel,yLabel])
3162            if dlg.Show():
3163                Title,xLabel,yLabel = dlg.GetValues()
3164            dlg.Destroy()
3165            Draw()
3166           
3167    def Draw():
3168        global Title,xLabel,yLabel
3169        Page.SetFocus()
3170        G2frame.G2plotNB.status.SetStatusText('press s to select X axis, t to change titles',1)
3171        Plot.clear()
3172        if G2frame.seqXaxis is not None:   
3173            xName,X,Xsig = Page.seqTableGet(G2frame.seqXaxis)
3174        else:
3175            X = np.arange(0,G2frame.SeqTable.GetNumberRows(),1)
3176            xName = 'Data sequence number'
3177        for col in Page.seqYaxisList:
3178            name,Y,sig = Page.seqTableGet(col)
3179            # deal with missing (None) values
3180            Xnew = []
3181            Ynew = []
3182            Ysnew = []
3183            for i in range(len(X)):
3184                if X[i] is None or Y[i] is None: continue
3185                Xnew.append(X[i])
3186                Ynew.append(Y[i])
3187                if sig: Ysnew.append(sig[i])
3188            if Ysnew:
3189                if G2frame.seqReverse and not G2frame.seqXaxis:
3190                    Ynew = Ynew[::-1]
3191                    Ysnew = Ysnew[::-1]
3192                Plot.errorbar(Xnew,Ynew,yerr=Ysnew,label=name)
3193            else:
3194                if G2frame.seqReverse and not G2frame.seqXaxis:
3195                    Ynew = Ynew[::-1]
3196                Plot.plot(Xnew,Ynew)
3197                Plot.plot(Xnew,Ynew,'o',label=name)
3198        if Page.fitvals: # TODO: deal with fitting of None values
3199            if G2frame.seqReverse and not G2frame.seqXaxis:
3200                Page.fitvals = Page.fitvals[::-1]
3201            Plot.plot(X,Page.fitvals,label='Fit')
3202           
3203        Plot.legend(loc='best')
3204        print Title,xLabel,yLabel
3205        if Title:
3206            Plot.set_title(Title)
3207        else:
3208            Plot.set_title('')
3209        if xLabel:
3210            Plot.set_xlabel(xLabel)
3211        else:
3212            Plot.set_xlabel(xName)
3213        if yLabel:
3214            Plot.set_ylabel(yLabel)
3215        else:
3216            Plot.set_ylabel('Parameter values')
3217        Page.canvas.draw()           
3218           
3219    G2frame.seqXselect = SelectX
3220    try:
3221        G2frame.seqXaxis
3222    except:
3223        G2frame.seqXaxis = None
3224
3225    if fitnum is None:
3226        label = 'Sequential refinement'
3227    else:
3228        label = 'Parametric fit #'+str(fitnum+1)
3229    try:
3230        plotNum = G2frame.G2plotNB.plotList.index(label)
3231        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3232        Page.figure.clf()
3233        Plot = Page.figure.gca()
3234        if not Page.IsShown():
3235            Page.Show()
3236    except ValueError:
3237        Plot = G2frame.G2plotNB.addMpl(label).gca()
3238        plotNum = G2frame.G2plotNB.plotList.index(label)
3239        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3240        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
3241        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3242    Page.Choice = ['s - select x-axis','t - change titles',]
3243    Page.keyPress = OnKeyPress
3244    Page.seqYaxisList = ColumnList
3245    Page.seqTableGet = TableGet
3246    Page.fitvals = fitvals
3247       
3248    Draw()
3249    G2frame.G2plotNB.nb.SetSelection(plotNum) # raises plot tab
3250               
3251################################################################################
3252##### PlotExposedImage & PlotImage
3253################################################################################
3254           
3255def PlotExposedImage(G2frame,newPlot=False,event=None):
3256    '''General access module for 2D image plotting
3257    '''
3258    plotNo = G2frame.G2plotNB.nb.GetSelection()
3259    if G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
3260        PlotImage(G2frame,newPlot,event,newImage=True)
3261    elif G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
3262        PlotIntegration(G2frame,newPlot,event)
3263
3264def OnStartMask(G2frame):
3265    '''Initiate the start of a Frame or Polygon map
3266
3267    :param wx.Frame G2frame: The main GSAS-II tree "window"
3268    :param str eventkey: a single letter ('f' or 'p') that
3269      determines what type of mask is created.   
3270    '''
3271    Masks = G2frame.PatternTree.GetItemPyData(
3272        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
3273    if G2frame.MaskKey == 'f':
3274        Masks['Frames'] = []
3275    elif G2frame.MaskKey == 'p':
3276        Masks['Polygons'].append([])
3277    elif G2frame.MaskKey == 's':
3278        Masks['Points'].append([])
3279    elif G2frame.MaskKey == 'a':
3280        Masks['Arcs'].append([])
3281    elif G2frame.MaskKey == 'r':
3282        Masks['Rings'].append([])
3283    G2imG.UpdateMasks(G2frame,Masks)
3284    PlotImage(G2frame,newImage=True)
3285   
3286def OnStartNewDzero(G2frame):
3287    '''Initiate the start of adding a new d-zero to a strain data set
3288
3289    :param wx.Frame G2frame: The main GSAS-II tree "window"
3290    :param str eventkey: a single letter ('a') that
3291      triggers the addition of a d-zero.   
3292    '''
3293    G2frame.dataFrame.GetStatusBar().SetStatusText('Add strain ring active - LB pick d-zero value',0)
3294    G2frame.PickId = G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain')
3295    data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
3296    return data
3297
3298def PlotImage(G2frame,newPlot=False,event=None,newImage=True):
3299    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
3300    masks, etc.
3301    '''
3302    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
3303    import numpy.ma as ma
3304    Dsp = lambda tth,wave: wave/(2.*npsind(tth/2.))
3305    global Data,Masks,StrSta
3306    colors=['b','g','r','c','m','k']
3307    Data = G2frame.PatternTree.GetItemPyData(
3308        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3309# patch
3310    if 'invert_x' not in Data:
3311        Data['invert_x'] = False
3312        Data['invert_y'] = True
3313# end patch
3314    Masks = G2frame.PatternTree.GetItemPyData(
3315        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
3316    try:    #may be absent
3317        StrSta = G2frame.PatternTree.GetItemPyData(
3318            G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain'))
3319    except TypeError:   #is missing
3320        StrSta = {}
3321
3322    def OnImMotion(event):
3323        Page.canvas.SetToolTipString('')
3324        sizexy = Data['size']
3325        if event.xdata and event.ydata and len(G2frame.ImageZ):                 #avoid out of frame errors
3326            Page.canvas.SetToolTipString('%8.2f %8.2fmm'%(event.xdata,event.ydata))
3327            Page.canvas.SetCursor(wx.CROSS_CURSOR)
3328            item = G2frame.itemPicked
3329            pixelSize = Data['pixelSize']
3330            scalex = 1000./pixelSize[0]
3331            scaley = 1000./pixelSize[1]
3332            if item and G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Image Controls':
3333                if 'Text' in str(item):
3334                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
3335                else:
3336                    xcent,ycent = Data['center']
3337                    xpos = event.xdata-xcent
3338                    ypos = event.ydata-ycent
3339                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
3340                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
3341                        Page.canvas.SetToolTipString('%6d deg'%(azm))
3342                    elif 'line1' in  str(item) or 'line2' in str(item):
3343                        Page.canvas.SetToolTipString('%8.3f deg'%(tth))                           
3344            else:
3345                xpos = event.xdata
3346                ypos = event.ydata
3347                xpix = xpos*scalex
3348                ypix = ypos*scaley
3349                Int = 0
3350                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
3351                    Int = G2frame.ImageZ[ypix][xpix]
3352                tth,azm,D,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
3353                Q = 2.*math.pi/dsp
3354                fields = ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)]
3355                if G2frame.MaskKey in ['p','f']:
3356                    fields[1] = 'Polygon/frame mask pick - LB next point, RB close polygon'
3357                elif G2frame.StrainKey:
3358                    fields[0] = 'd-zero pick active'
3359                G2frame.G2plotNB.status.SetFields(fields)
3360
3361    def OnImPlotKeyPress(event):
3362        try:
3363            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
3364        except TypeError:
3365            return
3366        if PickName == 'Masks':
3367            if event.key in ['l','p','f','s','a','r']:
3368                G2frame.MaskKey = event.key
3369                OnStartMask(G2frame)
3370                PlotImage(G2frame,newPlot=False)
3371               
3372        elif PickName == 'Stress/Strain':
3373            if event.key in ['a',]:
3374                G2frame.StrainKey = event.key
3375                StrSta = OnStartNewDzero(G2frame)
3376                PlotImage(G2frame,newPlot=False)
3377               
3378        elif PickName == 'Image Controls':
3379            if event.key in ['c',]:
3380                Xpos = event.xdata
3381                if not Xpos:            #got point out of frame
3382                    return
3383                Ypos = event.ydata
3384                dlg = wx.MessageDialog(G2frame,'Are you sure you want to change the center?',
3385                    'Center change',style=wx.OK|wx.CANCEL)
3386                try:
3387                    if dlg.ShowModal() == wx.ID_OK:
3388                        print 'move center to: ',Xpos,Ypos
3389                        Data['center'] = [Xpos,Ypos]
3390                        G2imG.UpdateImageControls(G2frame,Data,Masks)
3391                        PlotImage(G2frame,newPlot=False)
3392                finally:
3393                    dlg.Destroy()
3394                return
3395            elif event.key == 'l':
3396                G2frame.logPlot = not G2frame.logPlot
3397            elif event.key in ['x',]:
3398                Data['invert_x'] = not Data['invert_x']
3399            elif event.key in ['y',]:
3400                Data['invert_y'] = not Data['invert_y']
3401            PlotImage(G2frame,newPlot=True)
3402           
3403    def OnKeyBox(event):
3404        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):
3405            event.key = cb.GetValue()[0]
3406            cb.SetValue(' key press')
3407            if event.key in ['l','s','a','r','p','x','y']:
3408                wx.CallAfter(OnImPlotKeyPress,event)
3409        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
3410                       
3411    def OnImPick(event):
3412        if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']:
3413            return
3414        if G2frame.itemPicked is not None: return
3415        G2frame.itemPicked = event.artist
3416        G2frame.mousePicked = event.mouseevent
3417       
3418    def OnImRelease(event):
3419        try:
3420            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
3421        except TypeError:
3422            return
3423        if PickName not in ['Image Controls','Masks','Stress/Strain']:
3424            return
3425        pixelSize = Data['pixelSize']
3426        scalex = 1000./pixelSize[0]
3427        scaley = 1000./pixelSize[1]
3428#        pixLimit = Data['pixLimit']    #can be too tight
3429        pixLimit = 20       #this makes the search box 40x40 pixels
3430        if G2frame.itemPicked is None and PickName == 'Image Controls' and len(G2frame.ImageZ):
3431            Xpos = event.xdata
3432            if not (Xpos and G2frame.ifGetRing):                   #got point out of frame
3433                return
3434            Ypos = event.ydata
3435            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
3436                if event.button == 1:
3437                    Xpix = Xpos*scalex
3438                    Ypix = Ypos*scaley
3439                    xpos,ypos,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
3440                    if I and J:
3441                        xpos += .5                              #shift to pixel center
3442                        ypos += .5
3443                        xpos /= scalex                          #convert to mm
3444                        ypos /= scaley
3445                        Data['ring'].append([xpos,ypos])
3446                elif event.button == 3:
3447                    G2frame.dataFrame.GetStatusBar().SetStatusText('Calibrating...',0)
3448                    if G2img.ImageCalibrate(G2frame,Data):
3449                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check',0)
3450                        print 'Calibration successful'
3451                    else:
3452                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose',0)
3453                        print 'Calibration failed'
3454                    G2frame.ifGetRing = False
3455                    G2imG.UpdateImageControls(G2frame,Data,Masks)
3456                    return
3457                PlotImage(G2frame,newImage=False)
3458            return
3459        elif G2frame.MaskKey and PickName == 'Masks':
3460            Xpos,Ypos = [event.xdata,event.ydata]
3461            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3462                return
3463            if G2frame.MaskKey == 's' and event.button == 1:
3464                Masks['Points'][-1] = [Xpos,Ypos,1.]
3465                G2frame.MaskKey = ''               
3466            elif G2frame.MaskKey == 'r' and event.button == 1:
3467                tth = G2img.GetTth(Xpos,Ypos,Data)
3468                Masks['Rings'][-1] = [tth,0.1]
3469                G2frame.MaskKey = ''               
3470            elif G2frame.MaskKey == 'a' and event.button == 1:
3471                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
3472                azm = int(azm)               
3473                Masks['Arcs'][-1] = [tth,[azm-5,azm+5],0.1]
3474                G2frame.MaskKey = ''               
3475            elif G2frame.MaskKey =='p':
3476                polygon = Masks['Polygons'][-1]
3477                if len(polygon) > 2 and event.button == 3:
3478                    x0,y0 = polygon[0]
3479                    polygon.append([x0,y0])
3480                    G2frame.MaskKey = ''
3481                    G2frame.G2plotNB.status.SetFields(['','Polygon closed'])
3482                else:
3483                    G2frame.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(Xpos,Ypos)])
3484                    polygon.append([Xpos,Ypos])
3485            elif G2frame.MaskKey =='f':
3486                frame = Masks['Frames']
3487                if len(frame) > 2 and event.button == 3:
3488                    x0,y0 = frame[0]
3489                    frame.append([x0,y0])
3490                    G2frame.MaskKey = ''
3491                    G2frame.G2plotNB.status.SetFields(['','Frame closed'])
3492                else:
3493                    G2frame.G2plotNB.status.SetFields(['','New frame point: %.1f,%.1f'%(Xpos,Ypos)])
3494                    frame.append([Xpos,Ypos])
3495            G2imG.UpdateMasks(G2frame,Masks)
3496            PlotImage(G2frame,newImage=False)
3497        elif PickName == 'Stress/Strain' and G2frame.StrainKey:
3498            Xpos,Ypos = [event.xdata,event.ydata]
3499            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3500                return
3501            dsp = G2img.GetDsp(Xpos,Ypos,Data)
3502            StrSta['d-zero'].append({'Dset':dsp,'Dcalc':0.0,'pixLimit':10,'cutoff':1.0,
3503                'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]],'Emat':[1.0,1.0,1.0]})
3504            R,r = G2img.MakeStrStaRing(StrSta['d-zero'][-1],G2frame.ImageZ,Data)
3505            if not len(R):
3506                del StrSta['d-zero'][-1]
3507                G2frame.ErrorDialog('Strain peak selection','WARNING - No points found for this ring selection')
3508            StrSta['d-zero'] = G2mth.sortArray(StrSta['d-zero'],'Dset',reverse=True)
3509            G2frame.StrainKey = ''
3510            G2imG.UpdateStressStrain(G2frame,StrSta)
3511            PlotStrain(G2frame,StrSta)
3512            PlotImage(G2frame,newPlot=False)           
3513        else:
3514            Xpos,Ypos = [event.xdata,event.ydata]
3515            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
3516                return
3517            if G2frame.ifGetRing:                          #delete a calibration ring pick
3518                xypos = [Xpos,Ypos]
3519                rings = Data['ring']
3520                for ring in rings:
3521                    if np.allclose(ring,xypos,.01,0):
3522                        rings.remove(ring)
3523            else:
3524                tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data)[:3]
3525                itemPicked = str(G2frame.itemPicked)
3526                if 'Line2D' in itemPicked and PickName == 'Image Controls':
3527                    if 'line1' in itemPicked:
3528                        Data['IOtth'][0] = max(tth,0.001)
3529                    elif 'line2' in itemPicked:
3530                        Data['IOtth'][1] = tth
3531                    elif 'line3' in itemPicked:
3532                        Data['LRazimuth'][0] = int(azm)
3533                    elif 'line4' in itemPicked and not Data['fullIntegrate']:
3534                        Data['LRazimuth'][1] = int(azm)
3535                   
3536                    Data['LRazimuth'][0] %= 360
3537                    Data['LRazimuth'][1] %= 360
3538                    if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
3539                        Data['LRazimuth'][1] += 360                       
3540                    if Data['fullIntegrate']:
3541                        Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
3542                       
3543                    if  Data['IOtth'][0] > Data['IOtth'][1]:
3544                        Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
3545                       
3546                    G2frame.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
3547                    G2frame.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
3548                    G2frame.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
3549                    G2frame.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
3550                elif 'Circle' in itemPicked and PickName == 'Masks':
3551                    spots = Masks['Points']
3552                    newPos = itemPicked.split(')')[0].split('(')[2].split(',')
3553                    newPos = np.array([float(newPos[0]),float(newPos[1])])
3554                    for spot in spots:
3555                        if spot and np.allclose(np.array([spot[:2]]),newPos):
3556                            spot[:2] = Xpos,Ypos
3557                    G2imG.UpdateMasks(G2frame,Masks)
3558                elif 'Line2D' in itemPicked and PickName == 'Masks':
3559                    Obj = G2frame.itemPicked.findobj()
3560                    rings = Masks['Rings']
3561                    arcs = Masks['Arcs']
3562                    polygons = Masks['Polygons']
3563                    frame = Masks['Frames']
3564                    for ring in G2frame.ringList:
3565                        if Obj == ring[0]:
3566                            rN = ring[1]
3567                            if ring[2] == 'o':
3568                                rings[rN][0] = G2img.GetTth(Xpos,Ypos,Data)-rings[rN][1]/2.
3569                            else:
3570                                rings[rN][0] = G2img.GetTth(Xpos,Ypos,Data)+rings[rN][1]/2.
3571                    for arc in G2frame.arcList:
3572                        if Obj == arc[0]:
3573                            aN = arc[1]
3574                            if arc[2] == 'o':
3575                                arcs[aN][0] = G2img.GetTth(Xpos,Ypos,Data)-arcs[aN][2]/2
3576                            elif arc[2] == 'i':
3577                                arcs[aN][0] = G2img.GetTth(Xpos,Ypos,Data)+arcs[aN][2]/2
3578                            elif arc[2] == 'l':
3579                                arcs[aN][1][0] = int(G2img.GetAzm(Xpos,Ypos,Data))
3580                            else:
3581                                arcs[aN][1][1] = int(G2img.GetAzm(Xpos,Ypos,Data))
3582                    for poly in G2frame.polyList:   #merging points problem here?
3583                        if Obj == poly[0]:
3584                            ind = G2frame.itemPicked.contains(G2frame.mousePicked)[1]['ind'][0]
3585                            oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
3586                            pN = poly[1]
3587                            for i,xy in enumerate(polygons[pN]):
3588                                if np.allclose(np.array([xy]),oldPos,atol=1.0):
3589                                    if event.button == 1:
3590                                        polygons[pN][i] = Xpos,Ypos
3591                                    elif event.button == 3:
3592                                        polygons[pN].insert(i,[Xpos,Ypos])
3593                                        break
3594                    if frame:
3595                        oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
3596                        for i,xy in enumerate(frame):
3597                            if np.allclose(np.array([xy]),oldPos,atol=1.0):
3598                                if event.button == 1:
3599                                    frame[i] = Xpos,Ypos
3600                                elif event.button == 3:
3601                                    frame.insert(i,[Xpos,Ypos])
3602                                    break
3603                    G2imG.UpdateMasks(G2frame,Masks)
3604#                else:                  #keep for future debugging
3605#                    print str(G2frame.itemPicked),event.xdata,event.ydata,event.button
3606            PlotImage(G2frame,newImage=True)
3607            G2frame.itemPicked = None
3608           
3609    try:
3610        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
3611        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3612        if not newPlot:
3613            Plot = Page.figure.gca()          #get previous powder plot & get limits
3614            xylim = Plot.get_xlim(),Plot.get_ylim()
3615        if newImage:
3616            Page.figure.clf()
3617            Plot = Page.figure.gca()          #get a fresh plot after clf()
3618    except ValueError:
3619        Plot = G2frame.G2plotNB.addMpl('2D Powder Image').gca()
3620        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
3621        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3622        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
3623        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
3624        Page.canvas.mpl_connect('pick_event', OnImPick)
3625        Page.canvas.mpl_connect('button_release_event', OnImRelease)
3626        xylim = []
3627    Page.Choice = None
3628    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
3629        Page.SetFocus()
3630    Title = G2frame.PatternTree.GetItemText(G2frame.Image)[4:]
3631    G2frame.G2plotNB.status.DestroyChildren()
3632    if G2frame.logPlot:
3633        Title = 'log('+Title+')'
3634    Plot.set_title(Title)
3635    try:
3636        if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Image Controls',]:
3637            Page.Choice = (' key press','l: log(I) on','x: flip x','y: flip y',)
3638            if G2frame.logPlot:
3639                Page.Choice[1] = 'l: log(I) off'
3640            Page.keyPress = OnImPlotKeyPress
3641        elif G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Masks',]:
3642            Page.Choice = (' key press','l: log(I) on','s: spot mask','a: arc mask','r: ring mask',
3643                'p: polygon mask','f: frame mask',)
3644            if G2frame.logPlot:
3645                Page.Choice[1] = 'l: log(I) off'
3646            Page.keyPress = OnImPlotKeyPress
3647        elif G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Stress/Strain',]:
3648            Page.Choice = (' key press','a: add new ring',)
3649            Page.keyPress = OnImPlotKeyPress
3650    except TypeError:
3651        pass
3652    size,imagefile = G2frame.PatternTree.GetItemPyData(G2frame.Image)
3653    dark = Data.get('dark image',[0,''])
3654    if dark[0]:
3655        darkfile = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, 
3656            G2frame.root,dark[0]))[1]
3657    if imagefile != G2frame.oldImagefile:
3658        imagefile = G2IO.CheckImageFile(G2frame,imagefile)
3659        if not imagefile:
3660            G2frame.G2plotNB.Delete('2D Powder Image')
3661            return
3662        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
3663        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
3664        if dark[0]:
3665            darkImg = G2IO.GetImageData(G2frame,darkfile,imageOnly=True)
3666            G2frame.ImageZ += dark[1]*darkImg
3667        G2frame.oldImagefile = imagefile
3668
3669    imScale = 1
3670    if len(G2frame.ImageZ) > 1024:
3671        imScale = len(G2frame.ImageZ)/1024
3672    sizexy = Data['size']
3673    pixelSize = Data['pixelSize']
3674    scalex = 1000./pixelSize[0]
3675    scaley = 1000./pixelSize[1]
3676    Xmax = sizexy[0]*pixelSize[0]/1000.
3677    Ymax = sizexy[1]*pixelSize[1]/1000.
3678    xlim = (0,Xmax)
3679    ylim = (Ymax,0)
3680    Imin,Imax = Data['range'][1]
3681    acolor = mpl.cm.get_cmap(Data['color'])
3682    xcent,ycent = Data['center']
3683    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
3684    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
3685    #do threshold mask - "real" mask - others are just bondaries
3686    Zlim = Masks['Thresholds'][1]
3687    wx.BeginBusyCursor()
3688    try:
3689           
3690        if newImage:                   
3691            Imin,Imax = Data['range'][1]
3692            MA = ma.masked_greater(ma.masked_less(G2frame.ImageZ,Zlim[0]),Zlim[1])
3693            MaskA = ma.getmaskarray(MA)
3694            A = G2img.ImageCompress(MA,imScale)
3695            AM = G2img.ImageCompress(MaskA,imScale)
3696            if G2frame.logPlot:
3697                A = np.where(A>Imin,np.where(A<Imax,A,0),0)
3698                A = np.where(A>0,np.log(A),0)
3699                AM = np.where(AM>0,np.log(AM),0)
3700                Imin,Imax = [np.amin(A),np.amax(A)]
3701            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
3702                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Ymax,0])
3703            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
3704                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
3705   
3706        Plot.plot(xcent,ycent,'x')
3707        #G2frame.PatternTree.GetItemText(item)
3708        if Data['showLines']:
3709            LRAzim = Data['LRazimuth']                  #NB: integers
3710            Nazm = Data['outAzimuths']
3711            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
3712            AzmthOff = Data['azmthOff']
3713            IOtth = Data['IOtth']
3714            wave = Data['wavelength']
3715            dspI = wave/(2.0*sind(IOtth[0]/2.0))
3716            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
3717            dspO = wave/(2.0*sind(IOtth[1]/2.0))
3718            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
3719            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
3720            if ellI:
3721                xyI = []
3722                for azm in Azm:
3723                    xy = G2img.GetDetectorXY(dspI,azm,Data)
3724                    if np.any(xy):
3725                        xyI.append(xy)
3726                if len(xyI):
3727                    xyI = np.array(xyI)
3728                    arcxI,arcyI = xyI.T
3729                    Plot.plot(arcxI,arcyI,picker=3)
3730            if ellO:
3731                xyO = []
3732                for azm in Azm:
3733                    xy = G2img.GetDetectorXY(dspO,azm,Data)
3734                    if np.any(xy):
3735                        xyO.append(xy)
3736                if len(xyO):
3737                    xyO = np.array(xyO)
3738                    arcxO,arcyO = xyO.T               
3739                    Plot.plot(arcxO,arcyO,picker=3)
3740            if ellO and ellI:
3741                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
3742                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
3743            for i in range(Nazm):
3744                cake = LRAzim[0]+i*delAzm-AzmthOff
3745                if Data['centerAzm']:
3746                    cake += delAzm/2.
3747                ind = np.searchsorted(Azm,cake)
3748                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
3749                   
3750        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Image Controls':
3751            for xring,yring in Data['ring']:
3752                Plot.plot(xring,yring,'r+',picker=3)
3753            if Data['setRings']:
3754                N = 0
3755                for ring in Data['rings']:
3756                    xring,yring = np.array(ring).T[:2]
3757                    Plot.plot(xring,yring,'.',color=colors[N%6])
3758                    N += 1           
3759            for ellipse in Data['ellipses']:      #what about hyperbola?
3760                cent,phi,[width,height],col = ellipse
3761                if width > 0:       #ellipses
3762                    Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
3763                    Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
3764        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Stress/Strain':
3765            for N,ring in enumerate(StrSta['d-zero']):
3766                xring,yring = ring['ImxyObs']
3767                Plot.plot(xring,yring,colors[N%6]+'.')
3768        #masks - mask lines numbered after integration limit lines
3769        spots = Masks['Points']
3770        rings = Masks['Rings']
3771        arcs = Masks['Arcs']
3772        polygons = Masks['Polygons']
3773        if 'Frames' not in Masks:
3774            Masks['Frames'] = []
3775        frame = Masks['Frames']
3776        for spot in spots:
3777            if spot:
3778                x,y,d = spot
3779                Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
3780        G2frame.ringList = []
3781        for iring,ring in enumerate(rings):
3782            if ring:
3783                tth,thick = ring
3784                wave = Data['wavelength']
3785                xy1 = []
3786                xy2 = []
3787                Azm = np.linspace(0,362,181)
3788                for azm in Azm:
3789                    xy1.append(G2img.GetDetectorXY(Dsp(tth+thick/2.,wave),azm,Data))      #what about hyperbola
3790                    xy2.append(G2img.GetDetectorXY(Dsp(tth-thick/2.,wave),azm,Data))      #what about hyperbola
3791                x1,y1 = np.array(xy1).T
3792                x2,y2 = np.array(xy2).T
3793                G2frame.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
3794                G2frame.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
3795        G2frame.arcList = []
3796        for iarc,arc in enumerate(arcs):
3797            if arc:
3798                tth,azm,thick = arc           
3799                wave = Data['wavelength']
3800                xy1 = []
3801                xy2 = []
3802                aR = azm[0],azm[1],azm[1]-azm[0]
3803                if azm[1]-azm[0] > 180:
3804                    aR[2] /= 2
3805                Azm = np.linspace(aR[0],aR[1],aR[2])
3806                for azm in Azm:
3807                    xy1.append(G2img.GetDetectorXY(Dsp(tth+thick/2.,wave),azm,Data))      #what about hyperbola
3808                    xy2.append(G2img.GetDetectorXY(Dsp(tth-thick/2.,wave),azm,Data))      #what about hyperbola
3809                x1,y1 = np.array(xy1).T
3810                x2,y2 = np.array(xy2).T
3811                G2frame.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
3812                G2frame.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
3813                G2frame.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
3814                G2frame.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
3815        G2frame.polyList = []
3816        for ipoly,polygon in enumerate(polygons):
3817            if polygon:
3818                x,y = np.hsplit(np.array(polygon),2)
3819                G2frame.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
3820                Plot.plot(x,y,'r')           
3821        G2frame.frameList = []
3822        if frame:
3823            x,y = np.hsplit(np.array(frame),2)
3824            G2frame.frameList.append([Plot.plot(x,y,'g+',picker=10),0])
3825            Plot.plot(x,y,'g')           
3826        if newImage:
3827            colorBar = Page.figure.colorbar(Img)
3828        Plot.set_xlim(xlim)
3829        Plot.set_ylim(ylim)
3830        if Data['invert_x']:
3831            Plot.invert_xaxis()
3832        if Data['invert_y']:
3833            Plot.invert_yaxis()
3834        if not newPlot and xylim:
3835            Page.toolbar.push_current()
3836            Plot.set_xlim(xylim[0])
3837            Plot.set_ylim(xylim[1])
3838            xylim = []
3839            Page.toolbar.push_current()
3840            Page.toolbar.draw()
3841            # patch for wx 2.9 on Mac, to force a redraw
3842            i,j= wx.__version__.split('.')[0:2]
3843            if int(i)+int(j)/10. > 2.8 and 'wxOSX' in wx.PlatformInfo:
3844                Page.canvas.draw()
3845        else:
3846            Page.canvas.draw()
3847    finally:
3848        wx.EndBusyCursor()
3849       
3850################################################################################
3851##### PlotIntegration
3852################################################################################
3853           
3854def PlotIntegration(G2frame,newPlot=False,event=None):
3855    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
3856    '''
3857           
3858    def OnMotion(event):
3859        Page.canvas.SetToolTipString('')
3860        Page.canvas.SetCursor(wx.CROSS_CURSOR)
3861        azm = event.ydata
3862        tth = event.xdata
3863        if azm and tth:
3864            G2frame.G2plotNB.status.SetFields(\
3865                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
3866                               
3867    try:
3868        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
3869        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3870        if not newPlot:
3871            Plot = Page.figure.gca()          #get previous plot & get limits
3872            xylim = Plot.get_xlim(),Plot.get_ylim()
3873        Page.figure.clf()
3874        Plot = Page.figure.gca()          #get a fresh plot after clf()
3875       
3876    except ValueError:
3877        Plot = G2frame.G2plotNB.addMpl('2D Integration').gca()
3878        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
3879        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3880        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3881        Page.views = False
3882        view = False
3883    Page.Choice = None
3884    if not event:
3885        Page.SetFocus()
3886       
3887    Data = G2frame.PatternTree.GetItemPyData(
3888        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3889    image = G2frame.Integrate[0]
3890    xsc = G2frame.Integrate[1]
3891    ysc = G2frame.Integrate[2]
3892    Imin,Imax = Data['range'][1]
3893    acolor = mpl.cm.get_cmap(Data['color'])
3894    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
3895    Plot.set_ylabel('azimuth',fontsize=12)
3896    Plot.set_xlabel('2-theta',fontsize=12)
3897    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
3898        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
3899    colorBar = Page.figure.colorbar(Img)
3900#    if Data['ellipses']:           
3901#        for ellipse in Data['ellipses']:
3902#            x,y = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
3903#            tth,azm = G2img.GetTthAzm(x,y,Data)
3904##            azm = np.where(azm < 0.,azm+360,azm)
3905#            Plot.plot(tth,azm,'b,')
3906    if not newPlot:
3907        Page.toolbar.push_current()
3908        Plot.set_xlim(xylim[0])
3909        Plot.set_ylim(xylim[1])
3910        xylim = []
3911        Page.toolbar.push_current()
3912        Page.toolbar.draw()
3913    else:
3914        Page.canvas.draw()
3915               
3916################################################################################
3917##### PlotTRImage
3918################################################################################
3919           
3920def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):
3921    '''a test plot routine - not normally used
3922    ''' 
3923           
3924    def OnMotion(event):
3925        Page.canvas.SetToolTipString('')
3926        Page.canvas.SetCursor(wx.CROSS_CURSOR)
3927        azm = event.xdata
3928        tth = event.ydata
3929        if azm and tth:
3930            G2frame.G2plotNB.status.SetFields(\
3931                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
3932                               
3933    try:
3934        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
3935        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3936        if not newPlot:
3937            Plot = Page.figure.gca()          #get previous plot & get limits
3938            xylim = Plot.get_xlim(),Plot.get_ylim()
3939        Page.figure.clf()
3940        Plot = Page.figure.gca()          #get a fresh plot after clf()
3941       
3942    except ValueError:
3943        Plot = G2frame.G2plotNB.addMpl('2D Transformed Powder Image').gca()
3944        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
3945        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3946        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
3947        Page.views = False
3948        view = False
3949    Page.Choice = None
3950    Page.SetFocus()
3951       
3952    Data = G2frame.PatternTree.GetItemPyData(
3953        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
3954    Imin,Imax = Data['range'][1]
3955    step = (Imax-Imin)/5.
3956    V = np.arange(Imin,Imax,step)
3957    acolor = mpl.cm.get_cmap(Data['color'])
3958    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
3959    Plot.set_xlabel('azimuth',fontsize=12)
3960    Plot.set_ylabel('2-theta',fontsize=12)
3961    Plot.contour(tax,tay,taz,V,cmap=acolor)
3962    if Data['showLines']:
3963        IOtth = Data['IOtth']
3964        if Data['fullIntegrate']:
3965            LRAzim = [-180,180]
3966        else:
3967            LRAzim = Data['LRazimuth']                  #NB: integers
3968        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
3969        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
3970        if not Data['fullIntegrate']:
3971            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
3972            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
3973    if Data['setRings']:
3974        rings = np.concatenate((Data['rings']),axis=0)
3975        for xring,yring,dsp in rings:
3976            x,y = G2img.GetTthAzm(xring,yring,Data)
3977            Plot.plot(y,x,'r+')           
3978    if Data['ellipses']:           
3979        for ellipse in Data['ellipses']:
3980            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
3981            x,y = np.hsplit(ring,2)
3982            tth,azm = G2img.GetTthAzm(x,y,Data)
3983            Plot.plot(azm,tth,'b,')
3984    if not newPlot:
3985        Page.toolbar.push_current()
3986        Plot.set_xlim(xylim[0])
3987        Plot.set_ylim(xylim[1])
3988        xylim = []
3989        Page.toolbar.push_current()
3990        Page.toolbar.draw()
3991    else:
3992        Page.canvas.draw()
3993       
3994################################################################################
3995##### PlotStructure
3996################################################################################
3997           
3998def PlotStructure(G2frame,data,firstCall=False):
3999    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
4000    thermal motion ellipsoids and polyhedra
4001    '''
4002
4003    def FindPeaksBonds(XYZ):
4004        rFact = data['Drawing'].get('radiusFactor',0.85)    #data['Drawing'] could be empty!
4005        Bonds = [[] for x in XYZ]
4006        for i,xyz in enumerate(XYZ):
4007            Dx = XYZ-xyz
4008            dist = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=1))
4009            IndB = ma.nonzero(ma.masked_greater(dist,rFact*2.2))
4010            for j in IndB[0]:
4011                Bonds[i].append(Dx[j]/2.)
4012                Bonds[j].append(-Dx[j]/2.)
4013        return Bonds
4014
4015    # PlotStructure initialization here
4016    ForthirdPI = 4.0*math.pi/3.0
4017    generalData = data['General']
4018    cell = generalData['Cell'][1:7]
4019    Vol = generalData['Cell'][7:8][0]
4020    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
4021    Gmat,gmat = G2lat.cell2Gmat(cell)
4022    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
4023    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
4024    SGData = generalData['SGData']
4025    Mydir = generalData['Mydir']
4026    atomData = data['Atoms']
4027    mapPeaks = []
4028    drawingData = data['Drawing']
4029    if not drawingData:
4030        return          #nothing setup, nothing to draw   
4031    if 'Map Peaks' in data:
4032        mapPeaks = np.array(data['Map Peaks'])
4033        peakMax = 100.
4034        if len(mapPeaks):
4035            peakMax = np.max(mapPeaks.T[0])
4036    resRBData = data['RBModels'].get('Residue',[])
4037    vecRBData = data['RBModels'].get('Vector',[])
4038    rbAtmDict = {}
4039    for rbObj in resRBData+vecRBData:
4040        exclList = ['X' for i in range(len(rbObj['Ids']))]
4041        rbAtmDict.update(dict(zip(rbObj['Ids'],exclList)))
4042    testRBObj = data.get('testRBObj',{})
4043    rbObj = testRBObj.get('rbObj',{})
4044    MCSA = data.get('MCSA',{})
4045    mcsaModels = MCSA.get('Models',[])
4046    if len(mcsaModels) > 1:
4047        XYZs,Types = G2mth.UpdateMCSAxyz(Bmat,MCSA)
4048        mcsaXYZ = []
4049        mcsaTypes = []
4050        neqv = 0
4051        for xyz,atyp in zip(XYZs,Types):
4052            equiv = G2spc.GenAtom(xyz,SGData,All=True,Move=False)
4053            neqv = max(neqv,len(equiv))
4054            for item in equiv:
4055                mcsaXYZ.append(item[0]) 
4056                mcsaTypes.append(atyp)
4057        mcsaXYZ = np.array(mcsaXYZ)
4058        mcsaTypes = np.array(mcsaTypes)
4059        nuniq = mcsaXYZ.shape[0]/neqv
4060        mcsaXYZ = np.reshape(mcsaXYZ,(nuniq,neqv,3))
4061        mcsaTypes = np.reshape(mcsaTypes,(nuniq,neqv))
4062        cent = np.fix(np.sum(mcsaXYZ+2.,axis=0)/nuniq)-2
4063        cent[0] = [0,0,0]   #make sure 1st one isn't moved
4064        mcsaXYZ = np.swapaxes(mcsaXYZ,0,1)-cent[:,np.newaxis,:]
4065        mcsaTypes = np.swapaxes(mcsaTypes,0,1)
4066        mcsaXYZ = np.reshape(mcsaXYZ,(nuniq*neqv,3))
4067        mcsaTypes = np.reshape(mcsaTypes,(nuniq*neqv))                       
4068        mcsaBonds = FindPeaksBonds(mcsaXYZ)       
4069    drawAtoms = drawingData.get('Atoms',[])
4070    mapData = {}
4071    flipData = {}
4072    rhoXYZ = []
4073    showBonds = False
4074    if 'Map' in generalData:
4075        mapData = generalData['Map']
4076        showBonds = mapData.get('Show bonds',False)
4077    if 'Flip' in generalData:
4078        flipData = generalData['Flip']                       
4079    Wt = np.array([255,255,255])
4080    Rd = np.array([255,0,0])
4081    Gr = np.array([0,255,0])
4082    wxGreen = wx.Colour(0,255,0)
4083    Bl = np.array([0,0,255])
4084    Or = np.array([255,128,0])
4085    wxOrange = wx.Colour(255,128,0)
4086    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]])
4087    uEdges = np.array([
4088        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
4089        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
4090        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
4091    mD = 0.1
4092    mV = np.array([[[-mD,0,0],[mD,0,0]],[[0,-mD,0],[0,mD,0]],[[0,0,-mD],[0,0,mD]]])
4093    mapPeakVecs = np.inner(mV,Bmat)
4094
4095    backColor = np.array(list(drawingData['backColor'])+[0,])
4096    Bc = np.array(list(drawingData['backColor']))
4097    uColors = [Rd,Gr,Bl,Wt-Bc, Wt-Bc,Wt-Bc,Wt-Bc,Wt-Bc, Wt-Bc,Wt-Bc,Wt-Bc,Wt-Bc]
4098    altDown = False
4099    shiftDown = False
4100    ctrlDown = False
4101   
4102    def OnKeyBox(event):
4103#        Draw()                          #make sure plot is fresh!!
4104        mode = cb.GetValue()
4105        if mode in ['jpeg','bmp','tiff',]:
4106            try:
4107                import Image as Im
4108            except ImportError:
4109                try:
4110                    from PIL import Image as Im
4111                except ImportError:
4112                    print "PIL/pillow Image module not present. Cannot save images without this"
4113                    raise Exception("PIL/pillow Image module not found")
4114            Fname = os.path.join(Mydir,generalData['Name']+'.'+mode)
4115            size = Page.canvas.GetSize()
4116            glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
4117            if mode in ['jpeg',]:
4118                Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
4119                im = Im.new("RGBA", (size[0],size[1]))
4120            else:
4121                Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
4122                im = Im.new("RGB", (size[0],size[1]))
4123            im.fromstring(Pix)
4124            im.save(Fname,mode)
4125            cb.SetValue(' save as/key:')
4126            G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
4127        else:
4128            event.key = cb.GetValue()[0]
4129            cb.SetValue(' save as/key:')
4130            wx.CallAfter(OnKey,event)
4131        Page.canvas.SetFocus() # redirect the Focus from the button back to the plot
4132
4133    def OnKey(event):           #on key UP!!
4134#        Draw()                          #make sure plot is fresh!!
4135        try:
4136            keyCode = event.GetKeyCode()
4137            if keyCode > 255:
4138                keyCode = 0
4139            key = chr(keyCode)
4140        except AttributeError:       #if from OnKeyBox above
4141            key = str(event.key).upper()
4142        indx = drawingData['selectedAtoms']
4143        cx,ct = drawingData['atomPtrs'][:2]
4144        if key in ['C']:
4145            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
4146            drawingData['viewDir'] = [0,0,1]
4147            drawingData['oldxy'] = []
4148            V0 = np.array([0,0,1])
4149            V = np.inner(Amat,V0)
4150            V /= np.sqrt(np.sum(V**2))
4151            A = np.arccos(np.sum(V*V0))
4152            Q = G2mth.AV2Q(A,[0,1,0])
4153            drawingData['Quaternion'] = Q
4154            SetViewPointText(drawingData['viewPoint'][0])
4155            SetViewDirText(drawingData['viewDir'])
4156            Q = drawingData['Quaternion']
4157            G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4158        elif key in ['N']:
4159            drawAtoms = drawingData['Atoms']
4160            if not len(drawAtoms):      #no atoms
4161                return
4162            pI = drawingData['viewPoint'][1]
4163            if not len(pI):
4164                pI = [0,0]
4165            if indx:
4166                pI[0] = indx[pI[1]]
4167                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
4168                pI[1] += 1
4169                if pI[1] >= len(indx):
4170                    pI[1] = 0
4171            else:
4172                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
4173                pI[0] += 1
4174                if pI[0] >= len(drawAtoms):
4175                    pI[0] = 0
4176            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
4177            SetViewPointText(drawingData['viewPoint'][0])
4178            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
4179               
4180        elif key in ['P']:
4181            drawAtoms = drawingData['Atoms']
4182            if not len(drawAtoms):      #no atoms
4183                return
4184            pI = drawingData['viewPoint'][1]
4185            if not len(pI):
4186                pI = [0,0]
4187            if indx:
4188                pI[0] = indx[pI[1]]
4189                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
4190                pI[1] -= 1
4191                if pI[1] < 0:
4192                    pI[1] = len(indx)-1
4193            else:
4194                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
4195                pI[0] -= 1
4196                if pI[0] < 0:
4197                    pI[0] = len(drawAtoms)-1
4198            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
4199            SetViewPointText(drawingData['viewPoint'][0])           
4200            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
4201        elif key in ['U','D','L','R'] and mapData['Flip'] == True:
4202            dirDict = {'U':[0,1],'D':[0,-1],'L':[-1,0],'R':[1,0]}
4203            SetMapRoll(dirDict[key])
4204            if 'rho' in generalData.get('4DmapData',{}):
4205                Set4DMapRoll(dirDict[key])
4206            SetPeakRoll(dirDict[key])
4207            SetMapPeaksText(mapPeaks)
4208        Draw('key')
4209           
4210    def GetTruePosition(xy,Add=False):
4211        View = glGetIntegerv(GL_VIEWPORT)
4212        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
4213        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
4214        Zmax = 1.
4215        if Add:
4216            Indx = GetSelectedAtoms()
4217        if G2frame.dataDisplay.GetPageText(getSelection()) == 'Map peaks':
4218            for i,peak in enumerate(mapPeaks):
4219                x,y,z = peak[1:4]
4220                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
4221                XY = [int(X),int(View[3]-Y)]
4222                if np.allclose(xy,XY,atol=10) and Z < Zmax:
4223                    Zmax = Z
4224                    try:
4225                        Indx.remove(i)
4226                        ClearSelectedAtoms()
4227                        for id in Indx:
4228                            SetSelectedAtoms(id,Add)
4229                    except:
4230                        SetSelectedAtoms(i,Add)
4231        else:
4232            cx = drawingData['atomPtrs'][0]
4233            for i,atom in enumerate(drawAtoms):
4234                x,y,z = atom[cx:cx+3]
4235                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
4236                XY = [int(X),int(View[3]-Y)]
4237                if np.allclose(xy,XY,atol=10) and Z < Zmax:
4238                    Zmax = Z
4239                    try:
4240                        Indx.remove(i)
4241                        ClearSelectedAtoms()
4242                        for id in Indx:
4243                            SetSelectedAtoms(id,Add)
4244                    except:
4245                        SetSelectedAtoms(i,Add)
4246                                       
4247    def OnMouseDown(event):
4248        xy = event.GetPosition()
4249        if event.ShiftDown():
4250            if event.LeftIsDown():
4251                GetTruePosition(xy)
4252            elif event.RightIsDown():
4253                GetTruePosition(xy,True)
4254        else:
4255            drawingData['oldxy'] = list(xy)
4256       
4257    def OnMouseMove(event):
4258        if event.ShiftDown():           #don't want any inadvertant moves when picking
4259            return
4260        newxy = event.GetPosition()
4261                               
4262        if event.Dragging():
4263            if event.AltDown() and rbObj:
4264                if event.LeftIsDown():
4265                    SetRBRotation(newxy)
4266                    Q = rbObj['Orient'][0]
4267                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4268                elif event.RightIsDown():
4269                    SetRBTranslation(newxy)
4270                    Tx,Ty,Tz = rbObj['Orig'][0]
4271                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
4272                elif event.MiddleIsDown():
4273                    SetRBRotationZ(newxy)
4274                    Q = rbObj['Orient'][0]
4275                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4276                Draw('move')
4277            elif not event.ControlDown():
4278                if event.LeftIsDown():
4279                    SetRotation(newxy)
4280                    Q = drawingData['Quaternion']
4281                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4282                elif event.RightIsDown():
4283                    SetTranslation(newxy)
4284                    Tx,Ty,Tz = drawingData['viewPoint'][0]
4285                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
4286                elif event.MiddleIsDown():
4287                    SetRotationZ(newxy)
4288                    Q = drawingData['Quaternion']
4289                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
4290                Draw('move')
4291           
4292       
4293    def OnMouseWheel(event):
4294        if event.ShiftDown():
4295            return
4296        drawingData['cameraPos'] += event.GetWheelRotation()/24.
4297        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
4298        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(drawingData['cameraPos']),1)
4299        page = getSelection()
4300        if page:
4301            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
4302                G2frame.dataDisplay.cameraPosTxt.SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
4303                G2frame.dataDisplay.cameraSlider.SetValue(drawingData['cameraPos'])
4304            Draw('wheel')
4305       
4306    def getSelection():
4307        try:
4308            return G2frame.dataDisplay.GetSelection()
4309        except AttributeError:
4310            G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!',1)
4311            return 0
4312           
4313    def SetViewPointText(VP):
4314        page = getSelection()
4315        if page:
4316            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
4317                G2frame.dataDisplay.viewPoint.SetValue('%.3f %.3f %.3f'%(VP[0],VP[1],VP[2]))
4318               
4319    def SetRBOrigText():
4320        page = getSelection()
4321        if page:
4322            if G2frame.dataDisplay.GetPageText(page) == 'RB Models':
4323                for i,sizer in enumerate(testRBObj['Sizers']['Xsizers']):
4324                    sizer.SetValue('%8.5f'%(testRBObj['rbObj']['Orig'][0][i]))
4325                   
4326    def SetRBOrienText():
4327        page = getSelection