source: trunk/GSASIIplot.py @ 795

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

implement TOF input, peak search & fitting
auto peak search
convert instrument parms to dictionary
add charge flip on max rho

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