source: trunk/GSASIIplot.py @ 243

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

GSASIIpeak.py - added factorize for future use in finding good numbers for fft
GSASIIimgGUI.py - minor tweak to q-value formatting & make small angle name 'SASD'
GSASIIIO.py - major refactor of GetTifData?

  • Property svn:keywords set to Date Author Revision URL Id
File size: 79.2 KB
Line 
1#GSASII plotting routines
2########### SVN repository information ###################
3# $Date: 2011-02-08 21:47:07 +0000 (Tue, 08 Feb 2011) $
4# $Author: vondreele $
5# $Revision: 243 $
6# $URL: trunk/GSASIIplot.py $
7# $Id: GSASIIplot.py 243 2011-02-08 21:47: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    global Data,Masks
835    Data = self.PatternTree.GetItemPyData(
836        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
837    Masks = self.PatternTree.GetItemPyData(
838        G2gd.GetPatternTreeItemId(self,self.Image, 'Masks'))
839
840    def OnImMotion(event):
841        Page.canvas.SetToolTipString('')
842        sizexy = Data['size']
843        if event.xdata and event.ydata:                 #avoid out of frame errors
844            Page.canvas.SetCursor(wx.CROSS_CURSOR)
845            item = self.itemPicked
846            pixelSize = Data['pixelSize']
847            scalex = 1000./pixelSize[0]
848            scaley = 1000./pixelSize[1]
849            if item and self.PatternTree.GetItemText(self.PickId) == 'Image Controls':
850                if 'Text' in str(item):
851                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
852                else:
853                    xcent,ycent = Data['center']
854                    xpos = event.xdata-xcent
855                    ypos = event.ydata-ycent
856                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
857                        ang = int(atan2d(xpos,ypos))
858                        Page.canvas.SetToolTipString('%6d deg'%(ang))
859                    elif 'line1' in  str(item) or 'line2' in str(item):
860                        tth = G2img.GetTth(event.xdata,event.ydata,Data)
861                        Page.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
862            else:
863                xpos = event.xdata
864                ypos = event.ydata
865                xpix = xpos*scalex
866                ypix = ypos*scaley
867                Int = 0
868                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
869                    Int = self.ImageZ[ypix][xpix]
870                tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
871                Q = 2.*math.pi/dsp
872                if self.setPoly:
873                    self.G2plotNB.status.SetFields(['','Polygon mask pick - LB next point, RB close polygon'])
874                else:
875                    self.G2plotNB.status.SetFields(\
876                        ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)])
877
878    def OnImPlotKeyPress(event):
879        if self.PatternTree.GetItemText(self.PickId) == 'Masks':
880            Xpos = event.xdata
881            if not Xpos:            #got point out of frame
882                return
883            Ypos = event.ydata
884            if event.key == 's':
885                print 'spot mask @ ',Xpos,Ypos
886                Masks['Points'].append([Xpos,Ypos,1.])
887            elif event.key == 'r':
888                tth = G2img.GetTth(Xpos,Ypos,Data)
889                print 'ring mask @ ',Xpos,Ypos,tth
890                Masks['Rings'].append([tth,0.1])
891            elif event.key == 'a':
892                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
893                azm = int(azm)               
894                print 'arc mask @ ', Xpos,Ypos
895                Masks['Arcs'].append([tth,[azm-5,azm+5],0.1])
896            elif event.key == 'p':
897                self.setPoly = True
898                Masks['Polygons'].append([])
899                self.G2plotNB.status.SetFields(['','Polygon mask active - LB pick next point, RB close polygon'])
900            G2imG.UpdateMasks(self,Masks)
901        else:
902            Xpos = event.xdata
903            if not Xpos:            #got point out of frame
904                return
905            Ypos = event.ydata
906            if event.key == 'c':
907                print 'move center to: ',Xpos,Ypos
908                Data['center'] = [Xpos,Ypos]
909                G2imG.UpdateImageControls(self,Data,Masks)
910        PlotImage(self)
911           
912    def OnImPick(event):
913        if self.PatternTree.GetItemText(self.PickId) not in ['Image Controls','Masks']:
914            return
915        if self.setPoly:
916            polygon = Masks['Polygons'][-1]
917            xpos,ypos = event.mouseevent.xdata,event.mouseevent.ydata
918            if xpos and ypos:                       #point inside image
919                if len(polygon) > 2 and event.mouseevent.button == 3:
920                    x0,y0 = polygon[0]
921                    polygon.append([x0,y0])
922                    self.setPoly = False
923                    self.G2plotNB.status.SetFields(['','Polygon closed - RB drag a vertex to change shape'])
924                else:
925                    self.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(xpos,ypos)])
926                    polygon.append([xpos,ypos])
927                G2imG.UpdateMasks(self,Masks)
928        else:
929            if self.itemPicked is not None: return
930            self.itemPicked = event.artist
931            self.mousePicked = event.mouseevent
932       
933    def OnImRelease(event):
934        PickName = self.PatternTree.GetItemText(self.PickId)
935        if PickName not in ['Image Controls','Masks']:
936            return
937        pixelSize = Data['pixelSize']
938        scalex = 1000./pixelSize[0]
939        scaley = 1000./pixelSize[1]
940        pixLimit = Data['pixLimit']
941        if self.itemPicked is None and PickName == 'Image Controls':
942#            sizexy = Data['size']
943            Xpos = event.xdata
944            if not (Xpos and self.ifGetRing):                   #got point out of frame
945                return
946            Ypos = event.ydata
947            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
948                if event.button == 1:
949                    Xpix = Xpos*scalex
950                    Ypix = Ypos*scaley
951                    xpos,ypos,I,J = G2img.ImageLocalMax(self.ImageZ,pixLimit,Xpix,Ypix)
952                    if I and J:
953                        xpos += .5                              #shift to pixel center
954                        ypos += .5
955                        xpos /= scalex                          #convert to mm
956                        ypos /= scaley
957                        Data['ring'].append([xpos,ypos])
958                PlotImage(self)
959            return
960        else:
961            xpos = event.xdata
962            if xpos:                                        #avoid out of frame mouse position
963                ypos = event.ydata
964                if self.ifGetRing:
965                    xypos = [xpos,ypos]
966                    rings = Data['ring']
967                    for ring in rings:
968                        if np.allclose(ring,xypos,.01,0):
969                            rings.remove(ring)                                                                       
970                else:
971                    tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
972                    itemPicked = str(self.itemPicked)
973                    if 'Line2D' in itemPicked and PickName == 'Image Controls':
974                        if 'line1' in itemPicked:
975                            Data['IOtth'][0] = tth
976                        elif 'line2' in itemPicked:
977                            Data['IOtth'][1] = tth
978                        elif 'line3' in itemPicked and not Data['fullIntegrate']:
979                            Data['LRazimuth'][0] = int(azm)
980                        elif 'line4' in itemPicked and not Data['fullIntegrate']:
981                            Data['LRazimuth'][1] = int(azm)
982                           
983                        if Data['LRazimuth'][1] < Data['LRazimuth'][0]:
984                            Data['LRazimuth'][1] += 360
985                        if  Data['IOtth'][0] > Data['IOtth'][1]:
986                            Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
987                           
988                        self.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
989                        self.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
990                        self.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
991                        self.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
992                    elif 'Circle' in itemPicked and PickName == 'Masks':
993                        spots = Masks['Points']
994                        newPos = itemPicked.split(')')[0].split('(')[2].split(',')
995                        newPos = np.array([float(newPos[0]),float(newPos[1])])
996                        for spot in spots:
997                            if np.allclose(np.array([spot[:2]]),newPos):
998                                spot[:2] = xpos,ypos
999                        G2imG.UpdateMasks(self,Masks)
1000                    elif 'Line2D' in itemPicked and PickName == 'Masks':
1001                        Obj = self.itemPicked.findobj()
1002                        rings = Masks['Rings']
1003                        arcs = Masks['Arcs']
1004                        polygons = Masks['Polygons']
1005                        for ring in self.ringList:
1006                            if Obj == ring[0]:
1007                                rN = ring[1]
1008                                if ring[2] == 'o':
1009                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)-rings[rN][1]/2.
1010                                else:
1011                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)+rings[rN][1]/2.
1012                        for arc in self.arcList:
1013                            if Obj == arc[0]:
1014                                aN = arc[1]
1015                                if arc[2] == 'o':
1016                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)-arcs[aN][2]/2
1017                                elif arc[2] == 'i':
1018                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)+arcs[aN][2]/2
1019                                elif arc[2] == 'l':
1020                                    arcs[aN][1][0] = int(G2img.GetAzm(xpos,ypos,Data))
1021                                else:
1022                                    arcs[aN][1][1] = int(G2img.GetAzm(xpos,ypos,Data))
1023                        for poly in self.polyList:
1024                            if Obj == poly[0]:
1025                                ind = self.itemPicked.contains(self.mousePicked)[1]['ind'][0]
1026                                oldPos = np.array([self.mousePicked.xdata,self.mousePicked.ydata])
1027                                pN = poly[1]
1028                                for i,xy in enumerate(polygons[pN]):
1029                                    if np.allclose(np.array([xy]),oldPos,atol=1.0):
1030                                        polygons[pN][i] = xpos,ypos
1031                        G2imG.UpdateMasks(self,Masks)
1032#                    else:                  #keep for future debugging
1033#                        print str(self.itemPicked),event.xdata,event.ydata,event.button
1034                PlotImage(self)
1035            self.itemPicked = None
1036           
1037    try:
1038        plotNum = self.G2plotNB.plotList.index('2D Powder Image')
1039        Page = self.G2plotNB.nb.GetPage(plotNum)
1040        if not newPlot:
1041            Plot = Page.figure.gca()          #get previous powder plot & get limits
1042            xylim = Plot.get_xlim(),Plot.get_ylim()
1043        Page.figure.clf()
1044        Plot = Page.figure.gca()          #get a fresh plot after clf()
1045       
1046    except ValueError,error:
1047        Plot = self.G2plotNB.addMpl('2D Powder Image').gca()
1048        plotNum = self.G2plotNB.plotList.index('2D Powder Image')
1049        Page = self.G2plotNB.nb.GetPage(plotNum)
1050        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
1051        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
1052        Page.canvas.mpl_connect('pick_event', OnImPick)
1053        Page.canvas.mpl_connect('button_release_event', OnImRelease)
1054        xylim = []
1055    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
1056        Page.SetFocus()
1057    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1058    size,imagefile = self.PatternTree.GetItemPyData(self.Image)
1059    if imagefile != self.oldImagefile:
1060        imagefile = G2IO.CheckImageFile(self,imagefile)
1061        if not imagefile:
1062            self.G2plotNB.Delete('2D Powder Image')
1063            return
1064        self.PatternTree.SetItemPyData(self.Image,[size,imagefile])
1065        self.ImageZ = G2IO.GetImageData(self,imagefile,imageOnly=True)
1066        self.oldImagefile = imagefile
1067
1068    imScale = 1
1069    if len(self.ImageZ) > 1024:
1070        imScale = len(self.ImageZ)/1024
1071    sizexy = Data['size']
1072    pixelSize = Data['pixelSize']
1073    scalex = 1000./pixelSize[0]
1074    scaley = 1000./pixelSize[1]
1075    Xmax = sizexy[0]*pixelSize[0]/1000.
1076    Ymax = sizexy[1]*pixelSize[1]/1000.
1077    xlim = (0,Xmax)
1078    ylim = (Ymax,0)
1079    Imin,Imax = Data['range'][1]
1080    acolor = mpl.cm.get_cmap(Data['color'])
1081    xcent,ycent = Data['center']
1082    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
1083    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
1084    #do threshold mask - "real" mask - others are just bondaries
1085    Zlim = Masks['Thresholds'][1]
1086    wx.BeginBusyCursor()
1087    try:
1088        MA = ma.masked_greater(ma.masked_less(self.ImageZ,Zlim[0]),Zlim[1])
1089        MaskA = ma.getmaskarray(MA)
1090        A = G2img.ImageCompress(MA,imScale)
1091        AM = G2img.ImageCompress(MaskA,imScale)
1092       
1093        ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
1094            interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Xmax,0])
1095        Img = Plot.imshow(A,aspect='equal',cmap=acolor,
1096            interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
1097        if self.setPoly:
1098            Img.set_picker(True)
1099   
1100        Plot.plot(xcent,ycent,'x')
1101        if Data['showLines']:
1102            LRAzim = Data['LRazimuth']                  #NB: integers
1103            IOtth = Data['IOtth']
1104            wave = Data['wavelength']
1105            dspI = wave/(2.0*sind(IOtth[0]/2.0))
1106            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
1107            dspO = wave/(2.0*sind(IOtth[1]/2.0))
1108            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
1109            if Data['fullIntegrate']:
1110                Azm = np.array(range(0,361))
1111            else:
1112                Azm = np.array(range(LRAzim[0],LRAzim[1]+1))
1113            if ellI:
1114                xyI = []
1115                for azm in Azm:
1116                    xyI.append(G2img.GetDetectorXY(dspI,azm,Data))
1117                xyI = np.array(xyI)
1118                arcxI,arcyI = xyI.T
1119                Plot.plot(arcxI,arcyI,picker=3)
1120            if ellO:
1121                xyO = []
1122                for azm in Azm:
1123                    xyO.append(G2img.GetDetectorXY(dspO,azm,Data))
1124                xyO = np.array(xyO)
1125                arcxO,arcyO = xyO.T
1126                Plot.plot(arcxO,arcyO,picker=3)
1127            if ellO and ellI and not Data['fullIntegrate']:
1128                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
1129                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
1130        for xring,yring in Data['ring']:
1131            Plot.plot(xring,yring,'r+',picker=3)
1132        if Data['setRings']:
1133            rings = np.concatenate((Data['rings']),axis=0)
1134            for xring,yring,dsp in rings:
1135                Plot.plot(xring,yring,'r+')           
1136        for ellipse in Data['ellipses']:
1137            cent,phi,[width,height],col = ellipse
1138            Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
1139            Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
1140        #masks - mask lines numbered after integration limit lines
1141        spots = Masks['Points']
1142        rings = Masks['Rings']
1143        arcs = Masks['Arcs']
1144        polygons = Masks['Polygons']
1145        for x,y,d in spots:
1146            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
1147        self.ringList = []
1148        for iring,(tth,thick) in enumerate(rings):
1149            wave = Data['wavelength']
1150            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
1151            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
1152            self.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
1153            self.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
1154        self.arcList = []
1155        for iarc,(tth,azm,thick) in enumerate(arcs):           
1156            wave = Data['wavelength']
1157            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
1158            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
1159            self.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
1160            self.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
1161            self.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
1162            self.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
1163        self.polyList = []
1164        for ipoly,polygon in enumerate(polygons):
1165            x,y = np.hsplit(np.array(polygon),2)
1166            self.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
1167            Plot.plot(x,y,'r')           
1168        colorBar = Page.figure.colorbar(Img)
1169        Plot.set_xlim(xlim)
1170        Plot.set_ylim(ylim)
1171        if not newPlot and xylim:
1172            Page.toolbar.push_current()
1173            Plot.set_xlim(xylim[0])
1174            Plot.set_ylim(xylim[1])
1175            xylim = []
1176            Page.toolbar.push_current()
1177            Page.toolbar.draw()
1178        else:
1179            Page.canvas.draw()
1180    finally:
1181        wx.EndBusyCursor()
1182       
1183def PlotIntegration(self,newPlot=False,event=None):
1184           
1185    def OnMotion(event):
1186        Page.canvas.SetToolTipString('')
1187        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1188        azm = event.ydata
1189        tth = event.xdata
1190        if azm and tth:
1191            self.G2plotNB.status.SetFields(\
1192                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1193                               
1194    try:
1195        plotNum = self.G2plotNB.plotList.index('2D Integration')
1196        Page = self.G2plotNB.nb.GetPage(plotNum)
1197        if not newPlot:
1198            Plot = Page.figure.gca()          #get previous plot & get limits
1199            xylim = Plot.get_xlim(),Plot.get_ylim()
1200        Page.figure.clf()
1201        Plot = Page.figure.gca()          #get a fresh plot after clf()
1202       
1203    except ValueError,error:
1204        Plot = self.G2plotNB.addMpl('2D Integration').gca()
1205        plotNum = self.G2plotNB.plotList.index('2D Integration')
1206        Page = self.G2plotNB.nb.GetPage(plotNum)
1207        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1208        Page.views = False
1209        view = False
1210    if not event:
1211        Page.SetFocus()
1212       
1213    Data = self.PatternTree.GetItemPyData(
1214        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1215    image = self.Integrate[0]
1216    xsc = self.Integrate[1]
1217    ysc = self.Integrate[2]
1218    Imin,Imax = Data['range'][1]
1219    acolor = mpl.cm.get_cmap(Data['color'])
1220    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1221    Plot.set_ylabel('azimuth',fontsize=12)
1222    Plot.set_xlabel('2-theta',fontsize=12)
1223    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
1224        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
1225    colorBar = Page.figure.colorbar(Img)
1226    if Data['setRings']:
1227        rings = np.concatenate((Data['rings']),axis=0)
1228        for xring,yring,dsp in rings:
1229            x,y = G2img.GetTthAzm(xring,yring,Data)
1230            Plot.plot(x,y,'r+')
1231    if Data['ellipses']:           
1232        for ellipse in Data['ellipses']:
1233            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1234            x,y = np.hsplit(ring,2)
1235            tth,azm = G2img.GetTthAzm(x,y,Data)
1236            Plot.plot(tth,azm,'b,')
1237    if not newPlot:
1238        Page.toolbar.push_current()
1239        Plot.set_xlim(xylim[0])
1240        Plot.set_ylim(xylim[1])
1241        xylim = []
1242        Page.toolbar.push_current()
1243        Page.toolbar.draw()
1244    else:
1245        Page.canvas.draw()
1246               
1247def PlotTRImage(self,tax,tay,taz,newPlot=False):
1248    #a test plot routine - not normally used
1249           
1250    def OnMotion(event):
1251        Page.canvas.SetToolTipString('')
1252        Page.canvas.SetCursor(wx.CROSS_CURSOR)
1253        azm = event.xdata
1254        tth = event.ydata
1255        if azm and tth:
1256            self.G2plotNB.status.SetFields(\
1257                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
1258                               
1259    try:
1260        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1261        Page = self.G2plotNB.nb.GetPage(plotNum)
1262        if not newPlot:
1263            Plot = Page.figure.gca()          #get previous plot & get limits
1264            xylim = Plot.get_xlim(),Plot.get_ylim()
1265        Page.figure.clf()
1266        Plot = Page.figure.gca()          #get a fresh plot after clf()
1267       
1268    except ValueError,error:
1269        Plot = self.G2plotNB.addMpl('2D Transformed Powder Image').gca()
1270        plotNum = self.G2plotNB.plotList.index('2D Transformed Powder Image')
1271        Page = self.G2plotNB.nb.GetPage(plotNum)
1272        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1273        Page.views = False
1274        view = False
1275    Page.SetFocus()
1276       
1277    Data = self.PatternTree.GetItemPyData(
1278        G2gd.GetPatternTreeItemId(self,self.Image, 'Image Controls'))
1279    Imin,Imax = Data['range'][1]
1280    step = (Imax-Imin)/5.
1281    V = np.arange(Imin,Imax,step)
1282    acolor = mpl.cm.get_cmap(Data['color'])
1283    Plot.set_title(self.PatternTree.GetItemText(self.Image)[4:])
1284    Plot.set_xlabel('azimuth',fontsize=12)
1285    Plot.set_ylabel('2-theta',fontsize=12)
1286    Plot.contour(tax,tay,taz,V,cmap=acolor)
1287    if Data['showLines']:
1288        IOtth = Data['IOtth']
1289        if Data['fullIntegrate']:
1290            LRAzim = [-180,180]
1291        else:
1292            LRAzim = Data['LRazimuth']                  #NB: integers
1293        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
1294        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
1295        if not Data['fullIntegrate']:
1296            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
1297            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
1298    if Data['setRings']:
1299        rings = np.concatenate((Data['rings']),axis=0)
1300        for xring,yring,dsp in rings:
1301            x,y = G2img.GetTthAzm(xring,yring,Data)
1302            Plot.plot(y,x,'r+')           
1303    if Data['ellipses']:           
1304        for ellipse in Data['ellipses']:
1305            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
1306            x,y = np.hsplit(ring,2)
1307            tth,azm = G2img.GetTthAzm(x,y,Data)
1308            Plot.plot(azm,tth,'b,')
1309    if not newPlot:
1310        Page.toolbar.push_current()
1311        Plot.set_xlim(xylim[0])
1312        Plot.set_ylim(xylim[1])
1313        xylim = []
1314        Page.toolbar.push_current()
1315        Page.toolbar.draw()
1316    else:
1317        Page.canvas.draw()
1318       
1319def PlotStructure(self,data):
1320    generalData = data['General']
1321    cell = generalData['Cell'][1:7]
1322    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
1323    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1324    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1325    Mydir = generalData['Mydir']
1326    atomData = data['Atoms']
1327    drawingData = data['Drawing']
1328    drawAtoms = drawingData['Atoms']
1329    cx,ct,cs = drawingData['atomPtrs']
1330    Wt = [255,255,255]
1331    Rd = [255,0,0]
1332    Gr = [0,255,0]
1333    Bl = [0,0,255]
1334    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]])
1335    uEdges = np.array([
1336        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
1337        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
1338        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
1339    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
1340    altDown = False
1341    shiftDown = False
1342    ctrlDown = False
1343   
1344    def OnKeyBox(event):
1345        import Image
1346        Draw()                          #make sure plot is fresh!!
1347        mode = cb.GetValue()
1348        Fname = Mydir+'\\'+generalData['Name']+'.'+mode
1349        size = Page.canvas.GetSize()
1350        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
1351        if mode in ['jpeg',]:
1352            Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
1353            im = Image.new("RGBA", (size[0],size[1]))
1354        else:
1355            Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
1356            im = Image.new("RGB", (size[0],size[1]))
1357        im.fromstring(Pix)
1358        im.save(Fname,mode)
1359        cb.SetValue(' Save as:')
1360        self.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
1361   
1362    def GetTruePosition(xy):
1363        View = glGetIntegerv(GL_VIEWPORT)
1364        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
1365        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
1366        Zmax = 1.
1367        for i,atom in enumerate(drawAtoms):
1368            x,y,z = atom[cx:cx+3]
1369            X,Y,Z = gluProject(x,y,z,Model,Proj,View)
1370            XY = [int(X),int(View[3]-Y)]
1371            if np.allclose(xy,XY,atol=10) and Z < Zmax:
1372                Zmax = Z
1373                SetSelectedAtoms(i)
1374                   
1375    def OnMouseDown(event):
1376        xy = event.GetPosition()
1377        if event.ShiftDown():
1378            GetTruePosition(xy)
1379        else:
1380            drawingData['Rotation'][3] = xy
1381        Draw()
1382       
1383    def OnMouseMove(event):
1384        newxy = event.GetPosition()
1385        if event.ControlDown() and drawingData['showABC']:
1386            if event.LeftIsDown():
1387                ctrlDown = True
1388                SetTestRot(newxy)
1389            elif event.RightIsDown():
1390                SetTestPos(newxy)
1391            elif event.MiddleIsDown():
1392                SetTestRotZ(newxy)
1393               
1394               
1395        if event.Dragging() and not event.ControlDown():
1396            if event.LeftIsDown():
1397                SetRotation(newxy)
1398            elif event.RightIsDown():
1399                SetTranslation(newxy)
1400            elif event.MiddleIsDown():
1401                SetRotationZ(newxy)
1402        Draw()
1403       
1404    def OnMouseWheel(event):
1405        drawingData['cameraPos'] += event.GetWheelRotation()/24
1406        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
1407        page = self.dataDisplay.GetSelection()
1408        if page:
1409            if self.dataDisplay.GetPageText(page) == 'Draw Options':
1410                panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1411                names = [child.GetName() for child in panel]
1412                panel[names.index('cameraPos')].SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
1413                panel[names.index('cameraSlider')].SetValue(drawingData['cameraPos'])
1414        Draw()
1415           
1416    def SetViewPointText(VP):
1417        page = self.dataDisplay.GetSelection()
1418        if page:
1419            if self.dataDisplay.GetPageText(page) == 'Draw Options':
1420                panel = self.dataDisplay.GetPage(page).GetChildren()[0].GetChildren()
1421                names = [child.GetName() for child in panel]
1422                panel[names.index('viewPoint')].SetValue('%.3f, %.3f, %.3f'%(VP[0],VP[1],VP[2]))
1423           
1424    def ClearSelectedAtoms():
1425        page = self.dataDisplay.GetSelection()
1426        if page:
1427            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1428                self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
1429            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1430                self.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
1431                   
1432    def SetSelectedAtoms(ind):
1433        page = self.dataDisplay.GetSelection()
1434        if page:
1435            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1436                self.dataDisplay.GetPage(page).SelectRow(ind)      #this is the Atoms grid in Draw Atoms
1437            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1438                Id = drawAtoms[ind][-2]
1439                for i,atom in enumerate(atomData):
1440                    if atom[-1] == Id:
1441                        self.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
1442                 
1443    def GetSelectedAtoms():
1444        page = self.dataDisplay.GetSelection()
1445        Ind = []
1446        if page:
1447            if self.dataDisplay.GetPageText(page) == 'Draw Atoms':
1448                Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
1449            elif self.dataDisplay.GetPageText(page) == 'Atoms':
1450                Ind = self.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
1451        return Ind
1452                                       
1453    def OnKey(event):           #on key UP!!
1454        keyCode = event.GetKeyCode()
1455        if keyCode > 255:
1456            keyCode = 0
1457        key,xyz = chr(keyCode),event.GetPosition()
1458        indx = drawingData['selectedAtoms']
1459        if key in ['c','C']:
1460            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
1461            drawingData['testPos'] = [[-.1,-.1,-.1],[0.0,0.0,0.0],[0,0]]
1462            drawingData['Rotation'] = [0.0,0.0,0.0,[]]
1463            SetViewPointText(drawingData['viewPoint'][0])
1464        elif key in ['n','N']:
1465            drawAtoms = drawingData['Atoms']
1466            pI = drawingData['viewPoint'][1]
1467            if indx:
1468                pI[0] = indx[pI[1]]
1469                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1470                pI[1] += 1
1471                if pI[1] >= len(indx):
1472                    pI[1] = 0
1473            else:
1474                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1475                pI[0] += 1
1476                if pI[0] >= len(drawAtoms):
1477                    pI[0] = 0
1478            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1479            SetViewPointText(drawingData['viewPoint'][0])
1480               
1481        elif key in ['p','P']:
1482            drawAtoms = drawingData['Atoms']
1483            pI = drawingData['viewPoint'][1]
1484            if indx:
1485                pI[0] = indx[pI[1]]
1486                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
1487                pI[1] -= 1
1488                if pI[1] < 0:
1489                    pI[1] = len(indx)-1
1490            else:
1491                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
1492                pI[0] -= 1
1493                if pI[0] < 0:
1494                    pI[0] = len(drawAtoms)-1
1495            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
1496            SetViewPointText(drawingData['viewPoint'][0])           
1497        Draw()
1498           
1499    def SetBackground():
1500        R,G,B,A = Page.camera['backColor']
1501        glClearColor(R,G,B,A)
1502        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
1503       
1504    def SetLights():
1505        glEnable(GL_DEPTH_TEST)
1506        glShadeModel(GL_SMOOTH)
1507        glEnable(GL_LIGHTING)
1508        glEnable(GL_LIGHT0)
1509        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
1510        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
1511        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
1512       
1513    def SetTranslation(newxy):
1514        Tx,Ty,Tz = drawingData['viewPoint'][0]
1515        anglex,angley,anglez,oldxy = drawingData['Rotation']
1516        Rx = G2lat.rotdMat(anglex,0)
1517        Ry = G2lat.rotdMat(angley,1)
1518        Rz = G2lat.rotdMat(anglez,2)
1519        dxy = list(newxy-oldxy)+[0,]
1520        dxy = np.inner(Bmat,np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy))))
1521        Tx -= dxy[0]*0.01
1522        Ty += dxy[1]*0.01
1523        Tz -= dxy[2]*0.01
1524        drawingData['Rotation'][3] = newxy
1525        drawingData['viewPoint'][0] =  Tx,Ty,Tz
1526        SetViewPointText([Tx,Ty,Tz])
1527       
1528    def SetTestPos(newxy):
1529        Tx,Ty,Tz = drawingData['testPos'][0]
1530        anglex,angley,anglez,oldxy = drawingData['Rotation']
1531        Rx = G2lat.rotdMat(anglex,0)
1532        Ry = G2lat.rotdMat(angley,1)
1533        Rz = G2lat.rotdMat(anglez,2)
1534        dxy = list(newxy-oldxy)+[0,]
1535        dxy = np.inner(Rz,np.inner(Ry,np.inner(Rx,dxy)))
1536        Tx += dxy[0]*0.001
1537        Ty -= dxy[1]*0.001
1538        Tz += dxy[2]*0.001
1539        drawingData['Rotation'][3] = newxy
1540        drawingData['testPos'][0] =  Tx,Ty,Tz
1541       
1542    def SetTestRot(newxy):
1543        Txyz = np.array(drawingData['testPos'][0])
1544        oldxy = drawingData['testPos'][2]
1545        Ax,Ay,Az = drawingData['testPos'][1]
1546        Vxyz = np.array(drawingData['viewPoint'][0])
1547        Dxyz = np.inner(Amat,Txyz-Vxyz)
1548        dxy = list(newxy-oldxy)+[0,]
1549        Ax += dxy[1]*0.01
1550        Ay += dxy[0]*0.01
1551        Rx = G2lat.rotdMat(Ax,0)
1552        Ry = G2lat.rotdMat(Ay,1)
1553        Dxyz = np.inner(Ry,np.inner(Rx,Dxyz))       
1554        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
1555        drawingData['testPos'][1] = [Ax,Ay,Az]
1556        drawingData['testPos'][2] = newxy
1557        drawingData['testPos'][0] = Dxyz
1558       
1559    def SetTestRotZ(newxy):
1560        Txyz = np.array(drawingData['testPos'][0])
1561        oldxy = drawingData['testPos'][2]
1562        Ax,Ay,Az = drawingData['testPos'][1]
1563        Vxyz = np.array(drawingData['viewPoint'][0])
1564        Dxyz = np.inner(Amat,Txyz-Vxyz)       
1565        dxy = list(newxy-oldxy)+[0,]
1566        Az += (dxy[0]+dxy[1])*.01
1567        Rz = G2lat.rotdMat(Az,2)
1568        Dxyz = np.inner(Rz,Dxyz)       
1569        Dxyz = np.inner(Bmat,Dxyz)+Vxyz
1570        drawingData['testPos'][1] = [Ax,Ay,Az]
1571        drawingData['testPos'][2] = newxy
1572        drawingData['testPos'][0] = Dxyz
1573                             
1574    def SetRotation(newxy):       
1575        anglex,angley,anglez,oldxy = drawingData['Rotation']
1576        dxy = newxy-oldxy
1577        anglex += dxy[1]*.25
1578        angley += dxy[0]*.25
1579        oldxy = newxy
1580        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
1581       
1582    def SetRotationZ(newxy):                       
1583        anglex,angley,anglez,oldxy = drawingData['Rotation']
1584        dxy = newxy-oldxy
1585        anglez += (dxy[0]+dxy[1])*.25
1586        oldxy = newxy
1587        drawingData['Rotation'] = [anglex,angley,anglez,oldxy]
1588       
1589    def RenderBox():
1590        glEnable(GL_COLOR_MATERIAL)
1591        glLineWidth(2)
1592        glEnable(GL_BLEND)
1593        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
1594        glEnable(GL_LINE_SMOOTH)
1595        glBegin(GL_LINES)
1596        for line,color in zip(uEdges,uColors):
1597            glColor3ubv(color)
1598            glVertex3fv(line[0])
1599            glVertex3fv(line[1])
1600        glEnd()
1601        glColor4ubv([0,0,0,0])
1602        glDisable(GL_LINE_SMOOTH)
1603        glDisable(GL_BLEND)
1604        glDisable(GL_COLOR_MATERIAL)
1605       
1606    def RenderUnitVectors(x,y,z):
1607        xyz = np.array([x,y,z])
1608        glEnable(GL_COLOR_MATERIAL)
1609        glLineWidth(1)
1610        glPushMatrix()
1611        glTranslate(x,y,z)
1612        glScalef(1/cell[0],1/cell[1],1/cell[2])
1613        glBegin(GL_LINES)
1614        for line,color in zip(uEdges,uColors)[:3]:
1615            glColor3ubv(color)
1616            glVertex3fv(line[0])
1617            glVertex3fv(line[1])
1618        glEnd()
1619        glPopMatrix()
1620        glColor4ubv([0,0,0,0])
1621        glDisable(GL_COLOR_MATERIAL)
1622               
1623    def RenderSphere(x,y,z,radius,color):
1624        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1625        glPushMatrix()
1626        glTranslate(x,y,z)
1627        glMultMatrixf(B4mat.T)
1628        q = gluNewQuadric()
1629        gluSphere(q,radius,20,10)
1630        glPopMatrix()
1631       
1632    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
1633        s1,s2,s3 = E
1634        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1635        glPushMatrix()
1636        glTranslate(x,y,z)
1637        glMultMatrixf(B4mat.T)
1638        glMultMatrixf(R4.T)
1639        glEnable(GL_NORMALIZE)
1640        glScale(s1,s2,s3)
1641        q = gluNewQuadric()
1642        gluSphere(q,ellipseProb,20,10)
1643        glDisable(GL_NORMALIZE)
1644        glPopMatrix()
1645       
1646    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
1647        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1648        glPushMatrix()
1649        glTranslate(x,y,z)
1650        glMultMatrixf(B4mat.T)
1651        for bond in Bonds:
1652            glPushMatrix()
1653            Dx = np.inner(Amat,bond)
1654            Z = np.sqrt(np.sum(Dx**2))
1655            azm = atan2d(-Dx[1],-Dx[0])
1656            phi = acosd(Dx[2]/Z)
1657            glRotate(-azm,0,0,1)
1658            glRotate(phi,1,0,0)
1659            q = gluNewQuadric()
1660            gluCylinder(q,radius,radius,Z,slice,2)
1661            glPopMatrix()           
1662        glPopMatrix()
1663               
1664    def RenderLines(x,y,z,Bonds,color):
1665        xyz = np.array([x,y,z])
1666        glEnable(GL_COLOR_MATERIAL)
1667        glLineWidth(1)
1668        glColor3fv(color)
1669        glPushMatrix()
1670        glBegin(GL_LINES)
1671        for bond in Bonds:
1672            glVertex3fv(xyz)
1673            glVertex3fv(xyz+bond)
1674        glEnd()
1675        glColor4ubv([0,0,0,0])
1676        glPopMatrix()
1677        glDisable(GL_COLOR_MATERIAL)
1678       
1679    def RenderPolyhedra(x,y,z,Faces,color):
1680        glPushMatrix()
1681        glTranslate(x,y,z)
1682        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
1683        glShadeModel(GL_SMOOTH)
1684        glMultMatrixf(B4mat.T)
1685        for face,norm in Faces:
1686            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
1687            glFrontFace(GL_CW)
1688            glNormal3fv(norm)
1689            glBegin(GL_TRIANGLES)
1690            for vert in face:
1691                glVertex3fv(vert)
1692            glEnd()
1693        glPopMatrix()
1694       
1695    def RenderBackbone(Backbone,BackboneColor,radius):
1696        glPushMatrix()
1697        glMultMatrixf(B4mat.T)
1698        glEnable(GL_COLOR_MATERIAL)
1699        glShadeModel(GL_SMOOTH)
1700        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
1701        glePolyCylinder(Backbone,BackboneColor,radius)
1702        glPopMatrix()       
1703        glDisable(GL_COLOR_MATERIAL)
1704       
1705    def RenderLabel(x,y,z,label,r):       
1706        glPushMatrix()
1707        glTranslate(x,y,z)
1708        glMultMatrixf(B4mat.T)
1709        glDisable(GL_LIGHTING)
1710        glColor3f(0,1.,0)
1711        glRasterPos3f(r,r,r)
1712        for c in list(label):
1713            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
1714        glEnable(GL_LIGHTING)
1715        glPopMatrix()
1716                           
1717    def Draw():
1718        import numpy.linalg as nl
1719        Ind = GetSelectedAtoms()
1720        VS = np.array(Page.canvas.GetSize())
1721        aspect = float(VS[0])/float(VS[1])
1722        cPos = drawingData['cameraPos']
1723        Zclip = drawingData['Zclip']*cPos/200.
1724        anglex,angley,anglez = drawingData['Rotation'][:3]
1725        Tx,Ty,Tz = drawingData['viewPoint'][0]
1726        cx,ct,cs = drawingData['atomPtrs']
1727        bondR = drawingData['bondRadius']
1728        G,g = G2lat.cell2Gmat(cell)
1729        GS = G
1730        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
1731        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
1732        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
1733        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
1734       
1735        SetBackground()
1736        glInitNames()
1737        glPushName(0)
1738       
1739        glMatrixMode(GL_PROJECTION)
1740        glLoadIdentity()
1741        glViewport(0,0,VS[0],VS[1])
1742        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
1743        gluLookAt(0,0,cPos,0,0,0,0,1,0)
1744        SetLights()           
1745           
1746        glMatrixMode(GL_MODELVIEW)
1747        glLoadIdentity()
1748        glRotate(anglez,0,0,1)
1749        glRotate(anglex,cosd(anglez),-sind(anglez),0)
1750        glRotate(angley,sind(anglez),cosd(anglez),0)
1751        glMultMatrixf(A4mat.T)
1752        glTranslate(-Tx,-Ty,-Tz)
1753        if drawingData['unitCellBox']:
1754            RenderBox()
1755        if drawingData['showABC']:
1756#            try:            #temporary fix - not needed further?
1757#                x,y,z = drawingData['testPos'][0]
1758#            except TypeError:
1759#                x,y,z = drawingData['testPos']
1760            x,y,z = drawingData['testPos'][0]
1761            if altDown:
1762                self.G2plotNB.status.SetStatusText('moving test point %.4f,%.4f,%.4f'%(x,y,z),1)
1763            else:
1764                self.G2plotNB.status.SetStatusText('test point %.4f,%.4f,%.4f'%(x,y,z),1)           
1765            RenderUnitVectors(x,y,z)
1766        Backbone = []
1767        BackboneColor = []
1768        time0 = time.time()
1769        for iat,atom in enumerate(drawingData['Atoms']):
1770            x,y,z = atom[cx:cx+3]
1771            Bonds = atom[-2]
1772            Faces = atom[-1]
1773            try:
1774                atNum = generalData['AtomTypes'].index(atom[ct])
1775            except ValueError:
1776                atNum = -1
1777            CL = atom[cs+2]
1778            color = np.array(CL)/255.
1779            if iat in Ind:
1780                color = np.array(Gr)/255.
1781            radius = 0.5
1782            if atom[cs] != '':
1783                glLoadName(atom[-3])                   
1784            if 'balls' in atom[cs]:
1785                vdwScale = drawingData['vdwScale']
1786                ballScale = drawingData['ballScale']
1787                if atNum < 0:
1788                    radius = 0.2
1789                elif 'H' == atom[ct]:
1790                    if drawingData['showHydrogen']:
1791                        if 'vdW' in atom[cs] and atNum >= 0:
1792                            radius = vdwScale*generalData['vdWRadii'][atNum]
1793                        else:
1794                            radius = ballScale*drawingData['sizeH']
1795                    else:
1796                        radius = 0.0
1797                else:
1798                    if 'vdW' in atom[cs]:
1799                        radius = vdwScale*generalData['vdWRadii'][atNum]
1800                    else:
1801                        radius = ballScale*generalData['BondRadii'][atNum]
1802                RenderSphere(x,y,z,radius,color)
1803                if 'sticks' in atom[cs]:
1804                    RenderBonds(x,y,z,Bonds,bondR,color)
1805            elif 'ellipsoids' in atom[cs]:
1806                RenderBonds(x,y,z,Bonds,bondR,color)
1807                if atom[cs+3] == 'A':                   
1808                    Uij = atom[cs+5:cs+11]
1809                    U = np.multiply(G2spc.Uij2U(Uij),GS)
1810                    U = np.inner(Amat,np.inner(U,Amat).T)
1811                    E,R = nl.eigh(U)
1812                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
1813                    E = np.sqrt(E)
1814                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
1815                        pass
1816                    else:
1817                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,color)                   
1818                else:
1819                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
1820                        pass
1821                    else:
1822                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
1823                        RenderSphere(x,y,z,radius,color)
1824            elif 'lines' in atom[cs]:
1825                radius = 0.1
1826                RenderLines(x,y,z,Bonds,color)
1827#                RenderBonds(x,y,z,Bonds,0.05,color,6)
1828            elif atom[cs] == 'sticks':
1829                radius = 0.1
1830                RenderBonds(x,y,z,Bonds,bondR,color)
1831            elif atom[cs] == 'polyhedra':
1832                RenderPolyhedra(x,y,z,Faces,color)
1833            elif atom[cs] == 'backbone':
1834                if atom[ct-1].split()[0] in ['C','N']:
1835                    Backbone.append(list(np.inner(Amat,np.array([x,y,z]))))
1836                    BackboneColor.append(list(color))
1837                   
1838            if atom[cs+1] == 'type':
1839                RenderLabel(x,y,z,atom[ct],radius)
1840            elif atom[cs+1] == 'name':
1841                RenderLabel(x,y,z,atom[ct-1],radius)
1842            elif atom[cs+1] == 'number':
1843                RenderLabel(x,y,z,str(iat+1),radius)
1844            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
1845                RenderLabel(x,y,z,atom[ct-4],radius)
1846            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
1847                RenderLabel(x,y,z,atom[ct-3],radius)
1848            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
1849                RenderLabel(x,y,z,atom[ct-2],radius)
1850        if Backbone:
1851            RenderBackbone(Backbone,BackboneColor,bondR)
1852#        print time.time()-time0
1853        Page.canvas.SwapBuffers()
1854       
1855    def OnSize(event):
1856        Draw()
1857       
1858    def OnFocus(event):
1859        Draw()
1860       
1861    try:
1862        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
1863        Page = self.G2plotNB.nb.GetPage(plotNum)       
1864    except (ValueError,error):
1865        Plot = self.G2plotNB.addOgl(generalData['Name'])
1866        plotNum = self.G2plotNB.plotList.index(generalData['Name'])
1867        Page = self.G2plotNB.nb.GetPage(plotNum)
1868        Page.views = False
1869        view = False
1870        altDown = False
1871    Page.SetFocus()
1872    cb = wx.ComboBox(self.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,
1873        choices=(' save as:','jpeg','tiff','bmp'))
1874    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
1875    cb.SetValue(' save as:')
1876    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
1877    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
1878    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
1879    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
1880    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
1881    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
1882    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
1883    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
1884    Page.camera['position'] = drawingData['cameraPos']
1885    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
1886    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
1887    Page.canvas.SetCurrent()
1888    Draw()
1889       
Note: See TracBrowser for help on using the repository browser.