source: trunk/GSASIIplot.py @ 311

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

finish new indexing refinement stuff
more on texture plotting

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