source: branch/MPbranch/GSASIIplot.py

Last change on this file was 517, checked in by toby, 10 years ago

make sure that texture is initialized for multiprocessing

  • Property svn:keywords set to Date Author Revision URL Id
File size: 117.8 KB
Line 
1#GSASII plotting routines
2########### SVN repository information ###################
3# $Date: 2012-03-09 21:10:26 +0000 (Fri, 09 Mar 2012) $
4# $Author: vondreele $
5# $Revision: 517 $
6# $URL: branch/MPbranch/GSASIIplot.py $
7# $Id: GSASIIplot.py 517 2012-03-09 21:10:26Z 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 'plot weighting:',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','+: 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                    W2 = np.sqrt(xye[2])
671                    D *= W2-Ymax*G2frame.delOffset
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 PlotISFG(G2frame,newPlot=False,type=''):
752    ''' PLotting package for PDF analysis; displays I(q), S(q), F(q) and G(r) as single
753    or multiple plots with waterfall and contour plots as options
754    '''
755    if not type:
756        type = G2frame.G2plotNB.plotList[G2frame.G2plotNB.nb.GetSelection()]
757    if type not in ['I(Q)','S(Q)','F(Q)','G(R)']:
758        return
759    superMinusOne = unichr(0xaf)+unichr(0xb9)
760   
761    def OnPlotKeyPress(event):
762        newPlot = False
763        if event.key == 'u':
764            if G2frame.Contour:
765                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
766            elif G2frame.Offset[0] < 100.:
767                G2frame.Offset[0] += 1.
768        elif event.key == 'd':
769            if G2frame.Contour:
770                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
771            elif G2frame.Offset[0] > 0.:
772                G2frame.Offset[0] -= 1.
773        elif event.key == 'l':
774            G2frame.Offset[1] -= 1.
775        elif event.key == 'r':
776            G2frame.Offset[1] += 1.
777        elif event.key == 'o':
778            G2frame.Offset = [0,0]
779        elif event.key == 'c':
780            newPlot = True
781            if G2frame.Contour:
782                G2frame.Contour = False
783            else:
784                G2frame.Contour = True
785                G2frame.SinglePlot = False
786                G2frame.Offset = [0.,0.]
787        elif event.key == 's':
788            if G2frame.Contour:
789                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
790                choice.sort()
791                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
792                if dlg.ShowModal() == wx.ID_OK:
793                    sel = dlg.GetSelection()
794                    G2frame.ContourColor = choice[sel]
795                else:
796                    G2frame.ContourColor = 'Paired'
797                dlg.Destroy()
798            else:               
799                if G2frame.SinglePlot:
800                    G2frame.SinglePlot = False
801                else:
802                    G2frame.SinglePlot = True
803        elif event.key == 'i':                  #for smoothing contour plot
804            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
805               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
806               'mitchell','sinc','lanczos']
807            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
808            if dlg.ShowModal() == wx.ID_OK:
809                sel = dlg.GetSelection()
810                G2frame.Interpolate = choice[sel]
811            else:
812                G2frame.Interpolate = 'nearest'
813            dlg.Destroy()
814        elif event.key == 't' and not G2frame.Contour:
815            if G2frame.Legend:
816                G2frame.Legend = False
817            else:
818                G2frame.Legend = True
819           
820           
821        PlotISFG(G2frame,newPlot=newPlot,type=type)
822       
823    def OnKeyBox(event):
824        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index(type):
825            event.key = cb.GetValue()[0]
826            cb.SetValue(' key press')
827            OnPlotKeyPress(event)
828                       
829    def OnMotion(event):
830        xpos = event.xdata
831        if xpos:                                        #avoid out of frame mouse position
832            ypos = event.ydata
833            Page.canvas.SetCursor(wx.CROSS_CURSOR)
834            try:
835                if G2frame.Contour:
836                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA pattern ID =%5d'%(xpos,int(ypos)),1)
837                else:
838                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA %s =%.2f'%(xpos,type,ypos),1)                   
839            except TypeError:
840                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
841   
842    xylim = []
843    try:
844        plotNum = G2frame.G2plotNB.plotList.index(type)
845        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
846        if not newPlot:
847            Plot = Page.figure.gca()          #get previous plot & get limits
848            xylim = Plot.get_xlim(),Plot.get_ylim()
849        Page.figure.clf()
850        Plot = Page.figure.gca()
851    except ValueError:
852        newPlot = True
853        G2frame.Cmax = 1.0
854        Plot = G2frame.G2plotNB.addMpl(type).gca()
855        plotNum = G2frame.G2plotNB.plotList.index(type)
856        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
857        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
858        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
859   
860    Page.SetFocus()
861    G2frame.G2plotNB.status.DestroyChildren()
862    if G2frame.Contour:
863        Choice = (' key press','d: lower contour max','u: raise contour max',
864            'i: interpolation method','s: color scheme','c: contour off')
865    else:
866        Choice = (' key press','l: offset left','r: offset right','d: offset down','u: offset up',
867            'o: reset offset','t: toggle legend','c: contour on','s: toggle single plot')
868    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
869        choices=Choice)
870    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
871    cb.SetValue(' key press')
872    PatternId = G2frame.PatternId
873    PickId = G2frame.PickId
874    Plot.set_title(type)
875    if type == 'G(R)':
876        Plot.set_xlabel(r'$R,\AA$',fontsize=14)
877    else:
878        Plot.set_xlabel(r'$Q,\AA$'+superMinusOne,fontsize=14)
879    Plot.set_ylabel(r''+type,fontsize=14)
880    colors=['b','g','r','c','m','k']
881    name = G2frame.PatternTree.GetItemText(PatternId)[4:]
882    Pattern = []   
883    if G2frame.SinglePlot:
884        name = G2frame.PatternTree.GetItemText(PatternId)
885        name = type+name[4:]
886        Id = G2gd.GetPatternTreeItemId(G2frame,PatternId,name)
887        Pattern = G2frame.PatternTree.GetItemPyData(Id)
888        if Pattern:
889            Pattern.append(name)
890        PlotList = [Pattern,]
891    else:
892        PlotList = []
893        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
894        while item:
895            if 'PDF' in G2frame.PatternTree.GetItemText(item):
896                name = type+G2frame.PatternTree.GetItemText(item)[4:]
897                Id = G2gd.GetPatternTreeItemId(G2frame,item,name)
898                Pattern = G2frame.PatternTree.GetItemPyData(Id)
899                if Pattern:
900                    Pattern.append(name)
901                    PlotList.append(Pattern)
902            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
903    PDFdata = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'PDF Controls'))
904    numbDen = G2pwd.GetNumDensity(PDFdata['ElList'],PDFdata['Form Vol'])
905    Xb = [0.,10.]
906    Yb = [0.,-40.*np.pi*numbDen]
907    Ymax = 1.0
908    lenX = 0
909    for Pattern in PlotList:
910        xye = Pattern[1]
911        Ymax = max(Ymax,max(xye[1]))
912    offset = G2frame.Offset[0]*Ymax/100.0
913    if G2frame.Contour:
914        ContourZ = []
915        ContourY = []
916        Nseq = 0
917    for N,Pattern in enumerate(PlotList):
918        xye = Pattern[1]
919        if PickId:
920            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
921        X = xye[0]
922        if not lenX:
923            lenX = len(X)           
924        Y = xye[1]+offset*N
925        if G2frame.Contour:
926            if lenX == len(X):
927                ContourY.append(N)
928                ContourZ.append(Y)
929                ContourX = X
930                Nseq += 1
931                Plot.set_ylabel('Data sequence',fontsize=12)
932        else:
933            X = xye[0]+G2frame.Offset[1]*.005*N
934            if ifpicked:
935                Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
936                Page.canvas.SetToolTipString('')
937            else:
938                if G2frame.Legend:
939                    Plot.plot(X,Y,colors[N%6],picker=False,label='Azm:'+Pattern[2].split('=')[1])
940                else:
941                    Plot.plot(X,Y,colors[N%6],picker=False)
942            if type == 'G(R)':
943                Plot.plot(Xb,Yb,color='k',dashes=(5,5))
944            elif type == 'F(Q)':
945                Plot.axhline(0.,color=wx.BLACK)
946            elif type == 'S(Q)':
947                Plot.axhline(1.,color=wx.BLACK)
948    if G2frame.Contour:
949        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
950        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
951            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
952        Page.figure.colorbar(Img)
953    elif G2frame.Legend:
954        Plot.legend(loc='best')
955    if not newPlot:
956        Page.toolbar.push_current()
957        Plot.set_xlim(xylim[0])
958        Plot.set_ylim(xylim[1])
959        xylim = []
960        Page.toolbar.push_current()
961        Page.toolbar.draw()
962    else:
963        Page.canvas.draw()
964       
965def PlotXY(G2frame,XY,newPlot=False,type=''):
966    '''simple plot of xy data, used for diagnostic purposes
967    '''
968    def OnMotion(event):
969        xpos = event.xdata
970        if xpos:                                        #avoid out of frame mouse position
971            ypos = event.ydata
972            Page.canvas.SetCursor(wx.CROSS_CURSOR)
973            try:
974                G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,type,ypos),1)                   
975            except TypeError:
976                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
977
978    try:
979        plotNum = G2frame.G2plotNB.plotList.index(type)
980        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
981        if not newPlot:
982            Plot = Page.figure.gca()
983            xylim = Plot.get_xlim(),Plot.get_ylim()
984        Page.figure.clf()
985        Plot = Page.figure.gca()
986    except ValueError:
987        newPlot = True
988        Plot = G2frame.G2plotNB.addMpl(type).gca()
989        plotNum = G2frame.G2plotNB.plotList.index(type)
990        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
991        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
992   
993    Page.SetFocus()
994    G2frame.G2plotNB.status.DestroyChildren()
995    Plot.set_title(type)
996    Plot.set_xlabel(r'X',fontsize=14)
997    Plot.set_ylabel(r''+type,fontsize=14)
998    colors=['b','g','r','c','m','k']
999    Ymax = 1.0
1000    lenX = 0
1001    X,Y = XY[:2]
1002    Ymax = max(Ymax,max(Y))
1003    Plot.plot(X,Y,'k',picker=False)
1004    if not newPlot:
1005        Page.toolbar.push_current()
1006        Plot.set_xlim(xylim[0])
1007        Plot.set_ylim(xylim[1])
1008        xylim = []
1009        Page.toolbar.push_current()
1010        Page.toolbar.draw()
1011    else:
1012        Page.canvas.draw()
1013
1014def PlotPowderLines(G2frame):
1015    ''' plotting of powder lines (i.e. no powder pattern) as sticks
1016    '''
1017    global HKL
1018
1019    def OnMotion(event):
1020        xpos = event.xdata
1021        if xpos:                                        #avoid out of frame mouse position
1022            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1023            G2frame.G2plotNB.status.SetFields(['','2-theta =%9.3f '%(xpos,)])
1024            if G2frame.PickId and G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List']:
1025                found = []
1026                if len(HKL):
1027                    view = Page.toolbar._views.forward()[0][:2]
1028                    wid = view[1]-view[0]
1029                    found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
1030                if len(found):
1031                    h,k,l = found[0][:3] 
1032                    Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
1033                else:
1034                    Page.canvas.SetToolTipString('')
1035
1036    try:
1037        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
1038        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1039        Page.figure.clf()
1040        Plot = Page.figure.gca()
1041    except ValueError:
1042        Plot = G2frame.G2plotNB.addMpl('Powder Lines').gca()
1043        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
1044        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1045        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1046       
1047    Page.SetFocus()
1048    Plot.set_title('Powder Pattern Lines')
1049    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
1050    PickId = G2frame.PickId
1051    PatternId = G2frame.PatternId
1052    peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))
1053    for peak in peaks:
1054        Plot.axvline(peak[0],color='b')
1055    HKL = np.array(G2frame.HKL)
1056    for hkl in G2frame.HKL:
1057        Plot.axvline(hkl[5],color='r',dashes=(5,5))
1058    xmin = peaks[0][0]
1059    xmax = peaks[-1][0]
1060    delt = xmax-xmin
1061    xlim = [max(0,xmin-delt/20.),min(180.,xmax+delt/20.)]
1062    Plot.set_xlim(xlim)
1063    Page.canvas.draw()
1064    Page.toolbar.push_current()
1065
1066def PlotPeakWidths(G2frame):
1067    ''' Plotting of instrument broadening terms as function of 2-theta (future TOF)
1068    Seen when "Instrument Parameters" chosen from powder pattern data tree
1069    '''
1070    PatternId = G2frame.PatternId
1071    limitID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
1072    if limitID:
1073        limits = G2frame.PatternTree.GetItemPyData(limitID)
1074    else:
1075        return
1076    instParms = G2frame.PatternTree.GetItemPyData( \
1077        G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
1078    if instParms[0][0] in ['PXC','PNC']:
1079        lam = instParms[1][1]
1080        if len(instParms[1]) == 13:
1081            GU,GV,GW,LX,LY = instParms[0][6:11]
1082        else:
1083            GU,GV,GW,LX,LY = instParms[0][4:9]
1084    peakID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List')
1085    if peakID:
1086        peaks = G2frame.PatternTree.GetItemPyData(peakID)
1087    else:
1088        peaks = []
1089   
1090    try:
1091        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
1092        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1093        Page.figure.clf()
1094        Plot = Page.figure.gca()
1095    except ValueError:
1096        Plot = G2frame.G2plotNB.addMpl('Peak Widths').gca()
1097        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
1098        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1099    Page.SetFocus()
1100   
1101    Page.canvas.SetToolTipString('')
1102    colors=['b','g','r','c','m','k']
1103    Xmin,Xmax = limits[1]
1104    Xmin = min(0.5,max(Xmin,1))
1105    Xmin /= 2
1106    Xmax /= 2
1107    nPts = 100
1108    delt = (Xmax-Xmin)/nPts
1109    thetas = []
1110    for i in range(nPts):
1111        thetas.append(Xmin+i*delt)
1112    X = []
1113    Y = []
1114    Z = []
1115    W = []
1116    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
1117    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
1118    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.)
1119#    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
1120    for theta in thetas:
1121        X.append(4.0*math.pi*sind(theta)/lam)              #q
1122        s = sig(theta,GU,GV,GW)
1123        g = gam(theta,LX,LY)
1124        G = gamFW(g,s)
1125        Y.append(s/tand(theta))
1126        Z.append(g/tand(theta))
1127        W.append(G/tand(theta))
1128    Plot.set_title('Instrument and sample peak widths')
1129    Plot.set_ylabel(r'$\Delta q/q, \Delta d/d$',fontsize=14)
1130    Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
1131    Plot.plot(X,Y,color='r',label='Gaussian')
1132    Plot.plot(X,Z,color='g',label='Lorentzian')
1133    Plot.plot(X,W,color='b',label='G+L')
1134    X = []
1135    Y = []
1136    Z = []
1137    W = []
1138    V = []
1139    for peak in peaks:
1140        X.append(4.0*math.pi*sind(peak[0]/2.0)/lam)
1141        try:
1142            s = 1.17741*math.sqrt(peak[4])*math.pi/18000.
1143        except ValueError:
1144            s = 0.01
1145        g = peak[6]*math.pi/18000.
1146        G = gamFW(g,s)
1147        Y.append(s/tand(peak[0]/2.))
1148        Z.append(g/tand(peak[0]/2.))
1149        W.append(G/tand(peak[0]/2.))
1150    Plot.plot(X,Y,'+',color='r',label='G peak')
1151    Plot.plot(X,Z,'+',color='g',label='L peak')
1152    Plot.plot(X,W,'+',color='b',label='G+L peak')
1153    Plot.legend(loc='best')
1154    Page.canvas.draw()
1155   
1156def PlotSizeStrainPO(G2frame,data,Start=False):
1157    '''Plot 3D mustrain/size/preferred orientation figure. In this instance data is for a phase
1158    '''
1159    # Start is not used -- fix later
1160   
1161    PatternId = G2frame.PatternId
1162    generalData = data['General']
1163    SGData = generalData['SGData']
1164    SGLaue = SGData['SGLaue']
1165    G2lat.InitTexture()
1166    MuStrKeys = G2spc.MustrainNames(SGData)
1167    cell = generalData['Cell'][1:]
1168    A,B = G2lat.cell2AB(cell[:6])
1169    Vol = cell[6]
1170    useList = data['Histograms']
1171    phase = generalData['Name']
1172    plotType = generalData['Data plot type']
1173    plotDict = {'Mustrain':'Mustrain','Size':'Size','Preferred orientation':'Pref.Ori.'}
1174    for ptype in plotDict:
1175        G2frame.G2plotNB.Delete(ptype)       
1176
1177    for item in useList:
1178        if useList[item]['Show']:
1179            break
1180    else:
1181        return            #nothing to show!!
1182   
1183    numPlots = len(useList)
1184
1185    if plotType in ['Mustrain','Size']:
1186        Plot = mp3d.Axes3D(G2frame.G2plotNB.add3D(plotType))
1187    else:
1188        Plot = G2frame.G2plotNB.addMpl(plotType).gca()       
1189    plotNum = G2frame.G2plotNB.plotList.index(plotType)
1190    Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1191    Page.SetFocus()
1192    G2frame.G2plotNB.status.SetStatusText('Adjust frame size to get desired aspect ratio',1)
1193    if not Page.IsShown():
1194        Page.Show()
1195   
1196    for item in useList:
1197        if useList[item]['Show']:
1198            PHI = np.linspace(0.,360.,30,True)
1199            PSI = np.linspace(0.,180.,30,True)
1200            X = np.outer(npsind(PHI),npsind(PSI))
1201            Y = np.outer(npcosd(PHI),npsind(PSI))
1202            Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
1203            try:        #temp patch instead of 'mustrain' for old files with 'microstrain'
1204                coeff = useList[item][plotDict[plotType]]
1205            except KeyError:
1206                break
1207            if plotType in ['Mustrain','Size']:
1208                if coeff[0] == 'isotropic':
1209                    X *= coeff[1][0]
1210                    Y *= coeff[1][0]
1211                    Z *= coeff[1][0]                               
1212                elif coeff[0] == 'uniaxial':
1213                   
1214                    def uniaxCalc(xyz,iso,aniso,axes):
1215                        Z = np.array(axes)
1216                        cp = abs(np.dot(xyz,Z))
1217                        sp = np.sqrt(1.-cp**2)
1218                        R = iso*aniso/np.sqrt((iso*cp)**2+(aniso*sp)**2)
1219                        return R*xyz
1220                       
1221                    iso,aniso = coeff[1][:2]
1222                    axes = np.inner(A,np.array(coeff[3]))
1223                    axes /= nl.norm(axes)
1224                    Shkl = np.array(coeff[1])
1225                    XYZ = np.dstack((X,Y,Z))
1226                    XYZ = np.nan_to_num(np.apply_along_axis(uniaxCalc,2,XYZ,iso,aniso,axes))
1227                    X,Y,Z = np.dsplit(XYZ,3)
1228                    X = X[:,:,0]
1229                    Y = Y[:,:,0]
1230                    Z = Z[:,:,0]
1231               
1232                elif coeff[0] == 'ellipsoidal':
1233                   
1234                    def ellipseCalc(xyz,E,R):
1235                        XYZ = xyz*E.T
1236                        return np.inner(XYZ.T,R)
1237                       
1238                    S6 = coeff[4]
1239                    Sij = G2lat.U6toUij(S6)
1240                    E,R = nl.eigh(Sij)
1241                    XYZ = np.dstack((X,Y,Z))
1242                    XYZ = np.nan_to_num(np.apply_along_axis(ellipseCalc,2,XYZ,E,R))
1243                    X,Y,Z = np.dsplit(XYZ,3)
1244                    X = X[:,:,0]
1245                    Y = Y[:,:,0]
1246                    Z = Z[:,:,0]
1247                   
1248                elif coeff[0] == 'generalized':
1249                   
1250                    def genMustrain(xyz,SGData,A,Shkl):
1251                        uvw = np.inner(A.T,xyz)
1252                        Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
1253                        sum = np.sum(np.multiply(Shkl,Strm))
1254                        sum = np.where(sum > 0.01,sum,0.01)
1255                        sum = np.sqrt(sum)*math.pi/0.018      #centidegrees to radians!
1256                        return sum*xyz
1257                       
1258                    Shkl = np.array(coeff[4])
1259                    if np.any(Shkl):
1260                        XYZ = np.dstack((X,Y,Z))
1261                        XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,SGData,A,Shkl))
1262                        X,Y,Z = np.dsplit(XYZ,3)
1263                        X = X[:,:,0]
1264                        Y = Y[:,:,0]
1265                        Z = Z[:,:,0]
1266                           
1267                if np.any(X) and np.any(Y) and np.any(Z):
1268                    errFlags = np.seterr(all='ignore')
1269                    Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1)
1270                    np.seterr(all='ignore')
1271                    xyzlim = np.array([Plot.get_xlim3d(),Plot.get_ylim3d(),Plot.get_zlim3d()]).T
1272                    XYZlim = [min(xyzlim[0]),max(xyzlim[1])]
1273                    Plot.set_xlim3d(XYZlim)
1274                    Plot.set_ylim3d(XYZlim)
1275                    Plot.set_zlim3d(XYZlim)
1276                    Plot.set_aspect('equal')
1277                if plotType == 'Size':
1278                    Plot.set_title('Crystallite size for '+phase+'\n'+coeff[0]+' model')
1279                    Plot.set_xlabel(r'X, $\mu$m')
1280                    Plot.set_ylabel(r'Y, $\mu$m')
1281                    Plot.set_zlabel(r'Z, $\mu$m')
1282                else:   
1283                    Plot.set_title(r'$\mu$strain for '+phase+'\n'+coeff[0]+' model')
1284                    Plot.set_xlabel(r'X, $\mu$strain')
1285                    Plot.set_ylabel(r'Y, $\mu$strain')
1286                    Plot.set_zlabel(r'Z, $\mu$strain')
1287            else:
1288                h,k,l = generalData['POhkl']
1289                if coeff[0] == 'MD':
1290                    print 'March-Dollase preferred orientation plot'
1291               
1292                else:
1293                    PH = np.array(generalData['POhkl'])
1294                    phi,beta = G2lat.CrsAng(PH,cell[:6],SGData)
1295                    SHCoef = {}
1296                    for item in coeff[5]:
1297                        L,N = eval(item.strip('C'))
1298                        SHCoef['C%d,0,%d'%(L,N)] = coeff[5][item]                       
1299                    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1300                    X = np.linspace(0,90.0,26)
1301                    Y = G2lat.polfcal(ODFln,'0',X,0.0)
1302                    Plot.plot(X,Y,color='k',label=str(PH))
1303                    Plot.legend(loc='best')
1304                    Plot.set_title('Axial distribution for HKL='+str(PH)+' in '+phase)
1305                    Plot.set_xlabel(r'$\psi$',fontsize=16)
1306                    Plot.set_ylabel('MRD',fontsize=14)
1307    Page.canvas.draw()
1308   
1309def PlotTexture(G2frame,data,newPlot=False,Start=False):
1310    '''Pole figure, inverse pole figure, 3D pole distribution and 3D inverse pole distribution
1311    plotting.
1312    dict generalData contains all phase info needed which is in data
1313    '''
1314
1315    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1316    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1317    PatternId = G2frame.PatternId
1318    generalData = data['General']
1319    SGData = generalData['SGData']
1320    textureData = generalData['SH Texture']
1321    G2frame.G2plotNB.Delete('Texture')
1322    if not textureData['Order']:
1323        return                  #no plot!!
1324    SHData = generalData['SH Texture']
1325    SHCoef = SHData['SH Coeff'][1]
1326    cell = generalData['Cell'][1:7]
1327    Amat,Bmat = G2lat.cell2AB(cell)
1328    sq2 = 1.0/math.sqrt(2.0)
1329   
1330    def rp2xyz(r,p):
1331        z = npcosd(r)
1332        xy = np.sqrt(1.-z**2)
1333        return xy*npsind(p),xy*npcosd(p),z
1334           
1335    def OnMotion(event):
1336        SHData = data['General']['SH Texture']
1337        if event.xdata and event.ydata:                 #avoid out of frame errors
1338            xpos = event.xdata
1339            ypos = event.ydata
1340            if 'Inverse' in SHData['PlotType']:
1341                r = xpos**2+ypos**2
1342                if r <= 1.0:
1343                    if 'equal' in G2frame.Projection: 
1344                        r,p = 2.*npasind(np.sqrt(r)*sq2),npatan2d(ypos,xpos)
1345                    else:
1346                        r,p = 2.*npatand(np.sqrt(r)),npatan2d(ypos,xpos)
1347                    ipf = G2lat.invpolfcal(IODFln,SGData,np.array([r,]),np.array([p,]))
1348                    xyz = np.inner(Amat.T,np.array([rp2xyz(r,p)]))
1349                    y,x,z = list(xyz/np.max(np.abs(xyz)))
1350                   
1351                    G2frame.G2plotNB.status.SetFields(['',
1352                        'psi =%9.3f, beta =%9.3f, MRD =%9.3f hkl=%5.2f,%5.2f,%5.2f'%(r,p,ipf,x,y,z)])
1353                                   
1354            elif 'Axial' in SHData['PlotType']:
1355                pass
1356               
1357            else:                       #ordinary pole figure
1358                z = xpos**2+ypos**2
1359                if z <= 1.0:
1360                    z = np.sqrt(z)
1361                    if 'equal' in G2frame.Projection: 
1362                        r,p = 2.*npasind(z*sq2),npatan2d(ypos,xpos)
1363                    else:
1364                        r,p = 2.*npatand(z),npatan2d(ypos,xpos)
1365                    pf = G2lat.polfcal(ODFln,SamSym[textureData['Model']],np.array([r,]),np.array([p,]))
1366                    G2frame.G2plotNB.status.SetFields(['','phi =%9.3f, gam =%9.3f, MRD =%9.3f'%(r,p,pf)])
1367   
1368    if G2frame.Projection == '3D display' and 'Axial'  not in SHData['PlotType']:               
1369        Plot = mp3d.Axes3D(G2frame.G2plotNB.add3D('Texture'))
1370    else:
1371        Plot = G2frame.G2plotNB.addMpl('Texture').gca()
1372    plotNum = G2frame.G2plotNB.plotList.index('Texture')
1373    Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1374    if G2frame.Projection == '3D display':
1375        pass
1376    else:
1377        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1378
1379    Page.SetFocus()
1380    G2frame.G2plotNB.status.SetFields(['',''])   
1381    PH = np.array(SHData['PFhkl'])
1382    phi,beta = G2lat.CrsAng(PH,cell,SGData)
1383    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1384    PX = np.array(SHData['PFxyz'])
1385    gam = atan2d(PX[0],PX[1])
1386    xy = math.sqrt(PX[0]**2+PX[1]**2)
1387    xyz = math.sqrt(PX[0]**2+PX[1]**2+PX[2]**2)
1388    psi = asind(xy/xyz)
1389    IODFln = G2lat.Glnh(Start,SHCoef,psi,gam,SamSym[textureData['Model']])
1390    if 'Axial' in SHData['PlotType']:
1391        X = np.linspace(0,90.0,26)
1392        Y = G2lat.polfcal(ODFln,SamSym[textureData['Model']],X,0.0)
1393        Plot.plot(X,Y,color='k',label=str(SHData['PFhkl']))
1394        Plot.legend(loc='best')
1395        Plot.set_title('Axial distribution for HKL='+str(SHData['PFhkl']))
1396        Plot.set_xlabel(r'$\psi$',fontsize=16)
1397        Plot.set_ylabel('MRD',fontsize=14)
1398       
1399    else:       
1400        npts = 201
1401        if 'Inverse' in SHData['PlotType']:
1402            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
1403            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
1404            if 'equal' in G2frame.Projection:
1405                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
1406            else:
1407                R = np.where(R <= 1.,2.*npatand(R),0.0)
1408            Z = np.zeros_like(R)
1409            Z = G2lat.invpolfcal(IODFln,SGData,R,P)
1410            Z = np.reshape(Z,(npts,npts))
1411            CS = Plot.contour(Y,X,Z,aspect='equal')
1412            Plot.clabel(CS,fontsize=9,inline=1)
1413            try:
1414                Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
1415            except ValueError:
1416                pass
1417            if newPlot:
1418#                Page.figure.colorbar(Img)    #colorbar fails - crashes gsasii
1419                newPlot = False
1420            Plot.set_title('Inverse pole figure for XYZ='+str(SHData['PFxyz']))
1421            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
1422                       
1423        else:
1424            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
1425            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
1426            if 'equal' in G2frame.Projection:
1427                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
1428            else:
1429                R = np.where(R <= 1.,2.*npatand(R),0.0)
1430            Z = np.zeros_like(R)
1431            Z = G2lat.polfcal(ODFln,SamSym[textureData['Model']],R,P)
1432            Z = np.reshape(Z,(npts,npts))
1433            try:
1434                CS = Plot.contour(Y,X,Z,aspect='equal')
1435                Plot.clabel(CS,fontsize=9,inline=1)
1436            except ValueError:
1437                pass
1438            Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
1439            if newPlot:
1440#                Page.figure.colorbar(Img)    #colorbar fails - crashes gsasii
1441                newPlot = False
1442            Plot.set_title('Pole figure for HKL='+str(SHData['PFhkl']))
1443            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
1444    Page.canvas.draw()
1445
1446def PlotCovariance(G2frame,Data={}):
1447    if not Data:
1448        Data = G2frame.PatternTree.GetItemPyData(
1449            G2gd.GetPatternTreeItemId(G2frame,G2frame.root, 'Covariance'))
1450    if not Data:
1451        print 'No covariance matrix available'
1452        return
1453    varyList = Data['varyList']
1454    values = Data['variables']
1455    Xmax = len(varyList)
1456    covMatrix = Data['covMatrix']
1457    sig = np.sqrt(np.diag(covMatrix))
1458    xvar = np.outer(sig,np.ones_like(sig))
1459    covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
1460    title = ' for\n'+Data['title']
1461    newAtomDict = Data['newAtomDict']
1462
1463    def OnPlotKeyPress(event):
1464        newPlot = False
1465        if event.key == 's':
1466            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
1467            choice.sort()
1468            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
1469            if dlg.ShowModal() == wx.ID_OK:
1470                sel = dlg.GetSelection()
1471                G2frame.VcovColor = choice[sel]
1472            else:
1473                G2frame.VcovColor = 'RdYlGn'
1474            dlg.Destroy()
1475        PlotCovariance(G2frame)
1476
1477    def OnMotion(event):
1478        if event.button:
1479            ytics = imgAx.get_yticks()
1480            ytics = np.where(ytics<len(varyList),ytics,-1)
1481            ylabs = [np.where(0<=i ,varyList[int(i)],' ') for i in ytics]
1482            imgAx.set_yticklabels(ylabs)
1483           
1484        if event.xdata and event.ydata:                 #avoid out of frame errors
1485            xpos = int(event.xdata+.5)
1486            ypos = int(event.ydata+.5)
1487            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
1488                if xpos == ypos:
1489                    value = values[xpos]
1490                    name = varyList[xpos]
1491                    if varyList[xpos] in newAtomDict:
1492                        name,value = newAtomDict[name]                       
1493                    msg = '%s value = %.4g, esd = %.4g'%(name,value,sig[xpos])
1494                else:
1495                    msg = '%s - %s: %5.3f'%(varyList[xpos],varyList[ypos],covArray[xpos][ypos])
1496                Page.canvas.SetToolTipString(msg)
1497                G2frame.G2plotNB.status.SetFields(['Key: s to change colors',msg])
1498    try:
1499        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
1500        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1501        Page.figure.clf()
1502        Plot = Page.figure.gca()
1503        if not Page.IsShown():
1504            Page.Show()
1505    except ValueError:
1506        Plot = G2frame.G2plotNB.addMpl('Covariance').gca()
1507        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
1508        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1509        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1510        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1511
1512    Page.SetFocus()
1513    G2frame.G2plotNB.status.SetFields(['',''])   
1514    acolor = mpl.cm.get_cmap(G2frame.VcovColor)
1515    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower')
1516    imgAx = Img.get_axes()
1517    ytics = imgAx.get_yticks()
1518    ylabs = [varyList[int(i)] for i in ytics[:-1]]
1519    imgAx.set_yticklabels(ylabs)
1520    colorBar = Page.figure.colorbar(Img)
1521    Plot.set_title('V-Cov matrix'+title)
1522    Plot.set_xlabel('Variable number')
1523    Plot.set_ylabel('Variable name')
1524    Page.canvas.draw()
1525   
1526def PlotSeq(G2frame,SeqData,SeqSig,SeqNames,sampleParm):
1527   
1528    def OnKeyPress(event):
1529        if event.key == 's' and sampleParm:
1530            if G2frame.xAxis:
1531                G2frame.xAxis = False
1532            else:
1533                G2frame.xAxis = True
1534            Draw(False)
1535    try:
1536        plotNum = G2frame.G2plotNB.plotList.index('Sequential refinement')
1537        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1538        Page.figure.clf()
1539        Plot = Page.figure.gca()
1540        if not Page.IsShown():
1541            Page.Show()
1542    except ValueError:
1543        Plot = G2frame.G2plotNB.addMpl('Sequential refinement').gca()
1544        plotNum = G2frame.G2plotNB.plotList.index('Sequential refinement')
1545        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1546        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
1547        G2frame.xAxis = False
1548       
1549    def Draw(newPlot):
1550        Page.SetFocus()
1551        G2frame.G2plotNB.status.SetFields(['','press s to toggle x-axis = sample environment parameter'])
1552        if len(SeqData):
1553            Plot.clear()
1554            if G2frame.xAxis:   
1555                xName = sampleParm.keys()[0]
1556                X = sampleParm[xName]
1557            else:
1558                X = np.arange(0,len(SeqData[0]),1)
1559                xName = 'Data sequence number'
1560            for Y,sig,name in zip(SeqData,SeqSig,SeqNames):
1561                Plot.errorbar(X,Y,yerr=sig,label=name)       
1562            Plot.legend(loc='best')
1563            Plot.set_ylabel('Parameter values')
1564            Plot.set_xlabel(xName)
1565            Page.canvas.draw()           
1566    Draw(True)
1567           
1568def PlotExposedImage(G2frame,newPlot=False,event=None):
1569    '''General access module for 2D image plotting
1570    '''
1571    plotNo = G2frame.G2plotNB.nb.GetSelection()
1572    if G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
1573        PlotImage(G2frame,newPlot,event,newImage=True)
1574    elif G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
1575        PlotIntegration(G2frame,newPlot,event)
1576
1577def PlotImage(G2frame,newPlot=False,event=None,newImage=True):
1578    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
1579    masks, etc.
1580    '''
1581    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
1582    import numpy.ma as ma
1583    Dsp = lambda tth,wave: wave/(2.*sind(tth/2.))
1584    global Data,Masks
1585    colors=['b','g','r','c','m','k']
1586    Data = G2frame.PatternTree.GetItemPyData(
1587        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
1588    Masks = G2frame.PatternTree.GetItemPyData(
1589        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
1590
1591    def OnImMotion(event):
1592        Page.canvas.SetToolTipString('')
1593        sizexy = Data['size']
1594        if event.xdata and event.ydata:                 #avoid out of frame errors
1595            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1596            item = G2frame.itemPicked
1597            pixelSize = Data['pixelSize']
1598            scalex = 1000./pixelSize[0]
1599            scaley = 1000./pixelSize[1]
1600            if item and G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Image Controls':
1601                if 'Text' in str(item):
1602                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
1603                else:
1604                    xcent,ycent = Data['center']
1605                    xpos = event.xdata-xcent
1606                    ypos = event.ydata-ycent
1607                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
1608                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
1609                        Page.canvas.SetToolTipString('%6d deg'%(azm))
1610                    elif 'line1' in  str(item) or 'line2' in str(item):
1611                        Page.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
1612            else:
1613                xpos = event.xdata
1614                ypos = event.ydata
1615                xpix = xpos*scalex
1616                ypix = ypos*scaley
1617                Int = 0
1618                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
1619                    Int = G2frame.ImageZ[ypix][xpix]
1620                tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
1621                Q = 2.*math.pi/dsp
1622                if G2frame.setPoly:
1623                    G2frame.G2plotNB.status.SetFields(['','Polygon mask pick - LB next point, RB close polygon'])
1624                else:
1625                    G2frame.G2plotNB.status.SetFields(\
1626                        ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)])
1627
1628    def OnImPlotKeyPress(event):
1629        try:
1630            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
1631        except TypeError:
1632            return
1633        if PickName == 'Masks':
1634            Xpos = event.xdata
1635            if not Xpos:            #got point out of frame
1636                return
1637            Ypos = event.ydata
1638            if event.key == 's':
1639                Masks['Points'].append([Xpos,Ypos,1.])
1640            elif event.key == 'r':
1641                tth = G2img.GetTth(Xpos,Ypos,Data)
1642                Masks['Rings'].append([tth,0.1])
1643            elif event.key == 'a':
1644                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
1645                azm = int(azm)               
1646                Masks['Arcs'].append([tth,[azm-5,azm+5],0.1])
1647            elif event.key == 'p':
1648                G2frame.setPoly = True
1649                Masks['Polygons'].append([])
1650                G2frame.G2plotNB.status.SetFields(['','Polygon mask active - LB pick next point, RB close polygon'])
1651            G2imG.UpdateMasks(G2frame,Masks)
1652        elif PickName == 'Image Controls':
1653            if event.key == 'c':
1654                Xpos = event.xdata
1655                if not Xpos:            #got point out of frame
1656                    return
1657                Ypos = event.ydata
1658                dlg = wx.MessageDialog(G2frame,'Are you sure you want to change the center?',
1659                    'Center change',style=wx.OK|wx.CANCEL)
1660                try:
1661                    if dlg.ShowModal() == wx.ID_OK:
1662                        print 'move center to: ',Xpos,Ypos
1663                        Data['center'] = [Xpos,Ypos]
1664                        G2imG.UpdateImageControls(G2frame,Data,Masks)
1665                finally:
1666                    dlg.Destroy()
1667            elif event.key == 'l':
1668                if G2frame.logPlot:
1669                    G2frame.logPlot = False
1670                else:
1671                    G2frame.logPlot = True
1672        PlotImage(G2frame,newImage=True)
1673           
1674    def OnKeyBox(event):
1675        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):
1676            event.key = cb.GetValue()[0]
1677            cb.SetValue(' key press')
1678            if event.key in 'l':
1679                OnImPlotKeyPress(event)
1680                       
1681    def OnImPick(event):
1682        if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']:
1683            return
1684        if G2frame.setPoly:
1685            polygon = Masks['Polygons'][-1]
1686            xpos,ypos = event.mouseevent.xdata,event.mouseevent.ydata
1687            if xpos and ypos:                       #point inside image
1688                if len(polygon) > 2 and event.mouseevent.button == 3:
1689                    x0,y0 = polygon[0]
1690                    polygon.append([x0,y0])
1691                    G2frame.setPoly = False
1692                    G2frame.G2plotNB.status.SetFields(['','Polygon closed - RB drag a vertex to change shape'])
1693                else:
1694                    G2frame.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(xpos,ypos)])
1695                    polygon.append([xpos,ypos])
1696                G2imG.UpdateMasks(G2frame,Masks)
1697        else:
1698            if G2frame.itemPicked is not None: return
1699            G2frame.itemPicked = event.artist
1700            G2frame.mousePicked = event.mouseevent
1701       
1702    def OnImRelease(event):
1703        try:
1704            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
1705        except TypeError:
1706            return
1707        if PickName not in ['Image Controls','Masks']:
1708            return
1709        pixelSize = Data['pixelSize']
1710        scalex = 1000./pixelSize[0]
1711        scaley = 1000./pixelSize[1]
1712        pixLimit = Data['pixLimit']
1713        if G2frame.itemPicked is None and PickName == 'Image Controls':
1714#            sizexy = Data['size']
1715            Xpos = event.xdata
1716            if not (Xpos and G2frame.ifGetRing):                   #got point out of frame
1717                return
1718            Ypos = event.ydata
1719            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
1720                if event.button == 1:
1721                    Xpix = Xpos*scalex
1722                    Ypix = Ypos*scaley
1723                    xpos,ypos,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
1724                    if I and J:
1725                        xpos += .5                              #shift to pixel center
1726                        ypos += .5
1727                        xpos /= scalex                          #convert to mm
1728                        ypos /= scaley
1729                        Data['ring'].append([xpos,ypos])
1730                elif event.button == 3:
1731                    G2frame.dataFrame.GetStatusBar().SetStatusText('Calibrating...')
1732                    if G2img.ImageCalibrate(G2frame,Data):
1733                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check')
1734                        print 'Calibration successful'
1735                    else:
1736                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose')
1737                        print 'Calibration failed'
1738                    G2frame.ifGetRing = False
1739                    G2imG.UpdateImageControls(G2frame,Data,Masks)
1740                    return
1741                PlotImage(G2frame,newImage=False)
1742            return
1743        else:
1744            xpos = event.xdata
1745            if xpos:                                        #avoid out of frame mouse position
1746                ypos = event.ydata
1747                if G2frame.ifGetRing:                          #delete a calibration ring pick
1748                    xypos = [xpos,ypos]
1749                    rings = Data['ring']
1750                    for ring in rings:
1751                        if np.allclose(ring,xypos,.01,0):
1752                            rings.remove(ring)                                                                       
1753                else:
1754                    tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
1755                    itemPicked = str(G2frame.itemPicked)
1756                    if 'Line2D' in itemPicked and PickName == 'Image Controls':
1757                        if 'line1' in itemPicked:
1758                            Data['IOtth'][0] = max(tth,0.001)
1759                        elif 'line2' in itemPicked:
1760                            Data['IOtth'][1] = tth
1761                        elif 'line3' in itemPicked:
1762                            Data['LRazimuth'][0] = int(azm)
1763                            if Data['fullIntegrate']:
1764                                Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
1765                        elif 'line4' in itemPicked and not Data['fullIntegrate']:
1766                            Data['LRazimuth'][1] = int(azm)
1767                           
1768                        if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
1769                            Data['LRazimuth'][0] -= 360
1770                           
1771                        azLim = np.array(Data['LRazimuth'])   
1772                        if np.any(azLim>360):
1773                            azLim -= 360
1774                            Data['LRazimuth'] = list(azLim)
1775                           
1776                        if  Data['IOtth'][0] > Data['IOtth'][1]:
1777                            Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
1778                           
1779                        G2frame.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
1780                        G2frame.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
1781                        G2frame.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
1782                        G2frame.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
1783                    elif 'Circle' in itemPicked and PickName == 'Masks':
1784                        spots = Masks['Points']
1785                        newPos = itemPicked.split(')')[0].split('(')[2].split(',')
1786                        newPos = np.array([float(newPos[0]),float(newPos[1])])
1787                        for spot in spots:
1788                            if np.allclose(np.array([spot[:2]]),newPos):
1789                                spot[:2] = xpos,ypos
1790                        G2imG.UpdateMasks(G2frame,Masks)
1791                    elif 'Line2D' in itemPicked and PickName == 'Masks':
1792                        Obj = G2frame.itemPicked.findobj()
1793                        rings = Masks['Rings']
1794                        arcs = Masks['Arcs']
1795                        polygons = Masks['Polygons']
1796                        for ring in G2frame.ringList:
1797                            if Obj == ring[0]:
1798                                rN = ring[1]
1799                                if ring[2] == 'o':
1800                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)-rings[rN][1]/2.
1801                                else:
1802                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)+rings[rN][1]/2.
1803                        for arc in G2frame.arcList:
1804                            if Obj == arc[0]:
1805                                aN = arc[1]
1806                                if arc[2] == 'o':
1807                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)-arcs[aN][2]/2
1808                                elif arc[2] == 'i':
1809                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)+arcs[aN][2]/2
1810                                elif arc[2] == 'l':
1811                                    arcs[aN][1][0] = int(G2img.GetAzm(xpos,ypos,Data))
1812                                else:
1813                                    arcs[aN][1][1] = int(G2img.GetAzm(xpos,ypos,Data))
1814                        for poly in G2frame.polyList:
1815                            if Obj == poly[0]:
1816                                ind = G2frame.itemPicked.contains(G2frame.mousePicked)[1]['ind'][0]
1817                                oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
1818                                pN = poly[1]
1819                                for i,xy in enumerate(polygons[pN]):
1820                                    if np.allclose(np.array([xy]),oldPos,atol=1.0):
1821                                        polygons[pN][i] = xpos,ypos
1822                        G2imG.UpdateMasks(G2frame,Masks)
1823#                    else:                  #keep for future debugging
1824#                        print str(G2frame.itemPicked),event.xdata,event.ydata,event.button
1825                PlotImage(G2frame,newImage=True)
1826            G2frame.itemPicked = None
1827           
1828    try:
1829        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
1830        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1831        if not newPlot:
1832            Plot = Page.figure.gca()          #get previous powder plot & get limits
1833            xylim = Plot.get_xlim(),Plot.get_ylim()
1834        if newImage:
1835            Page.figure.clf()
1836            Plot = Page.figure.gca()          #get a fresh plot after clf()
1837       
1838    except ValueError:
1839        Plot = G2frame.G2plotNB.addMpl('2D Powder Image').gca()
1840        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
1841        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1842        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
1843        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
1844        Page.canvas.mpl_connect('pick_event', OnImPick)
1845        Page.canvas.mpl_connect('button_release_event', OnImRelease)
1846        xylim = []
1847    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
1848        Page.SetFocus()
1849    Title = G2frame.PatternTree.GetItemText(G2frame.Image)[4:]
1850    G2frame.G2plotNB.status.DestroyChildren()
1851    if G2frame.logPlot:
1852        Title = 'log('+Title+')'
1853    Plot.set_title(Title)
1854    try:
1855        if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Image Controls',]:
1856            Choice = (' key press','l: log(I) on',)
1857            if G2frame.logPlot:
1858                Choice[1] = 'l: log(I) off'
1859            cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
1860                choices=Choice)
1861            cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
1862            cb.SetValue(' key press')
1863    except TypeError:
1864        pass
1865    size,imagefile = G2frame.PatternTree.GetItemPyData(G2frame.Image)
1866    if imagefile != G2frame.oldImagefile:
1867        imagefile = G2IO.CheckImageFile(G2frame,imagefile)
1868        if not imagefile:
1869            G2frame.G2plotNB.Delete('2D Powder Image')
1870            return
1871        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
1872        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
1873#        print G2frame.ImageZ.shape,G2frame.ImageZ.size,Data['size'] #might be useful debugging line
1874        G2frame.oldImagefile = imagefile
1875
1876    imScale = 1
1877    if len(G2frame.ImageZ) > 1024:
1878        imScale = len(G2frame.ImageZ)/1024
1879    sizexy = Data['size']
1880    pixelSize = Data['pixelSize']
1881    scalex = 1000./pixelSize[0]
1882    scaley = 1000./pixelSize[1]
1883    Xmax = sizexy[0]*pixelSize[0]/1000.
1884    Ymax = sizexy[1]*pixelSize[1]/1000.
1885    xlim = (0,Xmax)
1886    ylim = (Ymax,0)
1887    Imin,Imax = Data['range'][1]
1888    acolor = mpl.cm.get_cmap(Data['color'])
1889    xcent,ycent = Data['center']
1890    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
1891    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
1892    #do threshold mask - "real" mask - others are just bondaries
1893    Zlim = Masks['Thresholds'][1]
1894    wx.BeginBusyCursor()
1895    try:
1896           
1897        if newImage:                   
1898            MA = ma.masked_greater(ma.masked_less(G2frame.ImageZ,Zlim[0]),Zlim[1])
1899            MaskA = ma.getmaskarray(MA)
1900            A = G2img.ImageCompress(MA,imScale)
1901            AM = G2img.ImageCompress(MaskA,imScale)
1902            if G2frame.logPlot:
1903                A = np.where(A>0,np.log(A),0)
1904                AM = np.where(AM>0,np.log(AM),0)
1905                Imin,Imax = [np.amin(A),np.amax(A)]
1906            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
1907                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Ymax,0])
1908            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
1909                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
1910            if G2frame.setPoly:
1911                Img.set_picker(True)
1912   
1913        Plot.plot(xcent,ycent,'x')
1914        if Data['showLines']:
1915            LRAzim = Data['LRazimuth']                  #NB: integers
1916            Nazm = Data['outAzimuths']
1917            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
1918            AzmthOff = Data['azmthOff']
1919            IOtth = Data['IOtth']
1920            wave = Data['wavelength']
1921            dspI = wave/(2.0*sind(IOtth[0]/2.0))
1922            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
1923            dspO = wave/(2.0*sind(IOtth[1]/2.0))
1924            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
1925            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
1926            if ellI:
1927                xyI = []
1928                for azm in Azm:
1929                    xyI.append(G2img.GetDetectorXY(dspI,azm-90.,Data))
1930                xyI = np.array(xyI)
1931                arcxI,arcyI = xyI.T
1932                Plot.plot(arcxI,arcyI,picker=3)
1933            if ellO:
1934                xyO = []
1935                for azm in Azm:
1936                    xyO.append(G2img.GetDetectorXY(dspO,azm-90.,Data))
1937                xyO = np.array(xyO)
1938                arcxO,arcyO = xyO.T
1939                Plot.plot(arcxO,arcyO,picker=3)
1940            if ellO and ellI:
1941                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
1942                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
1943            for i in range(Nazm):
1944                cake = LRAzim[0]+i*delAzm-AzmthOff
1945                ind = np.searchsorted(Azm,cake)
1946                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
1947                   
1948        for xring,yring in Data['ring']:
1949            Plot.plot(xring,yring,'r+',picker=3)
1950        if Data['setRings']:
1951#            rings = np.concatenate((Data['rings']),axis=0)
1952            N = 0
1953            for ring in Data['rings']:
1954                xring,yring = np.array(ring).T[:2]
1955                Plot.plot(xring,yring,'+',color=colors[N%6])
1956                N += 1           
1957        for ellipse in Data['ellipses']:
1958            cent,phi,[width,height],col = ellipse
1959            Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
1960            Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
1961        #masks - mask lines numbered after integration limit lines
1962        spots = Masks['Points']
1963        rings = Masks['Rings']
1964        arcs = Masks['Arcs']
1965        polygons = Masks['Polygons']
1966        for x,y,d in spots:
1967            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
1968        G2frame.ringList = []
1969        for iring,(tth,thick) in enumerate(rings):
1970            wave = Data['wavelength']
1971            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
1972            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
1973            G2frame.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
1974            G2frame.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
1975        G2frame.arcList = []
1976        for iarc,(tth,azm,thick) in enumerate(arcs):           
1977            wave = Data['wavelength']
1978            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
1979            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
1980            G2frame.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
1981            G2frame.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
1982            G2frame.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
1983            G2frame.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
1984        G2frame.polyList = []
1985        for ipoly,polygon in enumerate(polygons):
1986            x,y = np.hsplit(np.array(polygon),2)
1987            G2frame.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
1988            Plot.plot(x,y,'r')           
1989        if newImage:
1990            colorBar = Page.figure.colorbar(Img)
1991        Plot.set_xlim(xlim)
1992        Plot.set_ylim(ylim)
1993        Plot.invert_yaxis()
1994        if not newPlot and xylim:
1995            Page.toolbar.push_current()
1996            Plot.set_xlim(xylim[0])
1997            Plot.set_ylim(xylim[1])
1998            xylim = []
1999            Page.toolbar.push_current()
2000            Page.toolbar.draw()
2001        else:
2002            Page.canvas.draw()
2003    finally:
2004        wx.EndBusyCursor()
2005       
2006def PlotIntegration(G2frame,newPlot=False,event=None):
2007    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
2008    '''
2009           
2010    def OnMotion(event):
2011        Page.canvas.SetToolTipString('')
2012        Page.canvas.SetCursor(wx.CROSS_CURSOR)
2013        azm = event.ydata
2014        tth = event.xdata
2015        if azm and tth:
2016            G2frame.G2plotNB.status.SetFields(\
2017                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
2018                               
2019    try:
2020        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
2021        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2022        if not newPlot:
2023            Plot = Page.figure.gca()          #get previous plot & get limits
2024            xylim = Plot.get_xlim(),Plot.get_ylim()
2025        Page.figure.clf()
2026        Plot = Page.figure.gca()          #get a fresh plot after clf()
2027       
2028    except ValueError:
2029        Plot = G2frame.G2plotNB.addMpl('2D Integration').gca()
2030        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
2031        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2032        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2033        Page.views = False
2034        view = False
2035    if not event:
2036        Page.SetFocus()
2037       
2038    Data = G2frame.PatternTree.GetItemPyData(
2039        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2040    image = G2frame.Integrate[0]
2041    xsc = G2frame.Integrate[1]
2042    ysc = G2frame.Integrate[2]
2043    Imin,Imax = Data['range'][1]
2044    acolor = mpl.cm.get_cmap(Data['color'])
2045    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
2046    Plot.set_ylabel('azimuth',fontsize=12)
2047    Plot.set_xlabel('2-theta',fontsize=12)
2048    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
2049        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
2050    colorBar = Page.figure.colorbar(Img)
2051    if Data['ellipses']:           
2052        for ellipse in Data['ellipses']:
2053            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
2054            x,y = np.hsplit(ring,2)
2055            tth,azm = G2img.GetTthAzm(x,y,Data)
2056#            azm = np.where(azm < 0.,azm+360,azm)
2057            Plot.plot(tth,azm,'b,')
2058    if not newPlot:
2059        Page.toolbar.push_current()
2060        Plot.set_xlim(xylim[0])
2061        Plot.set_ylim(xylim[1])
2062        xylim = []
2063        Page.toolbar.push_current()
2064        Page.toolbar.draw()
2065    else:
2066        Page.canvas.draw()
2067               
2068def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):
2069    '''a test plot routine - not normally used
2070    ''' 
2071           
2072    def OnMotion(event):
2073        Page.canvas.SetToolTipString('')
2074        Page.canvas.SetCursor(wx.CROSS_CURSOR)
2075        azm = event.xdata
2076        tth = event.ydata
2077        if azm and tth:
2078            G2frame.G2plotNB.status.SetFields(\
2079                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
2080                               
2081    try:
2082        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
2083        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2084        if not newPlot:
2085            Plot = Page.figure.gca()          #get previous plot & get limits
2086            xylim = Plot.get_xlim(),Plot.get_ylim()
2087        Page.figure.clf()
2088        Plot = Page.figure.gca()          #get a fresh plot after clf()
2089       
2090    except ValueError:
2091        Plot = G2frame.G2plotNB.addMpl('2D Transformed Powder Image').gca()
2092        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
2093        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2094        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2095        Page.views = False
2096        view = False
2097    Page.SetFocus()
2098       
2099    Data = G2frame.PatternTree.GetItemPyData(
2100        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2101    Imin,Imax = Data['range'][1]
2102    step = (Imax-Imin)/5.
2103    V = np.arange(Imin,Imax,step)
2104    acolor = mpl.cm.get_cmap(Data['color'])
2105    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
2106    Plot.set_xlabel('azimuth',fontsize=12)
2107    Plot.set_ylabel('2-theta',fontsize=12)
2108    Plot.contour(tax,tay,taz,V,cmap=acolor)
2109    if Data['showLines']:
2110        IOtth = Data['IOtth']
2111        if Data['fullIntegrate']:
2112            LRAzim = [-180,180]
2113        else:
2114            LRAzim = Data['LRazimuth']                  #NB: integers
2115        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
2116        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
2117        if not Data['fullIntegrate']:
2118            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
2119            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
2120    if Data['setRings']:
2121        rings = np.concatenate((Data['rings']),axis=0)
2122        for xring,yring,dsp in rings:
2123            x,y = G2img.GetTthAzm(xring,yring,Data)
2124            Plot.plot(y,x,'r+')           
2125    if Data['ellipses']:           
2126        for ellipse in Data['ellipses']:
2127            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
2128            x,y = np.hsplit(ring,2)
2129            tth,azm = G2img.GetTthAzm(x,y,Data)
2130            Plot.plot(azm,tth,'b,')
2131    if not newPlot:
2132        Page.toolbar.push_current()
2133        Plot.set_xlim(xylim[0])
2134        Plot.set_ylim(xylim[1])
2135        xylim = []
2136        Page.toolbar.push_current()
2137        Page.toolbar.draw()
2138    else:
2139        Page.canvas.draw()
2140       
2141def PlotStructure(G2frame,data):
2142    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
2143    thermal motion ellipsoids and polyhedra
2144    '''
2145    generalData = data['General']
2146    cell = generalData['Cell'][1:7]
2147    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
2148    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2149    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2150    Mydir = generalData['Mydir']
2151    atomData = data['Atoms']
2152    drawingData = data['Drawing']
2153    drawAtoms = drawingData['Atoms']
2154    cx,ct,cs = drawingData['atomPtrs']
2155    Wt = [255,255,255]
2156    Rd = [255,0,0]
2157    Gr = [0,255,0]
2158    Bl = [0,0,255]
2159    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]])
2160    uEdges = np.array([
2161        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
2162        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
2163        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
2164    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
2165    altDown = False
2166    shiftDown = False
2167    ctrlDown = False
2168   
2169    def OnKeyBox(event):
2170        import Image
2171        Draw()                          #make sure plot is fresh!!
2172        mode = cb.GetValue()
2173        Fname = Mydir+'\\'+generalData['Name']+'.'+mode
2174        size = Page.canvas.GetSize()
2175        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
2176        if mode in ['jpeg',]:
2177            Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
2178            im = Image.new("RGBA", (size[0],size[1]))
2179        else:
2180            Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
2181            im = Image.new("RGB", (size[0],size[1]))
2182        im.fromstring(Pix)
2183        im.save(Fname,mode)
2184        cb.SetValue(' Save as:')
2185        G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
2186   
2187    def GetTruePosition(xy,Add=False):
2188        View = glGetIntegerv(GL_VIEWPORT)
2189        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
2190        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
2191        Zmax = 1.
2192        if Add:
2193            Indx = GetSelectedAtoms()
2194        for i,atom in enumerate(drawAtoms):
2195            x,y,z = atom[cx:cx+3]
2196            X,Y,Z = gluProject(x,y,z,Model,Proj,View)
2197            XY = [int(X),int(View[3]-Y)]
2198            if np.allclose(xy,XY,atol=10) and Z < Zmax:
2199                Zmax = Z
2200                try:
2201                    Indx.remove(i)
2202                    ClearSelectedAtoms()
2203                    for id in Indx:
2204                        SetSelectedAtoms(id,Add)
2205                except:
2206                    SetSelectedAtoms(i,Add)
2207                   
2208    def OnMouseDown(event):
2209        xy = event.GetPosition()
2210        if event.ShiftDown():
2211            if event.LeftIsDown():
2212                GetTruePosition(xy)
2213            elif event.RightIsDown():
2214                GetTruePosition(xy,True)
2215        else:
2216            drawingData['Rotation'][3] = list(xy)
2217        Draw()
2218       
2219    def OnMouseMove(event):
2220        if event.ShiftDown():
2221            return       
2222        newxy = event.GetPosition()
2223        page = getSelection()
2224        if event.ControlDown() and drawingData['showABC']:
2225            if event.LeftIsDown():
2226                SetTestRot(newxy)
2227            elif event.RightIsDown():
2228                SetTestPos(newxy)
2229            elif event.MiddleIsDown():
2230                SetTestRotZ(newxy)
2231            x,y,z = drawingData['testPos'][0]
2232            G2frame.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
2233            Draw()
2234                               
2235        if event.Dragging() and not event.ControlDown():
2236            if event.LeftIsDown():
2237                SetRotation(newxy)
2238                angX,angY,angZ = drawingData['Rotation'][:3]
2239                G2frame.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f ,%.2f'%(angX,angY,angZ),1)
2240            elif event.RightIsDown():
2241                SetTranslation(newxy)
2242                Tx,Ty,Tz = drawingData['viewPoint'][0]
2243                G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
2244            elif event.MiddleIsDown():
2245                SetRotationZ(newxy)
2246                angX,angY,angZ = drawingData['Rotation'][:3]
2247                G2frame.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f, %.2f'%(angX,angY,angZ),1)
2248            Draw()
2249       
2250    def OnMouseWheel(event):
2251        if event.ShiftDown():
2252            return
2253        drawingData['cameraPos'] += event.GetWheelRotation()/24
2254        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
2255        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(drawingData['cameraPos']),1)
2256        page = getSelection()
2257        if page:
2258            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2259                panel = G2frame.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
2260                names = [child.GetName() for child in panel]
2261                panel[names.index('cameraPos')].SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
2262                panel[names.index('cameraSlider')].SetValue(drawingData['cameraPos'])
2263            Draw()
2264       
2265    def getSelection():
2266        try:
2267            return G2frame.dataDisplay.GetSelection()
2268        except AttributeError:
2269            G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!')
2270            return 0
2271           
2272    def SetViewPointText(VP):
2273        page = getSelection()
2274        if page:
2275            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2276                panel = G2frame.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
2277                names = [child.GetName() for child in panel]
2278                panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
2279           
2280    def ClearSelectedAtoms():
2281        page = getSelection()
2282        if page:
2283            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2284                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
2285            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2286                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
2287                   
2288    def SetSelectedAtoms(ind,Add=False):
2289        page = getSelection()
2290        if page:
2291            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2292                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)      #this is the Atoms grid in Draw Atoms
2293            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2294                Id = drawAtoms[ind][-2]
2295                for i,atom in enumerate(atomData):
2296                    if atom[-1] == Id:
2297                        G2frame.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
2298                 
2299    def GetSelectedAtoms():
2300        page = getSelection()
2301        Ind = []
2302        if page:
2303            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2304                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
2305            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2306                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
2307        return Ind
2308                                       
2309    def OnKey(event):           #on key UP!!
2310        keyCode = event.GetKeyCode()
2311        if keyCode > 255:
2312            keyCode = 0
2313        key,xyz = chr(keyCode),event.GetPosition()
2314        indx = drawingData['selectedAtoms']
2315        if key in ['c','C']:
2316            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
2317            drawingData['testPos'] = [[-.1,-.1,-.1],[0.0,0.0,0.0],[0,0]]
2318            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
2319            SetViewPointText(drawingData['viewPoint'][0])
2320        elif key in ['n','N']:
2321            drawAtoms = drawingData['Atoms']
2322            pI = drawingData['viewPoint'][1]
2323            if indx:
2324                pI[0] = indx[pI[1]]
2325                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
2326                pI[1] += 1
2327                if pI[1] >= len(indx):
2328                    pI[1] = 0
2329            else:
2330                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
2331                pI[0] += 1
2332                if pI[0] >= len(drawAtoms):
2333                    pI[0] = 0
2334            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
2335            SetViewPointText(drawingData['viewPoint'][0])
2336            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
2337               
2338        elif key in ['p','P']:
2339            drawAtoms = drawingData['Atoms']
2340            pI = drawingData['viewPoint'][1]
2341            if indx:
2342                pI[0] = indx[pI[1]]
2343                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
2344                pI[1] -= 1
2345                if pI[1] < 0:
2346                    pI[1] = len(indx)-1
2347            else:
2348                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
2349                pI[0] -= 1
2350                if pI[0] < 0:
2351                    pI[0] = len(drawAtoms)-1
2352            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
2353            SetViewPointText(drawingData['viewPoint'][0])           
2354            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
2355        Draw()
2356           
2357    def SetBackground():
2358        R,G,B,A = Page.camera['backColor']
2359        glClearColor(R,G,B,A)
2360        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
2361       
2362    def SetLights():
2363        glEnable(GL_DEPTH_TEST)
2364        glShadeModel(GL_SMOOTH)
2365        glEnable(GL_LIGHTING)
2366        glEnable(GL_LIGHT0)
2367        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
2368        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
2369        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
2370       
2371    def SetTranslation(newxy):
2372        Tx,Ty,Tz = drawingData['viewPoint'][0]
2373        anglex,angley,anglez,oldxy = drawingData['Rotation']
2374        Rx = G2lat.rotdMat(anglex,0)
2375        Ry = G2lat.rotdMat(angley,1)
2376        Rz = G2lat.rotdMat(anglez,2)
2377        dxy = list(newxy-oldxy)+[0,]
2378        dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
2379        Tx -= dxy[0]*0.01
2380        Ty += dxy[1]*0.01
2381        Tz -= dxy[2]*0.01
2382        drawingData['Rotation'][3] = list(newxy)
2383        drawingData['viewPoint'][0] =  Tx,Ty,Tz
2384        SetViewPointText([Tx,Ty,Tz])
2385       
2386    def SetTestPos(newxy):
2387        Tx,Ty,Tz = drawingData['testPos'][0]
2388        anglex,angley,anglez,oldxy = drawingData['Rotation']
2389        Rx = G2lat.rotdMat(anglex,0)
2390        Ry = G2lat.rotdMat(angley,1)
2391        Rz = G2lat.rotdMat(anglez,2)
2392        dxy = list(newxy-oldxy)+[0,]
2393        dxy = np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy)))
2394        Tx += dxy[0]*0.001
2395        Ty -= dxy[1]*0.001
2396        Tz += dxy[2]*0.001
2397        drawingData['Rotation'][3] = list(newxy)
2398        drawingData['testPos'][0] =  Tx,Ty,Tz
2399       
2400    def SetTestRot(newxy):
2401        Txyz = np.array(drawingData['testPos'][0])
2402        oldxy = drawingData['testPos'][2]
2403        Ax,Ay,Az = drawingData['testPos'][1]
2404        Vxyz = np.array(drawingData['viewPoint'][0])
2405        Dxyz = np.inner(Amat,Txyz-Vxyz)
2406        dxy = list(newxy-oldxy)+[0,]
2407        Ax += dxy[1]*0.01
2408        Ay += dxy[0]*0.01
2409        Rx = G2lat.rotdMat(Ax,0)
2410        Ry = G2lat.rotdMat(Ay,1)
2411        Dxyz = np.inner(Ry,np.inner(Rx,Dxyz))       
2412        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
2413        drawingData['testPos'][1] = [Ax,Ay,Az]
2414        drawingData['testPos'][2] = newxy
2415        drawingData['testPos'][0] = Dxyz
2416       
2417    def SetTestRotZ(newxy):
2418        Txyz = np.array(drawingData['testPos'][0])
2419        oldxy = drawingData['testPos'][2]
2420        Ax,Ay,Az = drawingData['testPos'][1]
2421        Vxyz = np.array(drawingData['viewPoint'][0])
2422        Dxyz = np.inner(Amat,Txyz-Vxyz)       
2423        dxy = list(newxy-oldxy)+[0,]
2424        Az += (dxy[0]+dxy[1])*.01
2425        Rz = G2lat.rotdMat(Az,2)
2426        Dxyz = np.inner(Rz,Dxyz)       
2427        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
2428        drawingData['testPos'][1] = [Ax,Ay,Az]
2429        drawingData['testPos'][2] = newxy
2430        drawingData['testPos'][0] = Dxyz
2431                             
2432    def SetRotation(newxy):       
2433        anglex,angley,anglez,oldxy = drawingData['Rotation']
2434        dxy = newxy-oldxy
2435        anglex += dxy[1]*.25
2436        angley += dxy[0]*.25
2437        oldxy = newxy
2438        drawingData['Rotation'] = [anglex,angley,anglez,list(oldxy)]
2439       
2440    def SetRotationZ(newxy):                       
2441        anglex,angley,anglez,oldxy = drawingData['Rotation']
2442        dxy = newxy-oldxy
2443        anglez += (dxy[0]+dxy[1])*.25
2444        oldxy = newxy
2445        drawingData['Rotation'] = [anglex,angley,anglez,list(oldxy)]
2446       
2447    def RenderBox():
2448        glEnable(GL_COLOR_MATERIAL)
2449        glLineWidth(2)
2450        glEnable(GL_BLEND)
2451        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
2452        glEnable(GL_LINE_SMOOTH)
2453        glBegin(GL_LINES)
2454        for line,color in zip(uEdges,uColors):
2455            glColor3ubv(color)
2456            glVertex3fv(line[0])
2457            glVertex3fv(line[1])
2458        glEnd()
2459        glColor4ubv([0,0,0,0])
2460        glDisable(GL_LINE_SMOOTH)
2461        glDisable(GL_BLEND)
2462        glDisable(GL_COLOR_MATERIAL)
2463       
2464    def RenderUnitVectors(x,y,z):
2465        xyz = np.array([x,y,z])
2466        glEnable(GL_COLOR_MATERIAL)
2467        glLineWidth(1)
2468        glPushMatrix()
2469        glTranslate(x,y,z)
2470        glScalef(1/cell[0],1/cell[1],1/cell[2])
2471        glBegin(GL_LINES)
2472        for line,color in zip(uEdges,uColors)[:3]:
2473            glColor3ubv(color)
2474            glVertex3fv(line[0])
2475            glVertex3fv(line[1])
2476        glEnd()
2477        glPopMatrix()
2478        glColor4ubv([0,0,0,0])
2479        glDisable(GL_COLOR_MATERIAL)
2480               
2481    def RenderSphere(x,y,z,radius,color):
2482        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2483        glPushMatrix()
2484        glTranslate(x,y,z)
2485        glMultMatrixf(B4mat.T)
2486        q = gluNewQuadric()
2487        gluSphere(q,radius,20,10)
2488        glPopMatrix()
2489       
2490    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
2491        s1,s2,s3 = E
2492        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2493        glPushMatrix()
2494        glTranslate(x,y,z)
2495        glMultMatrixf(B4mat.T)
2496        glMultMatrixf(R4.T)
2497        glEnable(GL_NORMALIZE)
2498        glScale(s1,s2,s3)
2499        q = gluNewQuadric()
2500        gluSphere(q,ellipseProb,20,10)
2501        glDisable(GL_NORMALIZE)
2502        glPopMatrix()
2503       
2504    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
2505        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2506        glPushMatrix()
2507        glTranslate(x,y,z)
2508        glMultMatrixf(B4mat.T)
2509        for bond in Bonds:
2510            glPushMatrix()
2511            Dx = np.inner(Amat,bond)
2512            Z = np.sqrt(np.sum(Dx**2))
2513            azm = atan2d(-Dx[1],-Dx[0])
2514            phi = acosd(Dx[2]/Z)
2515            glRotate(-azm,0,0,1)
2516            glRotate(phi,1,0,0)
2517            q = gluNewQuadric()
2518            gluCylinder(q,radius,radius,Z,slice,2)
2519            glPopMatrix()           
2520        glPopMatrix()
2521               
2522    def RenderLines(x,y,z,Bonds,color):
2523        xyz = np.array([x,y,z])
2524        glEnable(GL_COLOR_MATERIAL)
2525        glLineWidth(1)
2526        glColor3fv(color)
2527        glPushMatrix()
2528        glBegin(GL_LINES)
2529        for bond in Bonds:
2530            glVertex3fv(xyz)
2531            glVertex3fv(xyz+bond)
2532        glEnd()
2533        glColor4ubv([0,0,0,0])
2534        glPopMatrix()
2535        glDisable(GL_COLOR_MATERIAL)
2536       
2537    def RenderPolyhedra(x,y,z,Faces,color):
2538        glPushMatrix()
2539        glTranslate(x,y,z)
2540        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2541        glShadeModel(GL_SMOOTH)
2542        glMultMatrixf(B4mat.T)
2543        for face,norm in Faces:
2544            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
2545            glFrontFace(GL_CW)
2546            glNormal3fv(norm)
2547            glBegin(GL_TRIANGLES)
2548            for vert in face:
2549                glVertex3fv(vert)
2550            glEnd()
2551        glPopMatrix()
2552       
2553    def RenderBackbone(Backbone,BackboneColor,radius):
2554        glPushMatrix()
2555        glMultMatrixf(B4mat.T)
2556        glEnable(GL_COLOR_MATERIAL)
2557        glShadeModel(GL_SMOOTH)
2558        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
2559        glePolyCylinder(Backbone,BackboneColor,radius)
2560        glPopMatrix()       
2561        glDisable(GL_COLOR_MATERIAL)
2562       
2563    def RenderLabel(x,y,z,label,r):       
2564        glPushMatrix()
2565        glTranslate(x,y,z)
2566        glMultMatrixf(B4mat.T)
2567        glDisable(GL_LIGHTING)
2568        glColor3f(0,1.,0)
2569        glRasterPos3f(r,r,r)
2570        for c in list(label):
2571            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
2572        glEnable(GL_LIGHTING)
2573        glPopMatrix()
2574                           
2575    def Draw():
2576        Ind = GetSelectedAtoms()
2577        VS = np.array(Page.canvas.GetSize())
2578        aspect = float(VS[0])/float(VS[1])
2579        cPos = drawingData['cameraPos']
2580        Zclip = drawingData['Zclip']*cPos/200.
2581        anglex,angley,anglez = drawingData['Rotation'][:3]
2582        Tx,Ty,Tz = drawingData['viewPoint'][0]
2583        cx,ct,cs = drawingData['atomPtrs']
2584        bondR = drawingData['bondRadius']
2585        G,g = G2lat.cell2Gmat(cell)
2586        GS = G
2587        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
2588        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
2589        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
2590        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
2591       
2592        SetBackground()
2593        glInitNames()
2594        glPushName(0)
2595       
2596        glMatrixMode(GL_PROJECTION)
2597        glLoadIdentity()
2598        glViewport(0,0,VS[0],VS[1])
2599        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
2600        gluLookAt(0,0,cPos,0,0,0,0,1,0)
2601        SetLights()           
2602           
2603        glMatrixMode(GL_MODELVIEW)
2604        glLoadIdentity()
2605        glRotate(anglez,0,0,1)
2606        glRotate(anglex,cosd(anglez),-sind(anglez),0)
2607        glRotate(angley,sind(anglez),cosd(anglez),0)
2608        glMultMatrixf(A4mat.T)
2609        glTranslate(-Tx,-Ty,-Tz)
2610        if drawingData['unitCellBox']:
2611            RenderBox()
2612        if drawingData['showABC']:
2613            x,y,z = drawingData['testPos'][0]
2614#            if altDown:
2615#                G2frame.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
2616#            else:
2617#                G2frame.G2plotNB.status.SetStatusText('test point %.4f,%.4f,%.4f'%(x,y,z),1)           
2618            RenderUnitVectors(x,y,z)
2619        Backbone = []
2620        BackboneColor = []
2621        time0 = time.time()
2622        for iat,atom in enumerate(drawingData['Atoms']):
2623            x,y,z = atom[cx:cx+3]
2624            Bonds = atom[-2]
2625            Faces = atom[-1]
2626            try:
2627                atNum = generalData['AtomTypes'].index(atom[ct])
2628            except ValueError:
2629                atNum = -1
2630            CL = atom[cs+2]
2631            color = np.array(CL)/255.
2632            if iat in Ind:
2633                color = np.array(Gr)/255.
2634            radius = 0.5
2635            if atom[cs] != '':
2636                glLoadName(atom[-3])                   
2637            if 'balls' in atom[cs]:
2638                vdwScale = drawingData['vdwScale']
2639                ballScale = drawingData['ballScale']
2640                if atNum < 0:
2641                    radius = 0.2
2642                elif 'H' == atom[ct]:
2643                    if drawingData['showHydrogen']:
2644                        if 'vdW' in atom[cs] and atNum >= 0:
2645                            radius = vdwScale*generalData['vdWRadii'][atNum]
2646                        else:
2647                            radius = ballScale*drawingData['sizeH']
2648                    else:
2649                        radius = 0.0
2650                else:
2651                    if 'vdW' in atom[cs]:
2652                        radius = vdwScale*generalData['vdWRadii'][atNum]
2653                    else:
2654                        radius = ballScale*generalData['BondRadii'][atNum]
2655                RenderSphere(x,y,z,radius,color)
2656                if 'sticks' in atom[cs]:
2657                    RenderBonds(x,y,z,Bonds,bondR,color)
2658            elif 'ellipsoids' in atom[cs]:
2659                RenderBonds(x,y,z,Bonds,bondR,color)
2660                if atom[cs+3] == 'A':                   
2661                    Uij = atom[cs+5:cs+11]
2662                    U = np.multiply(G2spc.Uij2U(Uij),GS)
2663                    U = np.inner(Amat,np.inner(U,Amat).T)
2664                    E,R = nl.eigh(U)
2665                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2666                    E = np.sqrt(E)
2667                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
2668                        pass
2669                    else:
2670                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,color)                   
2671                else:
2672                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
2673                        pass
2674                    else:
2675                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
2676                        RenderSphere(x,y,z,radius,color)
2677            elif 'lines' in atom[cs]:
2678                radius = 0.1
2679                RenderLines(x,y,z,Bonds,color)
2680#                RenderBonds(x,y,z,Bonds,0.05,color,6)
2681            elif atom[cs] == 'sticks':
2682                radius = 0.1
2683                RenderBonds(x,y,z,Bonds,bondR,color)
2684            elif atom[cs] == 'polyhedra':
2685                RenderPolyhedra(x,y,z,Faces,color)
2686            elif atom[cs] == 'backbone':
2687                if atom[ct-1].split()[0] in ['C','N']:
2688                    Backbone.append(list(np.inner(Amat,np.array([x,y,z]))))
2689                    BackboneColor.append(list(color))
2690                   
2691            if atom[cs+1] == 'type':
2692                RenderLabel(x,y,z,atom[ct],radius)
2693            elif atom[cs+1] == 'name':
2694                RenderLabel(x,y,z,atom[ct-1],radius)
2695            elif atom[cs+1] == 'number':
2696                RenderLabel(x,y,z,str(iat),radius)
2697            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
2698                RenderLabel(x,y,z,atom[ct-4],radius)
2699            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
2700                RenderLabel(x,y,z,atom[ct-3],radius)
2701            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
2702                RenderLabel(x,y,z,atom[ct-2],radius)
2703        if Backbone:
2704            RenderBackbone(Backbone,BackboneColor,bondR)
2705#        print time.time()-time0
2706        Page.canvas.SwapBuffers()
2707       
2708    def OnSize(event):
2709        Draw()
2710       
2711    def OnFocus(event):
2712        Draw()
2713       
2714    try:
2715        plotNum = G2frame.G2plotNB.plotList.index(generalData['Name'])
2716        Page = G2frame.G2plotNB.nb.GetPage(plotNum)       
2717    except ValueError:
2718        Plot = G2frame.G2plotNB.addOgl(generalData['Name'])
2719        plotNum = G2frame.G2plotNB.plotList.index(generalData['Name'])
2720        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2721        Page.views = False
2722        view = False
2723        altDown = False
2724    Page.SetFocus()
2725    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
2726        choices=(' save as:','jpeg','tiff','bmp'))
2727    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
2728    cb.SetValue(' save as:')
2729    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
2730    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
2731    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
2732    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
2733    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
2734    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
2735    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
2736    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
2737    Page.camera['position'] = drawingData['cameraPos']
2738    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
2739    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
2740    Page.canvas.SetCurrent()
2741    Draw()
2742       
Note: See TracBrowser for help on using the repository browser.