source: trunk/GSASIIplot.py @ 239

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

add 'any image file' to image file menu
add calibration skip & dmin to image data dictionary
fix to ellipse fitting
fix Pilatus reading - OK for 100K, not sure for 2M
now a image sizexy - 2 items for x & y sizes

  • Property svn:keywords set to Date Author Revision URL Id
File size: 80.6 KB
Line 
1#GSASII plotting routines
2########### SVN repository information ###################
3# $Date: 2011-01-13 19:34:07 +0000 (Thu, 13 Jan 2011) $
4# $Author: vondreele $
5# $Revision: 239 $
6# $URL: trunk/GSASIIplot.py $
7# $Id: GSASIIplot.py 239 2011-01-13 19:34:07Z vondreele $
8########### SVN repository information ###################
9import math
10import time
11import copy
12import os.path
13import numpy as np
14import numpy.linalg as nl
15import wx
16import wx.aui
17import wx.glcanvas
18import matplotlib as mpl
19import mpl_toolkits.mplot3d.axes3d as mp3d
20import GSASIIpath
21import GSASIIgrid as G2gd
22import GSASIIimage as G2img
23import 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        sizexy = Data['size']
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 <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
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#            sizexy = Data['size']
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    sizexy = Data['size']
1084    pixelSize = Data['pixelSize']
1085#    print sizexy,pixelSize,len(self.ImageZ),len(self.ImageZ[0])
1086    scalex = 1000./pixelSize[0]
1087    scaley = 1000./pixelSize[1]
1088    Xmax = sizexy[0]*pixelSize[0]/1000.
1089    Ymax = sizexy[1]*pixelSize[1]/1000.
1090    xlim = (-0.5,Xmax-.5)
1091    ylim = (Ymax-.5,-0.5,)
1092    Imin,Imax = Data['range'][1]
1093    acolor = mpl.cm.get_cmap(Data['color'])
1094    xcent,ycent = Data['center']
1095    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
1096    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
1097    #do threshold mask - "real" mask - others are just bondaries
1098    Zlim = Masks['Thresholds'][1]
1099    wx.BeginBusyCursor()
1100    try:
1101        MA = ma.masked_greater(ma.masked_less(self.ImageZ,Zlim[0]),Zlim[1])
1102        MaskA = ma.getmaskarray(MA)
1103        A = G2img.ImageCompress(MA,imScale)
1104        AM = G2img.ImageCompress(MaskA,imScale)
1105       
1106        ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
1107            interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Xmax,0])
1108        Img = Plot.imshow(A,aspect='equal',cmap=acolor,
1109            interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Xmax,0])
1110        if self.setPoly:
1111            Img.set_picker(True)
1112   
1113        Plot.plot(xcent,ycent,'x')
1114        if Data['showLines']:
1115            LRAzim = Data['LRazimuth']                  #NB: integers
1116            IOtth = Data['IOtth']
1117            wave = Data['wavelength']
1118            dspI = wave/(2.0*sind(IOtth[0]/2.0))
1119            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
1120            dspO = wave/(2.0*sind(IOtth[1]/2.0))
1121            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
1122            if Data['fullIntegrate']:
1123                Azm = np.array(range(0,361))
1124            else:
1125                Azm = np.array(range(LRAzim[0],LRAzim[1]+1))
1126            if ellI:
1127                xyI = []
1128                for azm in Azm:
1129                    xyI.append(G2img.GetDetectorXY(dspI,azm,Data))
1130                xyI = np.array(xyI)
1131                arcxI,arcyI = xyI.T
1132                Plot.plot(arcxI,arcyI,picker=3)
1133            if ellO:
1134                xyO = []
1135                for azm in Azm:
1136                    xyO.append(G2img.GetDetectorXY(dspO,azm,Data))
1137                xyO = np.array(xyO)
1138                arcxO,arcyO = xyO.T
1139                Plot.plot(arcxO,arcyO,picker=3)
1140            if ellO and ellI and not Data['fullIntegrate']:
1141                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
1142                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
1143        for xring,yring in Data['ring']:
1144            Plot.plot(xring,yring,'r+',picker=3)
1145        if Data['setRings']:
1146            rings = np.concatenate((Data['rings']),axis=0)
1147            for xring,yring,dsp in rings:
1148                Plot.plot(xring,yring,'r+')           
1149        for ellipse in Data['ellipses']:
1150            cent,phi,[width,height],col = ellipse
1151            Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
1152            Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
1153        #masks - mask lines numbered after integration limit lines
1154        spots = Masks['Points']
1155        rings = Masks['Rings']
1156        arcs = Masks['Arcs']
1157        polygons = Masks['Polygons']
1158        for x,y,d in spots:
1159            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
1160        self.ringList = []
1161        for iring,(tth,thick) in enumerate(rings):
1162            wave = Data['wavelength']
1163            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
1164            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
1165            self.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
1166            self.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
1167        self.arcList = []
1168        for iarc,(tth,azm,thick) in enumerate(arcs):           
1169            wave = Data['wavelength']
1170            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
1171            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
1172            self.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
1173            self.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
1174            self.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
1175            self.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
1176        self.polyList = []
1177        for ipoly,polygon in enumerate(polygons):
1178            x,y = np.hsplit(np.array(polygon),2)
1179            self.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
1180            Plot.plot(x,y,'r')           
1181        colorBar = Page.figure.colorbar(Img)
1182        Plot.set_xlim(xlim)
1183        Plot.set_ylim(ylim)
1184        if not newPlot and xylim:
1185            Page.toolbar.push_current()
1186            Plot.set_xlim(xylim[0])
1187            Plot.set_ylim(xylim[1])
1188            xylim = []
1189            Page.toolbar.push_current()
1190            Page.toolbar.draw()
1191        else:
1192            Page.canvas.draw()
1193    finally:
1194        wx.EndBusyCursor()
1195       
1196def PlotIntegration(self,newPlot=False,event=None):
1197           
1198    def OnMotion(event):
1199        Page.canvas.SetToolTipString('')
1200        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1201        azm = event.ydata
1202        tth = event.xdata
1203        if azm and tth:
1204            self.G2plotNB.status.SetFields(\
1205                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1206                               
1207    try:
1208        plotNum = self.G2plotNB.plotList.index('2D Integration')
1209        Page = self.G2plotNB.nb.GetPage(plotNum)
1210        if not newPlot:
1211            Plot = Page.figure.gca()          #get previous plot & get limits
1212            xylim = Plot.get_xlim(),Plot.get_ylim()
1213        Page.figure.clf()
1214        Plot = Page.figure.gca()          #get a fresh plot after clf()
1215       
1216    except ValueError,error:
1217        Plot = self.G2plotNB.addMpl('2D Integration').gca()
1218        plotNum = self.G2plotNB.plotList.index('2D Integration')
1219        Page = self.G2plotNB.nb.GetPage(plotNum)
1220        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1221        Page.views = False
1222        view = False
1223    if not event:
1224        Page.SetFocus()
1225       
1226    Data = self.PatternTree.GetItemPyData(
1227        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1228    image = self.Integrate[0]
1229    xsc = self.Integrate[1]
1230    ysc = self.Integrate[2]
1231    Imin,Imax = Data['range'][1]
1232    acolor = mpl.cm.get_cmap(Data['color'])
1233    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1234    Plot.set_ylabel('azimuth',fontsize=12)
1235    Plot.set_xlabel('2-theta',fontsize=12)
1236    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
1237        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
1238    colorBar = Page.figure.colorbar(Img)
1239    if Data['setRings']:
1240        rings = np.concatenate((Data['rings']),axis=0)
1241        for xring,yring,dsp in rings:
1242            x,y = G2img.GetTthAzm(xring,yring,Data)
1243            Plot.plot(x,y,'r+')
1244    if Data['ellipses']:           
1245        for ellipse in Data['ellipses']:
1246            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1247            x,y = np.hsplit(ring,2)
1248            tth,azm = G2img.GetTthAzm(x,y,Data)
1249            Plot.plot(tth,azm,'b,')
1250    if not newPlot:
1251        Page.toolbar.push_current()
1252        Plot.set_xlim(xylim[0])
1253        Plot.set_ylim(xylim[1])
1254        xylim = []
1255        Page.toolbar.push_current()
1256        Page.toolbar.draw()
1257    else:
1258        Page.canvas.draw()
1259               
1260def PlotTRImage(self,tax,tay,taz,newPlot=False):
1261    #a test plot routine - not normally used
1262           
1263    def OnMotion(event):
1264        Page.canvas.SetToolTipString('')
1265        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1266        azm = event.xdata
1267        tth = event.ydata
1268        if azm and tth:
1269            self.G2plotNB.status.SetFields(\
1270                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1271                               
1272    try:
1273        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1274        Page = self.G2plotNB.nb.GetPage(plotNum)
1275        if not newPlot:
1276            Plot = Page.figure.gca()          #get previous plot & get limits
1277            xylim = Plot.get_xlim(),Plot.get_ylim()
1278        Page.figure.clf()
1279        Plot = Page.figure.gca()          #get a fresh plot after clf()
1280       
1281    except ValueError,error:
1282        Plot = self.G2plotNB.addMpl('2D Transformed Powder Image').gca()
1283        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1284        Page = self.G2plotNB.nb.GetPage(plotNum)
1285        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1286        Page.views = False
1287        view = False
1288    Page.SetFocus()
1289       
1290    Data = self.PatternTree.GetItemPyData(
1291        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1292    Imin,Imax = Data['range'][1]
1293    step = (Imax-Imin)/5.
1294    V = np.arange(Imin,Imax,step)
1295    acolor = mpl.cm.get_cmap(Data['color'])
1296    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1297    Plot.set_xlabel('azimuth',fontsize=12)
1298    Plot.set_ylabel('2-theta',fontsize=12)
1299    Plot.contour(tax,tay,taz,V,cmap=acolor)
1300    if Data['showLines']:
1301        IOtth = Data['IOtth']
1302        if Data['fullIntegrate']:
1303            LRAzim = [-180,180]
1304        else:
1305            LRAzim = Data['LRazimuth']                  #NB: integers
1306        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
1307        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
1308        if not Data['fullIntegrate']:
1309            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
1310            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
1311    if Data['setRings']:
1312        rings = np.concatenate((Data['rings']),axis=0)
1313        for xring,yring,dsp in rings:
1314            x,y = G2img.GetTthAzm(xring,yring,Data)
1315            Plot.plot(y,x,'r+')           
1316    if Data['ellipses']:           
1317        for ellipse in Data['ellipses']:
1318            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1319            x,y = np.hsplit(ring,2)
1320            tth,azm = G2img.GetTthAzm(x,y,Data)
1321            Plot.plot(azm,tth,'b,')
1322    if not newPlot:
1323        Page.toolbar.push_current()
1324        Plot.set_xlim(xylim[0])
1325        Plot.set_ylim(xylim[1])
1326        xylim = []
1327        Page.toolbar.push_current()
1328        Page.toolbar.draw()
1329    else:
1330        Page.canvas.draw()
1331       
1332def PlotStructure(self,data):
1333    generalData = data['General']
1334    cell = generalData['Cell'][1:7]
1335    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
1336    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1337    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1338    Mydir = generalData['Mydir']
1339    atomData = data['Atoms']
1340    drawingData = data['Drawing']
1341    drawAtoms = drawingData['Atoms']
1342    cx,ct,cs = drawingData['atomPtrs']
1343    Wt = [255,255,255]
1344    Rd = [255,0,0]
1345    Gr = [0,255,0]
1346    Bl = [0,0,255]
1347    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]])
1348    uEdges = np.array([
1349        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
1350        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
1351        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
1352    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
1353    altDown = False
1354    shiftDown = False
1355    ctrlDown = False
1356   
1357    def OnKeyBox(event):
1358        import Image
1359        Draw()                          #make sure plot is fresh!!
1360        mode = cb.GetValue()
1361        Fname = Mydir+'\\'+generalData['Name']+'.'+mode
1362        size = Page.canvas.GetSize()
1363        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
1364        if mode in ['jpeg',]:
1365            Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
1366            im = Image.new("RGBA", (size[0],size[1]))
1367        else:
1368            Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
1369            im = Image.new("RGB", (size[0],size[1]))
1370        im.fromstring(Pix)
1371        im.save(Fname,mode)
1372        cb.SetValue(' Save as:')
1373        self.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
1374   
1375    def GetTruePosition(xy):
1376        View = glGetIntegerv(GL_VIEWPORT)
1377        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
1378        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
1379        Zmax = 1.
1380        for i,atom in enumerate(drawAtoms):
1381            x,y,z = atom[cx:cx+3]
1382            X,Y,Z = gluProject(x,y,z,Model,Proj,View)
1383            XY = [int(X),int(View[3]-Y)]
1384            if np.allclose(xy,XY,atol=10) and Z < Zmax:
1385                Zmax = Z
1386                SetSelectedAtoms(i)
1387                   
1388    def OnMouseDown(event):
1389        xy = event.GetPosition()
1390        if event.ShiftDown():
1391            GetTruePosition(xy)
1392        else:
1393            drawingData['Rotation'][3] = xy
1394        Draw()
1395       
1396    def OnMouseMove(event):
1397        newxy = event.GetPosition()
1398        if event.ControlDown() and drawingData['showABC']:
1399            if event.LeftIsDown():
1400                ctrlDown = True
1401                SetTestRot(newxy)
1402            elif event.RightIsDown():
1403                SetTestPos(newxy)
1404            elif event.MiddleIsDown():
1405                SetTestRotZ(newxy)
1406               
1407               
1408        if event.Dragging() and not event.ControlDown():
1409            if event.LeftIsDown():
1410                SetRotation(newxy)
1411            elif event.RightIsDown():
1412                SetTranslation(newxy)
1413            elif event.MiddleIsDown():
1414                SetRotationZ(newxy)
1415        Draw()
1416       
1417    def OnMouseWheel(event):
1418        drawingData['cameraPos'] += event.GetWheelRotation()/24
1419        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
1420        page = self.dataDisplay.GetSelection()
1421        if page:
1422            if self.dataDisplay.GetPageText(page) == 'Draw Options':
1423                panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1424                names = [child.GetName() for child in panel]
1425                panel[names.index('cameraPos')].SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
1426                panel[names.index('cameraSlider')].SetValue(drawingData['cameraPos'])
1427        Draw()
1428           
1429    def SetViewPointText(VP):
1430        page = self.dataDisplay.GetSelection()
1431        if page:
1432            if self.dataDisplay.GetPageText(page) == 'Draw Options':
1433                panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1434                names = [child.GetName() for child in panel]
1435                panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
1436           
1437    def ClearSelectedAtoms():
1438        page = self.dataDisplay.GetSelection()
1439        if page:
1440            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1441                self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
1442            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1443                self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
1444                   
1445    def SetSelectedAtoms(ind):
1446        page = self.dataDisplay.GetSelection()
1447        if page:
1448            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1449                self.dataDisplay.GetPage(page).SelectRow(ind)      #this is the Atoms grid in Draw Atoms
1450            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1451                Id = drawAtoms[ind][-2]
1452                for i,atom in enumerate(atomData):
1453                    if atom[-1] == Id:
1454                        self.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
1455                 
1456    def GetSelectedAtoms():
1457        page = self.dataDisplay.GetSelection()
1458        Ind = []
1459        if page:
1460            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1461                Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
1462            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1463                Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
1464        return Ind
1465                                       
1466    def OnKey(event):           #on key UP!!
1467        keyCode = event.GetKeyCode()
1468        if keyCode > 255:
1469            keyCode = 0
1470        key,xyz = chr(keyCode),event.GetPosition()
1471        indx = drawingData['selectedAtoms']
1472        if key in ['c','C']:
1473            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
1474            drawingData['testPos'] = [[-.1,-.1,-.1],[0.0,0.0,0.0],[0,0]]
1475            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
1476            SetViewPointText(drawingData['viewPoint'][0])
1477        elif key in ['n','N']:
1478            drawAtoms = drawingData['Atoms']
1479            pI = drawingData['viewPoint'][1]
1480            if indx:
1481                pI[0] = indx[pI[1]]
1482                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1483                pI[1] += 1
1484                if pI[1] >= len(indx):
1485                    pI[1] = 0
1486            else:
1487                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1488                pI[0] += 1
1489                if pI[0] >= len(drawAtoms):
1490                    pI[0] = 0
1491            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1492            SetViewPointText(drawingData['viewPoint'][0])
1493               
1494        elif key in ['p','P']:
1495            drawAtoms = drawingData['Atoms']
1496            pI = drawingData['viewPoint'][1]
1497            if indx:
1498                pI[0] = indx[pI[1]]
1499                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1500                pI[1] -= 1
1501                if pI[1] < 0:
1502                    pI[1] = len(indx)-1
1503            else:
1504                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1505                pI[0] -= 1
1506                if pI[0] < 0:
1507                    pI[0] = len(drawAtoms)-1
1508            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1509            SetViewPointText(drawingData['viewPoint'][0])           
1510        Draw()
1511           
1512    def SetBackground():
1513        R,G,B,A = Page.camera['backColor']
1514        glClearColor(R,G,B,A)
1515        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
1516       
1517    def SetLights():
1518        glEnable(GL_DEPTH_TEST)
1519        glShadeModel(GL_SMOOTH)
1520        glEnable(GL_LIGHTING)
1521        glEnable(GL_LIGHT0)
1522        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
1523        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
1524        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
1525       
1526    def SetTranslation(newxy):
1527        Tx,Ty,Tz = drawingData['viewPoint'][0]
1528        anglex,angley,anglez,oldxy = drawingData['Rotation']
1529        Rx = G2lat.rotdMat(anglex,0)
1530        Ry = G2lat.rotdMat(angley,1)
1531        Rz = G2lat.rotdMat(anglez,2)
1532        dxy = list(newxy-oldxy)+[0,]
1533        dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
1534        Tx -= dxy[0]*0.01
1535        Ty += dxy[1]*0.01
1536        Tz -= dxy[2]*0.01
1537        drawingData['Rotation'][3] = newxy
1538        drawingData['viewPoint'][0] =  Tx,Ty,Tz
1539        SetViewPointText([Tx,Ty,Tz])
1540       
1541    def SetTestPos(newxy):
1542        Tx,Ty,Tz = drawingData['testPos'][0]
1543        anglex,angley,anglez,oldxy = drawingData['Rotation']
1544        Rx = G2lat.rotdMat(anglex,0)
1545        Ry = G2lat.rotdMat(angley,1)
1546        Rz = G2lat.rotdMat(anglez,2)
1547        dxy = list(newxy-oldxy)+[0,]
1548        dxy = np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy)))
1549        Tx += dxy[0]*0.001
1550        Ty -= dxy[1]*0.001
1551        Tz += dxy[2]*0.001
1552        drawingData['Rotation'][3] = newxy
1553        drawingData['testPos'][0] =  Tx,Ty,Tz
1554       
1555    def SetTestRot(newxy):
1556        Txyz = np.array(drawingData['testPos'][0])
1557        oldxy = drawingData['testPos'][2]
1558        Ax,Ay,Az = drawingData['testPos'][1]
1559        Vxyz = np.array(drawingData['viewPoint'][0])
1560        Dxyz = np.inner(Amat,Txyz-Vxyz)
1561        dxy = list(newxy-oldxy)+[0,]
1562        Ax += dxy[1]*0.01
1563        Ay += dxy[0]*0.01
1564        Rx = G2lat.rotdMat(Ax,0)
1565        Ry = G2lat.rotdMat(Ay,1)
1566        Dxyz = np.inner(Ry,np.inner(Rx,Dxyz))       
1567        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
1568        drawingData['testPos'][1] = [Ax,Ay,Az]
1569        drawingData['testPos'][2] = newxy
1570        drawingData['testPos'][0] = Dxyz
1571       
1572    def SetTestRotZ(newxy):
1573        Txyz = np.array(drawingData['testPos'][0])
1574        oldxy = drawingData['testPos'][2]
1575        Ax,Ay,Az = drawingData['testPos'][1]
1576        Vxyz = np.array(drawingData['viewPoint'][0])
1577        Dxyz = np.inner(Amat,Txyz-Vxyz)       
1578        dxy = list(newxy-oldxy)+[0,]
1579        Az += (dxy[0]+dxy[1])*.01
1580        Rz = G2lat.rotdMat(Az,2)
1581        Dxyz = np.inner(Rz,Dxyz)       
1582        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
1583        drawingData['testPos'][1] = [Ax,Ay,Az]
1584        drawingData['testPos'][2] = newxy
1585        drawingData['testPos'][0] = Dxyz
1586                             
1587    def SetRotation(newxy):       
1588        anglex,angley,anglez,oldxy = drawingData['Rotation']
1589        dxy = newxy-oldxy
1590        anglex += dxy[1]*.25
1591        angley += dxy[0]*.25
1592        oldxy = newxy
1593        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
1594       
1595    def SetRotationZ(newxy):                       
1596        anglex,angley,anglez,oldxy = drawingData['Rotation']
1597        dxy = newxy-oldxy
1598        anglez += (dxy[0]+dxy[1])*.25
1599        oldxy = newxy
1600        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
1601       
1602    def RenderBox():
1603        glEnable(GL_COLOR_MATERIAL)
1604        glLineWidth(2)
1605        glEnable(GL_BLEND)
1606        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
1607        glEnable(GL_LINE_SMOOTH)
1608        glBegin(GL_LINES)
1609        for line,color in zip(uEdges,uColors):
1610            glColor3ubv(color)
1611            glVertex3fv(line[0])
1612            glVertex3fv(line[1])
1613        glEnd()
1614        glColor4ubv([0,0,0,0])
1615        glDisable(GL_LINE_SMOOTH)
1616        glDisable(GL_BLEND)
1617        glDisable(GL_COLOR_MATERIAL)
1618       
1619    def RenderUnitVectors(x,y,z):
1620        xyz = np.array([x,y,z])
1621        glEnable(GL_COLOR_MATERIAL)
1622        glLineWidth(1)
1623        glPushMatrix()
1624        glTranslate(x,y,z)
1625        glScalef(1/cell[0],1/cell[1],1/cell[2])
1626        glBegin(GL_LINES)
1627        for line,color in zip(uEdges,uColors)[:3]:
1628            glColor3ubv(color)
1629            glVertex3fv(line[0])
1630            glVertex3fv(line[1])
1631        glEnd()
1632        glPopMatrix()
1633        glColor4ubv([0,0,0,0])
1634        glDisable(GL_COLOR_MATERIAL)
1635               
1636    def RenderSphere(x,y,z,radius,color):
1637        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1638        glPushMatrix()
1639        glTranslate(x,y,z)
1640        glMultMatrixf(B4mat.T)
1641        q = gluNewQuadric()
1642        gluSphere(q,radius,20,10)
1643        glPopMatrix()
1644       
1645    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
1646        s1,s2,s3 = E
1647        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1648        glPushMatrix()
1649        glTranslate(x,y,z)
1650        glMultMatrixf(B4mat.T)
1651        glMultMatrixf(R4.T)
1652        glEnable(GL_NORMALIZE)
1653        glScale(s1,s2,s3)
1654        q = gluNewQuadric()
1655        gluSphere(q,ellipseProb,20,10)
1656        glDisable(GL_NORMALIZE)
1657        glPopMatrix()
1658       
1659    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
1660        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1661        glPushMatrix()
1662        glTranslate(x,y,z)
1663        glMultMatrixf(B4mat.T)
1664        for bond in Bonds:
1665            glPushMatrix()
1666            Dx = np.inner(Amat,bond)
1667            Z = np.sqrt(np.sum(Dx**2))
1668            azm = atan2d(-Dx[1],-Dx[0])
1669            phi = acosd(Dx[2]/Z)
1670            glRotate(-azm,0,0,1)
1671            glRotate(phi,1,0,0)
1672            q = gluNewQuadric()
1673            gluCylinder(q,radius,radius,Z,slice,2)
1674            glPopMatrix()           
1675        glPopMatrix()
1676               
1677    def RenderLines(x,y,z,Bonds,color):
1678        xyz = np.array([x,y,z])
1679        glEnable(GL_COLOR_MATERIAL)
1680        glLineWidth(1)
1681        glColor3fv(color)
1682        glPushMatrix()
1683        glBegin(GL_LINES)
1684        for bond in Bonds:
1685            glVertex3fv(xyz)
1686            glVertex3fv(xyz+bond)
1687        glEnd()
1688        glColor4ubv([0,0,0,0])
1689        glPopMatrix()
1690        glDisable(GL_COLOR_MATERIAL)
1691       
1692    def RenderPolyhedra(x,y,z,Faces,color):
1693        glPushMatrix()
1694        glTranslate(x,y,z)
1695        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1696        glShadeModel(GL_SMOOTH)
1697        glMultMatrixf(B4mat.T)
1698        for face,norm in Faces:
1699            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
1700            glFrontFace(GL_CW)
1701            glNormal3fv(norm)
1702            glBegin(GL_TRIANGLES)
1703            for vert in face:
1704                glVertex3fv(vert)
1705            glEnd()
1706        glPopMatrix()
1707       
1708    def RenderBackbone(Backbone,BackboneColor,radius):
1709        glPushMatrix()
1710        glMultMatrixf(B4mat.T)
1711        glEnable(GL_COLOR_MATERIAL)
1712        glShadeModel(GL_SMOOTH)
1713        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
1714        glePolyCylinder(Backbone,BackboneColor,radius)
1715        glPopMatrix()       
1716        glDisable(GL_COLOR_MATERIAL)
1717       
1718    def RenderLabel(x,y,z,label,r):       
1719        glPushMatrix()
1720        glTranslate(x,y,z)
1721        glMultMatrixf(B4mat.T)
1722        glDisable(GL_LIGHTING)
1723        glColor3f(0,1.,0)
1724        glRasterPos3f(r,r,r)
1725        for c in list(label):
1726            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
1727        glEnable(GL_LIGHTING)
1728        glPopMatrix()
1729                           
1730    def Draw():
1731        import numpy.linalg as nl
1732        Ind = GetSelectedAtoms()
1733        VS = np.array(Page.canvas.GetSize())
1734        aspect = float(VS[0])/float(VS[1])
1735        cPos = drawingData['cameraPos']
1736        Zclip = drawingData['Zclip']*cPos/200.
1737        anglex,angley,anglez = drawingData['Rotation'][:3]
1738        Tx,Ty,Tz = drawingData['viewPoint'][0]
1739        cx,ct,cs = drawingData['atomPtrs']
1740        bondR = drawingData['bondRadius']
1741        G,g = G2lat.cell2Gmat(cell)
1742        GS = G
1743        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
1744        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
1745        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
1746        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
1747       
1748        SetBackground()
1749        glInitNames()
1750        glPushName(0)
1751       
1752        glMatrixMode(GL_PROJECTION)
1753        glLoadIdentity()
1754        glViewport(0,0,VS[0],VS[1])
1755        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
1756        gluLookAt(0,0,cPos,0,0,0,0,1,0)
1757        SetLights()           
1758           
1759        glMatrixMode(GL_MODELVIEW)
1760        glLoadIdentity()
1761        glRotate(anglez,0,0,1)
1762        glRotate(anglex,cosd(anglez),-sind(anglez),0)
1763        glRotate(angley,sind(anglez),cosd(anglez),0)
1764        glMultMatrixf(A4mat.T)
1765        glTranslate(-Tx,-Ty,-Tz)
1766        if drawingData['unitCellBox']:
1767            RenderBox()
1768        if drawingData['showABC']:
1769#            try:            #temporary fix - not needed further?
1770#                x,y,z = drawingData['testPos'][0]
1771#            except TypeError:
1772#                x,y,z = drawingData['testPos']
1773            x,y,z = drawingData['testPos'][0]
1774            if altDown:
1775                self.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
1776            else:
1777                self.G2plotNB.status.SetStatusText('test point %.4f,%.4f,%.4f'%(x,y,z),1)           
1778            RenderUnitVectors(x,y,z)
1779        Backbone = []
1780        BackboneColor = []
1781        time0 = time.time()
1782        for iat,atom in enumerate(drawingData['Atoms']):
1783            x,y,z = atom[cx:cx+3]
1784            Bonds = atom[-1]
1785            try:
1786                atNum = generalData['AtomTypes'].index(atom[ct])
1787            except ValueError:
1788                atNum = -1
1789            CL = atom[cs+2]
1790            color = np.array(CL)/255.
1791            if iat in Ind:
1792                color = np.array(Gr)/255.
1793            radius = 0.5
1794            if atom[cs] != '':
1795                glLoadName(atom[-2])
1796            if 'balls' in atom[cs]:
1797                vdwScale = drawingData['vdwScale']
1798                ballScale = drawingData['ballScale']
1799                if atNum < 0:
1800                    radius = 0.2
1801                elif 'H' == atom[ct]:
1802                    if drawingData['showHydrogen']:
1803                        if 'vdW' in atom[cs] and atNum >= 0:
1804                            radius = vdwScale*generalData['vdWRadii'][atNum]
1805                        else:
1806                            radius = ballScale*drawingData['sizeH']
1807                    else:
1808                        radius = 0.0
1809                else:
1810                    if 'vdW' in atom[cs]:
1811                        radius = vdwScale*generalData['vdWRadii'][atNum]
1812                    else:
1813                        radius = ballScale*generalData['BondRadii'][atNum]
1814                RenderSphere(x,y,z,radius,color)
1815                if 'sticks' in atom[cs]:
1816                    RenderBonds(x,y,z,Bonds,bondR,color)
1817            elif 'ellipsoids' in atom[cs]:
1818                RenderBonds(x,y,z,Bonds,bondR,color)
1819                if atom[cs+3] == 'A':                   
1820                    Uij = atom[cs+5:cs+11]
1821                    U = np.multiply(G2spc.Uij2U(Uij),GS)
1822                    U = np.inner(Amat,np.inner(U,Amat).T)
1823                    E,R = nl.eigh(U)
1824                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1825                    E = np.sqrt(E)
1826                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
1827                        pass
1828                    else:
1829                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,color)                   
1830                else:
1831                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
1832                        pass
1833                    else:
1834                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
1835                        RenderSphere(x,y,z,radius,color)
1836            elif 'lines' in atom[cs]:
1837                radius = 0.1
1838                RenderLines(x,y,z,Bonds,color)
1839#                RenderBonds(x,y,z,Bonds,0.05,color,6)
1840            elif atom[cs] == 'sticks':
1841                radius = 0.1
1842                RenderBonds(x,y,z,Bonds,bondR,color)
1843            elif atom[cs] == 'polyhedra':
1844                if len(Bonds) > 2:
1845                    FaceGen = G2lat.uniqueCombinations(Bonds,3)     #N.B. this is a generator
1846                    Faces = []
1847                    for face in FaceGen:
1848                        vol = nl.det(face)
1849                        if abs(vol) > 1. or len(Bonds) == 3:
1850                            if vol < 0.:
1851                                face = [face[0],face[2],face[1]]
1852                            norm = np.cross(face[1]-face[0],face[2]-face[0])
1853                            norm /= np.sqrt(np.sum(norm**2))
1854                            Faces.append([face,norm])
1855                    RenderPolyhedra(x,y,z,Faces,color)
1856            elif atom[cs] == 'backbone':
1857                if atom[ct-1].split()[0] in ['C','N']:
1858                    Backbone.append(list(np.inner(Amat,np.array([x,y,z]))))
1859                    BackboneColor.append(list(color))
1860                   
1861            if atom[cs+1] == 'type':
1862                RenderLabel(x,y,z,atom[ct],radius)
1863            elif atom[cs+1] == 'name':
1864                RenderLabel(x,y,z,atom[ct-1],radius)
1865            elif atom[cs+1] == 'number':
1866                RenderLabel(x,y,z,str(iat+1),radius)
1867            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
1868                RenderLabel(x,y,z,atom[ct-4],radius)
1869            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
1870                RenderLabel(x,y,z,atom[ct-3],radius)
1871            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
1872                RenderLabel(x,y,z,atom[ct-2],radius)
1873        if Backbone:
1874            RenderBackbone(Backbone,BackboneColor,bondR)
1875#        print time.time()-time0
1876        Page.canvas.SwapBuffers()
1877       
1878    def OnSize(event):
1879        Draw()
1880       
1881    def OnFocus(event):
1882        Draw()
1883       
1884    try:
1885        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
1886        Page = self.G2plotNB.nb.GetPage(plotNum)       
1887    except (ValueError,error):
1888        Plot = self.G2plotNB.addOgl(generalData['Name'])
1889        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
1890        Page = self.G2plotNB.nb.GetPage(plotNum)
1891        Page.views = False
1892        view = False
1893        altDown = False
1894    Page.SetFocus()
1895    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
1896        choices=(' save as:','jpeg','tiff','bmp'))
1897    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
1898    cb.SetValue(' save as:')
1899    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
1900    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
1901    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
1902    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
1903    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
1904    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
1905    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
1906    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
1907    Page.camera['position'] = drawingData['cameraPos']
1908    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
1909    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
1910    Page.canvas.SetCurrent()
1911    Draw()
1912       
Note: See TracBrowser for help on using the repository browser.