source: trunk/GSASIIplot.py @ 1244

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

fix geometric correction in integrate - too many 1/cos(2-theta)
plot of size distribution from SASD
MaxEnt? size distribution in operation (some tuning/errors)

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