source: trunk/GSASIIplot.py @ 1401

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

trap keystrokes for powder plots when data window not on powder info.

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