source: trunk/GSASIIplot.py @ 1060

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

mods to MC/SA stuff - some speedup done

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