source: trunk/GSASIIplot.py @ 984

Last change on this file since 984 was 984, checked in by vondreele, 10 years ago

update Win2.7 g77 binaries
cov matrix now not modified by lam
no label plotting at all - fails
set MC/SA page after mc/sa run
fix 64 bit bug in switching to map peaks page

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