source: trunk/GSASIIplot.py @ 505

Last change on this file since 505 was 505, checked in by vondreele, 11 years ago

start on plotting "error analysis" for powder patterns

  • Property svn:keywords set to Date Author Revision URL Id
File size: 118.7 KB
Line 
1#GSASII plotting routines
2########### SVN repository information ###################
3# $Date: 2012-03-01 16:29:55 +0000 (Thu, 01 Mar 2012) $
4# $Author: vondreele $
5# $Revision: 505 $
6# $URL: trunk/GSASIIplot.py $
7# $Id: GSASIIplot.py 505 2012-03-01 16:29:55Z vondreele $
8########### SVN repository information ###################
9import math
10import time
11import copy
12import os.path
13import numpy as np
14import numpy.linalg as nl
15import wx
16import wx.aui
17import wx.glcanvas
18import matplotlib as mpl
19import mpl_toolkits.mplot3d.axes3d as mp3d
20import GSASIIpath
21import GSASIIgrid as G2gd
22import GSASIIimage as G2img
23import GSASIIpwd as G2pwd
24import GSASIIIO as G2IO
25import GSASIIpwdGUI as G2pdG
26import GSASIIimgGUI as G2imG
27import GSASIIphsGUI as G2phG
28import GSASIIlattice as G2lat
29import GSASIIspc as G2spc
30import pytexture as ptx
31from  OpenGL.GL import *
32from OpenGL.GLU import *
33from OpenGL.GLUT import *
34from OpenGL.GLE import *
35from matplotlib.backends.backend_wx import _load_bitmap
36from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as Canvas
37from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
38
39# useful degree trig functions
40sind = lambda x: math.sin(x*math.pi/180.)
41cosd = lambda x: math.cos(x*math.pi/180.)
42tand = lambda x: math.tan(x*math.pi/180.)
43asind = lambda x: 180.*math.asin(x)/math.pi
44acosd = lambda x: 180.*math.acos(x)/math.pi
45atan2d = lambda x,y: 180.*math.atan2(y,x)/math.pi
46atand = lambda x: 180.*math.atan(x)/math.pi
47# numpy versions
48npsind = lambda x: np.sin(x*np.pi/180.)
49npcosd = lambda x: np.cos(x*np.pi/180.)
50npacosd = lambda x: 180.*np.arccos(x)/np.pi
51npasind = lambda x: 180.*np.arcsin(x)/np.pi
52npatand = lambda x: 180.*np.arctan(x)/np.pi
53npatan2d = lambda x,y: 180.*np.arctan2(x,y)/np.pi
54   
55class G2PlotMpl(wx.Panel):   
56    def __init__(self,parent,id=-1,dpi=None,**kwargs):
57        wx.Panel.__init__(self,parent,id=id,**kwargs)
58        mpl.rcParams['legend.fontsize'] = 10
59        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(5,7))
60        self.canvas = Canvas(self,-1,self.figure)
61        self.toolbar = GSASIItoolbar(self.canvas)
62
63        self.toolbar.Realize()
64       
65        sizer=wx.BoxSizer(wx.VERTICAL)
66        sizer.Add(self.canvas,1,wx.EXPAND)
67        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
68        self.SetSizer(sizer)
69       
70class G2PlotOgl(wx.Panel):
71    def __init__(self,parent,id=-1,dpi=None,**kwargs):
72        self.figure = wx.Panel.__init__(self,parent,id=id,**kwargs)
73        self.canvas = wx.glcanvas.GLCanvas(self,-1,**kwargs)
74        self.camera = {}
75        sizer=wx.BoxSizer(wx.VERTICAL)
76        sizer.Add(self.canvas,1,wx.EXPAND)
77        self.SetSizer(sizer)
78       
79class G2Plot3D(wx.Panel):
80    def __init__(self,parent,id=-1,dpi=None,**kwargs):
81        wx.Panel.__init__(self,parent,id=id,**kwargs)
82        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(6,6))
83        self.canvas = Canvas(self,-1,self.figure)
84        self.toolbar = GSASIItoolbar(self.canvas)
85
86        self.toolbar.Realize()
87       
88        sizer=wx.BoxSizer(wx.VERTICAL)
89        sizer.Add(self.canvas,1,wx.EXPAND)
90        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
91        self.SetSizer(sizer)
92                             
93class G2PlotNoteBook(wx.Panel):
94    def __init__(self,parent,id=-1):
95        wx.Panel.__init__(self,parent,id=id)
96        #so one can't delete a plot page!!
97        self.nb = wx.aui.AuiNotebook(self, \
98            style=wx.aui.AUI_NB_DEFAULT_STYLE ^ wx.aui.AUI_NB_CLOSE_ON_ACTIVE_TAB)
99        sizer = wx.BoxSizer()
100        sizer.Add(self.nb,1,wx.EXPAND)
101        self.SetSizer(sizer)
102        self.status = parent.CreateStatusBar()
103        self.status.SetFieldsCount(2)
104        self.status.SetStatusWidths([150,-1])
105        self.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, self.OnPageChanged)
106       
107        self.plotList = []
108           
109    def addMpl(self,name=""):
110        page = G2PlotMpl(self.nb)
111        self.nb.AddPage(page,name)
112       
113        self.plotList.append(name)
114       
115        return page.figure
116       
117    def add3D(self,name=""):
118        page = G2Plot3D(self.nb)
119        self.nb.AddPage(page,name)
120       
121        self.plotList.append(name)
122       
123        return page.figure
124       
125    def addOgl(self,name=""):
126        page = G2PlotOgl(self.nb)
127        self.nb.AddPage(page,name)
128       
129        self.plotList.append(name)
130       
131        return page.figure
132       
133    def Delete(self,name):
134        try:
135            item = self.plotList.index(name)
136            del self.plotList[item]
137            self.nb.DeletePage(item)
138        except ValueError:          #no plot of this name - do nothing
139            return     
140               
141    def clear(self):
142        while self.nb.GetPageCount():
143            self.nb.DeletePage(0)
144        self.plotList = []
145        self.status.DestroyChildren()
146       
147    def Rename(self,oldName,newName):
148        try:
149            item = self.plotList.index(oldName)
150            self.plotList[item] = newName
151            self.nb.SetPageText(item,newName)
152        except ValueError:          #no plot of this name - do nothing
153            return     
154       
155    def OnPageChanged(self,event):       
156        if self.plotList:
157            self.status.SetStatusText('Better to select this from GSAS-II data tree',1)
158        self.status.DestroyChildren()                           #get rid of special stuff on status bar
159       
160class GSASIItoolbar(Toolbar):
161    ON_MPL_HELP = wx.NewId()
162    def __init__(self,plotCanvas):
163        Toolbar.__init__(self,plotCanvas)
164        POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6
165        self.DeleteToolByPos(POSITION_OF_CONFIGURE_SUBPLOTS_BTN)
166        help = os.path.join(os.path.split(__file__)[0],'help.ico')
167        self.AddSimpleTool(self.ON_MPL_HELP,_load_bitmap(help),'Help on','Show help on')
168        wx.EVT_TOOL(self,self.ON_MPL_HELP,self.OnHelp)
169    def OnHelp(self,event):
170        Page = self.GetParent().GetParent()
171        pageNo = Page.GetSelection()
172        bookmark = Page.GetPageText(pageNo)
173        bookmark = bookmark.strip(')').replace('(','_')
174        G2gd.ShowHelp(bookmark,self.TopLevelParent)
175           
176def PlotSngl(self,newPlot=False):
177    '''Single crystal structure factor plotting package - displays zone of reflections as rings proportional
178        to F, F**2, etc. as requested
179    '''
180    from matplotlib.patches import Circle
181    global HKL,HKLF
182
183    def OnSCMotion(event):
184        xpos = event.xdata
185        if xpos:
186            xpos = round(xpos)                                        #avoid out of frame mouse position
187            ypos = round(event.ydata)
188            zpos = Data['Layer']
189            if '100' in Data['Zone']:
190                HKLtxt = '(%3d,%3d,%3d)'%(zpos,xpos,ypos)
191            elif '010' in Data['Zone']:
192                HKLtxt = '(%3d,%3d,%3d)'%(xpos,zpos,ypos)
193            elif '001' in Data['Zone']:
194                HKLtxt = '(%3d,%3d,%3d)'%(xpos,ypos,zpos)
195            Page.canvas.SetToolTipString(HKLtxt)
196            self.G2plotNB.status.SetFields(['HKL = '+HKLtxt,''])
197               
198    def OnSCPick(event):
199        zpos = Data['Layer']
200        pos = event.artist.center
201        if '100' in Data['Zone']:
202            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(zpos,pos[0],pos[1]))
203            hkl = np.array([zpos,pos[0],pos[1]])
204        elif '010' in Data['Zone']:
205            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],zpos,pos[1]))
206            hkl = np.array([pos[0],zpos,pos[1]])
207        elif '001' in Data['Zone']:
208            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],pos[1],zpos))
209            hkl = np.array([pos[0],pos[1],zpos])
210        h,k,l = hkl
211        hklf = HKLF[np.where(np.all(HKL-hkl == [0,0,0],axis=1))]
212        if len(hklf):
213            Fosq,sig,Fcsq = hklf[0]
214            HKLtxt = '(%3d,%3d,%3d %.2f %.3f %.2f %.2f)'%(h,k,l,Fosq,sig,Fcsq,(Fosq-Fcsq)/(scale*sig))
215            self.G2plotNB.status.SetFields(['','HKL, Fosq, sig, Fcsq, delFsq/sig = '+HKLtxt])
216                                 
217    def OnSCKeyPress(event):
218        print event.key
219
220    try:
221        plotNum = self.G2plotNB.plotList.index('Structure Factors')
222        Page = self.G2plotNB.nb.GetPage(plotNum)
223        if not newPlot:
224            Plot = Page.figure.gca()          #get previous powder plot & get limits
225            xylim = Plot.get_xlim(),Plot.get_ylim()
226        Page.figure.clf()
227        Plot = Page.figure.gca()          #get a fresh plot after clf()
228    except ValueError:
229        Plot = self.G2plotNB.addMpl('Structure Factors').gca()
230        plotNum = self.G2plotNB.plotList.index('Structure Factors')
231        Page = self.G2plotNB.nb.GetPage(plotNum)
232#        Page.canvas.mpl_connect('key_press_event', OnSCKeyPress)
233        Page.canvas.mpl_connect('pick_event', OnSCPick)
234        Page.canvas.mpl_connect('motion_notify_event', OnSCMotion)
235    Page.SetFocus()
236   
237    Plot.set_aspect(aspect='equal')
238    HKLref = self.PatternTree.GetItemPyData(self.Sngl)
239    Data = self.PatternTree.GetItemPyData( \
240        G2gd.GetPatternTreeItemId(self,self.Sngl, 'HKL Plot Controls'))
241    Type = Data['Type']           
242    scale = Data['Scale']
243    HKLmax = Data['HKLmax']
244    HKLmin = Data['HKLmin']
245    FosqMax = Data['FoMax']
246    FoMax = math.sqrt(FosqMax)
247    ifFc = Data['ifFc']
248    xlabel = ['k, h=','h, k=','h, l=']
249    ylabel = ['l','l','k']
250    zones = ['100','010','001']
251    pzone = [[1,2],[0,2],[0,1]]
252    izone = zones.index(Data['Zone'])
253    Plot.set_title(self.PatternTree.GetItemText(self.Sngl)[5:])
254    HKL = []
255    HKLF = []
256    for H,Fosq,sig,Fcsq,x,x,x in HKLref:
257        HKL.append(H)
258        HKLF.append([Fosq,sig,Fcsq])
259        if H[izone] == Data['Layer']:
260            B = 0
261            if Type == 'Fosq':
262                A = scale*Fosq/FosqMax
263                B = scale*Fcsq/FosqMax
264                C = abs(A-B)
265            elif Type == 'Fo':
266                A = scale*math.sqrt(max(0,Fosq))/FoMax
267                B = scale*math.sqrt(max(0,Fcsq))/FoMax
268                C = abs(A-B)
269            elif Type == '|DFsq|/sig':
270                A = abs(Fosq-Fcsq)/(scale*sig)
271            elif Type == '|DFsq|>sig':
272                A = abs(Fosq-Fcsq)/(scale*sig)
273                if A < 1.0: A = 0                   
274            elif Type == '|DFsq|>3sig':
275                A = abs(Fosq-Fcsq)/(scale*sig)
276                if A < 3.0: A = 0                   
277            xy = (H[pzone[izone][0]],H[pzone[izone][1]])
278            if A > 0.0:
279                Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w',picker=3))
280            if B:
281                Plot.add_artist(Circle(xy,radius=B,ec='b',fc='w'))
282                radius = C
283                if radius > 0:
284                    if A > B:
285                        Plot.add_artist(Circle(xy,radius=radius,ec='g',fc='g'))
286                    else:                   
287                        Plot.add_artist(Circle(xy,radius=radius,ec='r',fc='r'))
288    HKL = np.array(HKL,dtype=np.int)
289    HKLF = np.array(HKLF)
290    Plot.set_xlabel(xlabel[izone]+str(Data['Layer']),fontsize=12)
291    Plot.set_ylabel(ylabel[izone],fontsize=12)
292    Plot.set_xlim((HKLmin[pzone[izone][0]],HKLmax[pzone[izone][0]]))
293    Plot.set_ylim((HKLmin[pzone[izone][1]],HKLmax[pzone[izone][1]]))
294    if not newPlot:
295        Page.toolbar.push_current()
296        Plot.set_xlim(xylim[0])
297        Plot.set_ylim(xylim[1])
298        xylim = []
299        Page.toolbar.push_current()
300        Page.toolbar.draw()
301    else:
302        Page.canvas.draw()
303       
304def PlotPatterns(G2frame,newPlot=False):
305    '''Powder pattern plotting package - displays single or multiple powder patterns as intensity vs
306    2-theta or q (future TOF). Can display multiple patterns as "waterfall plots" or contour plots. Log I
307    plotting available.
308    '''
309    global HKL
310   
311    def OnKeyBox(event):
312        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('Powder Patterns'):
313            event.key = cb.GetValue()[0]
314            cb.SetValue(' key press')
315            OnPlotKeyPress(event)
316                       
317    def OnPlotKeyPress(event):
318        newPlot = False
319        if event.key == 'w':
320            if G2frame.Weight:
321                G2frame.Weight = False
322            else:
323                G2frame.Weight = True
324            print 'error analysis:',G2frame.Weight
325        elif event.key == 'n':
326            if G2frame.Contour:
327                pass
328            else:
329                if G2frame.logPlot:
330                    G2frame.logPlot = False
331                else:
332                    G2frame.Offset[0] = 0
333                    G2frame.logPlot = True
334        elif event.key == 'u':
335            if G2frame.Contour:
336                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
337            elif G2frame.logPlot:
338                pass
339            elif G2frame.Offset[0] < 100.:
340                G2frame.Offset[0] += 1.
341        elif event.key == 'd':
342            if G2frame.Contour:
343                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
344            elif G2frame.logPlot:
345                pass
346            elif G2frame.Offset[0] > 0.:
347                G2frame.Offset[0] -= 1.
348        elif event.key == 'l':
349            G2frame.Offset[1] -= 1.
350        elif event.key == 'r':
351            G2frame.Offset[1] += 1.
352        elif event.key == 'o':
353            G2frame.Offset = [0,0]
354        elif event.key == 'c':
355            newPlot = True
356            if G2frame.Contour:
357                G2frame.Contour = False
358            else:
359                G2frame.Contour = True
360                G2frame.SinglePlot = False
361                G2frame.Offset = [0.,0.]
362        elif event.key == 'q':
363            newPlot = True
364            if G2frame.qPlot:
365                G2frame.qPlot = False
366            else:
367                G2frame.qPlot = True
368        elif event.key == 's':
369            if G2frame.Contour:
370                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
371                choice.sort()
372                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
373                if dlg.ShowModal() == wx.ID_OK:
374                    sel = dlg.GetSelection()
375                    G2frame.ContourColor = choice[sel]
376                else:
377                    G2frame.ContourColor = 'Paired'
378                dlg.Destroy()
379            else:               
380                if G2frame.SinglePlot:
381                    G2frame.SinglePlot = False
382                else:
383                    G2frame.SinglePlot = True
384        elif event.key == '+':
385            if G2frame.PickId:
386                G2frame.PickId = False
387        elif event.key == 'i':                  #for smoothing contour plot
388            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
389               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
390               'mitchell','sinc','lanczos']
391            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
392            if dlg.ShowModal() == wx.ID_OK:
393                sel = dlg.GetSelection()
394                G2frame.Interpolate = choice[sel]
395            else:
396                G2frame.Interpolate = 'nearest'
397            dlg.Destroy()
398           
399        PlotPatterns(G2frame,newPlot=newPlot)
400       
401    def OnMotion(event):
402        xpos = event.xdata
403        if xpos:                                        #avoid out of frame mouse position
404            ypos = event.ydata
405            Page.canvas.SetCursor(wx.CROSS_CURSOR)
406            try:
407                Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[1::2]
408                Parms = dict(zip(Names,Values))
409                try:
410                    wave = Parms['Lam']
411                except KeyError:
412                    wave = Parms['Lam1']
413                if G2frame.qPlot:
414                    xpos = 2.0*asind(xpos*wave/(4*math.pi))
415                dsp = 0.0
416                if abs(xpos) > 0.:                  #avoid possible singularity at beam center
417                    dsp = wave/(2.*sind(abs(xpos)/2.0))
418                if G2frame.Contour:
419                    G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f pattern ID =%5d'%(xpos,dsp,int(ypos)),1)
420                else:
421                    G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f Intensity =%9.1f'%(xpos,dsp,ypos),1)
422                if G2frame.itemPicked:
423                    Page.canvas.SetToolTipString('%9.3f'%(xpos))
424                if G2frame.PickId:
425                    found = []
426                    if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists']:
427                        if len(HKL):
428                            view = Page.toolbar._views.forward()[0][:2]
429                            wid = view[1]-view[0]
430                            found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
431                        if len(found):
432                            h,k,l = found[0][:3] 
433                            Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
434                        else:
435                            Page.canvas.SetToolTipString('')
436
437            except TypeError:
438                G2frame.G2plotNB.status.SetStatusText('Select PWDR powder pattern first',1)
439                                                   
440    def OnPick(event):
441        if G2frame.itemPicked is not None: return
442        PatternId = G2frame.PatternId
443        try:
444            Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[1::2]
445        except TypeError:
446            return
447        Parms = dict(zip(Names,Values))
448        try:
449            wave = Parms['Lam']
450        except KeyError:
451            wave = Parms['Lam1']
452        PickId = G2frame.PickId
453        pick = event.artist
454        mouse = event.mouseevent       
455        xpos = pick.get_xdata()
456        ypos = pick.get_ydata()
457        ind = event.ind
458        xy = list(zip(np.take(xpos,ind),np.take(ypos,ind))[0])
459        if G2frame.PatternTree.GetItemText(PickId) == 'Peak List':
460            if ind.all() != [0]:                                    #picked a data point
461                if 'C' in Parms['Type']:                            #CW data - TOF later in an elif
462                    ins = [Parms[x] for x in ['U','V','W','X','Y']]
463                    if G2frame.qPlot:                              #qplot - convert back to 2-theta
464                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
465                    sig = ins[0]*tand(xy[0]/2.0)**2+ins[1]*tand(xy[0]/2.0)+ins[2]
466                    gam = ins[3]/cosd(xy[0]/2.0)+ins[4]*tand(xy[0]/2.0)           
467                    data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
468                    XY = [xy[0],0, xy[1],1, sig,0, gam,0]       #default refine intensity 1st
469                data.append(XY)
470                G2pdG.UpdatePeakGrid(G2frame,data)
471                PlotPatterns(G2frame)
472            else:                                                   #picked a peak list line
473                G2frame.itemPicked = pick
474        elif G2frame.PatternTree.GetItemText(PickId) == 'Limits':
475            if ind.all() != [0]:                                    #picked a data point
476                LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
477                data = G2frame.PatternTree.GetItemPyData(LimitId)
478                if 'C' in Parms['Type']:                            #CW data - TOF later in an elif
479                    if G2frame.qPlot:                              #qplot - convert back to 2-theta
480                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
481                if mouse.button==1:
482                    data[1][0] = min(xy[0],data[1][1])
483                if mouse.button==3:
484                    data[1][1] = max(xy[0],data[1][0])
485                G2frame.PatternTree.SetItemPyData(LimitId,data)
486                G2pdG.UpdateLimitsGrid(G2frame,data)
487                PlotPatterns(G2frame)
488            else:                                                   #picked a limit line
489                G2frame.itemPicked = pick
490        elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists':
491            G2frame.itemPicked = pick
492            pick = str(pick)
493       
494    def OnRelease(event):
495        if G2frame.itemPicked is None: return
496        Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))[1::2]
497        Parms = dict(zip(Names,Values))
498        try:
499            wave = Parms['Lam']
500        except KeyError:
501            wave = Parms['Lam1']
502        xpos = event.xdata
503        PickId = G2frame.PickId
504        if G2frame.PatternTree.GetItemText(PickId) in ['Peak List','Limits'] and xpos:
505            lines = []
506            for line in G2frame.Lines: lines.append(line.get_xdata()[0])
507            lineNo = lines.index(G2frame.itemPicked.get_xdata()[0])
508            if  lineNo in [0,1]:
509                LimitId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Limits')
510                data = G2frame.PatternTree.GetItemPyData(LimitId)
511                if G2frame.qPlot:
512                    data[1][lineNo] = 2.0*asind(wave*xpos/(4*math.pi))
513                else:
514                    data[1][lineNo] = xpos
515                G2frame.PatternTree.SetItemPyData(LimitId,data)
516                if G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Limits':
517                    G2pdG.UpdateLimitsGrid(G2frame,data)
518            else:
519                PeakId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Peak List')
520                data = G2frame.PatternTree.GetItemPyData(PeakId)
521#                print 'peaks',xpos
522                if event.button == 3:
523                    del data[lineNo-2]
524                else:
525                    if G2frame.qPlot:
526                        data[lineNo-2][0] = 2.0*asind(wave*xpos/(4*math.pi))
527                    else:
528                        data[lineNo-2][0] = xpos
529                G2frame.PatternTree.SetItemPyData(PeakId,data)
530                G2pdG.UpdatePeakGrid(G2frame,data)
531        elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' and xpos:
532            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
533            pick = str(G2frame.itemPicked).split('(')[1].strip(')')
534            if 'Line' not in pick:       #avoid data points, etc.
535                num = Phases.keys().index(pick)
536                if num:
537                    G2frame.refDelt = -(event.ydata-G2frame.refOffset)/(num*Ymax)
538                else:       #1st row of refl ticks
539                    G2frame.refOffset = event.ydata
540        PlotPatterns(G2frame)
541        G2frame.itemPicked = None   
542
543    xylim = []
544    try:
545        plotNum = G2frame.G2plotNB.plotList.index('Powder Patterns')
546        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
547        if not newPlot:
548            Plot = Page.figure.gca()          #get previous powder plot & get limits
549            xylim = Plot.get_xlim(),Plot.get_ylim()
550        Page.figure.clf()
551        Plot = Page.figure.gca()          #get a fresh plot after clf()
552    except ValueError:
553        newPlot = True
554        G2frame.Cmax = 1.0
555        Plot = G2frame.G2plotNB.addMpl('Powder Patterns').gca()
556        plotNum = G2frame.G2plotNB.plotList.index('Powder Patterns')
557        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
558        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
559        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
560        Page.canvas.mpl_connect('pick_event', OnPick)
561        Page.canvas.mpl_connect('button_release_event', OnRelease)
562    Page.SetFocus()
563    G2frame.G2plotNB.status.DestroyChildren()
564    if G2frame.Contour:
565        Choice = (' key press','d: lower contour max','u: raise contour max',
566            'i: interpolation method','s: color scheme','c: contour off')
567    else:
568        if G2frame.logPlot:
569            Choice = (' key press','n: log(I) off','l: offset left','r: offset right',
570                'c: contour on','q: toggle q plot','s: toggle single plot','+: no selection')
571        else:
572            Choice = (' key press','l: offset left','r: offset right','d: offset down',
573                'u: offset up','o: reset offset','n: log(I) on','c: contour on',
574                'q: toggle q plot','s: toggle single plot','w: toggle error analysis','+: no selection')
575    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
576        choices=Choice)
577    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
578    cb.SetValue(' key press')
579   
580    PickId = G2frame.PickId
581    PatternId = G2frame.PatternId
582    colors=['b','g','r','c','m','k']
583    Lines = []
584    if G2frame.SinglePlot:
585        Pattern = G2frame.PatternTree.GetItemPyData(PatternId)
586        Pattern.append(G2frame.PatternTree.GetItemText(PatternId))
587        PlotList = [Pattern,]
588        ParmList = [G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
589            G2frame.PatternId, 'Instrument Parameters'))[1],]
590    else:       
591        PlotList = []
592        ParmList = []
593        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
594        while item:
595            if 'PWDR' in G2frame.PatternTree.GetItemText(item):
596                Pattern = G2frame.PatternTree.GetItemPyData(item)
597                if len(Pattern) < 3:                    # put name on end if needed
598                    Pattern.append(G2frame.PatternTree.GetItemText(item))
599                PlotList.append(Pattern)
600                ParmList.append(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
601                    item,'Instrument Parameters'))[1])
602            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)               
603    Ymax = 1.0
604    lenX = 0
605    if G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
606        Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
607        HKL = []
608        if Phases:
609            for peak in Phases[G2frame.RefList]:
610                HKL.append(peak[:6])
611            HKL = np.array(HKL)
612    else:
613        HKL = np.array(G2frame.HKL)
614    for Pattern in PlotList:
615        xye = Pattern[1]
616        Ymax = max(Ymax,max(xye[1]))
617    offset = G2frame.Offset[0]*Ymax/100.0
618    Title = 'Powder Patterns: '+os.path.split(G2frame.GSASprojectfile)[1]
619    if G2frame.logPlot:
620        Title = 'log('+Title+')'
621    Plot.set_title(Title)
622    if G2frame.qPlot:
623        Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
624    else:       
625        Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
626    Plot.set_ylabel('Intensity',fontsize=12)
627    if G2frame.Contour:
628        ContourZ = []
629        ContourY = []
630        Nseq = 0
631    if len(PlotList) < 2:
632        G2frame.Contour = False
633    for N,Pattern in enumerate(PlotList):
634        Parms = ParmList[N]
635        ifpicked = False
636        LimitId = 0
637        xye = np.array(Pattern[1])
638        if PickId:
639            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
640            LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
641        if G2frame.qPlot:
642            Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, Pattern[2])
643            X = 4*np.pi*npsind(xye[0]/2.0)/Parms[1]
644        else:
645            X = xye[0]
646        if not lenX:
647            lenX = len(X)           
648        Y = xye[1]+offset*N
649        if LimitId:
650            limits = np.array(G2frame.PatternTree.GetItemPyData(LimitId))
651            if G2frame.qPlot:
652                limits = 4*np.pi*npsind(limits/2.0)/Parms[1]
653            Lines.append(Plot.axvline(limits[1][0],color='g',dashes=(5,5),picker=3.))   
654            Lines.append(Plot.axvline(limits[1][1],color='r',dashes=(5,5),picker=3.))                   
655        if G2frame.Contour:
656           
657            if lenX == len(X):
658                ContourY.append(N)
659                ContourZ.append(Y)
660                ContourX = X
661                Nseq += 1
662                Plot.set_ylabel('Data sequence',fontsize=12)
663        else:
664            X += G2frame.Offset[1]*.005*N
665            if ifpicked:
666                Z = xye[3]+offset*N
667                W = xye[4]+offset*N
668                D = xye[5]-Ymax*G2frame.delOffset
669                if G2frame.Weight:
670                    DS = xye[5]*np.sqrt(xye[2])
671                    PlotDeltSig(G2frame,[X,DS])
672                if G2frame.logPlot:
673                    Plot.semilogy(X,Y,colors[N%6]+'+',picker=3.,clip_on=False,nonposy='mask')
674                    Plot.semilogy(X,Z,colors[(N+1)%6],picker=False,nonposy='mask')
675                    Plot.semilogy(X,W,colors[(N+2)%6],picker=False,nonposy='mask')
676                else:
677                    Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
678                    Plot.plot(X,Z,colors[(N+1)%6],picker=False)
679                    Plot.plot(X,W,colors[(N+2)%6],picker=False)
680                    Plot.plot(X,D,colors[(N+3)%6],picker=False)
681                    Plot.axhline(0.,color=wx.BLACK)
682                Page.canvas.SetToolTipString('')
683                if G2frame.PatternTree.GetItemText(PickId) == 'Peak List':
684                    tip = 'On data point: Pick peak - L or R MB. On line: L-move, R-delete'
685                    Page.canvas.SetToolTipString(tip)
686                    data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'))
687                    for item in data:
688                        if G2frame.qPlot:
689                            Lines.append(Plot.axvline(4*math.pi*sind(item[0]/2.)/Parms[1],color=colors[N%6],picker=2.))
690                        else:
691                            Lines.append(Plot.axvline(item[0],color=colors[N%6],picker=2.))
692                if G2frame.PatternTree.GetItemText(PickId) == 'Limits':
693                    tip = 'On data point: Lower limit - L MB; Upper limit - R MB. On limit: MB down to move'
694                    Page.canvas.SetToolTipString(tip)
695                    data = G2frame.LimitsTable.GetData()
696            else:
697                if G2frame.logPlot:
698                    Plot.semilogy(X,Y,colors[N%6],picker=False,nonposy='clip')
699                else:
700                    Plot.plot(X,Y,colors[N%6],picker=False)
701    if PickId and not G2frame.Contour:
702        Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))[1::2]
703        Parms = dict(zip(Names,Values))
704        try:
705            wave = Parms['Lam']
706        except KeyError:
707            wave = Parms['Lam1']
708        if G2frame.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
709            peaks = np.array((G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))))
710            for peak in peaks:
711                if G2frame.qPlot:
712                    Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
713                else:
714                    Plot.axvline(peak[0],color='b')
715            for hkl in G2frame.HKL:
716                if G2frame.qPlot:
717                    Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
718                else:
719                    Plot.axvline(hkl[5],color='r',dashes=(5,5))
720        elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
721            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
722            for pId,phase in enumerate(Phases):
723                peaks = Phases[phase]
724                peak = np.array([[peak[4],peak[5]] for peak in peaks])
725                pos = G2frame.refOffset-pId*Ymax*G2frame.refDelt*np.ones_like(peak)
726                if G2frame.qPlot:
727                    Plot.plot(2*np.pi/peak.T[0],pos,colors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
728                else:
729                    Plot.plot(peak.T[1],pos,colors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
730            handles,legends = Plot.get_legend_handles_labels()  #got double entries in the legends for some reason
731            Plot.legend(handles[::2],legends[::2],title='Phases',loc='best')    #skip every other one
732           
733    if G2frame.Contour:
734        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
735        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
736            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
737        Page.figure.colorbar(Img)
738    else:
739        G2frame.Lines = Lines
740    if not newPlot:
741        Page.toolbar.push_current()
742        Plot.set_xlim(xylim[0])
743        Plot.set_ylim(xylim[1])
744        xylim = []
745        Page.toolbar.push_current()
746        Page.toolbar.draw()
747    else:
748        Page.canvas.draw()
749    G2frame.Pwdr = True
750   
751def PlotDeltSig(G2frame,XDS):
752    try:
753        plotNum = G2frame.G2plotNB.plotList.index('Error analysis')
754        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
755        Page.figure.clf()
756        Plot = Page.figure.gca()          #get a fresh plot after clf()
757    except ValueError:
758        newPlot = True
759        G2frame.Cmax = 1.0
760        Plot = G2frame.G2plotNB.addMpl('Error analysis').gca()
761        plotNum = G2frame.G2plotNB.plotList.index('Error analysis')
762        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
763    Page.SetFocus()
764    G2frame.G2plotNB.status.DestroyChildren()
765    Plot.set_title('Delta/sigma')
766    X,DS = XDS
767    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
768    Plot.set_ylabel(r'$\mathsf{\Delta/\sigma}$',fontsize=14)
769    Plot.plot(X,DS,'k',picker=False)
770    Page.canvas.draw()
771   
772   
773def PlotISFG(G2frame,newPlot=False,type=''):
774    ''' PLotting package for PDF analysis; displays I(q), S(q), F(q) and G(r) as single
775    or multiple plots with waterfall and contour plots as options
776    '''
777    if not type:
778        type = G2frame.G2plotNB.plotList[G2frame.G2plotNB.nb.GetSelection()]
779    if type not in ['I(Q)','S(Q)','F(Q)','G(R)']:
780        return
781    superMinusOne = unichr(0xaf)+unichr(0xb9)
782   
783    def OnPlotKeyPress(event):
784        newPlot = False
785        if event.key == 'u':
786            if G2frame.Contour:
787                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
788            elif G2frame.Offset[0] < 100.:
789                G2frame.Offset[0] += 1.
790        elif event.key == 'd':
791            if G2frame.Contour:
792                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
793            elif G2frame.Offset[0] > 0.:
794                G2frame.Offset[0] -= 1.
795        elif event.key == 'l':
796            G2frame.Offset[1] -= 1.
797        elif event.key == 'r':
798            G2frame.Offset[1] += 1.
799        elif event.key == 'o':
800            G2frame.Offset = [0,0]
801        elif event.key == 'c':
802            newPlot = True
803            if G2frame.Contour:
804                G2frame.Contour = False
805            else:
806                G2frame.Contour = True
807                G2frame.SinglePlot = False
808                G2frame.Offset = [0.,0.]
809        elif event.key == 's':
810            if G2frame.Contour:
811                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
812                choice.sort()
813                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
814                if dlg.ShowModal() == wx.ID_OK:
815                    sel = dlg.GetSelection()
816                    G2frame.ContourColor = choice[sel]
817                else:
818                    G2frame.ContourColor = 'Paired'
819                dlg.Destroy()
820            else:               
821                if G2frame.SinglePlot:
822                    G2frame.SinglePlot = False
823                else:
824                    G2frame.SinglePlot = True
825        elif event.key == 'i':                  #for smoothing contour plot
826            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
827               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
828               'mitchell','sinc','lanczos']
829            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
830            if dlg.ShowModal() == wx.ID_OK:
831                sel = dlg.GetSelection()
832                G2frame.Interpolate = choice[sel]
833            else:
834                G2frame.Interpolate = 'nearest'
835            dlg.Destroy()
836        elif event.key == 't' and not G2frame.Contour:
837            if G2frame.Legend:
838                G2frame.Legend = False
839            else:
840                G2frame.Legend = True
841           
842           
843        PlotISFG(G2frame,newPlot=newPlot,type=type)
844       
845    def OnKeyBox(event):
846        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index(type):
847            event.key = cb.GetValue()[0]
848            cb.SetValue(' key press')
849            OnPlotKeyPress(event)
850                       
851    def OnMotion(event):
852        xpos = event.xdata
853        if xpos:                                        #avoid out of frame mouse position
854            ypos = event.ydata
855            Page.canvas.SetCursor(wx.CROSS_CURSOR)
856            try:
857                if G2frame.Contour:
858                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA pattern ID =%5d'%(xpos,int(ypos)),1)
859                else:
860                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA %s =%.2f'%(xpos,type,ypos),1)                   
861            except TypeError:
862                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
863   
864    xylim = []
865    try:
866        plotNum = G2frame.G2plotNB.plotList.index(type)
867        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
868        if not newPlot:
869            Plot = Page.figure.gca()          #get previous plot & get limits
870            xylim = Plot.get_xlim(),Plot.get_ylim()
871        Page.figure.clf()
872        Plot = Page.figure.gca()
873    except ValueError:
874        newPlot = True
875        G2frame.Cmax = 1.0
876        Plot = G2frame.G2plotNB.addMpl(type).gca()
877        plotNum = G2frame.G2plotNB.plotList.index(type)
878        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
879        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
880        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
881   
882    Page.SetFocus()
883    G2frame.G2plotNB.status.DestroyChildren()
884    if G2frame.Contour:
885        Choice = (' key press','d: lower contour max','u: raise contour max',
886            'i: interpolation method','s: color scheme','c: contour off')
887    else:
888        Choice = (' key press','l: offset left','r: offset right','d: offset down','u: offset up',
889            'o: reset offset','t: toggle legend','c: contour on','s: toggle single plot')
890    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
891        choices=Choice)
892    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
893    cb.SetValue(' key press')
894    PatternId = G2frame.PatternId
895    PickId = G2frame.PickId
896    Plot.set_title(type)
897    if type == 'G(R)':
898        Plot.set_xlabel(r'$R,\AA$',fontsize=14)
899    else:
900        Plot.set_xlabel(r'$Q,\AA$'+superMinusOne,fontsize=14)
901    Plot.set_ylabel(r''+type,fontsize=14)
902    colors=['b','g','r','c','m','k']
903    name = G2frame.PatternTree.GetItemText(PatternId)[4:]
904    Pattern = []   
905    if G2frame.SinglePlot:
906        name = G2frame.PatternTree.GetItemText(PatternId)
907        name = type+name[4:]
908        Id = G2gd.GetPatternTreeItemId(G2frame,PatternId,name)
909        Pattern = G2frame.PatternTree.GetItemPyData(Id)
910        if Pattern:
911            Pattern.append(name)
912        PlotList = [Pattern,]
913    else:
914        PlotList = []
915        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
916        while item:
917            if 'PDF' in G2frame.PatternTree.GetItemText(item):
918                name = type+G2frame.PatternTree.GetItemText(item)[4:]
919                Id = G2gd.GetPatternTreeItemId(G2frame,item,name)
920                Pattern = G2frame.PatternTree.GetItemPyData(Id)
921                if Pattern:
922                    Pattern.append(name)
923                    PlotList.append(Pattern)
924            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
925    PDFdata = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'PDF Controls'))
926    numbDen = G2pwd.GetNumDensity(PDFdata['ElList'],PDFdata['Form Vol'])
927    Xb = [0.,10.]
928    Yb = [0.,-40.*np.pi*numbDen]
929    Ymax = 1.0
930    lenX = 0
931    for Pattern in PlotList:
932        xye = Pattern[1]
933        Ymax = max(Ymax,max(xye[1]))
934    offset = G2frame.Offset[0]*Ymax/100.0
935    if G2frame.Contour:
936        ContourZ = []
937        ContourY = []
938        Nseq = 0
939    for N,Pattern in enumerate(PlotList):
940        xye = Pattern[1]
941        if PickId:
942            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
943        X = xye[0]
944        if not lenX:
945            lenX = len(X)           
946        Y = xye[1]+offset*N
947        if G2frame.Contour:
948            if lenX == len(X):
949                ContourY.append(N)
950                ContourZ.append(Y)
951                ContourX = X
952                Nseq += 1
953                Plot.set_ylabel('Data sequence',fontsize=12)
954        else:
955            X = xye[0]+G2frame.Offset[1]*.005*N
956            if ifpicked:
957                Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
958                Page.canvas.SetToolTipString('')
959            else:
960                if G2frame.Legend:
961                    Plot.plot(X,Y,colors[N%6],picker=False,label='Azm:'+Pattern[2].split('=')[1])
962                else:
963                    Plot.plot(X,Y,colors[N%6],picker=False)
964            if type == 'G(R)':
965                Plot.plot(Xb,Yb,color='k',dashes=(5,5))
966            elif type == 'F(Q)':
967                Plot.axhline(0.,color=wx.BLACK)
968            elif type == 'S(Q)':
969                Plot.axhline(1.,color=wx.BLACK)
970    if G2frame.Contour:
971        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
972        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
973            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
974        Page.figure.colorbar(Img)
975    elif G2frame.Legend:
976        Plot.legend(loc='best')
977    if not newPlot:
978        Page.toolbar.push_current()
979        Plot.set_xlim(xylim[0])
980        Plot.set_ylim(xylim[1])
981        xylim = []
982        Page.toolbar.push_current()
983        Page.toolbar.draw()
984    else:
985        Page.canvas.draw()
986       
987def PlotXY(G2frame,XY,newPlot=False,type=''):
988    '''simple plot of xy data, used for diagnostic purposes
989    '''
990    def OnMotion(event):
991        xpos = event.xdata
992        if xpos:                                        #avoid out of frame mouse position
993            ypos = event.ydata
994            Page.canvas.SetCursor(wx.CROSS_CURSOR)
995            try:
996                G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,type,ypos),1)                   
997            except TypeError:
998                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
999
1000    try:
1001        plotNum = G2frame.G2plotNB.plotList.index(type)
1002        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1003        if not newPlot:
1004            Plot = Page.figure.gca()
1005            xylim = Plot.get_xlim(),Plot.get_ylim()
1006        Page.figure.clf()
1007        Plot = Page.figure.gca()
1008    except ValueError:
1009        newPlot = True
1010        Plot = G2frame.G2plotNB.addMpl(type).gca()
1011        plotNum = G2frame.G2plotNB.plotList.index(type)
1012        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1013        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1014   
1015    Page.SetFocus()
1016    G2frame.G2plotNB.status.DestroyChildren()
1017    Plot.set_title(type)
1018    Plot.set_xlabel(r'X',fontsize=14)
1019    Plot.set_ylabel(r''+type,fontsize=14)
1020    colors=['b','g','r','c','m','k']
1021    Ymax = 1.0
1022    lenX = 0
1023    X,Y = XY[:2]
1024    Ymax = max(Ymax,max(Y))
1025    Plot.plot(X,Y,'k',picker=False)
1026    if not newPlot:
1027        Page.toolbar.push_current()
1028        Plot.set_xlim(xylim[0])
1029        Plot.set_ylim(xylim[1])
1030        xylim = []
1031        Page.toolbar.push_current()
1032        Page.toolbar.draw()
1033    else:
1034        Page.canvas.draw()
1035
1036def PlotPowderLines(G2frame):
1037    ''' plotting of powder lines (i.e. no powder pattern) as sticks
1038    '''
1039    global HKL
1040
1041    def OnMotion(event):
1042        xpos = event.xdata
1043        if xpos:                                        #avoid out of frame mouse position
1044            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1045            G2frame.G2plotNB.status.SetFields(['','2-theta =%9.3f '%(xpos,)])
1046            if G2frame.PickId and G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List']:
1047                found = []
1048                if len(HKL):
1049                    view = Page.toolbar._views.forward()[0][:2]
1050                    wid = view[1]-view[0]
1051                    found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
1052                if len(found):
1053                    h,k,l = found[0][:3] 
1054                    Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
1055                else:
1056                    Page.canvas.SetToolTipString('')
1057
1058    try:
1059        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
1060        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1061        Page.figure.clf()
1062        Plot = Page.figure.gca()
1063    except ValueError:
1064        Plot = G2frame.G2plotNB.addMpl('Powder Lines').gca()
1065        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
1066        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1067        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1068       
1069    Page.SetFocus()
1070    Plot.set_title('Powder Pattern Lines')
1071    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
1072    PickId = G2frame.PickId
1073    PatternId = G2frame.PatternId
1074    peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))
1075    for peak in peaks:
1076        Plot.axvline(peak[0],color='b')
1077    HKL = np.array(G2frame.HKL)
1078    for hkl in G2frame.HKL:
1079        Plot.axvline(hkl[5],color='r',dashes=(5,5))
1080    xmin = peaks[0][0]
1081    xmax = peaks[-1][0]
1082    delt = xmax-xmin
1083    xlim = [max(0,xmin-delt/20.),min(180.,xmax+delt/20.)]
1084    Plot.set_xlim(xlim)
1085    Page.canvas.draw()
1086    Page.toolbar.push_current()
1087
1088def PlotPeakWidths(G2frame):
1089    ''' Plotting of instrument broadening terms as function of 2-theta (future TOF)
1090    Seen when "Instrument Parameters" chosen from powder pattern data tree
1091    '''
1092    PatternId = G2frame.PatternId
1093    limitID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
1094    if limitID:
1095        limits = G2frame.PatternTree.GetItemPyData(limitID)
1096    else:
1097        return
1098    instParms = G2frame.PatternTree.GetItemPyData( \
1099        G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
1100    if instParms[0][0] in ['PXC','PNC']:
1101        lam = instParms[1][1]
1102        if len(instParms[1]) == 13:
1103            GU,GV,GW,LX,LY = instParms[0][6:11]
1104        else:
1105            GU,GV,GW,LX,LY = instParms[0][4:9]
1106    peakID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List')
1107    if peakID:
1108        peaks = G2frame.PatternTree.GetItemPyData(peakID)
1109    else:
1110        peaks = []
1111   
1112    try:
1113        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
1114        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1115        Page.figure.clf()
1116        Plot = Page.figure.gca()
1117    except ValueError:
1118        Plot = G2frame.G2plotNB.addMpl('Peak Widths').gca()
1119        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
1120        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1121    Page.SetFocus()
1122   
1123    Page.canvas.SetToolTipString('')
1124    colors=['b','g','r','c','m','k']
1125    Xmin,Xmax = limits[1]
1126    Xmin = min(0.5,max(Xmin,1))
1127    Xmin /= 2
1128    Xmax /= 2
1129    nPts = 100
1130    delt = (Xmax-Xmin)/nPts
1131    thetas = []
1132    for i in range(nPts):
1133        thetas.append(Xmin+i*delt)
1134    X = []
1135    Y = []
1136    Z = []
1137    W = []
1138    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
1139    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
1140    gamFW = lambda s,g: math.exp(math.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.)
1141#    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
1142    for theta in thetas:
1143        X.append(4.0*math.pi*sind(theta)/lam)              #q
1144        s = sig(theta,GU,GV,GW)
1145        g = gam(theta,LX,LY)
1146        G = gamFW(g,s)
1147        Y.append(s/tand(theta))
1148        Z.append(g/tand(theta))
1149        W.append(G/tand(theta))
1150    Plot.set_title('Instrument and sample peak widths')
1151    Plot.set_ylabel(r'$\Delta q/q, \Delta d/d$',fontsize=14)
1152    Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
1153    Plot.plot(X,Y,color='r',label='Gaussian')
1154    Plot.plot(X,Z,color='g',label='Lorentzian')
1155    Plot.plot(X,W,color='b',label='G+L')
1156    X = []
1157    Y = []
1158    Z = []
1159    W = []
1160    V = []
1161    for peak in peaks:
1162        X.append(4.0*math.pi*sind(peak[0]/2.0)/lam)
1163        try:
1164            s = 1.17741*math.sqrt(peak[4])*math.pi/18000.
1165        except ValueError:
1166            s = 0.01
1167        g = peak[6]*math.pi/18000.
1168        G = gamFW(g,s)
1169        Y.append(s/tand(peak[0]/2.))
1170        Z.append(g/tand(peak[0]/2.))
1171        W.append(G/tand(peak[0]/2.))
1172    Plot.plot(X,Y,'+',color='r',label='G peak')
1173    Plot.plot(X,Z,'+',color='g',label='L peak')
1174    Plot.plot(X,W,'+',color='b',label='G+L peak')
1175    Plot.legend(loc='best')
1176    Page.canvas.draw()
1177   
1178def PlotSizeStrainPO(G2frame,data,Start=False):
1179    '''Plot 3D mustrain/size/preferred orientation figure. In this instance data is for a phase
1180    '''
1181   
1182    PatternId = G2frame.PatternId
1183    generalData = data['General']
1184    SGData = generalData['SGData']
1185    SGLaue = SGData['SGLaue']
1186    if Start:                   #initialize the spherical harmonics qlmn arrays
1187        ptx.pyqlmninit()
1188        Start = False
1189#    MuStrKeys = G2spc.MustrainNames(SGData)
1190    cell = generalData['Cell'][1:]
1191    A,B = G2lat.cell2AB(cell[:6])
1192    Vol = cell[6]
1193    useList = data['Histograms']
1194    phase = generalData['Name']
1195    plotType = generalData['Data plot type']
1196    plotDict = {'Mustrain':'Mustrain','Size':'Size','Preferred orientation':'Pref.Ori.'}
1197    for ptype in plotDict:
1198        G2frame.G2plotNB.Delete(ptype)       
1199
1200    for item in useList:
1201        if useList[item]['Show']:
1202            break
1203    else:
1204        return            #nothing to show!!
1205   
1206    numPlots = len(useList)
1207
1208    if plotType in ['Mustrain','Size']:
1209        Plot = mp3d.Axes3D(G2frame.G2plotNB.add3D(plotType))
1210    else:
1211        Plot = G2frame.G2plotNB.addMpl(plotType).gca()       
1212    plotNum = G2frame.G2plotNB.plotList.index(plotType)
1213    Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1214    Page.SetFocus()
1215    G2frame.G2plotNB.status.SetStatusText('Adjust frame size to get desired aspect ratio',1)
1216    if not Page.IsShown():
1217        Page.Show()
1218   
1219    for item in useList:
1220        if useList[item]['Show']:
1221            PHI = np.linspace(0.,360.,30,True)
1222            PSI = np.linspace(0.,180.,30,True)
1223            X = np.outer(npsind(PHI),npsind(PSI))
1224            Y = np.outer(npcosd(PHI),npsind(PSI))
1225            Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
1226            try:        #temp patch instead of 'mustrain' for old files with 'microstrain'
1227                coeff = useList[item][plotDict[plotType]]
1228            except KeyError:
1229                break
1230            if plotType in ['Mustrain','Size']:
1231                if coeff[0] == 'isotropic':
1232                    X *= coeff[1][0]
1233                    Y *= coeff[1][0]
1234                    Z *= coeff[1][0]                               
1235                elif coeff[0] == 'uniaxial':
1236                   
1237                    def uniaxCalc(xyz,iso,aniso,axes):
1238                        Z = np.array(axes)
1239                        cp = abs(np.dot(xyz,Z))
1240                        sp = np.sqrt(1.-cp**2)
1241                        R = iso*aniso/np.sqrt((iso*cp)**2+(aniso*sp)**2)
1242                        return R*xyz
1243                       
1244                    iso,aniso = coeff[1][:2]
1245                    axes = np.inner(A,np.array(coeff[3]))
1246                    axes /= nl.norm(axes)
1247                    Shkl = np.array(coeff[1])
1248                    XYZ = np.dstack((X,Y,Z))
1249                    XYZ = np.nan_to_num(np.apply_along_axis(uniaxCalc,2,XYZ,iso,aniso,axes))
1250                    X,Y,Z = np.dsplit(XYZ,3)
1251                    X = X[:,:,0]
1252                    Y = Y[:,:,0]
1253                    Z = Z[:,:,0]
1254               
1255                elif coeff[0] == 'ellipsoidal':
1256                   
1257                    def ellipseCalc(xyz,E,R):
1258                        XYZ = xyz*E.T
1259                        return np.inner(XYZ.T,R)
1260                       
1261                    S6 = coeff[4]
1262                    Sij = G2lat.U6toUij(S6)
1263                    E,R = nl.eigh(Sij)
1264                    XYZ = np.dstack((X,Y,Z))
1265                    XYZ = np.nan_to_num(np.apply_along_axis(ellipseCalc,2,XYZ,E,R))
1266                    X,Y,Z = np.dsplit(XYZ,3)
1267                    X = X[:,:,0]
1268                    Y = Y[:,:,0]
1269                    Z = Z[:,:,0]
1270                   
1271                elif coeff[0] == 'generalized':
1272                   
1273                    def genMustrain(xyz,SGData,A,Shkl):
1274                        uvw = np.inner(A.T,xyz)
1275                        Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
1276                        sum = np.sum(np.multiply(Shkl,Strm))
1277                        sum = np.where(sum > 0.01,sum,0.01)
1278                        sum = np.sqrt(sum)*math.pi/0.018      #centidegrees to radians!
1279                        return sum*xyz
1280                       
1281                    Shkl = np.array(coeff[4])
1282                    if np.any(Shkl):
1283                        XYZ = np.dstack((X,Y,Z))
1284                        XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,SGData,A,Shkl))
1285                        X,Y,Z = np.dsplit(XYZ,3)
1286                        X = X[:,:,0]
1287                        Y = Y[:,:,0]
1288                        Z = Z[:,:,0]
1289                           
1290                if np.any(X) and np.any(Y) and np.any(Z):
1291                    errFlags = np.seterr(all='ignore')
1292                    Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1)
1293                    np.seterr(all='ignore')
1294                    xyzlim = np.array([Plot.get_xlim3d(),Plot.get_ylim3d(),Plot.get_zlim3d()]).T
1295                    XYZlim = [min(xyzlim[0]),max(xyzlim[1])]
1296                    Plot.set_xlim3d(XYZlim)
1297                    Plot.set_ylim3d(XYZlim)
1298                    Plot.set_zlim3d(XYZlim)
1299                    Plot.set_aspect('equal')
1300                if plotType == 'Size':
1301                    Plot.set_title('Crystallite size for '+phase+'\n'+coeff[0]+' model')
1302                    Plot.set_xlabel(r'X, $\mu$m')
1303                    Plot.set_ylabel(r'Y, $\mu$m')
1304                    Plot.set_zlabel(r'Z, $\mu$m')
1305                else:   
1306                    Plot.set_title(r'$\mu$strain for '+phase+'\n'+coeff[0]+' model')
1307                    Plot.set_xlabel(r'X, $\mu$strain')
1308                    Plot.set_ylabel(r'Y, $\mu$strain')
1309                    Plot.set_zlabel(r'Z, $\mu$strain')
1310            else:
1311                h,k,l = generalData['POhkl']
1312                if coeff[0] == 'MD':
1313                    print 'March-Dollase preferred orientation plot'
1314               
1315                else:
1316                    PH = np.array(generalData['POhkl'])
1317                    phi,beta = G2lat.CrsAng(PH,cell[:6],SGData)
1318                    SHCoef = {}
1319                    for item in coeff[5]:
1320                        L,N = eval(item.strip('C'))
1321                        SHCoef['C%d,0,%d'%(L,N)] = coeff[5][item]                       
1322                    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1323                    X = np.linspace(0,90.0,26)
1324                    Y = G2lat.polfcal(ODFln,'0',X,0.0)
1325                    Plot.plot(X,Y,color='k',label=str(PH))
1326                    Plot.legend(loc='best')
1327                    Plot.set_title('Axial distribution for HKL='+str(PH)+' in '+phase)
1328                    Plot.set_xlabel(r'$\psi$',fontsize=16)
1329                    Plot.set_ylabel('MRD',fontsize=14)
1330    Page.canvas.draw()
1331   
1332def PlotTexture(G2frame,data,newPlot=False,Start=False):
1333    '''Pole figure, inverse pole figure, 3D pole distribution and 3D inverse pole distribution
1334    plotting.
1335    dict generalData contains all phase info needed which is in data
1336    '''
1337
1338    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1339    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1340    PatternId = G2frame.PatternId
1341    generalData = data['General']
1342    SGData = generalData['SGData']
1343    textureData = generalData['SH Texture']
1344    G2frame.G2plotNB.Delete('Texture')
1345    if not textureData['Order']:
1346        return                  #no plot!!
1347    SHData = generalData['SH Texture']
1348    SHCoef = SHData['SH Coeff'][1]
1349    cell = generalData['Cell'][1:7]
1350    Amat,Bmat = G2lat.cell2AB(cell)
1351    sq2 = 1.0/math.sqrt(2.0)
1352   
1353    def rp2xyz(r,p):
1354        z = npcosd(r)
1355        xy = np.sqrt(1.-z**2)
1356        return xy*npsind(p),xy*npcosd(p),z
1357           
1358    def OnMotion(event):
1359        SHData = data['General']['SH Texture']
1360        if event.xdata and event.ydata:                 #avoid out of frame errors
1361            xpos = event.xdata
1362            ypos = event.ydata
1363            if 'Inverse' in SHData['PlotType']:
1364                r = xpos**2+ypos**2
1365                if r <= 1.0:
1366                    if 'equal' in G2frame.Projection: 
1367                        r,p = 2.*npasind(np.sqrt(r)*sq2),npatan2d(ypos,xpos)
1368                    else:
1369                        r,p = 2.*npatand(np.sqrt(r)),npatan2d(ypos,xpos)
1370                    ipf = G2lat.invpolfcal(IODFln,SGData,np.array([r,]),np.array([p,]))
1371                    xyz = np.inner(Amat.T,np.array([rp2xyz(r,p)]))
1372                    y,x,z = list(xyz/np.max(np.abs(xyz)))
1373                   
1374                    G2frame.G2plotNB.status.SetFields(['',
1375                        'psi =%9.3f, beta =%9.3f, MRD =%9.3f hkl=%5.2f,%5.2f,%5.2f'%(r,p,ipf,x,y,z)])
1376                                   
1377            elif 'Axial' in SHData['PlotType']:
1378                pass
1379               
1380            else:                       #ordinary pole figure
1381                z = xpos**2+ypos**2
1382                if z <= 1.0:
1383                    z = np.sqrt(z)
1384                    if 'equal' in G2frame.Projection: 
1385                        r,p = 2.*npasind(z*sq2),npatan2d(ypos,xpos)
1386                    else:
1387                        r,p = 2.*npatand(z),npatan2d(ypos,xpos)
1388                    pf = G2lat.polfcal(ODFln,SamSym[textureData['Model']],np.array([r,]),np.array([p,]))
1389                    G2frame.G2plotNB.status.SetFields(['','phi =%9.3f, gam =%9.3f, MRD =%9.3f'%(r,p,pf)])
1390   
1391    if G2frame.Projection == '3D display' and 'Axial'  not in SHData['PlotType']:               
1392        Plot = mp3d.Axes3D(G2frame.G2plotNB.add3D('Texture'))
1393    else:
1394        Plot = G2frame.G2plotNB.addMpl('Texture').gca()
1395    plotNum = G2frame.G2plotNB.plotList.index('Texture')
1396    Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1397    if G2frame.Projection == '3D display':
1398        pass
1399    else:
1400        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1401
1402    Page.SetFocus()
1403    G2frame.G2plotNB.status.SetFields(['',''])   
1404    PH = np.array(SHData['PFhkl'])
1405    phi,beta = G2lat.CrsAng(PH,cell,SGData)
1406    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1407    PX = np.array(SHData['PFxyz'])
1408    gam = atan2d(PX[0],PX[1])
1409    xy = math.sqrt(PX[0]**2+PX[1]**2)
1410    xyz = math.sqrt(PX[0]**2+PX[1]**2+PX[2]**2)
1411    psi = asind(xy/xyz)
1412    IODFln = G2lat.Glnh(Start,SHCoef,psi,gam,SamSym[textureData['Model']])
1413    if 'Axial' in SHData['PlotType']:
1414        X = np.linspace(0,90.0,26)
1415        Y = G2lat.polfcal(ODFln,SamSym[textureData['Model']],X,0.0)
1416        Plot.plot(X,Y,color='k',label=str(SHData['PFhkl']))
1417        Plot.legend(loc='best')
1418        Plot.set_title('Axial distribution for HKL='+str(SHData['PFhkl']))
1419        Plot.set_xlabel(r'$\psi$',fontsize=16)
1420        Plot.set_ylabel('MRD',fontsize=14)
1421       
1422    else:       
1423        npts = 201
1424        if 'Inverse' in SHData['PlotType']:
1425            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
1426            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
1427            if 'equal' in G2frame.Projection:
1428                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
1429            else:
1430                R = np.where(R <= 1.,2.*npatand(R),0.0)
1431            Z = np.zeros_like(R)
1432            Z = G2lat.invpolfcal(IODFln,SGData,R,P)
1433            Z = np.reshape(Z,(npts,npts))
1434            CS = Plot.contour(Y,X,Z,aspect='equal')
1435            Plot.clabel(CS,fontsize=9,inline=1)
1436            try:
1437                Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
1438            except ValueError:
1439                pass
1440            if newPlot:
1441#                Page.figure.colorbar(Img)    #colorbar fails - crashes gsasii
1442                newPlot = False
1443            Plot.set_title('Inverse pole figure for XYZ='+str(SHData['PFxyz']))
1444            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
1445                       
1446        else:
1447            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
1448            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
1449            if 'equal' in G2frame.Projection:
1450                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
1451            else:
1452                R = np.where(R <= 1.,2.*npatand(R),0.0)
1453            Z = np.zeros_like(R)
1454            Z = G2lat.polfcal(ODFln,SamSym[textureData['Model']],R,P)
1455            Z = np.reshape(Z,(npts,npts))
1456            try:
1457                CS = Plot.contour(Y,X,Z,aspect='equal')
1458                Plot.clabel(CS,fontsize=9,inline=1)
1459            except ValueError:
1460                pass
1461            Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
1462            if newPlot:
1463#                Page.figure.colorbar(Img)    #colorbar fails - crashes gsasii
1464                newPlot = False
1465            Plot.set_title('Pole figure for HKL='+str(SHData['PFhkl']))
1466            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
1467    Page.canvas.draw()
1468
1469def PlotCovariance(G2frame,Data={}):
1470    if not Data:
1471        Data = G2frame.PatternTree.GetItemPyData(
1472            G2gd.GetPatternTreeItemId(G2frame,G2frame.root, 'Covariance'))
1473    if not Data:
1474        print 'No covariance matrix available'
1475        return
1476    varyList = Data['varyList']
1477    values = Data['variables']
1478    Xmax = len(varyList)
1479    covMatrix = Data['covMatrix']
1480    sig = np.sqrt(np.diag(covMatrix))
1481    xvar = np.outer(sig,np.ones_like(sig))
1482    covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
1483    title = ' for\n'+Data['title']
1484    newAtomDict = Data['newAtomDict']
1485
1486    def OnPlotKeyPress(event):
1487        newPlot = False
1488        if event.key == 's':
1489            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
1490            choice.sort()
1491            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
1492            if dlg.ShowModal() == wx.ID_OK:
1493                sel = dlg.GetSelection()
1494                G2frame.VcovColor = choice[sel]
1495            else:
1496                G2frame.VcovColor = 'RdYlGn'
1497            dlg.Destroy()
1498        PlotCovariance(G2frame)
1499
1500    def OnMotion(event):
1501        if event.button:
1502            ytics = imgAx.get_yticks()
1503            ytics = np.where(ytics<len(varyList),ytics,-1)
1504            ylabs = [np.where(0<=i ,varyList[int(i)],' ') for i in ytics]
1505            imgAx.set_yticklabels(ylabs)
1506           
1507        if event.xdata and event.ydata:                 #avoid out of frame errors
1508            xpos = int(event.xdata+.5)
1509            ypos = int(event.ydata+.5)
1510            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
1511                if xpos == ypos:
1512                    value = values[xpos]
1513                    name = varyList[xpos]
1514                    if varyList[xpos] in newAtomDict:
1515                        name,value = newAtomDict[name]                       
1516                    msg = '%s value = %.4g, esd = %.4g'%(name,value,sig[xpos])
1517                else:
1518                    msg = '%s - %s: %5.3f'%(varyList[xpos],varyList[ypos],covArray[xpos][ypos])
1519                Page.canvas.SetToolTipString(msg)
1520                G2frame.G2plotNB.status.SetFields(['Key: s to change colors',msg])
1521    try:
1522        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
1523        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1524        Page.figure.clf()
1525        Plot = Page.figure.gca()
1526        if not Page.IsShown():
1527            Page.Show()
1528    except ValueError:
1529        Plot = G2frame.G2plotNB.addMpl('Covariance').gca()
1530        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
1531        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1532        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1533        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1534
1535    Page.SetFocus()
1536    G2frame.G2plotNB.status.SetFields(['',''])   
1537    acolor = mpl.cm.get_cmap(G2frame.VcovColor)
1538    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower')
1539    imgAx = Img.get_axes()
1540    ytics = imgAx.get_yticks()
1541    ylabs = [varyList[int(i)] for i in ytics[:-1]]
1542    imgAx.set_yticklabels(ylabs)
1543    colorBar = Page.figure.colorbar(Img)
1544    Plot.set_title('V-Cov matrix'+title)
1545    Plot.set_xlabel('Variable number')
1546    Plot.set_ylabel('Variable name')
1547    Page.canvas.draw()
1548   
1549def PlotSeq(G2frame,SeqData,SeqSig,SeqNames,sampleParm):
1550   
1551    def OnKeyPress(event):
1552        if event.key == 's' and sampleParm:
1553            if G2frame.xAxis:
1554                G2frame.xAxis = False
1555            else:
1556                G2frame.xAxis = True
1557            Draw(False)
1558    try:
1559        plotNum = G2frame.G2plotNB.plotList.index('Sequential refinement')
1560        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1561        Page.figure.clf()
1562        Plot = Page.figure.gca()
1563        if not Page.IsShown():
1564            Page.Show()
1565    except ValueError:
1566        Plot = G2frame.G2plotNB.addMpl('Sequential refinement').gca()
1567        plotNum = G2frame.G2plotNB.plotList.index('Sequential refinement')
1568        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1569        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
1570        G2frame.xAxis = False
1571       
1572    def Draw(newPlot):
1573        Page.SetFocus()
1574        G2frame.G2plotNB.status.SetFields(['','press s to toggle x-axis = sample environment parameter'])
1575        if len(SeqData):
1576            Plot.clear()
1577            if G2frame.xAxis:   
1578                xName = sampleParm.keys()[0]
1579                X = sampleParm[xName]
1580            else:
1581                X = np.arange(0,len(SeqData[0]),1)
1582                xName = 'Data sequence number'
1583            for Y,sig,name in zip(SeqData,SeqSig,SeqNames):
1584                Plot.errorbar(X,Y,yerr=sig,label=name)       
1585            Plot.legend(loc='best')
1586            Plot.set_ylabel('Parameter values')
1587            Plot.set_xlabel(xName)
1588            Page.canvas.draw()           
1589    Draw(True)
1590           
1591def PlotExposedImage(G2frame,newPlot=False,event=None):
1592    '''General access module for 2D image plotting
1593    '''
1594    plotNo = G2frame.G2plotNB.nb.GetSelection()
1595    if G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
1596        PlotImage(G2frame,newPlot,event,newImage=True)
1597    elif G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
1598        PlotIntegration(G2frame,newPlot,event)
1599
1600def PlotImage(G2frame,newPlot=False,event=None,newImage=True):
1601    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
1602    masks, etc.
1603    '''
1604    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
1605    import numpy.ma as ma
1606    Dsp = lambda tth,wave: wave/(2.*sind(tth/2.))
1607    global Data,Masks
1608    colors=['b','g','r','c','m','k']
1609    Data = G2frame.PatternTree.GetItemPyData(
1610        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
1611    Masks = G2frame.PatternTree.GetItemPyData(
1612        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
1613
1614    def OnImMotion(event):
1615        Page.canvas.SetToolTipString('')
1616        sizexy = Data['size']
1617        if event.xdata and event.ydata:                 #avoid out of frame errors
1618            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1619            item = G2frame.itemPicked
1620            pixelSize = Data['pixelSize']
1621            scalex = 1000./pixelSize[0]
1622            scaley = 1000./pixelSize[1]
1623            if item and G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Image Controls':
1624                if 'Text' in str(item):
1625                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
1626                else:
1627                    xcent,ycent = Data['center']
1628                    xpos = event.xdata-xcent
1629                    ypos = event.ydata-ycent
1630                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
1631                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
1632                        Page.canvas.SetToolTipString('%6d deg'%(azm))
1633                    elif 'line1' in  str(item) or 'line2' in str(item):
1634                        Page.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
1635            else:
1636                xpos = event.xdata
1637                ypos = event.ydata
1638                xpix = xpos*scalex
1639                ypix = ypos*scaley
1640                Int = 0
1641                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
1642                    Int = G2frame.ImageZ[ypix][xpix]
1643                tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
1644                Q = 2.*math.pi/dsp
1645                if G2frame.setPoly:
1646                    G2frame.G2plotNB.status.SetFields(['','Polygon mask pick - LB next point, RB close polygon'])
1647                else:
1648                    G2frame.G2plotNB.status.SetFields(\
1649                        ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)])
1650
1651    def OnImPlotKeyPress(event):
1652        try:
1653            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
1654        except TypeError:
1655            return
1656        if PickName == 'Masks':
1657            Xpos = event.xdata
1658            if not Xpos:            #got point out of frame
1659                return
1660            Ypos = event.ydata
1661            if event.key == 's':
1662                Masks['Points'].append([Xpos,Ypos,1.])
1663            elif event.key == 'r':
1664                tth = G2img.GetTth(Xpos,Ypos,Data)
1665                Masks['Rings'].append([tth,0.1])
1666            elif event.key == 'a':
1667                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
1668                azm = int(azm)               
1669                Masks['Arcs'].append([tth,[azm-5,azm+5],0.1])
1670            elif event.key == 'p':
1671                G2frame.setPoly = True
1672                Masks['Polygons'].append([])
1673                G2frame.G2plotNB.status.SetFields(['','Polygon mask active - LB pick next point, RB close polygon'])
1674            G2imG.UpdateMasks(G2frame,Masks)
1675        elif PickName == 'Image Controls':
1676            if event.key == 'c':
1677                Xpos = event.xdata
1678                if not Xpos:            #got point out of frame
1679                    return
1680                Ypos = event.ydata
1681                dlg = wx.MessageDialog(G2frame,'Are you sure you want to change the center?',
1682                    'Center change',style=wx.OK|wx.CANCEL)
1683                try:
1684                    if dlg.ShowModal() == wx.ID_OK:
1685                        print 'move center to: ',Xpos,Ypos
1686                        Data['center'] = [Xpos,Ypos]
1687                        G2imG.UpdateImageControls(G2frame,Data,Masks)
1688                finally:
1689                    dlg.Destroy()
1690            elif event.key == 'l':
1691                if G2frame.logPlot:
1692                    G2frame.logPlot = False
1693                else:
1694                    G2frame.logPlot = True
1695        PlotImage(G2frame,newImage=True)
1696           
1697    def OnKeyBox(event):
1698        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):
1699            event.key = cb.GetValue()[0]
1700            cb.SetValue(' key press')
1701            if event.key in 'l':
1702                OnImPlotKeyPress(event)
1703                       
1704    def OnImPick(event):
1705        if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']:
1706            return
1707        if G2frame.setPoly:
1708            polygon = Masks['Polygons'][-1]
1709            xpos,ypos = event.mouseevent.xdata,event.mouseevent.ydata
1710            if xpos and ypos:                       #point inside image
1711                if len(polygon) > 2 and event.mouseevent.button == 3:
1712                    x0,y0 = polygon[0]
1713                    polygon.append([x0,y0])
1714                    G2frame.setPoly = False
1715                    G2frame.G2plotNB.status.SetFields(['','Polygon closed - RB drag a vertex to change shape'])
1716                else:
1717                    G2frame.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(xpos,ypos)])
1718                    polygon.append([xpos,ypos])
1719                G2imG.UpdateMasks(G2frame,Masks)
1720        else:
1721            if G2frame.itemPicked is not None: return
1722            G2frame.itemPicked = event.artist
1723            G2frame.mousePicked = event.mouseevent
1724       
1725    def OnImRelease(event):
1726        try:
1727            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
1728        except TypeError:
1729            return
1730        if PickName not in ['Image Controls','Masks']:
1731            return
1732        pixelSize = Data['pixelSize']
1733        scalex = 1000./pixelSize[0]
1734        scaley = 1000./pixelSize[1]
1735        pixLimit = Data['pixLimit']
1736        if G2frame.itemPicked is None and PickName == 'Image Controls':
1737#            sizexy = Data['size']
1738            Xpos = event.xdata
1739            if not (Xpos and G2frame.ifGetRing):                   #got point out of frame
1740                return
1741            Ypos = event.ydata
1742            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
1743                if event.button == 1:
1744                    Xpix = Xpos*scalex
1745                    Ypix = Ypos*scaley
1746                    xpos,ypos,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
1747                    if I and J:
1748                        xpos += .5                              #shift to pixel center
1749                        ypos += .5
1750                        xpos /= scalex                          #convert to mm
1751                        ypos /= scaley
1752                        Data['ring'].append([xpos,ypos])
1753                elif event.button == 3:
1754                    G2frame.dataFrame.GetStatusBar().SetStatusText('Calibrating...')
1755                    if G2img.ImageCalibrate(G2frame,Data):
1756                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check')
1757                        print 'Calibration successful'
1758                    else:
1759                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose')
1760                        print 'Calibration failed'
1761                    G2frame.ifGetRing = False
1762                    G2imG.UpdateImageControls(G2frame,Data,Masks)
1763                    return
1764                PlotImage(G2frame,newImage=False)
1765            return
1766        else:
1767            xpos = event.xdata
1768            if xpos:                                        #avoid out of frame mouse position
1769                ypos = event.ydata
1770                if G2frame.ifGetRing:                          #delete a calibration ring pick
1771                    xypos = [xpos,ypos]
1772                    rings = Data['ring']
1773                    for ring in rings:
1774                        if np.allclose(ring,xypos,.01,0):
1775                            rings.remove(ring)                                                                       
1776                else:
1777                    tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
1778                    itemPicked = str(G2frame.itemPicked)
1779                    if 'Line2D' in itemPicked and PickName == 'Image Controls':
1780                        if 'line1' in itemPicked:
1781                            Data['IOtth'][0] = max(tth,0.001)
1782                        elif 'line2' in itemPicked:
1783                            Data['IOtth'][1] = tth
1784                        elif 'line3' in itemPicked:
1785                            Data['LRazimuth'][0] = int(azm)
1786                            if Data['fullIntegrate']:
1787                                Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
1788                        elif 'line4' in itemPicked and not Data['fullIntegrate']:
1789                            Data['LRazimuth'][1] = int(azm)
1790                           
1791                        if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
1792                            Data['LRazimuth'][0] -= 360
1793                           
1794                        azLim = np.array(Data['LRazimuth'])   
1795                        if np.any(azLim>360):
1796                            azLim -= 360
1797                            Data['LRazimuth'] = list(azLim)
1798                           
1799                        if  Data['IOtth'][0] > Data['IOtth'][1]:
1800                            Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
1801                           
1802                        G2frame.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
1803                        G2frame.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
1804                        G2frame.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
1805                        G2frame.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
1806                    elif 'Circle' in itemPicked and PickName == 'Masks':
1807                        spots = Masks['Points']
1808                        newPos = itemPicked.split(')')[0].split('(')[2].split(',')
1809                        newPos = np.array([float(newPos[0]),float(newPos[1])])
1810                        for spot in spots:
1811                            if np.allclose(np.array([spot[:2]]),newPos):
1812                                spot[:2] = xpos,ypos
1813                        G2imG.UpdateMasks(G2frame,Masks)
1814                    elif 'Line2D' in itemPicked and PickName == 'Masks':
1815                        Obj = G2frame.itemPicked.findobj()
1816                        rings = Masks['Rings']
1817                        arcs = Masks['Arcs']
1818                        polygons = Masks['Polygons']
1819                        for ring in G2frame.ringList:
1820                            if Obj == ring[0]:
1821                                rN = ring[1]
1822                                if ring[2] == 'o':
1823                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)-rings[rN][1]/2.
1824                                else:
1825                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)+rings[rN][1]/2.
1826                        for arc in G2frame.arcList:
1827                            if Obj == arc[0]:
1828                                aN = arc[1]
1829                                if arc[2] == 'o':
1830                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)-arcs[aN][2]/2
1831                                elif arc[2] == 'i':
1832                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)+arcs[aN][2]/2
1833                                elif arc[2] == 'l':
1834                                    arcs[aN][1][0] = int(G2img.GetAzm(xpos,ypos,Data))
1835                                else:
1836                                    arcs[aN][1][1] = int(G2img.GetAzm(xpos,ypos,Data))
1837                        for poly in G2frame.polyList:
1838                            if Obj == poly[0]:
1839                                ind = G2frame.itemPicked.contains(G2frame.mousePicked)[1]['ind'][0]
1840                                oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
1841                                pN = poly[1]
1842                                for i,xy in enumerate(polygons[pN]):
1843                                    if np.allclose(np.array([xy]),oldPos,atol=1.0):
1844                                        polygons[pN][i] = xpos,ypos
1845                        G2imG.UpdateMasks(G2frame,Masks)
1846#                    else:                  #keep for future debugging
1847#                        print str(G2frame.itemPicked),event.xdata,event.ydata,event.button
1848                PlotImage(G2frame,newImage=True)
1849            G2frame.itemPicked = None
1850           
1851    try:
1852        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
1853        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1854        if not newPlot:
1855            Plot = Page.figure.gca()          #get previous powder plot & get limits
1856            xylim = Plot.get_xlim(),Plot.get_ylim()
1857        if newImage:
1858            Page.figure.clf()
1859            Plot = Page.figure.gca()          #get a fresh plot after clf()
1860       
1861    except ValueError:
1862        Plot = G2frame.G2plotNB.addMpl('2D Powder Image').gca()
1863        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
1864        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1865        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
1866        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
1867        Page.canvas.mpl_connect('pick_event', OnImPick)
1868        Page.canvas.mpl_connect('button_release_event', OnImRelease)
1869        xylim = []
1870    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
1871        Page.SetFocus()
1872    Title = G2frame.PatternTree.GetItemText(G2frame.Image)[4:]
1873    G2frame.G2plotNB.status.DestroyChildren()
1874    if G2frame.logPlot:
1875        Title = 'log('+Title+')'
1876    Plot.set_title(Title)
1877    try:
1878        if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Image Controls',]:
1879            Choice = (' key press','l: log(I) on',)
1880            if G2frame.logPlot:
1881                Choice[1] = 'l: log(I) off'
1882            cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
1883                choices=Choice)
1884            cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
1885            cb.SetValue(' key press')
1886    except TypeError:
1887        pass
1888    size,imagefile = G2frame.PatternTree.GetItemPyData(G2frame.Image)
1889    if imagefile != G2frame.oldImagefile:
1890        imagefile = G2IO.CheckImageFile(G2frame,imagefile)
1891        if not imagefile:
1892            G2frame.G2plotNB.Delete('2D Powder Image')
1893            return
1894        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
1895        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
1896#        print G2frame.ImageZ.shape,G2frame.ImageZ.size,Data['size'] #might be useful debugging line
1897        G2frame.oldImagefile = imagefile
1898
1899    imScale = 1
1900    if len(G2frame.ImageZ) > 1024:
1901        imScale = len(G2frame.ImageZ)/1024
1902    sizexy = Data['size']
1903    pixelSize = Data['pixelSize']
1904    scalex = 1000./pixelSize[0]
1905    scaley = 1000./pixelSize[1]
1906    Xmax = sizexy[0]*pixelSize[0]/1000.
1907    Ymax = sizexy[1]*pixelSize[1]/1000.
1908    xlim = (0,Xmax)
1909    ylim = (Ymax,0)
1910    Imin,Imax = Data['range'][1]
1911    acolor = mpl.cm.get_cmap(Data['color'])
1912    xcent,ycent = Data['center']
1913    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
1914    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
1915    #do threshold mask - "real" mask - others are just bondaries
1916    Zlim = Masks['Thresholds'][1]
1917    wx.BeginBusyCursor()
1918    try:
1919           
1920        if newImage:                   
1921            MA = ma.masked_greater(ma.masked_less(G2frame.ImageZ,Zlim[0]),Zlim[1])
1922            MaskA = ma.getmaskarray(MA)
1923            A = G2img.ImageCompress(MA,imScale)
1924            AM = G2img.ImageCompress(MaskA,imScale)
1925            if G2frame.logPlot:
1926                A = np.where(A>0,np.log(A),0)
1927                AM = np.where(AM>0,np.log(AM),0)
1928                Imin,Imax = [np.amin(A),np.amax(A)]
1929            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
1930                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Ymax,0])
1931            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
1932                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
1933            if G2frame.setPoly:
1934                Img.set_picker(True)
1935   
1936        Plot.plot(xcent,ycent,'x')
1937        if Data['showLines']:
1938            LRAzim = Data['LRazimuth']                  #NB: integers
1939            Nazm = Data['outAzimuths']
1940            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
1941            AzmthOff = Data['azmthOff']
1942            IOtth = Data['IOtth']
1943            wave = Data['wavelength']
1944            dspI = wave/(2.0*sind(IOtth[0]/2.0))
1945            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
1946            dspO = wave/(2.0*sind(IOtth[1]/2.0))
1947            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
1948            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
1949            if ellI:
1950                xyI = []
1951                for azm in Azm:
1952                    xyI.append(G2img.GetDetectorXY(dspI,azm-90.,Data))
1953                xyI = np.array(xyI)
1954                arcxI,arcyI = xyI.T
1955                Plot.plot(arcxI,arcyI,picker=3)
1956            if ellO:
1957                xyO = []
1958                for azm in Azm:
1959                    xyO.append(G2img.GetDetectorXY(dspO,azm-90.,Data))
1960                xyO = np.array(xyO)
1961                arcxO,arcyO = xyO.T
1962                Plot.plot(arcxO,arcyO,picker=3)
1963            if ellO and ellI:
1964                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
1965                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
1966            for i in range(Nazm):
1967                cake = LRAzim[0]+i*delAzm-AzmthOff
1968                ind = np.searchsorted(Azm,cake)
1969                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
1970                   
1971        for xring,yring in Data['ring']:
1972            Plot.plot(xring,yring,'r+',picker=3)
1973        if Data['setRings']:
1974#            rings = np.concatenate((Data['rings']),axis=0)
1975            N = 0
1976            for ring in Data['rings']:
1977                xring,yring = np.array(ring).T[:2]
1978                Plot.plot(xring,yring,'+',color=colors[N%6])
1979                N += 1           
1980        for ellipse in Data['ellipses']:
1981            cent,phi,[width,height],col = ellipse
1982            Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
1983            Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
1984        #masks - mask lines numbered after integration limit lines
1985        spots = Masks['Points']
1986        rings = Masks['Rings']
1987        arcs = Masks['Arcs']
1988        polygons = Masks['Polygons']
1989        for x,y,d in spots:
1990            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
1991        G2frame.ringList = []
1992        for iring,(tth,thick) in enumerate(rings):
1993            wave = Data['wavelength']
1994            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
1995            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
1996            G2frame.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
1997            G2frame.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
1998        G2frame.arcList = []
1999        for iarc,(tth,azm,thick) in enumerate(arcs):           
2000            wave = Data['wavelength']
2001            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
2002            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
2003            G2frame.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
2004            G2frame.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
2005            G2frame.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
2006            G2frame.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
2007        G2frame.polyList = []
2008        for ipoly,polygon in enumerate(polygons):
2009            x,y = np.hsplit(np.array(polygon),2)
2010            G2frame.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
2011            Plot.plot(x,y,'r')           
2012        if newImage:
2013            colorBar = Page.figure.colorbar(Img)
2014        Plot.set_xlim(xlim)
2015        Plot.set_ylim(ylim)
2016        Plot.invert_yaxis()
2017        if not newPlot and xylim:
2018            Page.toolbar.push_current()
2019            Plot.set_xlim(xylim[0])
2020            Plot.set_ylim(xylim[1])
2021            xylim = []
2022            Page.toolbar.push_current()
2023            Page.toolbar.draw()
2024        else:
2025            Page.canvas.draw()
2026    finally:
2027        wx.EndBusyCursor()
2028       
2029def PlotIntegration(G2frame,newPlot=False,event=None):
2030    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
2031    '''
2032           
2033    def OnMotion(event):
2034        Page.canvas.SetToolTipString('')
2035        Page.canvas.SetCursor(wx.CROSS_CURSOR)
2036        azm = event.ydata
2037        tth = event.xdata
2038        if azm and tth:
2039            G2frame.G2plotNB.status.SetFields(\
2040                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
2041                               
2042    try:
2043        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
2044        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2045        if not newPlot:
2046            Plot = Page.figure.gca()          #get previous plot & get limits
2047            xylim = Plot.get_xlim(),Plot.get_ylim()
2048        Page.figure.clf()
2049        Plot = Page.figure.gca()          #get a fresh plot after clf()
2050       
2051    except ValueError:
2052        Plot = G2frame.G2plotNB.addMpl('2D Integration').gca()
2053        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
2054        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2055        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2056        Page.views = False
2057        view = False
2058    if not event:
2059        Page.SetFocus()
2060       
2061    Data = G2frame.PatternTree.GetItemPyData(
2062        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2063    image = G2frame.Integrate[0]
2064    xsc = G2frame.Integrate[1]
2065    ysc = G2frame.Integrate[2]
2066    Imin,Imax = Data['range'][1]
2067    acolor = mpl.cm.get_cmap(Data['color'])
2068    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
2069    Plot.set_ylabel('azimuth',fontsize=12)
2070    Plot.set_xlabel('2-theta',fontsize=12)
2071    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
2072        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
2073    colorBar = Page.figure.colorbar(Img)
2074    if Data['ellipses']:           
2075        for ellipse in Data['ellipses']:
2076            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
2077            x,y = np.hsplit(ring,2)
2078            tth,azm = G2img.GetTthAzm(x,y,Data)
2079#            azm = np.where(azm < 0.,azm+360,azm)
2080            Plot.plot(tth,azm,'b,')
2081    if not newPlot:
2082        Page.toolbar.push_current()
2083        Plot.set_xlim(xylim[0])
2084        Plot.set_ylim(xylim[1])
2085        xylim = []
2086        Page.toolbar.push_current()
2087        Page.toolbar.draw()
2088    else:
2089        Page.canvas.draw()
2090               
2091def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):
2092    '''a test plot routine - not normally used
2093    ''' 
2094           
2095    def OnMotion(event):
2096        Page.canvas.SetToolTipString('')
2097        Page.canvas.SetCursor(wx.CROSS_CURSOR)
2098        azm = event.xdata
2099        tth = event.ydata
2100        if azm and tth:
2101            G2frame.G2plotNB.status.SetFields(\
2102                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
2103                               
2104    try:
2105        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
2106        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2107        if not newPlot:
2108            Plot = Page.figure.gca()          #get previous plot & get limits
2109            xylim = Plot.get_xlim(),Plot.get_ylim()
2110        Page.figure.clf()
2111        Plot = Page.figure.gca()          #get a fresh plot after clf()
2112       
2113    except ValueError:
2114        Plot = G2frame.G2plotNB.addMpl('2D Transformed Powder Image').gca()
2115        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
2116        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2117        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2118        Page.views = False
2119        view = False
2120    Page.SetFocus()
2121       
2122    Data = G2frame.PatternTree.GetItemPyData(
2123        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2124    Imin,Imax = Data['range'][1]
2125    step = (Imax-Imin)/5.
2126    V = np.arange(Imin,Imax,step)
2127    acolor = mpl.cm.get_cmap(Data['color'])
2128    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
2129    Plot.set_xlabel('azimuth',fontsize=12)
2130    Plot.set_ylabel('2-theta',fontsize=12)
2131    Plot.contour(tax,tay,taz,V,cmap=acolor)
2132    if Data['showLines']:
2133        IOtth = Data['IOtth']
2134        if Data['fullIntegrate']:
2135            LRAzim = [-180,180]
2136        else:
2137            LRAzim = Data['LRazimuth']                  #NB: integers
2138        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
2139        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
2140        if not Data['fullIntegrate']:
2141            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
2142            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
2143    if Data['setRings']:
2144        rings = np.concatenate((Data['rings']),axis=0)
2145        for xring,yring,dsp in rings:
2146            x,y = G2img.GetTthAzm(xring,yring,Data)
2147            Plot.plot(y,x,'r+')           
2148    if Data['ellipses']:           
2149        for ellipse in Data['ellipses']:
2150            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
2151            x,y = np.hsplit(ring,2)
2152            tth,azm = G2img.GetTthAzm(x,y,Data)
2153            Plot.plot(azm,tth,'b,')
2154    if not newPlot:
2155        Page.toolbar.push_current()
2156        Plot.set_xlim(xylim[0])
2157        Plot.set_ylim(xylim[1])
2158        xylim = []
2159        Page.toolbar.push_current()
2160        Page.toolbar.draw()
2161    else:
2162        Page.canvas.draw()
2163       
2164def PlotStructure(G2frame,data):
2165    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
2166    thermal motion ellipsoids and polyhedra
2167    '''
2168    generalData = data['General']
2169    cell = generalData['Cell'][1:7]
2170    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
2171    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2172    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2173    Mydir = generalData['Mydir']
2174    atomData = data['Atoms']
2175    drawingData = data['Drawing']
2176    drawAtoms = drawingData['Atoms']
2177    cx,ct,cs = drawingData['atomPtrs']
2178    Wt = [255,255,255]
2179    Rd = [255,0,0]
2180    Gr = [0,255,0]
2181    Bl = [0,0,255]
2182    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]])
2183    uEdges = np.array([
2184        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
2185        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
2186        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
2187    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
2188    altDown = False
2189    shiftDown = False
2190    ctrlDown = False
2191   
2192    def OnKeyBox(event):
2193        import Image
2194        Draw()                          #make sure plot is fresh!!
2195        mode = cb.GetValue()
2196        Fname = Mydir+'\\'+generalData['Name']+'.'+mode
2197        size = Page.canvas.GetSize()
2198        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
2199        if mode in ['jpeg',]:
2200            Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
2201            im = Image.new("RGBA", (size[0],size[1]))
2202        else:
2203            Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
2204            im = Image.new("RGB", (size[0],size[1]))
2205        im.fromstring(Pix)
2206        im.save(Fname,mode)
2207        cb.SetValue(' Save as:')
2208        G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
2209   
2210    def GetTruePosition(xy,Add=False):
2211        View = glGetIntegerv(GL_VIEWPORT)
2212        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
2213        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
2214        Zmax = 1.
2215        if Add:
2216            Indx = GetSelectedAtoms()
2217        for i,atom in enumerate(drawAtoms):
2218            x,y,z = atom[cx:cx+3]
2219            X,Y,Z = gluProject(x,y,z,Model,Proj,View)
2220            XY = [int(X),int(View[3]-Y)]
2221            if np.allclose(xy,XY,atol=10) and Z < Zmax:
2222                Zmax = Z
2223                try:
2224                    Indx.remove(i)
2225                    ClearSelectedAtoms()
2226                    for id in Indx:
2227                        SetSelectedAtoms(id,Add)
2228                except:
2229                    SetSelectedAtoms(i,Add)
2230                   
2231    def OnMouseDown(event):
2232        xy = event.GetPosition()
2233        if event.ShiftDown():
2234            if event.LeftIsDown():
2235                GetTruePosition(xy)
2236            elif event.RightIsDown():
2237                GetTruePosition(xy,True)
2238        else:
2239            drawingData['Rotation'][3] = list(xy)
2240        Draw()
2241       
2242    def OnMouseMove(event):
2243        if event.ShiftDown():
2244            return       
2245        newxy = event.GetPosition()
2246        page = getSelection()
2247        if event.ControlDown() and drawingData['showABC']:
2248            if event.LeftIsDown():
2249                SetTestRot(newxy)
2250            elif event.RightIsDown():
2251                SetTestPos(newxy)
2252            elif event.MiddleIsDown():
2253                SetTestRotZ(newxy)
2254            x,y,z = drawingData['testPos'][0]
2255            G2frame.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
2256            Draw()
2257                               
2258        if event.Dragging() and not event.ControlDown():
2259            if event.LeftIsDown():
2260                SetRotation(newxy)
2261                angX,angY,angZ = drawingData['Rotation'][:3]
2262                G2frame.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f ,%.2f'%(angX,angY,angZ),1)
2263            elif event.RightIsDown():
2264                SetTranslation(newxy)
2265                Tx,Ty,Tz = drawingData['viewPoint'][0]
2266                G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
2267            elif event.MiddleIsDown():
2268                SetRotationZ(newxy)
2269                angX,angY,angZ = drawingData['Rotation'][:3]
2270                G2frame.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f, %.2f'%(angX,angY,angZ),1)
2271            Draw()
2272       
2273    def OnMouseWheel(event):
2274        if event.ShiftDown():
2275            return
2276        drawingData['cameraPos'] += event.GetWheelRotation()/24
2277        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
2278        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(drawingData['cameraPos']),1)
2279        page = getSelection()
2280        if page:
2281            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2282                panel = G2frame.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
2283                names = [child.GetName() for child in panel]
2284                panel[names.index('cameraPos')].SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
2285                panel[names.index('cameraSlider')].SetValue(drawingData['cameraPos'])
2286            Draw()
2287       
2288    def getSelection():
2289        try:
2290            return G2frame.dataDisplay.GetSelection()
2291        except AttributeError:
2292            G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!')
2293            return 0
2294           
2295    def SetViewPointText(VP):
2296        page = getSelection()
2297        if page:
2298            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2299                panel = G2frame.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
2300                names = [child.GetName() for child in panel]
2301                panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
2302           
2303    def ClearSelectedAtoms():
2304        page = getSelection()
2305        if page:
2306            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2307                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
2308            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2309                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
2310                   
2311    def SetSelectedAtoms(ind,Add=False):
2312        page = getSelection()
2313        if page:
2314            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2315                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)      #this is the Atoms grid in Draw Atoms
2316            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2317                Id = drawAtoms[ind][-2]
2318                for i,atom in enumerate(atomData):
2319                    if atom[-1] == Id:
2320                        G2frame.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
2321                 
2322    def GetSelectedAtoms():
2323        page = getSelection()
2324        Ind = []
2325        if page:
2326            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2327                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
2328            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2329                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
2330        return Ind
2331                                       
2332    def OnKey(event):           #on key UP!!
2333        keyCode = event.GetKeyCode()
2334        if keyCode > 255:
2335            keyCode = 0
2336        key,xyz = chr(keyCode),event.GetPosition()
2337        indx = drawingData['selectedAtoms']
2338        if key in ['c','C']:
2339            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
2340            drawingData['testPos'] = [[-.1,-.1,-.1],[0.0,0.0,0.0],[0,0]]
2341            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
2342            SetViewPointText(drawingData['viewPoint'][0])
2343        elif key in ['n','N']:
2344            drawAtoms = drawingData['Atoms']
2345            pI = drawingData['viewPoint'][1]
2346            if indx:
2347                pI[0] = indx[pI[1]]
2348                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
2349                pI[1] += 1
2350                if pI[1] >= len(indx):
2351                    pI[1] = 0
2352            else:
2353                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
2354                pI[0] += 1
2355                if pI[0] >= len(drawAtoms):
2356                    pI[0] = 0
2357            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
2358            SetViewPointText(drawingData['viewPoint'][0])
2359            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
2360               
2361        elif key in ['p','P']:
2362            drawAtoms = drawingData['Atoms']
2363            pI = drawingData['viewPoint'][1]
2364            if indx:
2365                pI[0] = indx[pI[1]]
2366                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
2367                pI[1] -= 1
2368                if pI[1] < 0:
2369                    pI[1] = len(indx)-1
2370            else:
2371                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
2372                pI[0] -= 1
2373                if pI[0] < 0:
2374                    pI[0] = len(drawAtoms)-1
2375            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
2376            SetViewPointText(drawingData['viewPoint'][0])           
2377            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
2378        Draw()
2379           
2380    def SetBackground():
2381        R,G,B,A = Page.camera['backColor']
2382        glClearColor(R,G,B,A)
2383        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
2384       
2385    def SetLights():
2386        glEnable(GL_DEPTH_TEST)
2387        glShadeModel(GL_SMOOTH)
2388        glEnable(GL_LIGHTING)
2389        glEnable(GL_LIGHT0)
2390        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
2391        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
2392        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
2393       
2394    def SetTranslation(newxy):
2395        Tx,Ty,Tz = drawingData['viewPoint'][0]
2396        anglex,angley,anglez,oldxy = drawingData['Rotation']
2397        Rx = G2lat.rotdMat(anglex,0)
2398        Ry = G2lat.rotdMat(angley,1)
2399        Rz = G2lat.rotdMat(anglez,2)
2400        dxy = list(newxy-oldxy)+[0,]
2401        dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
2402        Tx -= dxy[0]*0.01
2403        Ty += dxy[1]*0.01
2404        Tz -= dxy[2]*0.01
2405        drawingData['Rotation'][3] = list(newxy)
2406        drawingData['viewPoint'][0] =  Tx,Ty,Tz
2407        SetViewPointText([Tx,Ty,Tz])
2408       
2409    def SetTestPos(newxy):
2410        Tx,Ty,Tz = drawingData['testPos'][0]
2411        anglex,angley,anglez,oldxy = drawingData['Rotation']
2412        Rx = G2lat.rotdMat(anglex,0)
2413        Ry = G2lat.rotdMat(angley,1)
2414        Rz = G2lat.rotdMat(anglez,2)
2415        dxy = list(newxy-oldxy)+[0,]
2416        dxy = np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy)))
2417        Tx += dxy[0]*0.001
2418        Ty -= dxy[1]*0.001
2419        Tz += dxy[2]*0.001
2420        drawingData['Rotation'][3] = list(newxy)
2421        drawingData['testPos'][0] =  Tx,Ty,Tz
2422       
2423    def SetTestRot(newxy):
2424        Txyz = np.array(drawingData['testPos'][0])
2425        oldxy = drawingData['testPos'][2]
2426        Ax,Ay,Az = drawingData['testPos'][1]
2427        Vxyz = np.array(drawingData['viewPoint'][0])
2428        Dxyz = np.inner(Amat,Txyz-Vxyz)
2429        dxy = list(newxy-oldxy)+[0,]
2430        Ax += dxy[1]*0.01
2431        Ay += dxy[0]*0.01
2432        Rx = G2lat.rotdMat(Ax,0)
2433        Ry = G2lat.rotdMat(Ay,1)
2434        Dxyz = np.inner(Ry,np.inner(Rx,Dxyz))       
2435        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
2436        drawingData['testPos'][1] = [Ax,Ay,Az]
2437        drawingData['testPos'][2] = newxy
2438        drawingData['testPos'][0] = Dxyz
2439       
2440    def SetTestRotZ(newxy):
2441        Txyz = np.array(drawingData['testPos'][0])
2442        oldxy = drawingData['testPos'][2]
2443        Ax,Ay,Az = drawingData['testPos'][1]
2444        Vxyz = np.array(drawingData['viewPoint'][0])
2445        Dxyz = np.inner(Amat,Txyz-Vxyz)       
2446        dxy = list(newxy-oldxy)+[0,]
2447        Az += (dxy[0]+dxy[1])*.01
2448        Rz = G2lat.rotdMat(Az,2)
2449        Dxyz = np.inner(Rz,Dxyz)       
2450        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
2451        drawingData['testPos'][1] = [Ax,Ay,Az]
2452        drawingData['testPos'][2] = newxy
2453        drawingData['testPos'][0] = Dxyz
2454                             
2455    def SetRotation(newxy):       
2456        anglex,angley,anglez,oldxy = drawingData['Rotation']
2457        dxy = newxy-oldxy
2458        anglex += dxy[1]*.25
2459        angley += dxy[0]*.25
2460        oldxy = newxy
2461        drawingData['Rotation'] = [anglex,angley,anglez,list(oldxy)]
2462       
2463    def SetRotationZ(newxy):                       
2464        anglex,angley,anglez,oldxy = drawingData['Rotation']
2465        dxy = newxy-oldxy
2466        anglez += (dxy[0]+dxy[1])*.25
2467        oldxy = newxy
2468        drawingData['Rotation'] = [anglex,angley,anglez,list(oldxy)]
2469       
2470    def RenderBox():
2471        glEnable(GL_COLOR_MATERIAL)
2472        glLineWidth(2)
2473        glEnable(GL_BLEND)
2474        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
2475        glEnable(GL_LINE_SMOOTH)
2476        glBegin(GL_LINES)
2477        for line,color in zip(uEdges,uColors):
2478            glColor3ubv(color)
2479            glVertex3fv(line[0])
2480            glVertex3fv(line[1])
2481        glEnd()
2482        glColor4ubv([0,0,0,0])
2483        glDisable(GL_LINE_SMOOTH)
2484        glDisable(GL_BLEND)
2485        glDisable(GL_COLOR_MATERIAL)
2486       
2487    def RenderUnitVectors(x,y,z):
2488        xyz = np.array([x,y,z])
2489        glEnable(GL_COLOR_MATERIAL)
2490        glLineWidth(1)
2491        glPushMatrix()
2492        glTranslate(x,y,z)
2493        glScalef(1/cell[0],1/cell[1],1/cell[2])
2494        glBegin(GL_LINES)
2495        for line,color in zip(uEdges,uColors)[:3]:
2496            glColor3ubv(color)
2497            glVertex3fv(line[0])
2498            glVertex3fv(line[1])
2499        glEnd()
2500        glPopMatrix()
2501        glColor4ubv([0,0,0,0])
2502        glDisable(GL_COLOR_MATERIAL)
2503               
2504    def RenderSphere(x,y,z,radius,color):
2505        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2506        glPushMatrix()
2507        glTranslate(x,y,z)
2508        glMultMatrixf(B4mat.T)
2509        q = gluNewQuadric()
2510        gluSphere(q,radius,20,10)
2511        glPopMatrix()
2512       
2513    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
2514        s1,s2,s3 = E
2515        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2516        glPushMatrix()
2517        glTranslate(x,y,z)
2518        glMultMatrixf(B4mat.T)
2519        glMultMatrixf(R4.T)
2520        glEnable(GL_NORMALIZE)
2521        glScale(s1,s2,s3)
2522        q = gluNewQuadric()
2523        gluSphere(q,ellipseProb,20,10)
2524        glDisable(GL_NORMALIZE)
2525        glPopMatrix()
2526       
2527    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
2528        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2529        glPushMatrix()
2530        glTranslate(x,y,z)
2531        glMultMatrixf(B4mat.T)
2532        for bond in Bonds:
2533            glPushMatrix()
2534            Dx = np.inner(Amat,bond)
2535            Z = np.sqrt(np.sum(Dx**2))
2536            azm = atan2d(-Dx[1],-Dx[0])
2537            phi = acosd(Dx[2]/Z)
2538            glRotate(-azm,0,0,1)
2539            glRotate(phi,1,0,0)
2540            q = gluNewQuadric()
2541            gluCylinder(q,radius,radius,Z,slice,2)
2542            glPopMatrix()           
2543        glPopMatrix()
2544               
2545    def RenderLines(x,y,z,Bonds,color):
2546        xyz = np.array([x,y,z])
2547        glEnable(GL_COLOR_MATERIAL)
2548        glLineWidth(1)
2549        glColor3fv(color)
2550        glPushMatrix()
2551        glBegin(GL_LINES)
2552        for bond in Bonds:
2553            glVertex3fv(xyz)
2554            glVertex3fv(xyz+bond)
2555        glEnd()
2556        glColor4ubv([0,0,0,0])
2557        glPopMatrix()
2558        glDisable(GL_COLOR_MATERIAL)
2559       
2560    def RenderPolyhedra(x,y,z,Faces,color):
2561        glPushMatrix()
2562        glTranslate(x,y,z)
2563        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2564        glShadeModel(GL_SMOOTH)
2565        glMultMatrixf(B4mat.T)
2566        for face,norm in Faces:
2567            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
2568            glFrontFace(GL_CW)
2569            glNormal3fv(norm)
2570            glBegin(GL_TRIANGLES)
2571            for vert in face:
2572                glVertex3fv(vert)
2573            glEnd()
2574        glPopMatrix()
2575       
2576    def RenderBackbone(Backbone,BackboneColor,radius):
2577        glPushMatrix()
2578        glMultMatrixf(B4mat.T)
2579        glEnable(GL_COLOR_MATERIAL)
2580        glShadeModel(GL_SMOOTH)
2581        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
2582        glePolyCylinder(Backbone,BackboneColor,radius)
2583        glPopMatrix()       
2584        glDisable(GL_COLOR_MATERIAL)
2585       
2586    def RenderLabel(x,y,z,label,r):       
2587        glPushMatrix()
2588        glTranslate(x,y,z)
2589        glMultMatrixf(B4mat.T)
2590        glDisable(GL_LIGHTING)
2591        glColor3f(0,1.,0)
2592        glRasterPos3f(r,r,r)
2593        for c in list(label):
2594            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
2595        glEnable(GL_LIGHTING)
2596        glPopMatrix()
2597                           
2598    def Draw():
2599        Ind = GetSelectedAtoms()
2600        VS = np.array(Page.canvas.GetSize())
2601        aspect = float(VS[0])/float(VS[1])
2602        cPos = drawingData['cameraPos']
2603        Zclip = drawingData['Zclip']*cPos/200.
2604        anglex,angley,anglez = drawingData['Rotation'][:3]
2605        Tx,Ty,Tz = drawingData['viewPoint'][0]
2606        cx,ct,cs = drawingData['atomPtrs']
2607        bondR = drawingData['bondRadius']
2608        G,g = G2lat.cell2Gmat(cell)
2609        GS = G
2610        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
2611        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
2612        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
2613        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
2614       
2615        SetBackground()
2616        glInitNames()
2617        glPushName(0)
2618       
2619        glMatrixMode(GL_PROJECTION)
2620        glLoadIdentity()
2621        glViewport(0,0,VS[0],VS[1])
2622        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
2623        gluLookAt(0,0,cPos,0,0,0,0,1,0)
2624        SetLights()           
2625           
2626        glMatrixMode(GL_MODELVIEW)
2627        glLoadIdentity()
2628        glRotate(anglez,0,0,1)
2629        glRotate(anglex,cosd(anglez),-sind(anglez),0)
2630        glRotate(angley,sind(anglez),cosd(anglez),0)
2631        glMultMatrixf(A4mat.T)
2632        glTranslate(-Tx,-Ty,-Tz)
2633        if drawingData['unitCellBox']:
2634            RenderBox()
2635        if drawingData['showABC']:
2636            x,y,z = drawingData['testPos'][0]
2637#            if altDown:
2638#                G2frame.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
2639#            else:
2640#                G2frame.G2plotNB.status.SetStatusText('test point %.4f,%.4f,%.4f'%(x,y,z),1)           
2641            RenderUnitVectors(x,y,z)
2642        Backbone = []
2643        BackboneColor = []
2644        time0 = time.time()
2645        for iat,atom in enumerate(drawingData['Atoms']):
2646            x,y,z = atom[cx:cx+3]
2647            Bonds = atom[-2]
2648            Faces = atom[-1]
2649            try:
2650                atNum = generalData['AtomTypes'].index(atom[ct])
2651            except ValueError:
2652                atNum = -1
2653            CL = atom[cs+2]
2654            color = np.array(CL)/255.
2655            if iat in Ind:
2656                color = np.array(Gr)/255.
2657            radius = 0.5
2658            if atom[cs] != '':
2659                glLoadName(atom[-3])                   
2660            if 'balls' in atom[cs]:
2661                vdwScale = drawingData['vdwScale']
2662                ballScale = drawingData['ballScale']
2663                if atNum < 0:
2664                    radius = 0.2
2665                elif 'H' == atom[ct]:
2666                    if drawingData['showHydrogen']:
2667                        if 'vdW' in atom[cs] and atNum >= 0:
2668                            radius = vdwScale*generalData['vdWRadii'][atNum]
2669                        else:
2670                            radius = ballScale*drawingData['sizeH']
2671                    else:
2672                        radius = 0.0
2673                else:
2674                    if 'vdW' in atom[cs]:
2675                        radius = vdwScale*generalData['vdWRadii'][atNum]
2676                    else:
2677                        radius = ballScale*generalData['BondRadii'][atNum]
2678                RenderSphere(x,y,z,radius,color)
2679                if 'sticks' in atom[cs]:
2680                    RenderBonds(x,y,z,Bonds,bondR,color)
2681            elif 'ellipsoids' in atom[cs]:
2682                RenderBonds(x,y,z,Bonds,bondR,color)
2683                if atom[cs+3] == 'A':                   
2684                    Uij = atom[cs+5:cs+11]
2685                    U = np.multiply(G2spc.Uij2U(Uij),GS)
2686                    U = np.inner(Amat,np.inner(U,Amat).T)
2687                    E,R = nl.eigh(U)
2688                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2689                    E = np.sqrt(E)
2690                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
2691                        pass
2692                    else:
2693                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,color)                   
2694                else:
2695                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
2696                        pass
2697                    else:
2698                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
2699                        RenderSphere(x,y,z,radius,color)
2700            elif 'lines' in atom[cs]:
2701                radius = 0.1
2702                RenderLines(x,y,z,Bonds,color)
2703#                RenderBonds(x,y,z,Bonds,0.05,color,6)
2704            elif atom[cs] == 'sticks':
2705                radius = 0.1
2706                RenderBonds(x,y,z,Bonds,bondR,color)
2707            elif atom[cs] == 'polyhedra':
2708                RenderPolyhedra(x,y,z,Faces,color)
2709            elif atom[cs] == 'backbone':
2710                if atom[ct-1].split()[0] in ['C','N']:
2711                    Backbone.append(list(np.inner(Amat,np.array([x,y,z]))))
2712                    BackboneColor.append(list(color))
2713                   
2714            if atom[cs+1] == 'type':
2715                RenderLabel(x,y,z,atom[ct],radius)
2716            elif atom[cs+1] == 'name':
2717                RenderLabel(x,y,z,atom[ct-1],radius)
2718            elif atom[cs+1] == 'number':
2719                RenderLabel(x,y,z,str(iat),radius)
2720            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
2721                RenderLabel(x,y,z,atom[ct-4],radius)
2722            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
2723                RenderLabel(x,y,z,atom[ct-3],radius)
2724            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
2725                RenderLabel(x,y,z,atom[ct-2],radius)
2726        if Backbone:
2727            RenderBackbone(Backbone,BackboneColor,bondR)
2728#        print time.time()-time0
2729        Page.canvas.SwapBuffers()
2730       
2731    def OnSize(event):
2732        Draw()
2733       
2734    def OnFocus(event):
2735        Draw()
2736       
2737    try:
2738        plotNum = G2frame.G2plotNB.plotList.index(generalData['Name'])
2739        Page = G2frame.G2plotNB.nb.GetPage(plotNum)       
2740    except ValueError:
2741        Plot = G2frame.G2plotNB.addOgl(generalData['Name'])
2742        plotNum = G2frame.G2plotNB.plotList.index(generalData['Name'])
2743        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2744        Page.views = False
2745        view = False
2746        altDown = False
2747    Page.SetFocus()
2748    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
2749        choices=(' save as:','jpeg','tiff','bmp'))
2750    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
2751    cb.SetValue(' save as:')
2752    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
2753    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
2754    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
2755    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
2756    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
2757    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
2758    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
2759    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
2760    Page.camera['position'] = drawingData['cameraPos']
2761    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
2762    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
2763    Page.canvas.SetCurrent()
2764    Draw()
2765       
Note: See TracBrowser for help on using the repository browser.