source: trunk/GSASIIplot.py @ 1087

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

simulated CW powder data now complete with Poisson errors

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