source: trunk/GSASIIplot.py @ 808

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

Add Babinet modification to form factors - works OK
Use atom Ids for restraints
add 'all' in atom constraint editing routines
New add routines for bond restraints (angle add not done yet)
FindAtomIndexByIDs, FillAtomLookUp?, GetAtomsByID, & GetAtomItemsById? all now in GSASIImath.py
removed from GSASIIphsGUI.py
change phase tick mark colors
get GSASIIstruct.py to recognize protein atoms & do a protein calc.

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