source: trunk/GSASIIplot.py @ 1252

Last change on this file since 1252 was 1252, checked in by vondreele, 9 years ago

unsuccessful implementation of IPG method (commented out) for size distribution
revisit later
plot SASD background on changes
document SASD shape form factors & volumes

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