source: trunk/GSASIIplot.py @ 953

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

fixes to RB refinement & MC/SA

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