source: trunk/GSASIIplot.py @ 1466

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

make extinction math use FcTsq? not Fcsq
fix 3D HKL point lighting
fix reflection listing for CW data

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