source: trunk/GSASIIplot.py @ 1107

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

make reflist np.array
fix bug - hfx+'Type' needs to be in parmDict

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