source: trunk/GSASIIplot.py @ 1017

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

implement excluded regions - 1st attempt

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