source: trunk/GSASIIplot.py @ 1268

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

WACV = wx.ALIGN_CENTER_VERTICAL in G2imgGUI
add more sasd modeling code

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