source: trunk/GSASIIplot.py @ 835

Last change on this file since 835 was 835, checked in by vondreele, 9 years ago

more work on rigid bodies
fixes for crash in selecting wavelengths

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