source: trunk/GSASIIplot.py @ 1148

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

major changes to 2-D detector calibration; now works for strongly tilted detectors & short sample to detector distances.
Distance is now defined as sample to detector plane. Previously it was sample to intercept of detector plane with incident beam (Bragg cone axis).
The "penetration" parameter is still suspect.

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