source: trunk/GSASIIplot.py @ 1077

Last change on this file since 1077 was 1077, checked in by toby, 8 years ago

cleanup plot & svn bugs; set missing keywords; CIF export done; update docs

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