source: trunk/GSASIIplot.py @ 792

Last change on this file since 792 was 792, checked in by vondreele, 10 years ago

change instrument parameters to dict from lists
chase down effects - got them all?
start on TOF; read TimeMap? style data - not complete

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