source: trunk/GSASIIplot.py @ 797

Last change on this file since 797 was 797, checked in by vondreele, 10 years ago

some fixes to TOF stuff
some fixes to SXC stuff

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