source: trunk/GSASIIplot.py @ 1403

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

allow renaming of multiple columns in seq refinement results
allow changing of plot titles for seq plots.
new dialog - MultiStringDialog?

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