source: trunk/GSASIIplot.py @ 1341

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

SASD smearing now works
some work with SASD plots

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