source: trunk/GSASIIplot.py @ 1559

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

use deepcopy for various copy operations in G2ddataGUI
use complementary colors to background for cell edges
trap superlattice in cubics
more fixes to cell indexing routines from possible use of superlattice
add HStrainVals to G2spc
use it (maybe) in G2strIO, Math
make arguments for GetReflPos? & GetReflPosDeriv? the same (use A not G)
add GetDij? to G2strMath
topas xye file comments start with "'" sometimes

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