source: trunk/GSASIIplot.py @ 303

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

Allow reading of xye (Topas style) files
Begin implementation of spherical harmonics texture
Refactor indexing; replace cell refinement from peak positions

  • Property svn:keywords set to Date Author Revision URL Id
File size: 99.3 KB
Line 
1#GSASII plotting routines
2########### SVN repository information ###################
3# $Date: 2011-06-20 19:37:07 +0000 (Mon, 20 Jun 2011) $
4# $Author: vondreele $
5# $Revision: 303 $
6# $URL: trunk/GSASIIplot.py $
7# $Id: GSASIIplot.py 303 2011-06-20 19:37:07Z vondreele $
8########### SVN repository information ###################
9import math
10import time
11import copy
12import os.path
13import numpy as np
14import numpy.linalg as nl
15import wx
16import wx.aui
17import wx.glcanvas
18import matplotlib as mpl
19import mpl_toolkits.mplot3d.axes3d as mp3d
20import GSASIIpath
21import GSASIIgrid as G2gd
22import GSASIIimage as G2img
23import GSASIIpwd as G2pwd
24import GSASIIIO as G2IO
25import GSASIIpwdGUI as G2pdG
26import GSASIIimgGUI as G2imG
27import GSASIIphsGUI as G2phG
28import GSASIIlattice as G2lat
29import GSASIIspc as G2spc
30import pytexture as ptx
31from  OpenGL.GL import *
32from OpenGL.GLU import *
33from OpenGL.GLUT import *
34from OpenGL.GLE import *
35from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as Canvas
36from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
37
38# useful degree trig functions
39sind = lambda x: math.sin(x*math.pi/180.)
40cosd = lambda x: math.cos(x*math.pi/180.)
41tand = lambda x: math.tan(x*math.pi/180.)
42asind = lambda x: 180.*math.asin(x)/math.pi
43acosd = lambda x: 180.*math.acos(x)/math.pi
44atan2d = lambda x,y: 180.*math.atan2(y,x)/math.pi
45atand = lambda x: 180.*math.atan(x)/math.pi
46# numpy versions
47npsind = lambda x: np.sin(x*np.pi/180.)
48npcosd = lambda x: np.cos(x*np.pi/180.)
49npacosd = lambda x: 180.*np.arccos(x)/np.pi
50   
51class G2PlotMpl(wx.Panel):   
52    def __init__(self,parent,id=-1,dpi=None,**kwargs):
53        wx.Panel.__init__(self,parent,id=id,**kwargs)
54        mpl.rcParams['legend.fontsize'] = 8
55        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(5,7))
56        self.canvas = Canvas(self,-1,self.figure)
57        self.toolbar = Toolbar(self.canvas)
58
59        self.toolbar.Realize()
60       
61        sizer=wx.BoxSizer(wx.VERTICAL)
62        sizer.Add(self.canvas,1,wx.EXPAND)
63        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
64        self.SetSizer(sizer)
65       
66class G2PlotOgl(wx.Panel):
67    def __init__(self,parent,id=-1,dpi=None,**kwargs):
68        self.figure = wx.Panel.__init__(self,parent,id=id,**kwargs)
69        self.canvas = wx.glcanvas.GLCanvas(self,-1,**kwargs)
70        self.camera = {}
71        sizer=wx.BoxSizer(wx.VERTICAL)
72        sizer.Add(self.canvas,1,wx.EXPAND)
73        self.SetSizer(sizer)
74       
75class G2Plot3D(wx.Panel):
76    def __init__(self,parent,id=-1,dpi=None,**kwargs):
77        wx.Panel.__init__(self,parent,id=id,**kwargs)
78        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(6,6))
79        self.canvas = Canvas(self,-1,self.figure)
80        self.toolbar = Toolbar(self.canvas)
81
82        self.toolbar.Realize()
83       
84        sizer=wx.BoxSizer(wx.VERTICAL)
85        sizer.Add(self.canvas,1,wx.EXPAND)
86        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
87        self.SetSizer(sizer)
88                             
89class G2PlotNoteBook(wx.Panel):
90    def __init__(self,parent,id=-1):
91        wx.Panel.__init__(self,parent,id=id)
92        #so one can't delete a plot page!!
93        self.nb = wx.aui.AuiNotebook(self, \
94            style=wx.aui.AUI_NB_DEFAULT_STYLE ^ wx.aui.AUI_NB_CLOSE_ON_ACTIVE_TAB)
95        sizer = wx.BoxSizer()
96        sizer.Add(self.nb,1,wx.EXPAND)
97        self.SetSizer(sizer)
98        self.status = parent.CreateStatusBar()
99        self.status.SetFieldsCount(2)
100        self.status.SetStatusWidths([150,-1])
101        self.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, self.OnPageChanged)
102       
103        self.plotList = []
104           
105    def addMpl(self,name=""):
106        page = G2PlotMpl(self.nb)
107        self.nb.AddPage(page,name)
108       
109        self.plotList.append(name)
110       
111        return page.figure
112       
113    def add3D(self,name=""):
114        page = G2Plot3D(self.nb)
115        self.nb.AddPage(page,name)
116       
117        self.plotList.append(name)
118       
119        return page.figure
120       
121    def addOgl(self,name=""):
122        page = G2PlotOgl(self.nb)
123        self.nb.AddPage(page,name)
124       
125        self.plotList.append(name)
126       
127        return page.figure
128       
129    def Delete(self,name):
130        try:
131            item = self.plotList.index(name)
132            del self.plotList[item]
133            self.nb.DeletePage(item)
134        except ValueError:          #no plot of this name - do nothing
135            return     
136               
137    def clear(self):
138        while self.nb.GetPageCount():
139            self.nb.DeletePage(0)
140        self.plotList = []
141        self.status.DestroyChildren()
142       
143    def Rename(self,oldName,newName):
144        try:
145            item = self.plotList.index(oldName)
146            self.plotList[item] = newName
147            self.nb.SetPageText(item,newName)
148        except ValueError:          #no plot of this name - do nothing
149            return     
150       
151    def OnPageChanged(self,event):       
152        if self.plotList:
153            self.status.SetStatusText('Better to select this from GSAS-II data tree',1)
154        self.status.DestroyChildren()                           #get rid of special stuff on status bar
155       
156def PlotSngl(self,newPlot=False):
157    '''Single crystal structure factor plotting package - displays zone of reflections as rings proportional
158        to F, F**2, etc. as requested
159    '''
160    from matplotlib.patches import Circle
161    global HKL,HKLF
162
163    def OnSCMotion(event):
164        xpos = event.xdata
165        if xpos:
166            xpos = round(xpos)                                        #avoid out of frame mouse position
167            ypos = round(event.ydata)
168            zpos = Data['Layer']
169            if '100' in Data['Zone']:
170                HKLtxt = '(%3d,%3d,%3d)'%(zpos,xpos,ypos)
171            elif '010' in Data['Zone']:
172                HKLtxt = '(%3d,%3d,%3d)'%(xpos,zpos,ypos)
173            elif '001' in Data['Zone']:
174                HKLtxt = '(%3d,%3d,%3d)'%(xpos,ypos,zpos)
175            Page.canvas.SetToolTipString(HKLtxt)
176            self.G2plotNB.status.SetFields(['HKL = '+HKLtxt,''])
177               
178    def OnSCPick(event):
179        zpos = Data['Layer']
180        pos = event.artist.center
181        if '100' in Data['Zone']:
182            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(zpos,pos[0],pos[1]))
183            hkl = np.array([zpos,pos[0],pos[1]])
184        elif '010' in Data['Zone']:
185            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],zpos,pos[1]))
186            hkl = np.array([pos[0],zpos,pos[1]])
187        elif '001' in Data['Zone']:
188            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],pos[1],zpos))
189            hkl = np.array([pos[0],pos[1],zpos])
190        h,k,l = hkl
191        hklf = HKLF[np.where(np.all(HKL-hkl == [0,0,0],axis=1))]
192        if len(hklf):
193            Fosq,sig,Fcsq = hklf[0]
194            HKLtxt = '(%3d,%3d,%3d %.2f %.3f %.2f %.2f)'%(h,k,l,Fosq,sig,Fcsq,(Fosq-Fcsq)/(scale*sig))
195            self.G2plotNB.status.SetFields(['','HKL, Fosq, sig, Fcsq, delFsq/sig = '+HKLtxt])
196                                 
197    def OnSCKeyPress(event):
198        print event.key
199
200    try:
201        plotNum = self.G2plotNB.plotList.index('Structure Factors')
202        Page = self.G2plotNB.nb.GetPage(plotNum)
203        if not newPlot:
204            Plot = Page.figure.gca()          #get previous powder plot & get limits
205            xylim = Plot.get_xlim(),Plot.get_ylim()
206        Page.figure.clf()
207        Plot = Page.figure.gca()          #get a fresh plot after clf()
208    except ValueError:
209        Plot = self.G2plotNB.addMpl('Structure Factors').gca()
210        plotNum = self.G2plotNB.plotList.index('Structure Factors')
211        Page = self.G2plotNB.nb.GetPage(plotNum)
212#        Page.canvas.mpl_connect('key_press_event', OnSCKeyPress)
213        Page.canvas.mpl_connect('pick_event', OnSCPick)
214        Page.canvas.mpl_connect('motion_notify_event', OnSCMotion)
215    Page.SetFocus()
216   
217    Plot.set_aspect(aspect='equal')
218    HKLref = self.PatternTree.GetItemPyData(self.Sngl)
219    Data = self.PatternTree.GetItemPyData( \
220        G2gd.GetPatternTreeItemId(self,self.Sngl, 'HKL Plot Controls'))
221    Type = Data['Type']           
222    scale = Data['Scale']
223    HKLmax = Data['HKLmax']
224    HKLmin = Data['HKLmin']
225    FosqMax = Data['FoMax']
226    FoMax = math.sqrt(FosqMax)
227    ifFc = Data['ifFc']
228    xlabel = ['k, h=','h, k=','h, l=']
229    ylabel = ['l','l','k']
230    zones = ['100','010','001']
231    pzone = [[1,2],[0,2],[0,1]]
232    izone = zones.index(Data['Zone'])
233    Plot.set_title(self.PatternTree.GetItemText(self.Sngl)[5:])
234    HKL = []
235    HKLF = []
236    for H,Fosq,sig,Fcsq,x,x,x in HKLref:
237        HKL.append(H)
238        HKLF.append([Fosq,sig,Fcsq])
239        if H[izone] == Data['Layer']:
240            B = 0
241            if Type == 'Fosq':
242                A = scale*Fosq/FosqMax
243                B = scale*Fcsq/FosqMax
244                C = abs(A-B)
245            elif Type == 'Fo':
246                A = scale*math.sqrt(max(0,Fosq))/FoMax
247                B = scale*math.sqrt(max(0,Fcsq))/FoMax
248                C = abs(A-B)
249            elif Type == '|DFsq|/sig':
250                A = abs(Fosq-Fcsq)/(scale*sig)
251            elif Type == '|DFsq|>sig':
252                A = abs(Fosq-Fcsq)/(scale*sig)
253                if A < 1.0: A = 0                   
254            elif Type == '|DFsq|>3sig':
255                A = abs(Fosq-Fcsq)/(scale*sig)
256                if A < 3.0: A = 0                   
257            xy = (H[pzone[izone][0]],H[pzone[izone][1]])
258            if A > 0.0:
259                Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w',picker=3))
260            if B:
261                Plot.add_artist(Circle(xy,radius=B,ec='b',fc='w'))
262                radius = C
263                if radius > 0:
264                    if A > B:
265                        Plot.add_artist(Circle(xy,radius=radius,ec='g',fc='g'))
266                    else:                   
267                        Plot.add_artist(Circle(xy,radius=radius,ec='r',fc='r'))
268    HKL = np.array(HKL,dtype=np.int)
269    HKLF = np.array(HKLF)
270    Plot.set_xlabel(xlabel[izone]+str(Data['Layer']),fontsize=12)
271    Plot.set_ylabel(ylabel[izone],fontsize=12)
272    Plot.set_xlim((HKLmin[pzone[izone][0]],HKLmax[pzone[izone][0]]))
273    Plot.set_ylim((HKLmin[pzone[izone][1]],HKLmax[pzone[izone][1]]))
274    if not newPlot:
275        Page.toolbar.push_current()
276        Plot.set_xlim(xylim[0])
277        Plot.set_ylim(xylim[1])
278        xylim = []
279        Page.toolbar.push_current()
280        Page.toolbar.draw()
281    else:
282        Page.canvas.draw()
283       
284def PlotPatterns(self,newPlot=False):
285    '''Powder pattern plotting package - displays single or multiple powder patterns as intensity vs
286    2-theta or q (future TOF). Can display multiple patterns as "waterfall plots" or contour plots. Log I
287    plotting available.
288    '''
289    global HKL
290   
291    def OnPick(event):
292        if self.itemPicked is not None: return
293        PatternId = self.PatternId
294        try:
295            Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,self.PatternId, 'Instrument Parameters'))[1::2]
296        except TypeError:
297            return
298        Parms = dict(zip(Names,Values))
299        try:
300            wave = Parms['Lam']
301        except KeyError:
302            wave = Parms['Lam1']
303        PickId = self.PickId
304        pick = event.artist
305        mouse = event.mouseevent       
306        xpos = pick.get_xdata()
307        ypos = pick.get_ydata()
308        ind = event.ind
309        xy = list(zip(np.take(xpos,ind),np.take(ypos,ind))[0])
310        if self.PatternTree.GetItemText(PickId) == 'Peak List':
311            if ind.all() != [0]:                                    #picked a data point
312                if 'C' in Parms['Type']:                            #CW data - TOF later in an elif
313                    ins = [Parms[x] for x in ['U','V','W','X','Y']]
314                    if self.qPlot:                              #qplot - convert back to 2-theta
315                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
316                    sig = ins[0]*tand(xy[0]/2.0)**2+ins[1]*tand(xy[0]/2.0)+ins[2]
317                    gam = ins[3]/cosd(xy[0]/2.0)+ins[4]*tand(xy[0]/2.0)           
318                    data = self.PatternTree.GetItemPyData(self.PickId)
319                    XY = [xy[0],0, xy[1],1, sig,0, gam,0]       #default refine intensity 1st
320                data.append(XY)
321                G2pdG.UpdatePeakGrid(self,data)
322                PlotPatterns(self)
323            else:                                                   #picked a peak list line
324                self.itemPicked = pick
325        elif self.PatternTree.GetItemText(PickId) == 'Limits':
326            if ind.all() != [0]:                                    #picked a data point
327                LimitId = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
328                data = self.PatternTree.GetItemPyData(LimitId)
329                if 'C' in Parms['Type']:                            #CW data - TOF later in an elif
330                    if self.qPlot:                              #qplot - convert back to 2-theta
331                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
332                if mouse.button==1:
333                    data[1][0] = min(xy[0],data[1][1])
334                if mouse.button==3:
335                    data[1][1] = max(xy[0],data[1][0])
336                self.PatternTree.SetItemPyData(LimitId,data)
337                G2pdG.UpdateLimitsGrid(self,data)
338                PlotPatterns(self)
339            else:                                                   #picked a limit line
340                self.itemPicked = pick
341       
342    def OnPlotKeyPress(event):
343        newPlot = False
344        if event.key == 'w':
345            if self.Weight:
346                self.Weight = False
347            else:
348                self.Weight = True
349            print 'plot weighting:',self.Weight
350        elif event.key == 'n':
351            if self.Contour:
352                pass
353            else:
354                if self.logPlot:
355                    self.logPlot = False
356                else:
357                    self.Offset[0] = 0
358                    self.logPlot = True
359        elif event.key == 'u':
360            if self.Contour:
361                self.Cmax = min(1.0,self.Cmax*1.2)
362            elif self.logPlot:
363                pass
364            elif self.Offset[0] < 100.:
365                self.Offset[0] += 1.
366        elif event.key == 'd':
367            if self.Contour:
368                self.Cmax = max(0.0,self.Cmax*0.8)
369            elif self.logPlot:
370                pass
371            elif self.Offset[0] > 0.:
372                self.Offset[0] -= 1.
373        elif event.key == 'l':
374            self.Offset[1] -= 1.
375        elif event.key == 'r':
376            self.Offset[1] += 1.
377        elif event.key == 'o':
378            self.Offset = [0,0]
379        elif event.key == 'c':
380            newPlot = True
381            if self.Contour:
382                self.Contour = False
383            else:
384                self.Contour = True
385                self.SinglePlot = False
386                self.Offset = [0.,0.]
387        elif event.key == 'q':
388            newPlot = True
389            if self.qPlot:
390                self.qPlot = False
391            else:
392                self.qPlot = True
393        elif event.key == 's':
394            if self.Contour:
395                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
396                choice.sort()
397                dlg = wx.SingleChoiceDialog(self,'Select','Color scheme',choice)
398                if dlg.ShowModal() == wx.ID_OK:
399                    sel = dlg.GetSelection()
400                    self.ContourColor = choice[sel]
401                else:
402                    self.ContourColor = 'Paired'
403                dlg.Destroy()
404            else:               
405                if self.SinglePlot:
406                    self.SinglePlot = False
407                else:
408                    self.SinglePlot = True
409        elif event.key == '+':
410            if self.PickId:
411                self.PickId = False
412        elif event.key == 'i':                  #for smoothing contour plot
413            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
414               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
415               'mitchell','sinc','lanczos']
416            dlg = wx.SingleChoiceDialog(self,'Select','Interpolation',choice)
417            if dlg.ShowModal() == wx.ID_OK:
418                sel = dlg.GetSelection()
419                self.Interpolate = choice[sel]
420            else:
421                self.Interpolate = 'nearest'
422            dlg.Destroy()
423           
424        PlotPatterns(self,newPlot=newPlot)
425       
426    def OnKeyBox(event):
427        if self.G2plotNB.nb.GetSelection() == self.G2plotNB.plotList.index('Powder Patterns'):
428            event.key = cb.GetValue()[0]
429            cb.SetValue(' key press')
430            OnPlotKeyPress(event)
431                       
432    def OnMotion(event):
433        xpos = event.xdata
434        if xpos:                                        #avoid out of frame mouse position
435            ypos = event.ydata
436            Page.canvas.SetCursor(wx.CROSS_CURSOR)
437            try:
438                Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,self.PatternId, 'Instrument Parameters'))[1::2]
439                Parms = dict(zip(Names,Values))
440                try:
441                    wave = Parms['Lam']
442                except KeyError:
443                    wave = Parms['Lam1']
444                if self.qPlot:
445                    xpos = 2.0*asind(xpos*wave/(4*math.pi))
446                dsp = 0.0
447                if abs(xpos) > 0.:                  #avoid possible singularity at beam center
448                    dsp = wave/(2.*sind(abs(xpos)/2.0))
449                if self.Contour:
450                    self.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f pattern ID =%5d'%(xpos,dsp,int(ypos)),1)
451                else:
452                    self.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f Intensity =%9.1f'%(xpos,dsp,ypos),1)
453                if self.itemPicked:
454                    Page.canvas.SetToolTipString('%9.3f'%(xpos))
455                if self.PickId and self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List']:
456                    found = []
457                    if len(HKL):
458                        view = Page.toolbar._views.forward()[0][:2]
459                        wid = view[1]-view[0]
460                        found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
461                    if len(found):
462                        h,k,l = found[0][:3] 
463                        Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
464                    else:
465                        Page.canvas.SetToolTipString('')
466
467            except TypeError:
468                self.G2plotNB.status.SetStatusText('Select PWDR powder pattern first',1)
469                                                   
470    def OnRelease(event):
471        if self.itemPicked is None: return
472        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,self.PatternId, 'Instrument Parameters'))[1::2]
473        Parms = dict(zip(Names,Values))
474        try:
475            wave = Parms['Lam']
476        except KeyError:
477            wave = Parms['Lam1']
478        xpos = event.xdata
479        if xpos:                                        #avoid out of frame mouse position
480            lines = []
481            for line in self.Lines: lines.append(line.get_xdata()[0])
482            lineNo = lines.index(self.itemPicked.get_xdata()[0])
483            if  lineNo in [0,1]:
484                LimitId = G2gd.GetPatternTreeItemId(self,self.PatternId, 'Limits')
485                data = self.PatternTree.GetItemPyData(LimitId)
486#                print 'limits',xpos
487                if self.qPlot:
488                    data[1][lineNo] = 2.0*asind(wave*xpos/(4*math.pi))
489                else:
490                    data[1][lineNo] = xpos
491                self.PatternTree.SetItemPyData(LimitId,data)
492                if self.PatternTree.GetItemText(self.PickId) == 'Limits':
493                    G2pdG.UpdateLimitsGrid(self,data)
494            else:
495                PeakId = G2gd.GetPatternTreeItemId(self,self.PatternId, 'Peak List')
496                data = self.PatternTree.GetItemPyData(PeakId)
497#                print 'peaks',xpos
498                if event.button == 3:
499                    del data[lineNo-2]
500                else:
501                    if self.qPlot:
502                        data[lineNo-2][0] = 2.0*asind(wave*xpos/(4*math.pi))
503                    else:
504                        data[lineNo-2][0] = xpos
505                self.PatternTree.SetItemPyData(PeakId,data)
506                G2pdG.UpdatePeakGrid(self,data)
507        PlotPatterns(self)
508        self.itemPicked = None   
509
510    xylim = []
511    try:
512        plotNum = self.G2plotNB.plotList.index('Powder Patterns')
513        Page = self.G2plotNB.nb.GetPage(plotNum)
514        if not newPlot:
515            Plot = Page.figure.gca()          #get previous powder plot & get limits
516            xylim = Plot.get_xlim(),Plot.get_ylim()
517        Page.figure.clf()
518        Plot = Page.figure.gca()          #get a fresh plot after clf()
519    except ValueError:
520        newPlot = True
521        self.Cmax = 1.0
522        Plot = self.G2plotNB.addMpl('Powder Patterns').gca()
523        plotNum = self.G2plotNB.plotList.index('Powder Patterns')
524        Page = self.G2plotNB.nb.GetPage(plotNum)
525        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
526        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
527        Page.canvas.mpl_connect('pick_event', OnPick)
528        Page.canvas.mpl_connect('button_release_event', OnRelease)
529    Page.SetFocus()
530    self.G2plotNB.status.DestroyChildren()
531    if self.Contour:
532        Choice = (' key press','d: lower contour max','u: raise contour max',
533            'i: interpolation method','s: color scheme','c: contour off')
534    else:
535        if self.logPlot:
536            Choice = (' key press','n: log(I) off','l: offset left','r: offset right',
537                'c: contour on','q: toggle q plot','s: toggle single plot','+: no selection')
538        else:
539            Choice = (' key press','l: offset left','r: offset right','d: offset down',
540                'u: offset up','o: reset offset','n: log(I) on','c: contour on',
541                'q: toggle q plot','s: toggle single plot','+: no selection')
542    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
543        choices=Choice)
544    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
545    cb.SetValue(' key press')
546   
547    PickId = self.PickId
548    PatternId = self.PatternId
549    colors=['b','g','r','c','m','k']
550    Lines = []
551    if self.SinglePlot:
552        Pattern = self.PatternTree.GetItemPyData(PatternId)
553        Pattern.append(self.PatternTree.GetItemText(PatternId))
554        PlotList = [Pattern,]
555        ParmList = [self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,
556            self.PatternId, 'Instrument Parameters'))[1],]
557    else:       
558        PlotList = []
559        ParmList = []
560        item, cookie = self.PatternTree.GetFirstChild(self.root)
561        while item:
562            if 'PWDR' in self.PatternTree.GetItemText(item):
563                Pattern = self.PatternTree.GetItemPyData(item)
564                if len(Pattern) < 3:                    # put name on end if needed
565                    Pattern.append(self.PatternTree.GetItemText(item))
566                PlotList.append(Pattern)
567                ParmList.append(self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,
568                    item,'Instrument Parameters'))[1])
569            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)               
570    Ymax = 1.0
571    lenX = 0
572    HKL = np.array(self.HKL)
573    for Pattern in PlotList:
574        xye = Pattern[1]
575        Ymax = max(Ymax,max(xye[1]))
576    offset = self.Offset[0]*Ymax/100.0
577    Title = 'Powder Patterns: '+os.path.split(self.GSASprojectfile)[1]
578    if self.logPlot:
579        Title = 'log('+Title+')'
580    Plot.set_title(Title)
581    if self.qPlot:
582        Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
583    else:       
584        Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
585    Plot.set_ylabel('Intensity',fontsize=12)
586    if self.Contour:
587        ContourZ = []
588        ContourY = []
589        Nseq = 0
590    for N,Pattern in enumerate(PlotList):
591        Parms = ParmList[N]
592        ifpicked = False
593        LimitId = 0
594        xye = np.array(Pattern[1])
595        if PickId:
596            ifpicked = Pattern[2] == self.PatternTree.GetItemText(PatternId)
597            LimitId = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
598        if self.qPlot:
599            Id = G2gd.GetPatternTreeItemId(self,self.root, Pattern[2])
600            X = 4*np.pi*npsind(xye[0]/2.0)/Parms[1]
601        else:
602            X = xye[0]
603        if not lenX:
604            lenX = len(X)           
605        Y = xye[1]+offset*N
606        if LimitId:
607            limits = np.array(self.PatternTree.GetItemPyData(LimitId))
608            if self.qPlot:
609                limits = 4*np.pi*npsind(limits/2.0)/Parms[1]
610            Lines.append(Plot.axvline(limits[1][0],color='g',dashes=(5,5),picker=3.))   
611            Lines.append(Plot.axvline(limits[1][1],color='r',dashes=(5,5),picker=3.))                   
612        if self.Contour:
613           
614            if lenX == len(X):
615                ContourY.append(N)
616                ContourZ.append(Y)
617                ContourX = X
618                Nseq += 1
619                Plot.set_ylabel('Data sequence',fontsize=12)
620        else:
621            X += self.Offset[1]*.005*N
622            if ifpicked:
623                Z = xye[3]+offset*N
624                W = xye[4]+offset*N
625                D = xye[5]+offset*N-Ymax*.02
626                if self.Weight:
627                    W2 = np.sqrt(xye[2])
628                    D *= W2-Ymax*.02
629                if self.logPlot:
630                    Plot.semilogy(X,Y,colors[N%6]+'+',picker=3.,clip_on=False,nonposy='mask')
631                    Plot.semilogy(X,Z,colors[(N+1)%6],picker=False,nonposy='mask')
632                    Plot.semilogy(X,W,colors[(N+2)%6],picker=False,nonposy='mask')
633                else:
634                    Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
635                    Plot.plot(X,Z,colors[(N+1)%6],picker=False)
636                    Plot.plot(X,W,colors[(N+2)%6],picker=False)
637                    Plot.plot(X,D,colors[(N+3)%6],picker=False)
638                    Plot.axhline(0.,color=wx.BLACK)
639                Page.canvas.SetToolTipString('')
640                if self.PatternTree.GetItemText(PickId) == 'Peak List':
641                    tip = 'On data point: Pick peak - L or R MB. On line: L-move, R-delete'
642                    Page.canvas.SetToolTipString(tip)
643                    data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Peak List'))
644                    for item in data:
645                        if self.qPlot:
646                            Lines.append(Plot.axvline(4*math.pi*sind(item[0]/2.)/Parms[1],color=colors[N%6],picker=2.))
647                        else:
648                            Lines.append(Plot.axvline(item[0],color=colors[N%6],picker=2.))
649                if self.PatternTree.GetItemText(PickId) == 'Limits':
650                    tip = 'On data point: Lower limit - L MB; Upper limit - R MB. On limit: MB down to move'
651                    Page.canvas.SetToolTipString(tip)
652                    data = self.LimitsTable.GetData()
653            else:
654                if self.logPlot:
655                    Plot.semilogy(X,Y,colors[N%6],picker=False,nonposy='clip')
656                else:
657                    Plot.plot(X,Y,colors[N%6],picker=False)
658    if PickId and self.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
659        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))[1::2]
660        Parms = dict(zip(Names,Values))
661        try:
662            wave = Parms['Lam']
663        except KeyError:
664            wave = Parms['Lam1']
665        peaks = np.array((self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))))
666        for peak in peaks:
667            if self.qPlot:
668                Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
669            else:
670                Plot.axvline(peak[0],color='b')
671        for hkl in self.HKL:
672            if self.qPlot:
673                Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
674            else:
675                Plot.axvline(hkl[5],color='r',dashes=(5,5))
676    if self.Contour:
677        acolor = mpl.cm.get_cmap(self.ContourColor)
678        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*self.Cmax,interpolation=self.Interpolate, 
679            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
680        Page.figure.colorbar(Img)
681    else:
682        self.Lines = Lines
683    if not newPlot:
684        Page.toolbar.push_current()
685        Plot.set_xlim(xylim[0])
686        Plot.set_ylim(xylim[1])
687        xylim = []
688        Page.toolbar.push_current()
689        Page.toolbar.draw()
690    else:
691        Page.canvas.draw()
692    self.Pwdr = True
693   
694def PlotISFG(self,newPlot=False,type=''):
695    ''' PLotting package for PDF analysis; displays I(q), S(q), F(q) and G(r) as single
696    or multiple plots with waterfall and contour plots as options
697    '''
698    if not type:
699        type = self.G2plotNB.plotList[self.G2plotNB.nb.GetSelection()]
700    if type not in ['I(Q)','S(Q)','F(Q)','G(R)']:
701        return
702    superMinusOne = unichr(0xaf)+unichr(0xb9)
703   
704    def OnPlotKeyPress(event):
705        newPlot = False
706        if event.key == 'u':
707            if self.Contour:
708                self.Cmax = min(1.0,self.Cmax*1.2)
709            elif self.Offset[0] < 100.:
710                self.Offset[0] += 1.
711        elif event.key == 'd':
712            if self.Contour:
713                self.Cmax = max(0.0,self.Cmax*0.8)
714            elif self.Offset[0] > 0.:
715                self.Offset[0] -= 1.
716        elif event.key == 'l':
717            self.Offset[1] -= 1.
718        elif event.key == 'r':
719            self.Offset[1] += 1.
720        elif event.key == 'o':
721            self.Offset = [0,0]
722        elif event.key == 'c':
723            newPlot = True
724            if self.Contour:
725                self.Contour = False
726            else:
727                self.Contour = True
728                self.SinglePlot = False
729                self.Offset = [0.,0.]
730        elif event.key == 's':
731            if self.Contour:
732                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
733                choice.sort()
734                dlg = wx.SingleChoiceDialog(self,'Select','Color scheme',choice)
735                if dlg.ShowModal() == wx.ID_OK:
736                    sel = dlg.GetSelection()
737                    self.ContourColor = choice[sel]
738                else:
739                    self.ContourColor = 'Paired'
740                dlg.Destroy()
741            else:               
742                if self.SinglePlot:
743                    self.SinglePlot = False
744                else:
745                    self.SinglePlot = True
746        elif event.key == 'i':                  #for smoothing contour plot
747            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
748               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
749               'mitchell','sinc','lanczos']
750            dlg = wx.SingleChoiceDialog(self,'Select','Interpolation',choice)
751            if dlg.ShowModal() == wx.ID_OK:
752                sel = dlg.GetSelection()
753                self.Interpolate = choice[sel]
754            else:
755                self.Interpolate = 'nearest'
756            dlg.Destroy()
757        elif event.key == 't' and not self.Contour:
758            if self.Legend:
759                self.Legend = False
760            else:
761                self.Legend = True
762           
763           
764        PlotISFG(self,newPlot=newPlot,type=type)
765       
766    def OnKeyBox(event):
767        if self.G2plotNB.nb.GetSelection() == self.G2plotNB.plotList.index(type):
768            event.key = cb.GetValue()[0]
769            cb.SetValue(' key press')
770            OnPlotKeyPress(event)
771                       
772    def OnMotion(event):
773        xpos = event.xdata
774        if xpos:                                        #avoid out of frame mouse position
775            ypos = event.ydata
776            Page.canvas.SetCursor(wx.CROSS_CURSOR)
777            try:
778                if self.Contour:
779                    self.G2plotNB.status.SetStatusText('R =%.3fA pattern ID =%5d'%(xpos,int(ypos)),1)
780                else:
781                    self.G2plotNB.status.SetStatusText('R =%.3fA %s =%.2f'%(xpos,type,ypos),1)                   
782            except TypeError:
783                self.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
784   
785    xylim = []
786    try:
787        plotNum = self.G2plotNB.plotList.index(type)
788        Page = self.G2plotNB.nb.GetPage(plotNum)
789        if not newPlot:
790            Plot = Page.figure.gca()          #get previous plot & get limits
791            xylim = Plot.get_xlim(),Plot.get_ylim()
792        Page.figure.clf()
793        Plot = Page.figure.gca()
794    except ValueError:
795        newPlot = True
796        self.Cmax = 1.0
797        Plot = self.G2plotNB.addMpl(type).gca()
798        plotNum = self.G2plotNB.plotList.index(type)
799        Page = self.G2plotNB.nb.GetPage(plotNum)
800        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
801        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
802   
803    Page.SetFocus()
804    self.G2plotNB.status.DestroyChildren()
805    if self.Contour:
806        Choice = (' key press','d: lower contour max','u: raise contour max',
807            'i: interpolation method','s: color scheme','c: contour off')
808    else:
809        Choice = (' key press','l: offset left','r: offset right','d: offset down','u: offset up',
810            'o: reset offset','t: toggle legend','c: contour on','s: toggle single plot')
811    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
812        choices=Choice)
813    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
814    cb.SetValue(' key press')
815    PatternId = self.PatternId
816    PickId = self.PickId
817    Plot.set_title(type)
818    if type == 'G(R)':
819        Plot.set_xlabel(r'$R,\AA$',fontsize=14)
820    else:
821        Plot.set_xlabel(r'$Q,\AA$'+superMinusOne,fontsize=14)
822    Plot.set_ylabel(r''+type,fontsize=14)
823    colors=['b','g','r','c','m','k']
824    name = self.PatternTree.GetItemText(PatternId)[4:]
825    Pattern = []   
826    if self.SinglePlot:
827        name = self.PatternTree.GetItemText(PatternId)
828        name = type+name[4:]
829        Id = G2gd.GetPatternTreeItemId(self,PatternId,name)
830        Pattern = self.PatternTree.GetItemPyData(Id)
831        if Pattern:
832            Pattern.append(name)
833        PlotList = [Pattern,]
834    else:
835        PlotList = []
836        item, cookie = self.PatternTree.GetFirstChild(self.root)
837        while item:
838            if 'PDF' in self.PatternTree.GetItemText(item):
839                name = type+self.PatternTree.GetItemText(item)[4:]
840                Id = G2gd.GetPatternTreeItemId(self,item,name)
841                Pattern = self.PatternTree.GetItemPyData(Id)
842                if Pattern:
843                    Pattern.append(name)
844                    PlotList.append(Pattern)
845            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
846    PDFdata = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'PDF Controls'))
847    numbDen = G2pwd.GetNumDensity(PDFdata['ElList'],PDFdata['Form Vol'])
848    Xb = [0.,10.]
849    Yb = [0.,-40.*np.pi*numbDen]
850    Ymax = 1.0
851    lenX = 0
852    for Pattern in PlotList:
853        xye = Pattern[1]
854        Ymax = max(Ymax,max(xye[1]))
855    offset = self.Offset[0]*Ymax/100.0
856    if self.Contour:
857        ContourZ = []
858        ContourY = []
859        Nseq = 0
860    for N,Pattern in enumerate(PlotList):
861        xye = Pattern[1]
862        if PickId:
863            ifpicked = Pattern[2] == self.PatternTree.GetItemText(PatternId)
864        X = xye[0]
865        if not lenX:
866            lenX = len(X)           
867        Y = xye[1]+offset*N
868        if self.Contour:
869            if lenX == len(X):
870                ContourY.append(N)
871                ContourZ.append(Y)
872                ContourX = X
873                Nseq += 1
874                Plot.set_ylabel('Data sequence',fontsize=12)
875        else:
876            X = xye[0]+self.Offset[1]*.005*N
877            if ifpicked:
878                Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
879                Page.canvas.SetToolTipString('')
880            else:
881                if self.Legend:
882                    Plot.plot(X,Y,colors[N%6],picker=False,label='Azm:'+Pattern[2].split('=')[1])
883                else:
884                    Plot.plot(X,Y,colors[N%6],picker=False)
885            if type == 'G(R)':
886                Plot.plot(Xb,Yb,color='k',dashes=(5,5))
887            elif type == 'F(Q)':
888                Plot.axhline(0.,color=wx.BLACK)
889            elif type == 'S(Q)':
890                Plot.axhline(1.,color=wx.BLACK)
891    if self.Contour:
892        acolor = mpl.cm.get_cmap(self.ContourColor)
893        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*self.Cmax,interpolation=self.Interpolate, 
894            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
895        Page.figure.colorbar(Img)
896    elif self.Legend:
897        Plot.legend(loc='best')
898    if not newPlot:
899        Page.toolbar.push_current()
900        Plot.set_xlim(xylim[0])
901        Plot.set_ylim(xylim[1])
902        xylim = []
903        Page.toolbar.push_current()
904        Page.toolbar.draw()
905    else:
906        Page.canvas.draw()
907       
908def PlotXY(self,XY,newPlot=False,type=''):
909    '''simple plot of xy data, used for diagnostic purposes
910    '''
911    def OnMotion(event):
912        xpos = event.xdata
913        if xpos:                                        #avoid out of frame mouse position
914            ypos = event.ydata
915            Page.canvas.SetCursor(wx.CROSS_CURSOR)
916            try:
917                self.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,type,ypos),1)                   
918            except TypeError:
919                self.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
920
921    try:
922        plotNum = self.G2plotNB.plotList.index(type)
923        Page = self.G2plotNB.nb.GetPage(plotNum)
924        if not newPlot:
925            Plot = Page.figure.gca()
926            xylim = Plot.get_xlim(),Plot.get_ylim()
927        Page.figure.clf()
928        Plot = Page.figure.gca()
929    except ValueError:
930        newPlot = True
931        Plot = self.G2plotNB.addMpl(type).gca()
932        plotNum = self.G2plotNB.plotList.index(type)
933        Page = self.G2plotNB.nb.GetPage(plotNum)
934        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
935   
936    Page.SetFocus()
937    self.G2plotNB.status.DestroyChildren()
938    Plot.set_title(type)
939    Plot.set_xlabel(r'X',fontsize=14)
940    Plot.set_ylabel(r''+type,fontsize=14)
941    colors=['b','g','r','c','m','k']
942    Ymax = 1.0
943    lenX = 0
944    X,Y = XY[:2]
945    Ymax = max(Ymax,max(Y))
946    Plot.plot(X,Y,'k',picker=False)
947    if not newPlot:
948        Page.toolbar.push_current()
949        Plot.set_xlim(xylim[0])
950        Plot.set_ylim(xylim[1])
951        xylim = []
952        Page.toolbar.push_current()
953        Page.toolbar.draw()
954    else:
955        Page.canvas.draw()
956
957def PlotPowderLines(self):
958    ''' plotting of powder lines (i.e. no powder pattern) as sticks
959    '''
960    global HKL
961
962    def OnMotion(event):
963        xpos = event.xdata
964        if xpos:                                        #avoid out of frame mouse position
965            Page.canvas.SetCursor(wx.CROSS_CURSOR)
966            self.G2plotNB.status.SetFields(['','2-theta =%9.3f '%(xpos,)])
967            if self.PickId and self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List']:
968                found = []
969                if len(HKL):
970                    view = Page.toolbar._views.forward()[0][:2]
971                    wid = view[1]-view[0]
972                    found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
973                if len(found):
974                    h,k,l = found[0][:3] 
975                    Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
976                else:
977                    Page.canvas.SetToolTipString('')
978
979    try:
980        plotNum = self.G2plotNB.plotList.index('Powder Lines')
981        Page = self.G2plotNB.nb.GetPage(plotNum)
982        Page.figure.clf()
983        Plot = Page.figure.gca()
984    except ValueError:
985        Plot = self.G2plotNB.addMpl('Powder Lines').gca()
986        plotNum = self.G2plotNB.plotList.index('Powder Lines')
987        Page = self.G2plotNB.nb.GetPage(plotNum)
988        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
989       
990    Page.SetFocus()
991    Plot.set_title('Powder Pattern Lines')
992    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
993    PickId = self.PickId
994    PatternId = self.PatternId
995    peaks = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))
996    for peak in peaks:
997        Plot.axvline(peak[0],color='b')
998    HKL = np.array(self.HKL)
999    for hkl in self.HKL:
1000        Plot.axvline(hkl[5],color='r',dashes=(5,5))
1001    xmin = peaks[0][0]
1002    xmax = peaks[-1][0]
1003    delt = xmax-xmin
1004    xlim = [max(0,xmin-delt/20.),min(180.,xmax+delt/20.)]
1005    Plot.set_xlim(xlim)
1006    Page.canvas.draw()
1007    Page.toolbar.push_current()
1008
1009def PlotPeakWidths(self):
1010    ''' Plotting of instrument broadening terms as function of 2-theta (future TOF)
1011    Seen when "Instrument Parameters" chosen from powder pattern data tree
1012    '''
1013    PatternId = self.PatternId
1014    limitID = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
1015    if limitID:
1016        limits = self.PatternTree.GetItemPyData(limitID)
1017    else:
1018        return
1019    instParms = self.PatternTree.GetItemPyData( \
1020        G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))
1021    if instParms[0][0] == 'PXC':
1022        lam = instParms[1][1]
1023        if len(instParms[1]) == 13:
1024            GU,GV,GW,LX,LY = instParms[0][6:11]
1025        else:
1026            GU,GV,GW,LX,LY = instParms[0][4:9]
1027    peakID = G2gd.GetPatternTreeItemId(self,PatternId, 'Peak List')
1028    if peakID:
1029        peaks = self.PatternTree.GetItemPyData(peakID)
1030    else:
1031        peaks = []
1032   
1033    try:
1034        plotNum = self.G2plotNB.plotList.index('Peak Widths')
1035        Page = self.G2plotNB.nb.GetPage(plotNum)
1036        Page.figure.clf()
1037        Plot = Page.figure.gca()
1038    except ValueError:
1039        Plot = self.G2plotNB.addMpl('Peak Widths').gca()
1040        plotNum = self.G2plotNB.plotList.index('Peak Widths')
1041        Page = self.G2plotNB.nb.GetPage(plotNum)
1042    Page.SetFocus()
1043   
1044    Page.canvas.SetToolTipString('')
1045    colors=['b','g','r','c','m','k']
1046    Xmin,Xmax = limits[1]
1047    Xmin = min(0.5,max(Xmin,1))
1048    Xmin /= 2
1049    Xmax /= 2
1050    nPts = 100
1051    delt = (Xmax-Xmin)/nPts
1052    thetas = []
1053    for i in range(nPts):
1054        thetas.append(Xmin+i*delt)
1055    X = []
1056    Y = []
1057    Z = []
1058    W = []
1059    V = []
1060    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
1061    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
1062    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
1063    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*#Ubaldo Bafile - private communication
1064    for theta in thetas:
1065        X.append(4.0*math.pi*sind(theta)/lam)              #q
1066        s = sig(theta,GU,GV,GW)
1067        g = gam(theta,LX,LY)
1068        G = gamFW(g,s)
1069        H = gamFW2(g,s)
1070        Y.append(s/tand(theta))
1071        Z.append(g/tand(theta))
1072        W.append(G/tand(theta))
1073        V.append(H/tand(theta))
1074    Plot.set_title('Instrument and sample peak widths')
1075    Plot.set_ylabel(r'$\Delta q/q, \Delta d/d$',fontsize=14)
1076    Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
1077    Plot.plot(X,Y,color='r',label='Gaussian')
1078    Plot.plot(X,Z,color='g',label='Lorentzian')
1079    Plot.plot(X,W,color='b',label='G+L')
1080    Plot.plot(X,V,color='k',label='G+L2')
1081    X = []
1082    Y = []
1083    Z = []
1084    W = []
1085    V = []
1086    for peak in peaks:
1087        X.append(4.0*math.pi*sind(peak[0]/2.0)/lam)
1088        s = 1.17741*math.sqrt(peak[4])*math.pi/18000.
1089        g = peak[6]*math.pi/18000.
1090        G = gamFW(g,s)
1091        H = gamFW2(g,s)
1092        Y.append(s/tand(peak[0]/2.))
1093        Z.append(g/tand(peak[0]/2.))
1094        W.append(G/tand(peak[0]/2.))
1095        V.append(H/tand(peak[0]/2.))
1096    Plot.plot(X,Y,'+',color='r',label='G peak')
1097    Plot.plot(X,Z,'+',color='g',label='L peak')
1098    Plot.plot(X,W,'+',color='b',label='G+L peak')
1099    Plot.plot(X,V,'+',color='k',label='G+L2 peak')
1100    Plot.legend(loc='best')
1101    Page.canvas.draw()
1102   
1103def PlotStrain(self,data):
1104    '''Plot 3D microstrain figure. In this instance data is for a phase
1105    '''
1106    PatternId = self.PatternId
1107    generalData = data['General']
1108    SGData = generalData['SGData']
1109    MuStrKeys = G2spc.MustrainNames(SGData)
1110    cell = generalData['Cell'][1:]
1111    A,B = G2lat.cell2AB(cell[:6])
1112    Vol = cell[6]
1113    useList = data['Histograms']
1114    numPlots = len(useList)
1115   
1116    try:
1117        plotNum = self.G2plotNB.plotList.index('Microstrain')
1118        Page = self.G2plotNB.nb.GetPage(plotNum)
1119        Page.figure.clf()
1120        Plot = mp3d.Axes3D(Page.figure)
1121    except ValueError:
1122        Plot = mp3d.Axes3D(self.G2plotNB.add3D('Microstrain'))
1123        plotNum = self.G2plotNB.plotList.index('Microstrain')
1124        Page = self.G2plotNB.nb.GetPage(plotNum)
1125    Page.SetFocus()
1126    self.G2plotNB.status.SetStatusText('Adjust frame size to get desired aspect ratio',1)
1127   
1128    for item in useList:
1129        if useList[item]['Show']:
1130            muStrain = useList[item]['Mustrain']
1131            PHI = np.linspace(0.,360.,30,True)
1132            PSI = np.linspace(0.,180.,30,True)
1133            X = np.outer(npsind(PHI),npsind(PSI))
1134            Y = np.outer(npcosd(PHI),npsind(PSI))
1135            Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
1136            if muStrain[0] == 'isotropic':
1137                muiso = muStrain[1][0]*math.pi/0.018      #centidegrees to radians!
1138                X *= muiso
1139                Y *= muiso
1140                Z *= muiso               
1141               
1142            elif muStrain[0] == 'uniaxial':
1143                def uniaxMustrain(xyz,muiso,muaniso,axes):
1144                    cp = abs(np.dot(xyz,axes))
1145                    S = muiso+muaniso*cp
1146                    return S*xyz*math.pi/0.018      #centidegrees to radians!
1147                muiso,muaniso = muStrain[1][:2]
1148                axes = np.inner(A,np.array(muStrain[3]))
1149                axes /= nl.norm(axes)
1150                Shkl = np.array(muStrain[1])
1151                Shape = X.shape[0]
1152                XYZ = np.dstack((X,Y,Z))
1153                XYZ = np.nan_to_num(np.apply_along_axis(uniaxMustrain,2,XYZ,muiso,muaniso,axes))
1154                X,Y,Z = np.dsplit(XYZ,3)
1155                X = X[:,:,0]
1156                Y = Y[:,:,0]
1157                Z = Z[:,:,0]
1158               
1159            elif muStrain[0] == 'generalized':
1160                def genMustrain(xyz,SGData,A,Shkl):
1161                    uvw = np.inner(A.T,xyz)
1162                    Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
1163                    sum = np.sqrt(np.sum(np.multiply(Shkl,Strm)))*math.pi/0.018      #centidegrees to radians!
1164                    return sum*xyz
1165                Shkl = np.array(muStrain[4])
1166                if np.any(Shkl):
1167                    Shape = X.shape[0]
1168                    XYZ = np.dstack((X,Y,Z))
1169                    XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,SGData,A,Shkl))
1170                    X,Y,Z = np.dsplit(XYZ,3)
1171                    X = X[:,:,0]
1172                    Y = Y[:,:,0]
1173                    Z = Z[:,:,0]
1174                   
1175            if np.any(X) and np.any(Y) and np.any(Z):
1176                Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g')
1177               
1178            Plot.set_xlabel('X')
1179            Plot.set_ylabel('Y')
1180            Plot.set_zlabel('Z')
1181    Page.canvas.draw()
1182   
1183def PlotSphHarmTexture(self,generalData):
1184    '''Pole figure, inverse pole figure(?), 3D pole distribution and 3D inverse pole distribution(?)
1185    plotting; Need way to select 
1186    pole figure or pole distribution to be displayed - do in key enter menu
1187    dict generalData contains all phase info needed
1188    '''
1189    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1190    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1191    SGData = generalData['SGData']
1192    textureData = generalData['SH Texture']
1193    print 'Texture plot'
1194    SHData = generalData['SH Texture']
1195    SHCoef = SHData['SH Coeff'][1]
1196    cell = generalData['Cell'][1:7]
1197    PH = np.array(SHData['PFhkl']) 
1198    phi,beta = G2lat.CrsAng(PH,cell,SGData)
1199    Start = True
1200    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1201
1202           
1203def PlotExposedImage(self,newPlot=False,event=None):
1204    '''General access module for 2D image plotting
1205    '''
1206    plotNo = self.G2plotNB.nb.GetSelection()
1207    if self.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
1208        PlotImage(self,newPlot,event,newImage=True)
1209    elif self.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
1210        PlotIntegration(self,newPlot,event)
1211
1212def PlotImage(self,newPlot=False,event=None,newImage=True):
1213    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
1214    masks, etc.
1215    '''
1216    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
1217    import numpy.ma as ma
1218    Dsp = lambda tth,wave: wave/(2.*sind(tth/2.))
1219    global Data,Masks
1220    colors=['b','g','r','c','m','k']
1221    Data = self.PatternTree.GetItemPyData(
1222        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1223    Masks = self.PatternTree.GetItemPyData(
1224        G2gd.GetPatternTreeItemId(self,self.Image, 'Masks'))
1225
1226    def OnImMotion(event):
1227        Page.canvas.SetToolTipString('')
1228        sizexy = Data['size']
1229        if event.xdata and event.ydata:                 #avoid out of frame errors
1230            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1231            item = self.itemPicked
1232            pixelSize = Data['pixelSize']
1233            scalex = 1000./pixelSize[0]
1234            scaley = 1000./pixelSize[1]
1235            if item and self.PatternTree.GetItemText(self.PickId) == 'Image Controls':
1236                if 'Text' in str(item):
1237                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
1238                else:
1239                    xcent,ycent = Data['center']
1240                    xpos = event.xdata-xcent
1241                    ypos = event.ydata-ycent
1242                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
1243                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
1244                        Page.canvas.SetToolTipString('%6d deg'%(azm))
1245                    elif 'line1' in  str(item) or 'line2' in str(item):
1246                        Page.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
1247            else:
1248                xpos = event.xdata
1249                ypos = event.ydata
1250                xpix = xpos*scalex
1251                ypix = ypos*scaley
1252                Int = 0
1253                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
1254                    Int = self.ImageZ[ypix][xpix]
1255                tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
1256                Q = 2.*math.pi/dsp
1257                if self.setPoly:
1258                    self.G2plotNB.status.SetFields(['','Polygon mask pick - LB next point, RB close polygon'])
1259                else:
1260                    self.G2plotNB.status.SetFields(\
1261                        ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)])
1262
1263    def OnImPlotKeyPress(event):
1264        try:
1265            PickName = self.PatternTree.GetItemText(self.PickId)
1266        except TypeError:
1267            return
1268        if PickName == 'Masks':
1269            Xpos = event.xdata
1270            if not Xpos:            #got point out of frame
1271                return
1272            Ypos = event.ydata
1273            if event.key == 's':
1274                Masks['Points'].append([Xpos,Ypos,1.])
1275            elif event.key == 'r':
1276                tth = G2img.GetTth(Xpos,Ypos,Data)
1277                Masks['Rings'].append([tth,0.1])
1278            elif event.key == 'a':
1279                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
1280                azm = int(azm)               
1281                Masks['Arcs'].append([tth,[azm-5,azm+5],0.1])
1282            elif event.key == 'p':
1283                self.setPoly = True
1284                Masks['Polygons'].append([])
1285                self.G2plotNB.status.SetFields(['','Polygon mask active - LB pick next point, RB close polygon'])
1286            G2imG.UpdateMasks(self,Masks)
1287        elif PickName == 'Image Controls':
1288            if event.key == 'c':
1289                Xpos = event.xdata
1290                if not Xpos:            #got point out of frame
1291                    return
1292                Ypos = event.ydata
1293                dlg = wx.MessageDialog(self,'Are you sure you want to change the center?',
1294                    'Center change',style=wx.OK|wx.CANCEL)
1295                try:
1296                    if dlg.ShowModal() == wx.ID_OK:
1297                        print 'move center to: ',Xpos,Ypos
1298                        Data['center'] = [Xpos,Ypos]
1299                        G2imG.UpdateImageControls(self,Data,Masks)
1300                finally:
1301                    dlg.Destroy()
1302            elif event.key == 'l':
1303                if self.logPlot:
1304                    self.logPlot = False
1305                else:
1306                    self.logPlot = True
1307        PlotImage(self,newImage=True)
1308           
1309    def OnKeyBox(event):
1310        if self.G2plotNB.nb.GetSelection() == self.G2plotNB.plotList.index('2D Powder Image'):
1311            event.key = cb.GetValue()[0]
1312            cb.SetValue(' key press')
1313            if event.key in 'l':
1314                OnImPlotKeyPress(event)
1315                       
1316    def OnImPick(event):
1317        if self.PatternTree.GetItemText(self.PickId) not in ['Image Controls','Masks']:
1318            return
1319        if self.setPoly:
1320            polygon = Masks['Polygons'][-1]
1321            xpos,ypos = event.mouseevent.xdata,event.mouseevent.ydata
1322            if xpos and ypos:                       #point inside image
1323                if len(polygon) > 2 and event.mouseevent.button == 3:
1324                    x0,y0 = polygon[0]
1325                    polygon.append([x0,y0])
1326                    self.setPoly = False
1327                    self.G2plotNB.status.SetFields(['','Polygon closed - RB drag a vertex to change shape'])
1328                else:
1329                    self.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(xpos,ypos)])
1330                    polygon.append([xpos,ypos])
1331                G2imG.UpdateMasks(self,Masks)
1332        else:
1333            if self.itemPicked is not None: return
1334            self.itemPicked = event.artist
1335            self.mousePicked = event.mouseevent
1336       
1337    def OnImRelease(event):
1338        try:
1339            PickName = self.PatternTree.GetItemText(self.PickId)
1340        except TypeError:
1341            return
1342        if PickName not in ['Image Controls','Masks']:
1343            return
1344        pixelSize = Data['pixelSize']
1345        scalex = 1000./pixelSize[0]
1346        scaley = 1000./pixelSize[1]
1347        pixLimit = Data['pixLimit']
1348        if self.itemPicked is None and PickName == 'Image Controls':
1349#            sizexy = Data['size']
1350            Xpos = event.xdata
1351            if not (Xpos and self.ifGetRing):                   #got point out of frame
1352                return
1353            Ypos = event.ydata
1354            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
1355                if event.button == 1:
1356                    Xpix = Xpos*scalex
1357                    Ypix = Ypos*scaley
1358                    xpos,ypos,I,J = G2img.ImageLocalMax(self.ImageZ,pixLimit,Xpix,Ypix)
1359                    if I and J:
1360                        xpos += .5                              #shift to pixel center
1361                        ypos += .5
1362                        xpos /= scalex                          #convert to mm
1363                        ypos /= scaley
1364                        Data['ring'].append([xpos,ypos])
1365                elif event.button == 3:
1366                    self.dataFrame.GetStatusBar().SetStatusText('Calibrating...')
1367                    if G2img.ImageCalibrate(self,Data):
1368                        self.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check')
1369                        print 'Calibration successful'
1370                    else:
1371                        self.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose')
1372                        print 'Calibration failed'
1373                    self.ifGetRing = False
1374                    G2imG.UpdateImageControls(self,Data,Masks)
1375                    return
1376                PlotImage(self,newImage=False)
1377            return
1378        else:
1379            xpos = event.xdata
1380            if xpos:                                        #avoid out of frame mouse position
1381                ypos = event.ydata
1382                if self.ifGetRing:                          #delete a calibration ring pick
1383                    xypos = [xpos,ypos]
1384                    rings = Data['ring']
1385                    for ring in rings:
1386                        if np.allclose(ring,xypos,.01,0):
1387                            rings.remove(ring)                                                                       
1388                else:
1389                    tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
1390                    itemPicked = str(self.itemPicked)
1391                    if 'Line2D' in itemPicked and PickName == 'Image Controls':
1392                        if 'line1' in itemPicked:
1393                            Data['IOtth'][0] = max(tth,0.001)
1394                        elif 'line2' in itemPicked:
1395                            Data['IOtth'][1] = tth
1396                        elif 'line3' in itemPicked:
1397                            Data['LRazimuth'][0] = int(azm)
1398                            if Data['fullIntegrate']:
1399                                Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
1400                        elif 'line4' in itemPicked and not Data['fullIntegrate']:
1401                            Data['LRazimuth'][1] = int(azm)
1402                           
1403                        if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
1404                            Data['LRazimuth'][0] -= 360
1405                           
1406                        azLim = np.array(Data['LRazimuth'])   
1407                        if np.any(azLim>360):
1408                            azLim -= 360
1409                            Data['LRazimuth'] = list(azLim)
1410                           
1411                        if  Data['IOtth'][0] > Data['IOtth'][1]:
1412                            Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
1413                           
1414                        self.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
1415                        self.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
1416                        self.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
1417                        self.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
1418                    elif 'Circle' in itemPicked and PickName == 'Masks':
1419                        spots = Masks['Points']
1420                        newPos = itemPicked.split(')')[0].split('(')[2].split(',')
1421                        newPos = np.array([float(newPos[0]),float(newPos[1])])
1422                        for spot in spots:
1423                            if np.allclose(np.array([spot[:2]]),newPos):
1424                                spot[:2] = xpos,ypos
1425                        G2imG.UpdateMasks(self,Masks)
1426                    elif 'Line2D' in itemPicked and PickName == 'Masks':
1427                        Obj = self.itemPicked.findobj()
1428                        rings = Masks['Rings']
1429                        arcs = Masks['Arcs']
1430                        polygons = Masks['Polygons']
1431                        for ring in self.ringList:
1432                            if Obj == ring[0]:
1433                                rN = ring[1]
1434                                if ring[2] == 'o':
1435                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)-rings[rN][1]/2.
1436                                else:
1437                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)+rings[rN][1]/2.
1438                        for arc in self.arcList:
1439                            if Obj == arc[0]:
1440                                aN = arc[1]
1441                                if arc[2] == 'o':
1442                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)-arcs[aN][2]/2
1443                                elif arc[2] == 'i':
1444                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)+arcs[aN][2]/2
1445                                elif arc[2] == 'l':
1446                                    arcs[aN][1][0] = int(G2img.GetAzm(xpos,ypos,Data))
1447                                else:
1448                                    arcs[aN][1][1] = int(G2img.GetAzm(xpos,ypos,Data))
1449                        for poly in self.polyList:
1450                            if Obj == poly[0]:
1451                                ind = self.itemPicked.contains(self.mousePicked)[1]['ind'][0]
1452                                oldPos = np.array([self.mousePicked.xdata,self.mousePicked.ydata])
1453                                pN = poly[1]
1454                                for i,xy in enumerate(polygons[pN]):
1455                                    if np.allclose(np.array([xy]),oldPos,atol=1.0):
1456                                        polygons[pN][i] = xpos,ypos
1457                        G2imG.UpdateMasks(self,Masks)
1458#                    else:                  #keep for future debugging
1459#                        print str(self.itemPicked),event.xdata,event.ydata,event.button
1460                PlotImage(self,newImage=True)
1461            self.itemPicked = None
1462           
1463    try:
1464        plotNum = self.G2plotNB.plotList.index('2D Powder Image')
1465        Page = self.G2plotNB.nb.GetPage(plotNum)
1466        if not newPlot:
1467            Plot = Page.figure.gca()          #get previous powder plot & get limits
1468            xylim = Plot.get_xlim(),Plot.get_ylim()
1469        if newImage:
1470            Page.figure.clf()
1471            Plot = Page.figure.gca()          #get a fresh plot after clf()
1472       
1473    except ValueError:
1474        Plot = self.G2plotNB.addMpl('2D Powder Image').gca()
1475        plotNum = self.G2plotNB.plotList.index('2D Powder Image')
1476        Page = self.G2plotNB.nb.GetPage(plotNum)
1477        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
1478        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
1479        Page.canvas.mpl_connect('pick_event', OnImPick)
1480        Page.canvas.mpl_connect('button_release_event', OnImRelease)
1481        xylim = []
1482    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
1483        Page.SetFocus()
1484    Title = self.PatternTree.GetItemText(self.Image)[4:]
1485    self.G2plotNB.status.DestroyChildren()
1486    if self.logPlot:
1487        Title = 'log('+Title+')'
1488    Plot.set_title(Title)
1489    try:
1490        if self.PatternTree.GetItemText(self.PickId) in ['Image Controls',]:
1491            if self.logPlot:
1492                Choice = (' key press','l: log(I) off')
1493            else:
1494                Choice = (' key press','l: log(I) on')
1495            cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
1496                choices=Choice)
1497            cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
1498            cb.SetValue(' key press')
1499    except TypeError:
1500        pass
1501    size,imagefile = self.PatternTree.GetItemPyData(self.Image)
1502    if imagefile != self.oldImagefile:
1503        imagefile = G2IO.CheckImageFile(self,imagefile)
1504        if not imagefile:
1505            self.G2plotNB.Delete('2D Powder Image')
1506            return
1507        self.PatternTree.SetItemPyData(self.Image,[size,imagefile])
1508        self.ImageZ = G2IO.GetImageData(self,imagefile,imageOnly=True)
1509#        print self.ImageZ.shape,self.ImageZ.size,Data['size'] #might be useful debugging line
1510        self.oldImagefile = imagefile
1511
1512    imScale = 1
1513    if len(self.ImageZ) > 1024:
1514        imScale = len(self.ImageZ)/1024
1515    sizexy = Data['size']
1516    pixelSize = Data['pixelSize']
1517    scalex = 1000./pixelSize[0]
1518    scaley = 1000./pixelSize[1]
1519    Xmax = sizexy[0]*pixelSize[0]/1000.
1520    Ymax = sizexy[1]*pixelSize[1]/1000.
1521    xlim = (0,Xmax)
1522    ylim = (Ymax,0)
1523    Imin,Imax = Data['range'][1]
1524    acolor = mpl.cm.get_cmap(Data['color'])
1525    xcent,ycent = Data['center']
1526    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
1527    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
1528    #do threshold mask - "real" mask - others are just bondaries
1529    Zlim = Masks['Thresholds'][1]
1530    wx.BeginBusyCursor()
1531    try:
1532           
1533        if newImage:                   
1534            MA = ma.masked_greater(ma.masked_less(self.ImageZ,Zlim[0]),Zlim[1])
1535            MaskA = ma.getmaskarray(MA)
1536            A = G2img.ImageCompress(MA,imScale)
1537            AM = G2img.ImageCompress(MaskA,imScale)
1538            if self.logPlot:
1539                A = np.log(A)
1540                AM = np.log(AM)
1541                Imin,Imax = [np.amin(A),np.amax(A)]
1542            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
1543                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Xmax,0])
1544            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
1545                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
1546            if self.setPoly:
1547                Img.set_picker(True)
1548   
1549        Plot.plot(xcent,ycent,'x')
1550        if Data['showLines']:
1551            LRAzim = Data['LRazimuth']                  #NB: integers
1552            Nazm = Data['outAzimuths']
1553            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
1554            AzmthOff = Data['azmthOff']
1555            IOtth = Data['IOtth']
1556            wave = Data['wavelength']
1557            dspI = wave/(2.0*sind(IOtth[0]/2.0))
1558            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
1559            dspO = wave/(2.0*sind(IOtth[1]/2.0))
1560            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
1561            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
1562            if ellI:
1563                xyI = []
1564                for azm in Azm:
1565                    xyI.append(G2img.GetDetectorXY(dspI,azm-90.,Data))
1566                xyI = np.array(xyI)
1567                arcxI,arcyI = xyI.T
1568                Plot.plot(arcxI,arcyI,picker=3)
1569            if ellO:
1570                xyO = []
1571                for azm in Azm:
1572                    xyO.append(G2img.GetDetectorXY(dspO,azm-90.,Data))
1573                xyO = np.array(xyO)
1574                arcxO,arcyO = xyO.T
1575                Plot.plot(arcxO,arcyO,picker=3)
1576            if ellO and ellI:
1577                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
1578                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
1579            for i in range(Nazm):
1580                cake = LRAzim[0]+i*delAzm
1581                ind = np.searchsorted(Azm,cake)
1582                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
1583                   
1584        for xring,yring in Data['ring']:
1585            Plot.plot(xring,yring,'r+',picker=3)
1586        if Data['setRings']:
1587#            rings = np.concatenate((Data['rings']),axis=0)
1588            N = 0
1589            for ring in Data['rings']:
1590                xring,yring = np.array(ring).T[:2]
1591                Plot.plot(xring,yring,'+',color=colors[N%6])
1592                N += 1           
1593        for ellipse in Data['ellipses']:
1594            cent,phi,[width,height],col = ellipse
1595            Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
1596            Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
1597        #masks - mask lines numbered after integration limit lines
1598        spots = Masks['Points']
1599        rings = Masks['Rings']
1600        arcs = Masks['Arcs']
1601        polygons = Masks['Polygons']
1602        for x,y,d in spots:
1603            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
1604        self.ringList = []
1605        for iring,(tth,thick) in enumerate(rings):
1606            wave = Data['wavelength']
1607            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
1608            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
1609            self.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
1610            self.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
1611        self.arcList = []
1612        for iarc,(tth,azm,thick) in enumerate(arcs):           
1613            wave = Data['wavelength']
1614            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
1615            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
1616            self.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
1617            self.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
1618            self.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
1619            self.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
1620        self.polyList = []
1621        for ipoly,polygon in enumerate(polygons):
1622            x,y = np.hsplit(np.array(polygon),2)
1623            self.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
1624            Plot.plot(x,y,'r')           
1625        if newImage:
1626            colorBar = Page.figure.colorbar(Img)
1627        Plot.set_xlim(xlim)
1628        Plot.set_ylim(ylim)
1629        if not newPlot and xylim:
1630            Page.toolbar.push_current()
1631            Plot.set_xlim(xylim[0])
1632            Plot.set_ylim(xylim[1])
1633            xylim = []
1634            Page.toolbar.push_current()
1635            Page.toolbar.draw()
1636        else:
1637            Page.canvas.draw()
1638    finally:
1639        wx.EndBusyCursor()
1640       
1641def PlotIntegration(self,newPlot=False,event=None):
1642    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
1643    '''
1644           
1645    def OnMotion(event):
1646        Page.canvas.SetToolTipString('')
1647        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1648        azm = event.ydata
1649        tth = event.xdata
1650        if azm and tth:
1651            self.G2plotNB.status.SetFields(\
1652                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1653                               
1654    try:
1655        plotNum = self.G2plotNB.plotList.index('2D Integration')
1656        Page = self.G2plotNB.nb.GetPage(plotNum)
1657        if not newPlot:
1658            Plot = Page.figure.gca()          #get previous plot & get limits
1659            xylim = Plot.get_xlim(),Plot.get_ylim()
1660        Page.figure.clf()
1661        Plot = Page.figure.gca()          #get a fresh plot after clf()
1662       
1663    except ValueError:
1664        Plot = self.G2plotNB.addMpl('2D Integration').gca()
1665        plotNum = self.G2plotNB.plotList.index('2D Integration')
1666        Page = self.G2plotNB.nb.GetPage(plotNum)
1667        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1668        Page.views = False
1669        view = False
1670    if not event:
1671        Page.SetFocus()
1672       
1673    Data = self.PatternTree.GetItemPyData(
1674        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1675    image = self.Integrate[0]
1676    xsc = self.Integrate[1]
1677    ysc = self.Integrate[2]
1678    Imin,Imax = Data['range'][1]
1679    acolor = mpl.cm.get_cmap(Data['color'])
1680    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1681    Plot.set_ylabel('azimuth',fontsize=12)
1682    Plot.set_xlabel('2-theta',fontsize=12)
1683    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
1684        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
1685    colorBar = Page.figure.colorbar(Img)
1686    if Data['setRings'] and Data['rings']:
1687        rings = np.concatenate((Data['rings']),axis=0)
1688        for xring,yring,dsp in rings:
1689            x,y = G2img.GetTthAzm(xring,yring,Data)
1690            Plot.plot(x,y,'r+')
1691    if Data['ellipses']:           
1692        for ellipse in Data['ellipses']:
1693            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1694            x,y = np.hsplit(ring,2)
1695            tth,azm = G2img.GetTthAzm(x,y,Data)
1696            azm = np.where(azm < 0.,azm+360,azm)
1697            Plot.plot(tth,azm,'b,')
1698    if not newPlot:
1699        Page.toolbar.push_current()
1700        Plot.set_xlim(xylim[0])
1701        Plot.set_ylim(xylim[1])
1702        xylim = []
1703        Page.toolbar.push_current()
1704        Page.toolbar.draw()
1705    else:
1706        Page.canvas.draw()
1707               
1708def PlotTRImage(self,tax,tay,taz,newPlot=False):
1709    '''a test plot routine - not normally used
1710    ''' 
1711           
1712    def OnMotion(event):
1713        Page.canvas.SetToolTipString('')
1714        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1715        azm = event.xdata
1716        tth = event.ydata
1717        if azm and tth:
1718            self.G2plotNB.status.SetFields(\
1719                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1720                               
1721    try:
1722        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1723        Page = self.G2plotNB.nb.GetPage(plotNum)
1724        if not newPlot:
1725            Plot = Page.figure.gca()          #get previous plot & get limits
1726            xylim = Plot.get_xlim(),Plot.get_ylim()
1727        Page.figure.clf()
1728        Plot = Page.figure.gca()          #get a fresh plot after clf()
1729       
1730    except ValueError:
1731        Plot = self.G2plotNB.addMpl('2D Transformed Powder Image').gca()
1732        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1733        Page = self.G2plotNB.nb.GetPage(plotNum)
1734        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1735        Page.views = False
1736        view = False
1737    Page.SetFocus()
1738       
1739    Data = self.PatternTree.GetItemPyData(
1740        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1741    Imin,Imax = Data['range'][1]
1742    step = (Imax-Imin)/5.
1743    V = np.arange(Imin,Imax,step)
1744    acolor = mpl.cm.get_cmap(Data['color'])
1745    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1746    Plot.set_xlabel('azimuth',fontsize=12)
1747    Plot.set_ylabel('2-theta',fontsize=12)
1748    Plot.contour(tax,tay,taz,V,cmap=acolor)
1749    if Data['showLines']:
1750        IOtth = Data['IOtth']
1751        if Data['fullIntegrate']:
1752            LRAzim = [-180,180]
1753        else:
1754            LRAzim = Data['LRazimuth']                  #NB: integers
1755        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
1756        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
1757        if not Data['fullIntegrate']:
1758            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
1759            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
1760    if Data['setRings']:
1761        rings = np.concatenate((Data['rings']),axis=0)
1762        for xring,yring,dsp in rings:
1763            x,y = G2img.GetTthAzm(xring,yring,Data)
1764            Plot.plot(y,x,'r+')           
1765    if Data['ellipses']:           
1766        for ellipse in Data['ellipses']:
1767            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1768            x,y = np.hsplit(ring,2)
1769            tth,azm = G2img.GetTthAzm(x,y,Data)
1770            Plot.plot(azm,tth,'b,')
1771    if not newPlot:
1772        Page.toolbar.push_current()
1773        Plot.set_xlim(xylim[0])
1774        Plot.set_ylim(xylim[1])
1775        xylim = []
1776        Page.toolbar.push_current()
1777        Page.toolbar.draw()
1778    else:
1779        Page.canvas.draw()
1780       
1781def PlotStructure(self,data):
1782    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
1783    thermal motion ellipsoids and polyhedra
1784    '''
1785    generalData = data['General']
1786    cell = generalData['Cell'][1:7]
1787    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
1788    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1789    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1790    Mydir = generalData['Mydir']
1791    atomData = data['Atoms']
1792    drawingData = data['Drawing']
1793    drawAtoms = drawingData['Atoms']
1794    cx,ct,cs = drawingData['atomPtrs']
1795    Wt = [255,255,255]
1796    Rd = [255,0,0]
1797    Gr = [0,255,0]
1798    Bl = [0,0,255]
1799    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]])
1800    uEdges = np.array([
1801        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
1802        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
1803        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
1804    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
1805    altDown = False
1806    shiftDown = False
1807    ctrlDown = False
1808   
1809    def OnKeyBox(event):
1810        import Image
1811        Draw()                          #make sure plot is fresh!!
1812        mode = cb.GetValue()
1813        Fname = Mydir+'\\'+generalData['Name']+'.'+mode
1814        size = Page.canvas.GetSize()
1815        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
1816        if mode in ['jpeg',]:
1817            Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
1818            im = Image.new("RGBA", (size[0],size[1]))
1819        else:
1820            Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
1821            im = Image.new("RGB", (size[0],size[1]))
1822        im.fromstring(Pix)
1823        im.save(Fname,mode)
1824        cb.SetValue(' Save as:')
1825        self.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
1826   
1827    def GetTruePosition(xy):
1828        View = glGetIntegerv(GL_VIEWPORT)
1829        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
1830        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
1831        Zmax = 1.
1832        for i,atom in enumerate(drawAtoms):
1833            x,y,z = atom[cx:cx+3]
1834            X,Y,Z = gluProject(x,y,z,Model,Proj,View)
1835            XY = [int(X),int(View[3]-Y)]
1836            if np.allclose(xy,XY,atol=10) and Z < Zmax:
1837                Zmax = Z
1838                SetSelectedAtoms(i)
1839                   
1840    def OnMouseDown(event):
1841        xy = event.GetPosition()
1842        if event.ShiftDown():
1843            GetTruePosition(xy)
1844        else:
1845            drawingData['Rotation'][3] = xy
1846        Draw()
1847       
1848    def OnMouseMove(event):
1849        newxy = event.GetPosition()
1850        page = getSelection()
1851        if event.ControlDown() and drawingData['showABC']:
1852            if event.LeftIsDown():
1853                SetTestRot(newxy)
1854            elif event.RightIsDown():
1855                SetTestPos(newxy)
1856            elif event.MiddleIsDown():
1857                SetTestRotZ(newxy)
1858            x,y,z = drawingData['testPos'][0]
1859            self.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
1860               
1861               
1862        if event.Dragging() and not event.ControlDown():
1863            if event.LeftIsDown():
1864                SetRotation(newxy)
1865                angX,angY,angZ = drawingData['Rotation'][:3]
1866                self.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f ,%.2f'%(angX,angY,angZ),1)
1867            elif event.RightIsDown():
1868                SetTranslation(newxy)
1869                Tx,Ty,Tz = drawingData['viewPoint'][0]
1870                self.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
1871            elif event.MiddleIsDown():
1872                SetRotationZ(newxy)
1873                angX,angY,angZ = drawingData['Rotation'][:3]
1874                self.G2plotNB.status.SetStatusText('New rotation: %.2f, %.2f, %.2f'%(angX,angY,angZ),1)
1875        Draw()
1876       
1877    def OnMouseWheel(event):
1878        drawingData['cameraPos'] += event.GetWheelRotation()/24
1879        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
1880        self.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(drawingData['cameraPos']),1)
1881        page = getSelection()
1882        if page:
1883            if self.dataDisplay.GetPageText(page) == 'Draw Options':
1884                panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1885                names = [child.GetName() for child in panel]
1886                panel[names.index('cameraPos')].SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
1887                panel[names.index('cameraSlider')].SetValue(drawingData['cameraPos'])
1888        Draw()
1889       
1890    def getSelection():
1891        try:
1892            return self.dataDisplay.GetSelection()
1893        except AttributeError:
1894            print self.dataDisplay.GetLabel()
1895            self.G2plotNB.status.SetStatusText('Select this from Phase data window!')
1896            return 0
1897           
1898    def SetViewPointText(VP):
1899        page = getSelection()
1900        if page:
1901            if self.dataDisplay.GetPageText(page) == 'Draw Options':
1902                panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1903                names = [child.GetName() for child in panel]
1904                panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
1905           
1906    def ClearSelectedAtoms():
1907        page = getSelection()
1908        if page:
1909            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1910                self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
1911            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1912                self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
1913                   
1914    def SetSelectedAtoms(ind):
1915        page = getSelection()
1916        if page:
1917            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1918                self.dataDisplay.GetPage(page).SelectRow(ind)      #this is the Atoms grid in Draw Atoms
1919            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1920                Id = drawAtoms[ind][-2]
1921                for i,atom in enumerate(atomData):
1922                    if atom[-1] == Id:
1923                        self.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
1924                 
1925    def GetSelectedAtoms():
1926        page = getSelection()
1927        Ind = []
1928        if page:
1929            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1930                Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
1931            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1932                Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
1933        return Ind
1934                                       
1935    def OnKey(event):           #on key UP!!
1936        keyCode = event.GetKeyCode()
1937        if keyCode > 255:
1938            keyCode = 0
1939        key,xyz = chr(keyCode),event.GetPosition()
1940        indx = drawingData['selectedAtoms']
1941        if key in ['c','C']:
1942            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
1943            drawingData['testPos'] = [[-.1,-.1,-.1],[0.0,0.0,0.0],[0,0]]
1944            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
1945            SetViewPointText(drawingData['viewPoint'][0])
1946        elif key in ['n','N']:
1947            drawAtoms = drawingData['Atoms']
1948            pI = drawingData['viewPoint'][1]
1949            if indx:
1950                pI[0] = indx[pI[1]]
1951                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1952                pI[1] += 1
1953                if pI[1] >= len(indx):
1954                    pI[1] = 0
1955            else:
1956                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1957                pI[0] += 1
1958                if pI[0] >= len(drawAtoms):
1959                    pI[0] = 0
1960            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1961            SetViewPointText(drawingData['viewPoint'][0])
1962               
1963        elif key in ['p','P']:
1964            drawAtoms = drawingData['Atoms']
1965            pI = drawingData['viewPoint'][1]
1966            if indx:
1967                pI[0] = indx[pI[1]]
1968                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1969                pI[1] -= 1
1970                if pI[1] < 0:
1971                    pI[1] = len(indx)-1
1972            else:
1973                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1974                pI[0] -= 1
1975                if pI[0] < 0:
1976                    pI[0] = len(drawAtoms)-1
1977            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1978            SetViewPointText(drawingData['viewPoint'][0])           
1979        Draw()
1980           
1981    def SetBackground():
1982        R,G,B,A = Page.camera['backColor']
1983        glClearColor(R,G,B,A)
1984        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
1985       
1986    def SetLights():
1987        glEnable(GL_DEPTH_TEST)
1988        glShadeModel(GL_SMOOTH)
1989        glEnable(GL_LIGHTING)
1990        glEnable(GL_LIGHT0)
1991        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
1992        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
1993        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
1994       
1995    def SetTranslation(newxy):
1996        Tx,Ty,Tz = drawingData['viewPoint'][0]
1997        anglex,angley,anglez,oldxy = drawingData['Rotation']
1998        Rx = G2lat.rotdMat(anglex,0)
1999        Ry = G2lat.rotdMat(angley,1)
2000        Rz = G2lat.rotdMat(anglez,2)
2001        dxy = list(newxy-oldxy)+[0,]
2002        dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
2003        Tx -= dxy[0]*0.01
2004        Ty += dxy[1]*0.01
2005        Tz -= dxy[2]*0.01
2006        drawingData['Rotation'][3] = newxy
2007        drawingData['viewPoint'][0] =  Tx,Ty,Tz
2008        SetViewPointText([Tx,Ty,Tz])
2009       
2010    def SetTestPos(newxy):
2011        Tx,Ty,Tz = drawingData['testPos'][0]
2012        anglex,angley,anglez,oldxy = drawingData['Rotation']
2013        Rx = G2lat.rotdMat(anglex,0)
2014        Ry = G2lat.rotdMat(angley,1)
2015        Rz = G2lat.rotdMat(anglez,2)
2016        dxy = list(newxy-oldxy)+[0,]
2017        dxy = np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy)))
2018        Tx += dxy[0]*0.001
2019        Ty -= dxy[1]*0.001
2020        Tz += dxy[2]*0.001
2021        drawingData['Rotation'][3] = newxy
2022        drawingData['testPos'][0] =  Tx,Ty,Tz
2023       
2024    def SetTestRot(newxy):
2025        Txyz = np.array(drawingData['testPos'][0])
2026        oldxy = drawingData['testPos'][2]
2027        Ax,Ay,Az = drawingData['testPos'][1]
2028        Vxyz = np.array(drawingData['viewPoint'][0])
2029        Dxyz = np.inner(Amat,Txyz-Vxyz)
2030        dxy = list(newxy-oldxy)+[0,]
2031        Ax += dxy[1]*0.01
2032        Ay += dxy[0]*0.01
2033        Rx = G2lat.rotdMat(Ax,0)
2034        Ry = G2lat.rotdMat(Ay,1)
2035        Dxyz = np.inner(Ry,np.inner(Rx,Dxyz))       
2036        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
2037        drawingData['testPos'][1] = [Ax,Ay,Az]
2038        drawingData['testPos'][2] = newxy
2039        drawingData['testPos'][0] = Dxyz
2040       
2041    def SetTestRotZ(newxy):
2042        Txyz = np.array(drawingData['testPos'][0])
2043        oldxy = drawingData['testPos'][2]
2044        Ax,Ay,Az = drawingData['testPos'][1]
2045        Vxyz = np.array(drawingData['viewPoint'][0])
2046        Dxyz = np.inner(Amat,Txyz-Vxyz)       
2047        dxy = list(newxy-oldxy)+[0,]
2048        Az += (dxy[0]+dxy[1])*.01
2049        Rz = G2lat.rotdMat(Az,2)
2050        Dxyz = np.inner(Rz,Dxyz)       
2051        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
2052        drawingData['testPos'][1] = [Ax,Ay,Az]
2053        drawingData['testPos'][2] = newxy
2054        drawingData['testPos'][0] = Dxyz
2055                             
2056    def SetRotation(newxy):       
2057        anglex,angley,anglez,oldxy = drawingData['Rotation']
2058        dxy = newxy-oldxy
2059        anglex += dxy[1]*.25
2060        angley += dxy[0]*.25
2061        oldxy = newxy
2062        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
2063       
2064    def SetRotationZ(newxy):                       
2065        anglex,angley,anglez,oldxy = drawingData['Rotation']
2066        dxy = newxy-oldxy
2067        anglez += (dxy[0]+dxy[1])*.25
2068        oldxy = newxy
2069        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
2070       
2071    def RenderBox():
2072        glEnable(GL_COLOR_MATERIAL)
2073        glLineWidth(2)
2074        glEnable(GL_BLEND)
2075        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
2076        glEnable(GL_LINE_SMOOTH)
2077        glBegin(GL_LINES)
2078        for line,color in zip(uEdges,uColors):
2079            glColor3ubv(color)
2080            glVertex3fv(line[0])
2081            glVertex3fv(line[1])
2082        glEnd()
2083        glColor4ubv([0,0,0,0])
2084        glDisable(GL_LINE_SMOOTH)
2085        glDisable(GL_BLEND)
2086        glDisable(GL_COLOR_MATERIAL)
2087       
2088    def RenderUnitVectors(x,y,z):
2089        xyz = np.array([x,y,z])
2090        glEnable(GL_COLOR_MATERIAL)
2091        glLineWidth(1)
2092        glPushMatrix()
2093        glTranslate(x,y,z)
2094        glScalef(1/cell[0],1/cell[1],1/cell[2])
2095        glBegin(GL_LINES)
2096        for line,color in zip(uEdges,uColors)[:3]:
2097            glColor3ubv(color)
2098            glVertex3fv(line[0])
2099            glVertex3fv(line[1])
2100        glEnd()
2101        glPopMatrix()
2102        glColor4ubv([0,0,0,0])
2103        glDisable(GL_COLOR_MATERIAL)
2104               
2105    def RenderSphere(x,y,z,radius,color):
2106        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2107        glPushMatrix()
2108        glTranslate(x,y,z)
2109        glMultMatrixf(B4mat.T)
2110        q = gluNewQuadric()
2111        gluSphere(q,radius,20,10)
2112        glPopMatrix()
2113       
2114    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
2115        s1,s2,s3 = E
2116        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2117        glPushMatrix()
2118        glTranslate(x,y,z)
2119        glMultMatrixf(B4mat.T)
2120        glMultMatrixf(R4.T)
2121        glEnable(GL_NORMALIZE)
2122        glScale(s1,s2,s3)
2123        q = gluNewQuadric()
2124        gluSphere(q,ellipseProb,20,10)
2125        glDisable(GL_NORMALIZE)
2126        glPopMatrix()
2127       
2128    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
2129        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2130        glPushMatrix()
2131        glTranslate(x,y,z)
2132        glMultMatrixf(B4mat.T)
2133        for bond in Bonds:
2134            glPushMatrix()
2135            Dx = np.inner(Amat,bond)
2136            Z = np.sqrt(np.sum(Dx**2))
2137            azm = atan2d(-Dx[1],-Dx[0])
2138            phi = acosd(Dx[2]/Z)
2139            glRotate(-azm,0,0,1)
2140            glRotate(phi,1,0,0)
2141            q = gluNewQuadric()
2142            gluCylinder(q,radius,radius,Z,slice,2)
2143            glPopMatrix()           
2144        glPopMatrix()
2145               
2146    def RenderLines(x,y,z,Bonds,color):
2147        xyz = np.array([x,y,z])
2148        glEnable(GL_COLOR_MATERIAL)
2149        glLineWidth(1)
2150        glColor3fv(color)
2151        glPushMatrix()
2152        glBegin(GL_LINES)
2153        for bond in Bonds:
2154            glVertex3fv(xyz)
2155            glVertex3fv(xyz+bond)
2156        glEnd()
2157        glColor4ubv([0,0,0,0])
2158        glPopMatrix()
2159        glDisable(GL_COLOR_MATERIAL)
2160       
2161    def RenderPolyhedra(x,y,z,Faces,color):
2162        glPushMatrix()
2163        glTranslate(x,y,z)
2164        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
2165        glShadeModel(GL_SMOOTH)
2166        glMultMatrixf(B4mat.T)
2167        for face,norm in Faces:
2168            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
2169            glFrontFace(GL_CW)
2170            glNormal3fv(norm)
2171            glBegin(GL_TRIANGLES)
2172            for vert in face:
2173                glVertex3fv(vert)
2174            glEnd()
2175        glPopMatrix()
2176       
2177    def RenderBackbone(Backbone,BackboneColor,radius):
2178        glPushMatrix()
2179        glMultMatrixf(B4mat.T)
2180        glEnable(GL_COLOR_MATERIAL)
2181        glShadeModel(GL_SMOOTH)
2182        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
2183        glePolyCylinder(Backbone,BackboneColor,radius)
2184        glPopMatrix()       
2185        glDisable(GL_COLOR_MATERIAL)
2186       
2187    def RenderLabel(x,y,z,label,r):       
2188        glPushMatrix()
2189        glTranslate(x,y,z)
2190        glMultMatrixf(B4mat.T)
2191        glDisable(GL_LIGHTING)
2192        glColor3f(0,1.,0)
2193        glRasterPos3f(r,r,r)
2194        for c in list(label):
2195            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
2196        glEnable(GL_LIGHTING)
2197        glPopMatrix()
2198                           
2199    def Draw():
2200        import numpy.linalg as nl
2201        Ind = GetSelectedAtoms()
2202        VS = np.array(Page.canvas.GetSize())
2203        aspect = float(VS[0])/float(VS[1])
2204        cPos = drawingData['cameraPos']
2205        Zclip = drawingData['Zclip']*cPos/200.
2206        anglex,angley,anglez = drawingData['Rotation'][:3]
2207        Tx,Ty,Tz = drawingData['viewPoint'][0]
2208        cx,ct,cs = drawingData['atomPtrs']
2209        bondR = drawingData['bondRadius']
2210        G,g = G2lat.cell2Gmat(cell)
2211        GS = G
2212        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
2213        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
2214        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
2215        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
2216       
2217        SetBackground()
2218        glInitNames()
2219        glPushName(0)
2220       
2221        glMatrixMode(GL_PROJECTION)
2222        glLoadIdentity()
2223        glViewport(0,0,VS[0],VS[1])
2224        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
2225        gluLookAt(0,0,cPos,0,0,0,0,1,0)
2226        SetLights()           
2227           
2228        glMatrixMode(GL_MODELVIEW)
2229        glLoadIdentity()
2230        glRotate(anglez,0,0,1)
2231        glRotate(anglex,cosd(anglez),-sind(anglez),0)
2232        glRotate(angley,sind(anglez),cosd(anglez),0)
2233        glMultMatrixf(A4mat.T)
2234        glTranslate(-Tx,-Ty,-Tz)
2235        if drawingData['unitCellBox']:
2236            RenderBox()
2237        if drawingData['showABC']:
2238            x,y,z = drawingData['testPos'][0]
2239#            if altDown:
2240#                self.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
2241#            else:
2242#                self.G2plotNB.status.SetStatusText('test point %.4f,%.4f,%.4f'%(x,y,z),1)           
2243            RenderUnitVectors(x,y,z)
2244        Backbone = []
2245        BackboneColor = []
2246        time0 = time.time()
2247        for iat,atom in enumerate(drawingData['Atoms']):
2248            x,y,z = atom[cx:cx+3]
2249            Bonds = atom[-2]
2250            Faces = atom[-1]
2251            try:
2252                atNum = generalData['AtomTypes'].index(atom[ct])
2253            except ValueError:
2254                atNum = -1
2255            CL = atom[cs+2]
2256            color = np.array(CL)/255.
2257            if iat in Ind:
2258                color = np.array(Gr)/255.
2259            radius = 0.5
2260            if atom[cs] != '':
2261                glLoadName(atom[-3])                   
2262            if 'balls' in atom[cs]:
2263                vdwScale = drawingData['vdwScale']
2264                ballScale = drawingData['ballScale']
2265                if atNum < 0:
2266                    radius = 0.2
2267                elif 'H' == atom[ct]:
2268                    if drawingData['showHydrogen']:
2269                        if 'vdW' in atom[cs] and atNum >= 0:
2270                            radius = vdwScale*generalData['vdWRadii'][atNum]
2271                        else:
2272                            radius = ballScale*drawingData['sizeH']
2273                    else:
2274                        radius = 0.0
2275                else:
2276                    if 'vdW' in atom[cs]:
2277                        radius = vdwScale*generalData['vdWRadii'][atNum]
2278                    else:
2279                        radius = ballScale*generalData['BondRadii'][atNum]
2280                RenderSphere(x,y,z,radius,color)
2281                if 'sticks' in atom[cs]:
2282                    RenderBonds(x,y,z,Bonds,bondR,color)
2283            elif 'ellipsoids' in atom[cs]:
2284                RenderBonds(x,y,z,Bonds,bondR,color)
2285                if atom[cs+3] == 'A':                   
2286                    Uij = atom[cs+5:cs+11]
2287                    U = np.multiply(G2spc.Uij2U(Uij),GS)
2288                    U = np.inner(Amat,np.inner(U,Amat).T)
2289                    E,R = nl.eigh(U)
2290                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2291                    E = np.sqrt(E)
2292                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
2293                        pass
2294                    else:
2295                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,color)                   
2296                else:
2297                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
2298                        pass
2299                    else:
2300                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
2301                        RenderSphere(x,y,z,radius,color)
2302            elif 'lines' in atom[cs]:
2303                radius = 0.1
2304                RenderLines(x,y,z,Bonds,color)
2305#                RenderBonds(x,y,z,Bonds,0.05,color,6)
2306            elif atom[cs] == 'sticks':
2307                radius = 0.1
2308                RenderBonds(x,y,z,Bonds,bondR,color)
2309            elif atom[cs] == 'polyhedra':
2310                RenderPolyhedra(x,y,z,Faces,color)
2311            elif atom[cs] == 'backbone':
2312                if atom[ct-1].split()[0] in ['C','N']:
2313                    Backbone.append(list(np.inner(Amat,np.array([x,y,z]))))
2314                    BackboneColor.append(list(color))
2315                   
2316            if atom[cs+1] == 'type':
2317                RenderLabel(x,y,z,atom[ct],radius)
2318            elif atom[cs+1] == 'name':
2319                RenderLabel(x,y,z,atom[ct-1],radius)
2320            elif atom[cs+1] == 'number':
2321                RenderLabel(x,y,z,str(iat+1),radius)
2322            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
2323                RenderLabel(x,y,z,atom[ct-4],radius)
2324            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
2325                RenderLabel(x,y,z,atom[ct-3],radius)
2326            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
2327                RenderLabel(x,y,z,atom[ct-2],radius)
2328        if Backbone:
2329            RenderBackbone(Backbone,BackboneColor,bondR)
2330#        print time.time()-time0
2331        Page.canvas.SwapBuffers()
2332       
2333    def OnSize(event):
2334        Draw()
2335       
2336    def OnFocus(event):
2337        Draw()
2338       
2339    try:
2340        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
2341        Page = self.G2plotNB.nb.GetPage(plotNum)       
2342    except ValueError:
2343        Plot = self.G2plotNB.addOgl(generalData['Name'])
2344        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
2345        Page = self.G2plotNB.nb.GetPage(plotNum)
2346        Page.views = False
2347        view = False
2348        altDown = False
2349    Page.SetFocus()
2350    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
2351        choices=(' save as:','jpeg','tiff','bmp'))
2352    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
2353    cb.SetValue(' save as:')
2354    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
2355    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
2356    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
2357    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
2358    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
2359    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
2360    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
2361    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
2362    Page.camera['position'] = drawingData['cameraPos']
2363    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
2364    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
2365    Page.canvas.SetCurrent()
2366    Draw()
2367       
Note: See TracBrowser for help on using the repository browser.