source: trunk/GSASIIplot.py @ 1250

Last change on this file since 1250 was 1250, checked in by vondreele, 8 years ago

fix cursor position in SAS size plot
sky default = -3
fix background handling in SASD

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