source: trunk/GSASIIplot.py @ 1585

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

remove Pwr from FitPeaksZ, SSZ & T - wasn't used
put in derivatives for FitPeaksZ, etc. - works better
reorder arrows & make the expand/shrink symmetric about view point
remove zero & refine? from SS Unit Cells display - not refinable with modulation coeff.

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