source: trunk/GSASIIplot.py @ 1386

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

fix stress strain data references

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