source: trunk/GSASIIplot.py @ 855

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

RB input mostly finished; now has thermal motion models

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