source: trunk/GSASIIplot.py @ 1463

Last change on this file since 1463 was 1463, checked in by vondreele, 11 years ago

fix motion position in q & d plots

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