source: trunk/GSASIIplot.py @ 860

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

more fixes to rigid body stuff

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