source: trunk/GSASIIplot.py @ 1273

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

clean up SASD Model plotting
clean up form factor functions

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