source: trunk/GSASIIplot.py @ 1054

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

fix for RB plotting

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