source: trunk/GSASIIplot.py @ 1201

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

remove bind to leave window in ValidatedTxtCtrl?; interfered with plotting
implement plot of intensity*scale for small angle data
remove some dead code in G2pwdGUI

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