source: trunk/GSASIIplot.py @ 1370

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

plot labels cleanup

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