source: trunk/GSASIIplot.py @ 231

Last change on this file since 231 was 231, checked in by vondreele, 11 years ago

add top info line

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