source: trunk/GSASIIplot.py @ 180

Last change on this file since 180 was 180, checked in by vondreele, 13 years ago

fix image mask plots
File size: 74.2 KB
Line 
1import math
2import time
3import copy
4import numpy as np
5import numpy.linalg as nl
6import wx
7import wx.aui
8import wx.glcanvas
9import matplotlib as mpl
10import GSASIIpath
11import GSASIIgrid as G2gd
12import GSASIIimage as G2img
13import GSASIIIO as G2IO
14import GSASIIpwdGUI as G2pdG
15import GSASIIimgGUI as G2imG
16import GSASIIphsGUI as G2phG
17import GSASIIlattice as G2lat
18import GSASIIspc as G2spc
19from  OpenGL.GL import *
20from OpenGL.GLU import *
21from OpenGL.GLUT import *
22from OpenGL.GLE import *
23from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as Canvas
24from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
25
26# useful degree trig functions
27sind = lambda x: math.sin(x*math.pi/180.)
28cosd = lambda x: math.cos(x*math.pi/180.)
29tand = lambda x: math.tan(x*math.pi/180.)
30asind = lambda x: 180.*math.asin(x)/math.pi
31acosd = lambda x: 180.*math.acos(x)/math.pi
32atan2d = lambda x,y: 180.*math.atan2(y,x)/math.pi
33atand = lambda x: 180.*math.atan(x)/math.pi
34   
35class G2PlotMpl(wx.Panel):   
36    def __init__(self,parent,id=-1,dpi=None,**kwargs):
37        wx.Panel.__init__(self,parent,id=id,**kwargs)
38        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(5,7))
39        self.canvas = Canvas(self,-1,self.figure)
40        self.toolbar = Toolbar(self.canvas)
41
42        self.toolbar.Realize()
43       
44        sizer=wx.BoxSizer(wx.VERTICAL)
45        sizer.Add(self.canvas,1,wx.EXPAND)
46        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
47        self.SetSizer(sizer)
48       
49class G2PlotOgl(wx.Panel):
50    def __init__(self,parent,id=-1,dpi=None,**kwargs):
51        self.figure = wx.Panel.__init__(self,parent,id=id,**kwargs)
52        self.canvas = wx.glcanvas.GLCanvas(self,-1,**kwargs)
53        self.camera = {}
54        sizer=wx.BoxSizer(wx.VERTICAL)
55        sizer.Add(self.canvas,1,wx.EXPAND)
56        self.SetSizer(sizer)           
57               
58class G2PlotNoteBook(wx.Panel):
59    def __init__(self,parent,id=-1):
60        wx.Panel.__init__(self,parent,id=id)
61        #so one can't delete a plot page!!
62        self.nb = wx.aui.AuiNotebook(self, \
63            style=wx.aui.AUI_NB_DEFAULT_STYLE ^ wx.aui.AUI_NB_CLOSE_ON_ACTIVE_TAB)
64        sizer = wx.BoxSizer()
65        sizer.Add(self.nb,1,wx.EXPAND)
66        self.SetSizer(sizer)
67        self.status = parent.CreateStatusBar()
68        self.status.SetFieldsCount(2)
69        self.status.SetStatusWidths([125,-1])
70        self.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, self.OnPageChanged)
71       
72        self.plotList = []
73           
74    def addMpl(self,name=""):
75        page = G2PlotMpl(self.nb)
76        self.nb.AddPage(page,name)
77       
78        self.plotList.append(name)
79       
80        return page.figure
81       
82    def addOgl(self,name=""):
83        page = G2PlotOgl(self.nb)
84        self.nb.AddPage(page,name)
85       
86        self.plotList.append(name)
87       
88        return page.figure
89       
90    def Delete(self,name):
91        try:
92            item = self.plotList.index(name)
93            del self.plotList[item]
94            self.nb.DeletePage(item)
95        except ValueError:          #no plot of this name - do nothing
96            return     
97               
98    def clear(self):
99        while self.nb.GetPageCount():
100            self.nb.DeletePage(0)
101        self.plotList = []
102        self.status.DestroyChildren()
103       
104    def Rename(self,oldName,newName):
105        try:
106            item = self.plotList.index(oldName)
107            self.plotList[item] = newName
108            self.nb.SetPageText(item,newName)
109        except ValueError:          #no plot of this name - do nothing
110            return     
111       
112    def OnPageChanged(self,event):       
113        if self.plotList:
114            self.status.SetStatusText('Better to select this from GSAS-II data tree',1)
115        self.status.DestroyChildren()                           #get rid of special stuff on status bar
116       
117def PlotSngl(self,newPlot=False):
118    from matplotlib.patches import Circle
119    global HKL,HKLF
120
121    def OnSCMotion(event):
122        xpos = event.xdata
123        if xpos:
124            xpos = round(xpos)                                        #avoid out of frame mouse position
125            ypos = round(event.ydata)
126            zpos = Data['Layer']
127            if '100' in Data['Zone']:
128                HKLtxt = '(%3d,%3d,%3d)'%(zpos,xpos,ypos)
129            elif '010' in Data['Zone']:
130                HKLtxt = '(%3d,%3d,%3d)'%(xpos,zpos,ypos)
131            elif '001' in Data['Zone']:
132                HKLtxt = '(%3d,%3d,%3d)'%(xpos,ypos,zpos)
133            Page.canvas.SetToolTipString(HKLtxt)
134            self.G2plotNB.status.SetFields(['HKL = '+HKLtxt,''])
135               
136    def OnSCPick(event):
137        zpos = Data['Layer']
138        pos = event.artist.center
139        if '100' in Data['Zone']:
140            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(zpos,pos[0],pos[1]))
141            hkl = np.array([zpos,pos[0],pos[1]])
142        elif '010' in Data['Zone']:
143            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],zpos,pos[1]))
144            hkl = np.array([pos[0],zpos,pos[1]])
145        elif '001' in Data['Zone']:
146            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],pos[1],zpos))
147            hkl = np.array([pos[0],pos[1],zpos])
148        h,k,l = hkl
149        hklf = HKLF[np.where(np.all(HKL-hkl == [0,0,0],axis=1))]
150        if len(hklf):
151            Fosq,sig,Fcsq = hklf[0]
152            HKLtxt = '(%3d,%3d,%3d %.2f %.3f %.2f %.2f)'%(h,k,l,Fosq,sig,Fcsq,(Fosq-Fcsq)/(scale*sig))
153            self.G2plotNB.status.SetFields(['','HKL, Fosq, sig, Fcsq, delFsq/sig = '+HKLtxt])
154                                 
155    def OnSCKeyPress(event):
156        print event.key
157
158    try:
159        plotNum = self.G2plotNB.plotList.index('Structure Factors')
160        Page = self.G2plotNB.nb.GetPage(plotNum)
161        if not newPlot:
162            Plot = Page.figure.gca()          #get previous powder plot & get limits
163            xylim = Plot.get_xlim(),Plot.get_ylim()
164        Page.figure.clf()
165        Plot = Page.figure.gca()          #get a fresh plot after clf()
166    except ValueError,error:
167        Plot = self.G2plotNB.addMpl('Structure Factors').gca()
168        plotNum = self.G2plotNB.plotList.index('Structure Factors')
169        Page = self.G2plotNB.nb.GetPage(plotNum)
170#        Page.canvas.mpl_connect('key_press_event', OnSCKeyPress)
171        Page.canvas.mpl_connect('pick_event', OnSCPick)
172        Page.canvas.mpl_connect('motion_notify_event', OnSCMotion)
173    Page.SetFocus()
174   
175    Plot.set_aspect(aspect='equal')
176    HKLref = self.PatternTree.GetItemPyData(self.Sngl)
177    Data = self.PatternTree.GetItemPyData( \
178        G2gd.GetPatternTreeItemId(self,self.Sngl, 'HKL Plot Controls'))
179    Type = Data['Type']           
180    scale = Data['Scale']
181    HKLmax = Data['HKLmax']
182    HKLmin = Data['HKLmin']
183    FosqMax = Data['FoMax']
184    FoMax = math.sqrt(FosqMax)
185    ifFc = Data['ifFc']
186    xlabel = ['k, h=','h, k=','h, l=']
187    ylabel = ['l','l','k']
188    zones = ['100','010','001']
189    pzone = [[1,2],[0,2],[0,1]]
190    izone = zones.index(Data['Zone'])
191    Plot.set_title(self.PatternTree.GetItemText(self.Sngl)[5:])
192    HKL = []
193    HKLF = []
194    for H,Fosq,sig,Fcsq,x,x,x in HKLref:
195        HKL.append(H)
196        HKLF.append([Fosq,sig,Fcsq])
197        if H[izone] == Data['Layer']:
198            B = 0
199            if Type == 'Fosq':
200                A = scale*Fosq/FosqMax
201                B = scale*Fcsq/FosqMax
202                C = abs(A-B)
203            elif Type == 'Fo':
204                A = scale*math.sqrt(max(0,Fosq))/FoMax
205                B = scale*math.sqrt(max(0,Fcsq))/FoMax
206                C = abs(A-B)
207            elif Type == '|DFsq|/sig':
208                A = abs(Fosq-Fcsq)/(scale*sig)
209            elif Type == '|DFsq|>sig':
210                A = abs(Fosq-Fcsq)/(scale*sig)
211                if A < 1.0: A = 0                   
212            elif Type == '|DFsq|>3sig':
213                A = abs(Fosq-Fcsq)/(scale*sig)
214                if A < 3.0: A = 0                   
215            xy = (H[pzone[izone][0]],H[pzone[izone][1]])
216            if A > 0.0:
217                Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w',picker=3))
218            if B:
219                Plot.add_artist(Circle(xy,radius=B,ec='b',fc='w'))
220                radius = C
221                if radius > 0:
222                    if A > B:
223                        Plot.add_artist(Circle(xy,radius=radius,ec='g',fc='g'))
224                    else:                   
225                        Plot.add_artist(Circle(xy,radius=radius,ec='r',fc='r'))
226    HKL = np.array(HKL,dtype=np.int)
227    HKLF = np.array(HKLF)
228    Plot.set_xlabel(xlabel[izone]+str(Data['Layer']),fontsize=12)
229    Plot.set_ylabel(ylabel[izone],fontsize=12)
230    Plot.set_xlim((HKLmin[pzone[izone][0]],HKLmax[pzone[izone][0]]))
231    Plot.set_ylim((HKLmin[pzone[izone][1]],HKLmax[pzone[izone][1]]))
232    if not newPlot:
233        Page.toolbar.push_current()
234        Plot.set_xlim(xylim[0])
235        Plot.set_ylim(xylim[1])
236        xylim = []
237        Page.toolbar.push_current()
238        Page.toolbar.draw()
239    else:
240        Page.canvas.draw()
241       
242def PlotPatterns(self,newPlot=False):
243    global HKL
244   
245    def OnPick(event):
246        if self.itemPicked is not None: return
247        PatternId = self.PatternId
248        PickId = self.PickId
249        pick = event.artist
250        mouse = event.mouseevent
251        xpos = pick.get_xdata()
252        ypos = pick.get_ydata()
253        ind = event.ind
254        view = Page.toolbar._views.forward()
255        if view and 'line2' in str(pick):           #apply offset only for picked powder pattern points
256            ind += np.searchsorted(xye[0],view[0][0])
257        xy = zip(xpos[ind],ypos[ind])[0]
258        if self.PatternTree.GetItemText(PickId) == 'Peak List':
259            if ind.all() != [0]:                                    #picked a data point
260                inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))
261                if len(inst[1]) == 11:
262                    ins = inst[1][4:10]
263                else:
264                    ins = inst[1][6:12]   
265                sig = ins[0]*tand(xy[0]/2.0)**2+ins[1]*tand(xy[0]/2.0)+ins[2]
266                gam = ins[3]/cosd(xy[0]/2.0)+ins[4]*tand(xy[0]/2.0)           
267                data = self.PatternTree.GetItemPyData(self.PickId)
268                XY = [xy[0],0, xy[1],1, sig,0, gam,0]       #default refine intensity 1st   
269                data.append(XY)
270                G2pdG.UpdatePeakGrid(self,data)
271                PlotPatterns(self)
272            else:                                                   #picked a peak list line
273                self.itemPicked = pick
274        elif self.PatternTree.GetItemText(PickId) == 'Limits':
275            if ind.all() != [0]:                                    #picked a data point
276                LimitId = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
277                data = self.PatternTree.GetItemPyData(LimitId)
278                if mouse.button==1:
279                    data[1][0] = min(xy[0],data[1][1])
280                if mouse.button==3:
281                    data[1][1] = max(xy[0],data[1][0])
282                self.PatternTree.SetItemPyData(LimitId,data)
283                G2pdG.UpdateLimitsGrid(self,data)
284                PlotPatterns(self)
285            else:                                                   #picked a limit line
286                self.itemPicked = pick
287       
288    def OnPlotKeyPress(event):
289        newPlot = False
290        if event.key == 'w':
291            if self.Weight:
292                self.Weight = False
293            else:
294                self.Weight = True
295            print 'plot weighting:',self.Weight
296        elif event.key == 'u':
297            if self.Contour:
298                self.Cmax = min(1.0,self.Cmax*1.2)
299            elif self.Offset < 100.:
300                self.Offset += 1.
301        elif event.key == 'd':
302            if self.Contour:
303                self.Cmax = max(0.0,self.Cmax*0.8)
304            elif self.Offset > 0.:
305                self.Offset -= 1.
306        elif event.key == 'c':
307            newPlot = True
308            if self.Contour:
309                self.Contour = False
310            else:
311                self.Contour = True
312                self.SinglePlot = False
313                self.Offset = 0
314        elif event.key == 's':
315            if self.Contour:
316                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
317                dlg = wx.SingleChoiceDialog(self,'Select','Color scheme',choice)
318                if dlg.ShowModal() == wx.ID_OK:
319                    sel = dlg.GetSelection()
320                    self.ContourColor = choice[sel]
321                else:
322                    self.ContourColor = 'Paired'
323                dlg.Destroy()
324            else:               
325                if self.SinglePlot:
326                    self.SinglePlot = False
327                else:
328                    self.SinglePlot = True
329        elif event.key == '+':
330            if self.PickId:
331                self.PickId = False
332        elif event.key == 'i':
333            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
334               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
335               'mitchell','sinc','lanczos']
336            dlg = wx.SingleChoiceDialog(self,'Select','Interpolation',choice)
337            if dlg.ShowModal() == wx.ID_OK:
338                sel = dlg.GetSelection()
339                self.Interpolate = choice[sel]
340            else:
341                self.Interpolate = 'nearest'
342            dlg.Destroy()
343           
344        PlotPatterns(self,newPlot=newPlot)
345       
346    def OnKeyBox(event):
347        if self.G2plotNB.nb.GetSelection() == self.G2plotNB.plotList.index('Powder Patterns'):
348            event.key = cb.GetValue()[0]
349            cb.SetValue(' key press')
350            OnPlotKeyPress(event)
351                       
352    def OnMotion(event):
353        xpos = event.xdata
354        if xpos:                                        #avoid out of frame mouse position
355            ypos = event.ydata
356            Page.canvas.SetCursor(wx.CROSS_CURSOR)
357            try:
358                wave = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
359                    self.PatternId, 'Instrument Parameters'))[0][1]
360                dsp = 0.0
361                if abs(xpos) > 0.:                  #avoid possible singularity at beam center
362                    dsp = wave/(2.*sind(abs(xpos)/2.0))
363                if self.Contour:
364                    self.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f pattern ID =%5d'%(xpos,dsp,int(ypos)),1)
365                else:
366                    self.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f Intensity =%9.1f'%(xpos,dsp,ypos),1)
367                if self.itemPicked:
368                    Page.canvas.SetToolTipString('%9.3f'%(xpos))
369                if self.PickId and self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List']:
370                    found = []
371                    if len(HKL):
372                        view = Page.toolbar._views.forward()[0][:2]
373                        wid = view[1]-view[0]
374                        found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
375                    if len(found):
376                        h,k,l = found[0][:3] 
377                        Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
378                    else:
379                        Page.canvas.SetToolTipString('')
380
381            except TypeError:
382                self.G2plotNB.status.SetStatusText('Select PWDR powder pattern first',1)
383                                                   
384    def OnRelease(event):
385        if self.itemPicked is None: return
386        xpos = event.xdata
387        if xpos:                                        #avoid out of frame mouse position
388            lines = []
389            for line in self.Lines: lines.append(line.get_xdata()[0])
390            lineNo = lines.index(self.itemPicked.get_xdata()[0])
391            if  lineNo in [0,1]:
392                LimitId = G2gd.GetPatternTreeItemId(self,self.PatternId, 'Limits')
393                data = self.PatternTree.GetItemPyData(LimitId)
394#                print 'limits',xpos
395                data[1][lineNo] = xpos
396                self.PatternTree.SetItemPyData(LimitId,data)
397                if self.PatternTree.GetItemText(self.PickId) == 'Limits':
398                    G2pdG.UpdateLimitsGrid(self,data)
399            else:
400                PeakId = G2gd.GetPatternTreeItemId(self,self.PatternId, 'Peak List')
401                data = self.PatternTree.GetItemPyData(PeakId)
402#                print 'peaks',xpos
403                data[lineNo-2][0] = xpos
404                self.PatternTree.SetItemPyData(PeakId,data)
405                G2pdG.UpdatePeakGrid(self,data)
406        PlotPatterns(self)
407        self.itemPicked = None   
408
409    xylim = []
410    try:
411        plotNum = self.G2plotNB.plotList.index('Powder Patterns')
412        Page = self.G2plotNB.nb.GetPage(plotNum)
413        if not newPlot:
414            Plot = Page.figure.gca()          #get previous powder plot & get limits
415            xylim = Plot.get_xlim(),Plot.get_ylim()
416        Page.figure.clf()
417        Plot = Page.figure.gca()          #get a fresh plot after clf()
418    except ValueError,error:
419        newPlot = True
420        self.Cmax = 1.0
421        Plot = self.G2plotNB.addMpl('Powder Patterns').gca()
422        plotNum = self.G2plotNB.plotList.index('Powder Patterns')
423        Page = self.G2plotNB.nb.GetPage(plotNum)
424        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
425        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
426        Page.canvas.mpl_connect('pick_event', OnPick)
427        Page.canvas.mpl_connect('button_release_event', OnRelease)
428    Page.SetFocus()
429    self.G2plotNB.status.DestroyChildren()
430    if self.Contour:
431        Choice = (' key press','d: lower contour max','u: raise contour max',
432            'i: interpolation method','s: color scheme','c: contour off')
433    else:
434        Choice = (' key press','d: offset down','u: offset up','c: contour on','s: toggle single plot','+: no selection')
435    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
436        choices=Choice)
437    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
438    cb.SetValue(' key press')
439   
440    PickId = self.PickId
441    PatternId = self.PatternId
442    colors=['b','g','r','c','m','k']
443    Lines = []
444    if self.SinglePlot:
445        Pattern = self.PatternTree.GetItemPyData(self.PatternId)
446        Pattern.append(self.PatternTree.GetItemText(self.PatternId))
447        PlotList = [Pattern,]
448    else:       
449        PlotList = []
450        item, cookie = self.PatternTree.GetFirstChild(self.root)
451        while item:
452            if 'PWDR' in self.PatternTree.GetItemText(item):
453                Pattern = self.PatternTree.GetItemPyData(item)
454                if len(Pattern) < 3:                    # put name on end if needed
455                    Pattern.append(self.PatternTree.GetItemText(item))
456                PlotList.append(Pattern)
457            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)               
458    Ymax = 1.0
459    lenX = 0
460    HKL = np.array(self.HKL)
461    for Pattern in PlotList:
462        xye = Pattern[1]
463        Ymax = max(Ymax,max(xye[1]))
464    offset = self.Offset*Ymax/100.0
465    Plot.set_title('Powder Patterns')
466    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
467    Plot.set_ylabel('Intensity',fontsize=12)
468    if self.Contour:
469        ContourZ = []
470        ContourY = []
471        Nseq = 0
472    for N,Pattern in enumerate(PlotList):
473        ifpicked = False
474        LimitId = 0
475        xye = np.array(Pattern[1])
476        if PickId:
477            ifpicked = Pattern[2] == self.PatternTree.GetItemText(PatternId)
478            LimitId = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
479        X = xye[0]
480        if not lenX:
481            lenX = len(X)
482        Y = xye[1]+offset*N
483        if LimitId:
484            limits = self.PatternTree.GetItemPyData(LimitId)
485            Lines.append(Plot.axvline(limits[1][0],color='g',dashes=(5,5),picker=3.))   
486            Lines.append(Plot.axvline(limits[1][1],color='r',dashes=(5,5),picker=3.))                   
487        if self.Contour:
488            if lenX == len(X):
489                ContourY.append(N)
490                ContourZ.append(Y)
491                ContourX = X
492                Nseq += 1
493                Plot.set_ylabel('Data sequence',fontsize=12)
494        else:
495            if ifpicked:
496                Z = xye[3]+offset*N
497                W = xye[4]+offset*N
498                D = xye[5]+offset*N-Ymax*.02
499                if self.Weight:
500                    W2 = np.sqrt(xye[2])
501                    D *= W2-Ymax*.02
502                Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
503                Plot.plot(X,Z,colors[(N+1)%6],picker=False)
504                Plot.plot(X,W,colors[(N+2)%6],picker=False)
505                Plot.plot(X,D,colors[(N+3)%6],picker=False)
506                Plot.axhline(0.,color=wx.BLACK)
507                Page.canvas.SetToolTipString('')
508                if self.PatternTree.GetItemText(PickId) == 'Peak List':
509                    tip = 'On data point: Pick peak - L or R MB.On line: MB down to move'
510                    Page.canvas.SetToolTipString(tip)
511                    data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Peak List'))
512                    for item in data:
513                        Lines.append(Plot.axvline(item[0],color=colors[N%6],picker=2.))
514                if self.PatternTree.GetItemText(PickId) == 'Limits':
515                    tip = 'On data point: Lower limit - L MB; Upper limit - R MB. On limit: MB down to move'
516                    Page.canvas.SetToolTipString(tip)
517                    data = self.LimitsTable.GetData()
518            else:
519                Plot.plot(X,Y,colors[N%6],picker=False)
520    if PickId and self.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
521        peaks = np.array((self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))))
522        for peak in peaks:
523            Plot.axvline(peak[0],color='b')
524        for hkl in self.HKL:
525            Plot.axvline(hkl[5],color='r',dashes=(5,5))
526    if self.Contour:
527        acolor = mpl.cm.get_cmap(self.ContourColor)
528        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*self.Cmax,interpolation=self.Interpolate, 
529            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto')
530        Page.figure.colorbar(Img)
531    else:
532        self.Lines = Lines
533    if not newPlot:
534        Page.toolbar.push_current()
535        Plot.set_xlim(xylim[0])
536        Plot.set_ylim(xylim[1])
537        xylim = []
538        Page.toolbar.push_current()
539        Page.toolbar.draw()
540    else:
541        Page.canvas.draw()
542    self.Pwdr = True
543
544def PlotPowderLines(self):
545    global HKL
546
547    def OnMotion(event):
548        xpos = event.xdata
549        if xpos:                                        #avoid out of frame mouse position
550            Page.canvas.SetCursor(wx.CROSS_CURSOR)
551            self.G2plotNB.status.SetFields(['','2-theta =%9.3f '%(xpos,)])
552            if self.PickId and self.PatternTree.GetItemText(self.PickId) in ['Index Peak List','Unit Cells List']:
553                found = []
554                if len(HKL):
555                    view = Page.toolbar._views.forward()[0][:2]
556                    wid = view[1]-view[0]
557                    found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
558                if len(found):
559                    h,k,l = found[0][:3] 
560                    Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
561                else:
562                    Page.canvas.SetToolTipString('')
563
564    try:
565        plotNum = self.G2plotNB.plotList.index('Powder Lines')
566        Page = self.G2plotNB.nb.GetPage(plotNum)
567        Page.figure.clf()
568        Plot = Page.figure.gca()
569    except ValueError,error:
570        Plot = self.G2plotNB.addMpl('Powder Lines').gca()
571        plotNum = self.G2plotNB.plotList.index('Powder Lines')
572        Page = self.G2plotNB.nb.GetPage(plotNum)
573        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
574       
575    Page.SetFocus()
576    Plot.set_title('Powder Pattern Lines')
577    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
578    PickId = self.PickId
579    PatternId = self.PatternId
580    peaks = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,PatternId, 'Index Peak List'))
581    for peak in peaks:
582        Plot.axvline(peak[0],color='b')
583    HKL = np.array(self.HKL)
584    for hkl in self.HKL:
585        Plot.axvline(hkl[5],color='r',dashes=(5,5))
586    xmin = peaks[0][0]
587    xmax = peaks[-1][0]
588    delt = xmax-xmin
589    xlim = [max(0,xmin-delt/20.),min(180.,xmax+delt/20.)]
590    Plot.set_xlim(xlim)
591    Page.canvas.draw()
592
593def PlotPeakWidths(self):
594    PatternId = self.PatternId
595    limitID = G2gd.GetPatternTreeItemId(self,PatternId, 'Limits')
596    if limitID:
597        limits = self.PatternTree.GetItemPyData(limitID)
598    else:
599        return
600    instParms = self.PatternTree.GetItemPyData( \
601        G2gd.GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))
602    if instParms[0][0] == 'PXC':
603        lam = instParms[1][1]
604        if len(instParms[1]) == 13:
605            GU,GV,GW,LX,LY = instParms[0][6:11]
606        else:
607            GU,GV,GW,LX,LY = instParms[0][4:9]
608    peakID = G2gd.GetPatternTreeItemId(self,PatternId, 'Peak List')
609    if peakID:
610        peaks = self.PatternTree.GetItemPyData(peakID)
611    else:
612        peaks = []
613   
614    try:
615        plotNum = self.G2plotNB.plotList.index('Peak Widths')
616        Page = self.G2plotNB.nb.GetPage(plotNum)
617        Page.figure.clf()
618        Plot = Page.figure.gca()
619    except ValueError,error:
620        Plot = self.G2plotNB.addMpl('Peak Widths').gca()
621        plotNum = self.G2plotNB.plotList.index('Peak Widths')
622        Page = self.G2plotNB.nb.GetPage(plotNum)
623    Page.SetFocus()
624   
625    Page.canvas.SetToolTipString('')
626    colors=['b','g','r','c','m','k']
627    Xmin,Xmax = limits[1]
628    Xmin = min(0.5,max(Xmin,1))
629    Xmin /= 2
630    Xmax /= 2
631    nPts = 100
632    delt = (Xmax-Xmin)/nPts
633    thetas = []
634    for i in range(nPts):
635        thetas.append(Xmin+i*delt)
636    X = []
637    Y = []
638    Z = []
639    W = []
640    V = []
641    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
642    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
643    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.)
644    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
645    for theta in thetas:
646        X.append(4.0*math.pi*sind(theta)/lam)              #q
647        s = sig(theta,GU,GV,GW)
648        g = gam(theta,LX,LY)
649        G = gamFW(g,s)
650        H = gamFW2(g,s)
651        Y.append(s/tand(theta))
652        Z.append(g/tand(theta))
653        W.append(G/tand(theta))
654        V.append(H/tand(theta))
655    Plot.set_title('Instrument and sample peak widths')
656    Plot.set_ylabel(r'$\Delta q/q, \Delta d/d$',fontsize=14)
657    Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
658    Plot.plot(X,Y,color='r',label='Gaussian')
659    Plot.plot(X,Z,color='g',label='Lorentzian')
660    Plot.plot(X,W,color='b',label='G+L')
661    Plot.plot(X,V,color='k',label='G+L2')
662    X = []
663    Y = []
664    Z = []
665    W = []
666    V = []
667    for peak in peaks:
668        X.append(4.0*math.pi*sind(peak[0]/2.0)/lam)
669        s = 1.17741*math.sqrt(peak[4])*math.pi/18000.
670        g = peak[6]*math.pi/18000.
671        G = gamFW(g,s)
672        H = gamFW2(g,s)
673        Y.append(s/tand(peak[0]/2.))
674        Z.append(g/tand(peak[0]/2.))
675        W.append(G/tand(peak[0]/2.))
676        V.append(H/tand(peak[0]/2.))
677    Plot.plot(X,Y,'+',color='r',label='G peak')
678    Plot.plot(X,Z,'+',color='g',label='L peak')
679    Plot.plot(X,W,'+',color='b',label='G+L peak')
680    Plot.plot(X,V,'+',color='k',label='G+L2 peak')
681    Plot.legend(loc='best')
682    Page.canvas.draw()
683           
684def PlotExposedImage(self,newPlot=False,event=None):
685    plotNo = self.G2plotNB.nb.GetSelection()
686    if self.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
687        PlotImage(self,newPlot,event)
688    elif self.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
689        PlotIntegration(self,newPlot,event)
690
691def PlotImage(self,newPlot=False,event=None):
692    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
693    import numpy.ma as ma
694    Dsp = lambda tth,wave: wave/(2.*sind(tth/2.))
695
696    def OnImMotion(event):
697        Data = self.PatternTree.GetItemPyData(
698            G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
699        Page.canvas.SetToolTipString('')
700        size = len(self.ImageZ)
701        if event.xdata and event.ydata:                 #avoid out of frame errors
702            Page.canvas.SetCursor(wx.CROSS_CURSOR)
703            item = self.itemPicked
704            pixelSize = Data['pixelSize']
705            scalex = 1000./pixelSize[0]
706            scaley = 1000./pixelSize[1]                   
707            if item and self.PatternTree.GetItemText(self.PickId) == 'Image Controls':
708                if 'Text' in str(item):
709                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
710                else:
711                    xcent,ycent = Data['center']
712                    xpos = event.xdata-xcent
713                    ypos = event.ydata-ycent
714                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
715                        ang = int(atan2d(xpos,ypos))
716                        Page.canvas.SetToolTipString('%6d deg'%(ang))
717                    elif 'line1' in  str(item) or 'line2' in str(item):
718                        tth = G2img.GetTth(event.xdata,event.ydata,Data)
719                        Page.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
720            else:
721                xpos = event.xdata
722                ypos = event.ydata
723                xpix = xpos*scalex
724                ypix = ypos*scaley
725                Int = 0
726                if (0 <= xpix <= size) and (0 <= ypix <= size):
727                    Int = self.ImageZ[ypix][xpix]
728                tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
729                Q = 2.*math.pi/dsp
730                self.G2plotNB.status.SetFields(\
731                    ['','Detector 2-th =%9.2fdeg, dsp =%9.3fA, Q = %6.3fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)])
732
733    def OnImPlotKeyPress(event):
734        if self.PatternTree.GetItemText(self.PickId) == 'Masks':
735            Data = self.PatternTree.GetItemPyData(
736                G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
737            Masks = self.PatternTree.GetItemPyData(
738                G2gd.GetPatternTreeItemId(self,self.Image, 'Masks'))
739            Xpos = event.xdata
740            if not Xpos:            #got point out of frame
741                return
742            Ypos = event.ydata
743            if event.key == 's':
744                print 'spot mask @ ',Xpos,Ypos
745                Masks['Points'].append([Xpos,Ypos,1.])
746            elif event.key == 'r':
747                tth = G2img.GetTth(Xpos,Ypos,Data)
748                print 'ring mask @ ',Xpos,Ypos,tth
749                Masks['Rings'].append([tth,0.1])
750            elif event.key == 'a':
751                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
752                azm = int(azm)               
753                print 'arc mask @ ', Xpos,Ypos
754                Masks['Arcs'].append([tth,[azm-5,azm+5],0.1])
755            elif event.key == 'p':
756                self.setPoly = True
757                Masks['Polygons'].append([])
758                self.G2plotNB.status.SetFields(['','Polygon mask active - LB pick points, RB closepolygon'])
759            G2imG.UpdateMasks(self,Masks)
760        PlotImage(self)
761           
762    def OnImPick(event):
763        if self.PatternTree.GetItemText(self.PickId) not in ['Image Controls','Masks']:
764            return
765        if self.setPoly:
766            Masks = self.PatternTree.GetItemPyData(
767                G2gd.GetPatternTreeItemId(self,self.Image, 'Masks'))
768            polygon = Masks['Polygons'][-1]
769            xpos,ypos = event.mouseevent.xdata,event.mouseevent.ydata
770            if xpos and ypos:                       #point inside image
771                if len(polygon) > 2 and event.mouseevent.button == 3:
772                    x0,y0 = polygon[0]
773                    polygon.append([x0,y0])
774                    self.setPoly = False
775                    self.G2plotNB.status.SetFields(['','Polygon closed: %.1f,%.1f'%(x0,y0)])
776                else:
777                    self.G2plotNB.status.SetFields(['','Add polygon point: %.1f,%.1f'%(xpos,ypos)])
778                    polygon.append([xpos,ypos])
779                G2imG.UpdateMasks(self,Masks)
780        else:
781            if self.itemPicked is not None: return
782            self.itemPicked = event.artist
783            self.mousePicked = event.mouseevent
784       
785    def OnImRelease(event):
786        PickName = self.PatternTree.GetItemText(self.PickId)
787        if PickName not in ['Image Controls','Masks']:
788            return
789        Data = self.PatternTree.GetItemPyData(
790            G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
791        Masks = self.PatternTree.GetItemPyData(
792            G2gd.GetPatternTreeItemId(self,self.Image, 'Masks'))
793        pixelSize = Data['pixelSize']
794        scalex = 1000./pixelSize[0]
795        scaley = 1000./pixelSize[1]
796        pixLimit = Data['pixLimit']
797        if self.itemPicked is None and PickName == 'Image Controls':
798            size = len(self.ImageZ)
799            Xpos = event.xdata
800            if not (Xpos and self.ifGetRing):                   #got point out of frame
801                return
802            Ypos = event.ydata
803            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
804                if event.button == 1:
805                    Xpix = Xpos*scalex
806                    Ypix = Ypos*scaley
807                    xpos,ypos,I,J = G2img.ImageLocalMax(self.ImageZ,pixLimit,Xpix,Ypix)
808                    if I and J:
809                        xpos += .5                              #shift to pixel center
810                        ypos += .5
811                        xpos /= scalex                          #convert to mm
812                        ypos /= scaley
813                        Data['ring'].append([xpos,ypos])
814                PlotImage(self)
815            return
816        else:
817            xpos = event.xdata
818            if xpos:                                        #avoid out of frame mouse position
819                ypos = event.ydata
820                if self.ifGetRing:
821                    xypos = [xpos,ypos]
822                    rings = Data['ring']
823                    for ring in rings:
824                        if np.allclose(ring,xypos,.01,0):
825                            rings.remove(ring)                                                                       
826                else:
827                    tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
828                    itemPicked = str(self.itemPicked)
829                    if 'Line2D' in itemPicked and PickName == 'Image Controls':
830                        print int(itemPicked.split('_line')[1].strip(')'))
831                        if 'line1' in itemPicked:
832                            Data['IOtth'][0] = tth
833                        elif 'line2' in itemPicked:
834                            Data['IOtth'][1] = tth
835                        elif 'line3' in itemPicked and not Data['fullIntegrate']:
836                            Data['LRazimuth'][0] = int(azm)
837                        elif 'line4' in itemPicked and not Data['fullIntegrate']:
838                            Data['LRazimuth'][1] = int(azm)
839                           
840                        if Data['LRazimuth'][1] < Data['LRazimuth'][0]:
841                            Data['LRazimuth'][1] += 360
842                        if  Data['IOtth'][0] > Data['IOtth'][1]:
843                            Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
844                           
845                        self.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
846                        self.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
847                        self.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
848                        self.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
849                    elif 'Circle' in itemPicked and PickName == 'Masks':
850                        spots = Masks['Points']
851                        newPos = itemPicked.split(')')[0].split('(')[2].split(',')
852                        newPos = np.array([float(newPos[0]),float(newPos[1])])
853                        for spot in spots:
854                            if np.allclose(np.array([spot[:2]]),newPos):
855                                spot[:2] = xpos,ypos
856                        G2imG.UpdateMasks(self,Masks)
857                    elif 'Line2D' in itemPicked and PickName == 'Masks':
858                        Obj = self.itemPicked.findobj()
859                        rings = Masks['Rings']
860                        arcs = Masks['Arcs']
861                        polygons = Masks['Polygons']
862                        for ring in self.ringList:
863                            if Obj == ring[0]:
864                                rN = ring[1]
865                                if ring[2] == 'o':
866                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)-rings[rN][1]/2.
867                                else:
868                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)+rings[rN][1]/2.
869                        for arc in self.arcList:
870                            if Obj == arc[0]:
871                                aN = arc[1]
872                                if arc[2] == 'o':
873                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)-arcs[aN][2]/2
874                                elif arc[2] == 'i':
875                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)+arcs[aN][2]/2
876                                elif arc[2] == 'l':
877                                    arcs[aN][1][0] = int(G2img.GetAzm(xpos,ypos,Data))
878                                else:
879                                    arcs[aN][1][1] = int(G2img.GetAzm(xpos,ypos,Data))
880                        for poly in self.polyList:
881                            if Obj == poly[0]:
882                                ind = self.itemPicked.contains(self.mousePicked)[1]['ind'][0]
883                                oldPos = np.array([self.mousePicked.xdata,self.mousePicked.ydata])
884                                pN = poly[1]
885                                for i,xy in enumerate(polygons[pN]):
886                                    if np.allclose(np.array([xy]),oldPos,atol=1.0):
887                                        polygons[pN][i] = xpos,ypos
888                        G2imG.UpdateMasks(self,Masks)
889#                    else:                  #keep for future debugging
890#                        print str(self.itemPicked),event.xdata,event.ydata,event.button
891                PlotImage(self)
892            self.itemPicked = None
893           
894    try:
895        plotNum = self.G2plotNB.plotList.index('2D Powder Image')
896        Page = self.G2plotNB.nb.GetPage(plotNum)
897        if not newPlot:
898            Plot = Page.figure.gca()          #get previous powder plot & get limits
899            xylim = Plot.get_xlim(),Plot.get_ylim()
900        Page.figure.clf()
901        Plot = Page.figure.gca()          #get a fresh plot after clf()
902       
903    except ValueError,error:
904        Plot = self.G2plotNB.addMpl('2D Powder Image').gca()
905        plotNum = self.G2plotNB.plotList.index('2D Powder Image')
906        Page = self.G2plotNB.nb.GetPage(plotNum)
907        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
908        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
909        Page.canvas.mpl_connect('pick_event', OnImPick)
910        Page.canvas.mpl_connect('button_release_event', OnImRelease)
911    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
912        Page.SetFocus()
913    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
914    size,imagefile = self.PatternTree.GetItemPyData(self.Image)
915    if imagefile != self.oldImagefile:
916        self.ImageZ = G2IO.GetImageData(imagefile,imageOnly=True)
917        self.oldImagefile = imagefile
918    Data = self.PatternTree.GetItemPyData(
919        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
920    Masks = self.PatternTree.GetItemPyData(
921        G2gd.GetPatternTreeItemId(self,self.Image, 'Masks'))
922
923    imScale = 1
924    if len(self.ImageZ) > 1024:
925        imScale = len(self.ImageZ)/1024
926    pixelSize = Data['pixelSize']
927    scalex = 1000./pixelSize[0]
928    scaley = 1000./pixelSize[1]
929    xmax = len(self.ImageZ)
930    Xmax = len(self.ImageZ)*pixelSize[0]/1000.
931    xlim = (-0.5,Xmax-.5)
932    ylim = (Xmax-.5,-0.5,)
933    Imin,Imax = Data['range'][1]
934    acolor = mpl.cm.get_cmap(Data['color'])
935    xcent,ycent = Data['center']
936    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
937    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
938    #do threshold mask - "real" mask - others are just bondaries
939    Zlim = Masks['Thresholds'][1]
940    wx.BeginBusyCursor()
941    try:
942        MA = ma.masked_greater(ma.masked_less(self.ImageZ,Zlim[0]),Zlim[1])
943        MaskA = ma.getmaskarray(MA)
944        A = G2img.ImageCompress(MA,imScale)
945        AM = G2img.ImageCompress(MaskA,imScale)
946       
947        ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
948            interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Xmax,0])
949        Img = Plot.imshow(A,aspect='equal',cmap=acolor,
950            interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Xmax,0])
951        if self.setPoly:
952            Img.set_picker(True)
953   
954        Plot.plot(xcent,ycent,'x')
955        if Data['showLines']:
956            LRAzim = Data['LRazimuth']                  #NB: integers
957            IOtth = Data['IOtth']
958            wave = Data['wavelength']
959            dspI = wave/(2.0*sind(IOtth[0]/2.0))
960            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
961            dspO = wave/(2.0*sind(IOtth[1]/2.0))
962            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
963            if Data['fullIntegrate']:
964                Azm = np.array(range(0,361))
965            else:
966                Azm = np.array(range(LRAzim[0],LRAzim[1]+1))
967            if ellI:
968                xyI = []
969                for azm in Azm:
970                    xyI.append(G2img.GetDetectorXY(dspI,azm,Data))
971                xyI = np.array(xyI)
972                arcxI,arcyI = xyI.T
973                Plot.plot(arcxI,arcyI,picker=3)
974            if ellO:
975                xyO = []
976                for azm in Azm:
977                    xyO.append(G2img.GetDetectorXY(dspO,azm,Data))
978                xyO = np.array(xyO)
979                arcxO,arcyO = xyO.T
980                Plot.plot(arcxO,arcyO,picker=3)
981            if ellO and ellI and not Data['fullIntegrate']:
982                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
983                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
984        for xring,yring in Data['ring']:
985            Plot.plot(xring,yring,'r+',picker=3)
986        if Data['setRings']:
987            rings = np.concatenate((Data['rings']),axis=0)
988            for xring,yring,dsp in rings:
989                Plot.plot(xring,yring,'r+')           
990        for ellipse in Data['ellipses']:
991            cent,phi,[width,height],col = ellipse
992            Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
993            Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
994        #masks - mask lines numbered after integration limit lines
995        spots = Masks['Points']
996        rings = Masks['Rings']
997        arcs = Masks['Arcs']
998        polygons = Masks['Polygons']
999        for x,y,d in spots:
1000            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
1001        self.ringList = []
1002        for iring,(tth,thick) in enumerate(rings):
1003            wave = Data['wavelength']
1004            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
1005            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
1006            self.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
1007            self.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
1008        self.arcList = []
1009        for iarc,(tth,azm,thick) in enumerate(arcs):           
1010            wave = Data['wavelength']
1011            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
1012            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
1013            self.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
1014            self.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
1015            self.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
1016            self.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
1017        self.polyList = []
1018        for ipoly,polygon in enumerate(polygons):
1019            x,y = np.hsplit(np.array(polygon),2)
1020            self.polyList.append([Plot.plot(x,y,'r',picker=3),ipoly])           
1021        colorBar = Page.figure.colorbar(Img)
1022        Plot.set_xlim(xlim)
1023        Plot.set_ylim(ylim)
1024        if not newPlot:
1025            Page.toolbar.push_current()
1026            Plot.set_xlim(xylim[0])
1027            Plot.set_ylim(xylim[1])
1028            xylim = []
1029            Page.toolbar.push_current()
1030            Page.toolbar.draw()
1031        else:
1032            Page.canvas.draw()
1033    finally:
1034        wx.EndBusyCursor()
1035       
1036def PlotIntegration(self,newPlot=False,event=None):
1037           
1038    def OnMotion(event):
1039        Page.canvas.SetToolTipString('')
1040        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1041        azm = event.ydata
1042        tth = event.xdata
1043        if azm and tth:
1044            self.G2plotNB.status.SetFields(\
1045                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1046                               
1047    try:
1048        plotNum = self.G2plotNB.plotList.index('2D Integration')
1049        Page = self.G2plotNB.nb.GetPage(plotNum)
1050        if not newPlot:
1051            Plot = Page.figure.gca()          #get previous plot & get limits
1052            xylim = Plot.get_xlim(),Plot.get_ylim()
1053        Page.figure.clf()
1054        Plot = Page.figure.gca()          #get a fresh plot after clf()
1055       
1056    except ValueError,error:
1057        Plot = self.G2plotNB.addMpl('2D Integration').gca()
1058        plotNum = self.G2plotNB.plotList.index('2D Integration')
1059        Page = self.G2plotNB.nb.GetPage(plotNum)
1060        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1061        Page.views = False
1062        view = False
1063    if not event:
1064        Page.SetFocus()
1065       
1066    Data = self.PatternTree.GetItemPyData(
1067        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1068    image = self.Integrate[0]
1069    xsc = self.Integrate[1]
1070    ysc = self.Integrate[2]
1071    Imin,Imax = Data['range'][1]
1072    acolor = mpl.cm.get_cmap(Data['color'])
1073    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1074    Plot.set_ylabel('azimuth',fontsize=12)
1075    Plot.set_xlabel('2-theta',fontsize=12)
1076    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
1077        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
1078    colorBar = Page.figure.colorbar(Img)
1079    if Data['setRings']:
1080        rings = np.concatenate((Data['rings']),axis=0)
1081        for xring,yring,dsp in rings:
1082            x,y = G2img.GetTthAzm(xring,yring,Data)
1083            Plot.plot(x,y,'r+')
1084    if Data['ellipses']:           
1085        for ellipse in Data['ellipses']:
1086            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1087            x,y = np.hsplit(ring,2)
1088            tth,azm = G2img.GetTthAzm(x,y,Data)
1089            Plot.plot(tth,azm,'b,')
1090    if not newPlot:
1091        Page.toolbar.push_current()
1092        Plot.set_xlim(xylim[0])
1093        Plot.set_ylim(xylim[1])
1094        xylim = []
1095        Page.toolbar.push_current()
1096        Page.toolbar.draw()
1097    else:
1098        Page.canvas.draw()
1099               
1100def PlotTRImage(self,tax,tay,taz,newPlot=False):
1101    #a test plot routine - not normally used
1102           
1103    def OnMotion(event):
1104        Page.canvas.SetToolTipString('')
1105        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1106        azm = event.xdata
1107        tth = event.ydata
1108        if azm and tth:
1109            self.G2plotNB.status.SetFields(\
1110                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1111                               
1112    try:
1113        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1114        Page = self.G2plotNB.nb.GetPage(plotNum)
1115        if not newPlot:
1116            Plot = Page.figure.gca()          #get previous plot & get limits
1117            xylim = Plot.get_xlim(),Plot.get_ylim()
1118        Page.figure.clf()
1119        Plot = Page.figure.gca()          #get a fresh plot after clf()
1120       
1121    except ValueError,error:
1122        Plot = self.G2plotNB.addMpl('2D Transformed Powder Image').gca()
1123        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1124        Page = self.G2plotNB.nb.GetPage(plotNum)
1125        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1126        Page.views = False
1127        view = False
1128    Page.SetFocus()
1129       
1130    Data = self.PatternTree.GetItemPyData(
1131        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1132    Imin,Imax = Data['range'][1]
1133    step = (Imax-Imin)/5.
1134    V = np.arange(Imin,Imax,step)
1135    acolor = mpl.cm.get_cmap(Data['color'])
1136    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1137    Plot.set_xlabel('azimuth',fontsize=12)
1138    Plot.set_ylabel('2-theta',fontsize=12)
1139    Plot.contour(tax,tay,taz,V,cmap=acolor)
1140    if Data['showLines']:
1141        IOtth = Data['IOtth']
1142        if Data['fullIntegrate']:
1143            LRAzim = [-180,180]
1144        else:
1145            LRAzim = Data['LRazimuth']                  #NB: integers
1146        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
1147        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
1148        if not Data['fullIntegrate']:
1149            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
1150            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
1151    if Data['setRings']:
1152        rings = np.concatenate((Data['rings']),axis=0)
1153        for xring,yring,dsp in rings:
1154            x,y = G2img.GetTthAzm(xring,yring,Data)
1155            Plot.plot(y,x,'r+')           
1156    if Data['ellipses']:           
1157        for ellipse in Data['ellipses']:
1158            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1159            x,y = np.hsplit(ring,2)
1160            tth,azm = G2img.GetTthAzm(x,y,Data)
1161            Plot.plot(azm,tth,'b,')
1162    if not newPlot:
1163        Page.toolbar.push_current()
1164        Plot.set_xlim(xylim[0])
1165        Plot.set_ylim(xylim[1])
1166        xylim = []
1167        Page.toolbar.push_current()
1168        Page.toolbar.draw()
1169    else:
1170        Page.canvas.draw()
1171       
1172def PlotStructure(self,data):
1173    generalData = data['General']
1174    cell = generalData['Cell'][1:7]
1175    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
1176    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1177    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1178    Mydir = generalData['Mydir']
1179    atomData = data['Atoms']
1180    drawingData = data['Drawing']
1181    drawAtoms = drawingData['Atoms']
1182    cx,ct,cs = drawingData['atomPtrs']
1183    Wt = [255,255,255]
1184    Rd = [255,0,0]
1185    Gr = [0,255,0]
1186    Bl = [0,0,255]
1187    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]])
1188    uEdges = np.array([
1189        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
1190        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
1191        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
1192    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
1193    altDown = False
1194    shiftDown = False
1195    ctrlDown = False
1196   
1197    def OnKeyBox(event):
1198        import Image
1199        Draw()                          #make sure plot is fresh!!
1200        mode = cb.GetValue()
1201        Fname = Mydir+'\\'+generalData['Name']+'.'+mode
1202        size = Page.canvas.GetSize()
1203        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
1204        if mode in ['jpeg',]:
1205            Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
1206            im = Image.new("RGBA", (size[0],size[1]))
1207        else:
1208            Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
1209            im = Image.new("RGB", (size[0],size[1]))
1210        im.fromstring(Pix)
1211        im.save(Fname,mode)
1212        cb.SetValue(' Save as:')
1213        self.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
1214   
1215    def GetTruePosition(xy):
1216        View = glGetIntegerv(GL_VIEWPORT)
1217        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
1218        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
1219        Zmax = 1.
1220        for i,atom in enumerate(drawAtoms):
1221            x,y,z = atom[cx:cx+3]
1222            X,Y,Z = gluProject(x,y,z,Model,Proj,View)
1223            XY = [int(X),int(View[3]-Y)]
1224            if np.allclose(xy,XY,atol=10) and Z < Zmax:
1225                Zmax = Z
1226                SetSelectedAtoms(i)
1227                   
1228    def OnMouseDown(event):
1229        xy = event.GetPosition()
1230        if event.ShiftDown():
1231            GetTruePosition(xy)
1232        else:
1233            drawingData['Rotation'][3] = xy
1234        Draw()
1235       
1236    def OnMouseMove(event):
1237        newxy = event.GetPosition()
1238        if event.ControlDown() and drawingData['showABC']:
1239            if event.LeftIsDown():
1240                ctrlDown = True
1241                SetTestRot(newxy)
1242            elif event.RightIsDown():
1243                SetTestPos(newxy)
1244            elif event.MiddleIsDown():
1245                SetTestRotZ(newxy)
1246               
1247               
1248        if event.Dragging() and not event.ControlDown():
1249            if event.LeftIsDown():
1250                SetRotation(newxy)
1251            elif event.RightIsDown():
1252                SetTranslation(newxy)
1253            elif event.MiddleIsDown():
1254                SetRotationZ(newxy)
1255        Draw()
1256       
1257    def OnMouseWheel(event):
1258        drawingData['cameraPos'] += event.GetWheelRotation()/24
1259        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
1260        page = self.dataDisplay.GetSelection()
1261        if self.dataDisplay.GetPageText(page) == 'Draw Options':
1262            panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1263            names = [child.GetName() for child in panel]
1264            panel[names.index('cameraPos')].SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
1265            panel[names.index('cameraSlider')].SetValue(drawingData['cameraPos'])
1266        Draw()
1267           
1268    def SetViewPointText(VP):
1269        page = self.dataDisplay.GetSelection()
1270        if self.dataDisplay.GetPageText(page) == 'Draw Options':
1271            panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1272            names = [child.GetName() for child in panel]
1273            panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
1274           
1275    def ClearSelectedAtoms():
1276        page = self.dataDisplay.GetSelection()
1277        if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1278            self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
1279        elif self.dataDisplay.GetPageText(page) == 'Atoms':
1280            self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
1281                   
1282    def SetSelectedAtoms(ind):
1283        page = self.dataDisplay.GetSelection()
1284        if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1285            self.dataDisplay.GetPage(page).SelectRow(ind)      #this is the Atoms grid in Draw Atoms
1286        elif self.dataDisplay.GetPageText(page) == 'Atoms':
1287            Id = drawAtoms[ind][-2]
1288            for i,atom in enumerate(atomData):
1289                if atom[-1] == Id:
1290                    self.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
1291                 
1292    def GetSelectedAtoms():
1293        page = self.dataDisplay.GetSelection()
1294        Ind = []
1295        if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1296            Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
1297        elif self.dataDisplay.GetPageText(page) == 'Atoms':
1298            Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
1299        return Ind
1300                                       
1301    def OnKey(event):           #on key UP!!
1302        keyCode = event.GetKeyCode()
1303        if keyCode > 255:
1304            keyCode = 0
1305        key,xyz = chr(keyCode),event.GetPosition()
1306        indx = drawingData['selectedAtoms']
1307        if key in ['c','C']:
1308            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
1309            drawingData['testPos'] = [[-.1,-.1,-.1],[0.0,0.0,0.0],[0,0]]
1310            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
1311            SetViewPointText(drawingData['viewPoint'][0])
1312        elif key in ['n','N']:
1313            drawAtoms = drawingData['Atoms']
1314            pI = drawingData['viewPoint'][1]
1315            if indx:
1316                pI[0] = indx[pI[1]]
1317                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1318                pI[1] += 1
1319                if pI[1] >= len(indx):
1320                    pI[1] = 0
1321            else:
1322                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1323                pI[0] += 1
1324                if pI[0] >= len(drawAtoms):
1325                    pI[0] = 0
1326            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1327            SetViewPointText(drawingData['viewPoint'][0])
1328               
1329        elif key in ['p','P']:
1330            drawAtoms = drawingData['Atoms']
1331            pI = drawingData['viewPoint'][1]
1332            if indx:
1333                pI[0] = indx[pI[1]]
1334                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1335                pI[1] -= 1
1336                if pI[1] < 0:
1337                    pI[1] = len(indx)-1
1338            else:
1339                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1340                pI[0] -= 1
1341                if pI[0] < 0:
1342                    pI[0] = len(drawAtoms)-1
1343            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1344            SetViewPointText(drawingData['viewPoint'][0])           
1345        Draw()
1346           
1347    def SetBackground():
1348        R,G,B,A = Page.camera['backColor']
1349        glClearColor(R,G,B,A)
1350        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
1351       
1352    def SetLights():
1353        glEnable(GL_DEPTH_TEST)
1354        glShadeModel(GL_SMOOTH)
1355        glEnable(GL_LIGHTING)
1356        glEnable(GL_LIGHT0)
1357        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
1358        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
1359        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
1360       
1361    def SetTranslation(newxy):
1362        Tx,Ty,Tz = drawingData['viewPoint'][0]
1363        anglex,angley,anglez,oldxy = drawingData['Rotation']
1364        Rx = G2lat.rotdMat(anglex,0)
1365        Ry = G2lat.rotdMat(angley,1)
1366        Rz = G2lat.rotdMat(anglez,2)
1367        dxy = list(newxy-oldxy)+[0,]
1368        dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
1369        Tx -= dxy[0]*0.01
1370        Ty += dxy[1]*0.01
1371        Tz -= dxy[2]*0.01
1372        drawingData['Rotation'][3] = newxy
1373        drawingData['viewPoint'][0] =  Tx,Ty,Tz
1374        SetViewPointText([Tx,Ty,Tz])
1375       
1376    def SetTestPos(newxy):
1377        Tx,Ty,Tz = drawingData['testPos'][0]
1378        anglex,angley,anglez,oldxy = drawingData['Rotation']
1379        Rx = G2lat.rotdMat(anglex,0)
1380        Ry = G2lat.rotdMat(angley,1)
1381        Rz = G2lat.rotdMat(anglez,2)
1382        dxy = list(newxy-oldxy)+[0,]
1383        dxy = np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy)))
1384        Tx += dxy[0]*0.001
1385        Ty -= dxy[1]*0.001
1386        Tz += dxy[2]*0.001
1387        drawingData['Rotation'][3] = newxy
1388        drawingData['testPos'][0] =  Tx,Ty,Tz
1389       
1390    def SetTestRot(newxy):
1391        Txyz = np.array(drawingData['testPos'][0])
1392        oldxy = drawingData['testPos'][2]
1393        Ax,Ay,Az = drawingData['testPos'][1]
1394        Vxyz = np.array(drawingData['viewPoint'][0])
1395        Dxyz = np.inner(Amat,Txyz-Vxyz)
1396        dxy = list(newxy-oldxy)+[0,]
1397        Ax += dxy[1]*0.01
1398        Ay += dxy[0]*0.01
1399        Rx = G2lat.rotdMat(Ax,0)
1400        Ry = G2lat.rotdMat(Ay,1)
1401        Dxyz = np.inner(Ry,np.inner(Rx,Dxyz))       
1402        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
1403        drawingData['testPos'][1] = [Ax,Ay,Az]
1404        drawingData['testPos'][2] = newxy
1405        drawingData['testPos'][0] = Dxyz
1406       
1407    def SetTestRotZ(newxy):
1408        Txyz = np.array(drawingData['testPos'][0])
1409        oldxy = drawingData['testPos'][2]
1410        Ax,Ay,Az = drawingData['testPos'][1]
1411        Vxyz = np.array(drawingData['viewPoint'][0])
1412        Dxyz = np.inner(Amat,Txyz-Vxyz)       
1413        dxy = list(newxy-oldxy)+[0,]
1414        Az += (dxy[0]+dxy[1])*.01
1415        Rz = G2lat.rotdMat(Az,2)
1416        Dxyz = np.inner(Rz,Dxyz)       
1417        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
1418        drawingData['testPos'][1] = [Ax,Ay,Az]
1419        drawingData['testPos'][2] = newxy
1420        drawingData['testPos'][0] = Dxyz
1421                             
1422    def SetRotation(newxy):       
1423        anglex,angley,anglez,oldxy = drawingData['Rotation']
1424        dxy = newxy-oldxy
1425        anglex += dxy[1]*.25
1426        angley += dxy[0]*.25
1427        oldxy = newxy
1428        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
1429       
1430    def SetRotationZ(newxy):                       
1431        anglex,angley,anglez,oldxy = drawingData['Rotation']
1432        dxy = newxy-oldxy
1433        anglez += (dxy[0]+dxy[1])*.25
1434        oldxy = newxy
1435        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
1436       
1437    def RenderBox():
1438        glEnable(GL_COLOR_MATERIAL)
1439        glLineWidth(2)
1440        glEnable(GL_BLEND)
1441        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
1442        glEnable(GL_LINE_SMOOTH)
1443        glBegin(GL_LINES)
1444        for line,color in zip(uEdges,uColors):
1445            glColor3ubv(color)
1446            glVertex3fv(line[0])
1447            glVertex3fv(line[1])
1448        glEnd()
1449        glColor4ubv([0,0,0,0])
1450        glDisable(GL_LINE_SMOOTH)
1451        glDisable(GL_BLEND)
1452        glDisable(GL_COLOR_MATERIAL)
1453       
1454    def RenderUnitVectors(x,y,z):
1455        xyz = np.array([x,y,z])
1456        glEnable(GL_COLOR_MATERIAL)
1457        glLineWidth(1)
1458        glPushMatrix()
1459        glTranslate(x,y,z)
1460        glScalef(1/cell[0],1/cell[1],1/cell[2])
1461        glBegin(GL_LINES)
1462        for line,color in zip(uEdges,uColors)[:3]:
1463            glColor3ubv(color)
1464            glVertex3fv(line[0])
1465            glVertex3fv(line[1])
1466        glEnd()
1467        glPopMatrix()
1468        glColor4ubv([0,0,0,0])
1469        glDisable(GL_COLOR_MATERIAL)
1470               
1471    def RenderSphere(x,y,z,radius,color):
1472        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1473        glPushMatrix()
1474        glTranslate(x,y,z)
1475        glMultMatrixf(B4mat.T)
1476        q = gluNewQuadric()
1477        gluSphere(q,radius,20,10)
1478        glPopMatrix()
1479       
1480    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
1481        s1,s2,s3 = E
1482        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1483        glPushMatrix()
1484        glTranslate(x,y,z)
1485        glMultMatrixf(B4mat.T)
1486        glMultMatrixf(R4.T)
1487        glEnable(GL_NORMALIZE)
1488        glScale(s1,s2,s3)
1489        q = gluNewQuadric()
1490        gluSphere(q,ellipseProb,20,10)
1491        glDisable(GL_NORMALIZE)
1492        glPopMatrix()
1493       
1494    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
1495        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1496        glPushMatrix()
1497        glTranslate(x,y,z)
1498        glMultMatrixf(B4mat.T)
1499        for bond in Bonds:
1500            glPushMatrix()
1501            Dx = np.inner(Amat,bond)
1502            Z = np.sqrt(np.sum(Dx**2))
1503            azm = atan2d(-Dx[1],-Dx[0])
1504            phi = acosd(Dx[2]/Z)
1505            glRotate(-azm,0,0,1)
1506            glRotate(phi,1,0,0)
1507            q = gluNewQuadric()
1508            gluCylinder(q,radius,radius,Z,slice,2)
1509            glPopMatrix()           
1510        glPopMatrix()
1511               
1512    def RenderLines(x,y,z,Bonds,color):
1513        xyz = np.array([x,y,z])
1514        glEnable(GL_COLOR_MATERIAL)
1515        glLineWidth(1)
1516        glColor3fv(color)
1517        glPushMatrix()
1518        glBegin(GL_LINES)
1519        for bond in Bonds:
1520            glVertex3fv(xyz)
1521            glVertex3fv(xyz+bond)
1522        glEnd()
1523        glColor4ubv([0,0,0,0])
1524        glPopMatrix()
1525        glDisable(GL_COLOR_MATERIAL)
1526       
1527    def RenderPolyhedra(x,y,z,Faces,color):
1528        glPushMatrix()
1529        glTranslate(x,y,z)
1530        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1531        glShadeModel(GL_SMOOTH)
1532        glMultMatrixf(B4mat.T)
1533        for face,norm in Faces:
1534            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
1535            glFrontFace(GL_CW)
1536            glNormal3fv(norm)
1537            glBegin(GL_TRIANGLES)
1538            for vert in face:
1539                glVertex3fv(vert)
1540            glEnd()
1541        glPopMatrix()
1542       
1543    def RenderBackbone(Backbone,BackboneColor,radius):
1544        glPushMatrix()
1545        glMultMatrixf(B4mat.T)
1546        glEnable(GL_COLOR_MATERIAL)
1547        glShadeModel(GL_SMOOTH)
1548        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
1549        glePolyCylinder(Backbone,BackboneColor,radius)
1550        glPopMatrix()       
1551        glDisable(GL_COLOR_MATERIAL)
1552       
1553    def RenderLabel(x,y,z,label,r):       
1554        glPushMatrix()
1555        glTranslate(x,y,z)
1556        glMultMatrixf(B4mat.T)
1557        glDisable(GL_LIGHTING)
1558        glColor3f(0,1.,0)
1559        glRasterPos3f(r,r,r)
1560        for c in list(label):
1561            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
1562        glEnable(GL_LIGHTING)
1563        glPopMatrix()
1564                           
1565    def Draw():
1566        import numpy.linalg as nl
1567        Ind = GetSelectedAtoms()
1568        VS = np.array(Page.canvas.GetSize())
1569        aspect = float(VS[0])/float(VS[1])
1570        cPos = drawingData['cameraPos']
1571        Zclip = drawingData['Zclip']*cPos/200.
1572        anglex,angley,anglez = drawingData['Rotation'][:3]
1573        Tx,Ty,Tz = drawingData['viewPoint'][0]
1574        cx,ct,cs = drawingData['atomPtrs']
1575        bondR = drawingData['bondRadius']
1576        G,g = G2lat.cell2Gmat(cell)
1577        GS = G
1578        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
1579        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
1580        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
1581        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
1582       
1583        SetBackground()
1584        glInitNames()
1585        glPushName(0)
1586       
1587        glMatrixMode(GL_PROJECTION)
1588        glLoadIdentity()
1589        glViewport(0,0,VS[0],VS[1])
1590        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
1591        gluLookAt(0,0,cPos,0,0,0,0,1,0)
1592        SetLights()           
1593           
1594        glMatrixMode(GL_MODELVIEW)
1595        glLoadIdentity()
1596        glRotate(anglez,0,0,1)
1597        glRotate(anglex,cosd(anglez),-sind(anglez),0)
1598        glRotate(angley,sind(anglez),cosd(anglez),0)
1599        glMultMatrixf(A4mat.T)
1600        glTranslate(-Tx,-Ty,-Tz)
1601        if drawingData['unitCellBox']:
1602            RenderBox()
1603        if drawingData['showABC']:
1604#            try:            #temporary fix - not needed further?
1605#                x,y,z = drawingData['testPos'][0]
1606#            except TypeError:
1607#                x,y,z = drawingData['testPos']
1608            x,y,z = drawingData['testPos'][0]
1609            if altDown:
1610                self.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
1611            else:
1612                self.G2plotNB.status.SetStatusText('test point %.4f,%.4f,%.4f'%(x,y,z),1)           
1613            RenderUnitVectors(x,y,z)
1614        Backbone = []
1615        BackboneColor = []
1616        time0 = time.time()
1617        for iat,atom in enumerate(drawingData['Atoms']):
1618            x,y,z = atom[cx:cx+3]
1619            Bonds = atom[-1]
1620            try:
1621                atNum = generalData['AtomTypes'].index(atom[ct])
1622            except ValueError:
1623                atNum = -1
1624            CL = atom[cs+2]
1625            color = np.array(CL)/255.
1626            if iat in Ind:
1627                color = np.array(Gr)/255.
1628            radius = 0.5
1629            if atom[cs] != '':
1630                glLoadName(atom[-2])
1631            if 'balls' in atom[cs]:
1632                vdwScale = drawingData['vdwScale']
1633                ballScale = drawingData['ballScale']
1634                if atNum < 0:
1635                    radius = 0.2
1636                elif 'H' == atom[ct]:
1637                    if drawingData['showHydrogen']:
1638                        if 'vdW' in atom[cs] and atNum >= 0:
1639                            radius = vdwScale*generalData['vdWRadii'][atNum]
1640                        else:
1641                            radius = ballScale*drawingData['sizeH']
1642                    else:
1643                        radius = 0.0
1644                else:
1645                    if 'vdW' in atom[cs]:
1646                        radius = vdwScale*generalData['vdWRadii'][atNum]
1647                    else:
1648                        radius = ballScale*generalData['BondRadii'][atNum]
1649                RenderSphere(x,y,z,radius,color)
1650                if 'sticks' in atom[cs]:
1651                    RenderBonds(x,y,z,Bonds,bondR,color)
1652            elif 'ellipsoids' in atom[cs]:
1653                RenderBonds(x,y,z,Bonds,bondR,color)
1654                if atom[cs+3] == 'A':                   
1655                    Uij = atom[cs+5:cs+11]
1656                    U = np.multiply(G2spc.Uij2U(Uij),GS)
1657                    U = np.inner(Amat,np.inner(U,Amat).T)
1658                    E,R = nl.eigh(U)
1659                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1660                    E = np.sqrt(E)
1661                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
1662                        pass
1663                    else:
1664                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,color)                   
1665                else:
1666                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
1667                        pass
1668                    else:
1669                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
1670                        RenderSphere(x,y,z,radius,color)
1671            elif 'lines' in atom[cs]:
1672                radius = 0.1
1673                RenderLines(x,y,z,Bonds,color)
1674#                RenderBonds(x,y,z,Bonds,0.05,color,6)
1675            elif atom[cs] == 'sticks':
1676                radius = 0.1
1677                RenderBonds(x,y,z,Bonds,bondR,color)
1678            elif atom[cs] == 'polyhedra':
1679                if len(Bonds) > 2:
1680                    FaceGen = G2lat.uniqueCombinations(Bonds,3)     #N.B. this is a generator
1681                    Faces = []
1682                    for face in FaceGen:
1683                        vol = nl.det(face)
1684                        if abs(vol) > 1. or len(Bonds) == 3:
1685                            if vol < 0.:
1686                                face = [face[0],face[2],face[1]]
1687                            norm = np.cross(face[1]-face[0],face[2]-face[0])
1688                            norm /= np.sqrt(np.sum(norm**2))
1689                            Faces.append([face,norm])
1690                    RenderPolyhedra(x,y,z,Faces,color)
1691            elif atom[cs] == 'backbone':
1692                if atom[ct-1].split()[0] in ['C','N']:
1693                    Backbone.append(list(np.inner(Amat,np.array([x,y,z]))))
1694                    BackboneColor.append(list(color))
1695                   
1696            if atom[cs+1] == 'type':
1697                RenderLabel(x,y,z,atom[ct],radius)
1698            elif atom[cs+1] == 'name':
1699                RenderLabel(x,y,z,atom[ct-1],radius)
1700            elif atom[cs+1] == 'number':
1701                RenderLabel(x,y,z,str(iat+1),radius)
1702            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
1703                RenderLabel(x,y,z,atom[ct-4],radius)
1704            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
1705                RenderLabel(x,y,z,atom[ct-3],radius)
1706            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
1707                RenderLabel(x,y,z,atom[ct-2],radius)
1708        if Backbone:
1709            RenderBackbone(Backbone,BackboneColor,bondR)
1710#        print time.time()-time0
1711        Page.canvas.SwapBuffers()
1712       
1713    def OnSize(event):
1714        Draw()
1715       
1716    def OnFocus(event):
1717        Draw()
1718       
1719    try:
1720        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
1721        Page = self.G2plotNB.nb.GetPage(plotNum)       
1722    except ValueError,error:
1723        Plot = self.G2plotNB.addOgl(generalData['Name'])
1724        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
1725        Page = self.G2plotNB.nb.GetPage(plotNum)
1726        Page.views = False
1727        view = False
1728        altDown = False
1729    Page.SetFocus()
1730    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
1731        choices=(' save as:','jpeg','tiff','bmp'))
1732    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
1733    cb.SetValue(' save as:')
1734    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
1735    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
1736    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
1737    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
1738    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
1739    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
1740    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
1741    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
1742    Page.camera['position'] = drawingData['cameraPos']
1743    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
1744    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
1745    Page.canvas.SetCurrent()
1746    Draw()
1747       
Note: See TracBrowser for help on using the repository browser.