source: trunk/GSASIIplot.py @ 1396

Last change on this file since 1396 was 1396, checked in by toby, 8 years ago

fix bugs in multieq seq fits; print results better

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