source: trunk/GSASIIplot.py @ 1366

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

change 's' to 'm' - makes more sense

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