source: trunk/GSASIIplot.py @ 1368

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

implement sqrt(intensity) powder plots

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