source: trunk/GSASIIplot.py @ 1212

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

fix scaling of log images

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