source: trunk/GSASIIplot.py @ 1199

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

1st import of small angle data, create SASD data from image.
Work on proper plotting of SASD data
work on image GUI stuff for small angle data

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