source: trunk/GSASIIplot.py @ 951

Last change on this file since 951 was 951, checked in by vondreele, 9 years ago

more on MC/SA

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 158.4 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIplot: plotting routines*
4===============================
5
6'''
7########### SVN repository information ###################
8# $Date: 2013-06-14 19:38:20 +0000 (Fri, 14 Jun 2013) $
9# $Author: vondreele $
10# $Revision: 951 $
11# $URL: trunk/GSASIIplot.py $
12# $Id: GSASIIplot.py 951 2013-06-14 19:38:20Z vondreele $
13########### SVN repository information ###################
14import math
15import time
16import copy
17import sys
18import os.path
19import numpy as np
20import numpy.ma as ma
21import numpy.linalg as nl
22import wx
23import wx.aui
24import wx.glcanvas
25import matplotlib as mpl
26import mpl_toolkits.mplot3d.axes3d as mp3d
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 951 $")
29import GSASIIgrid as G2gd
30import GSASIIimage as G2img
31import GSASIIpwd as G2pwd
32import GSASIIIO as G2IO
33import GSASIIpwdGUI as G2pdG
34import GSASIIimgGUI as G2imG
35import GSASIIphsGUI as G2phG
36import GSASIIlattice as G2lat
37import GSASIIspc as G2spc
38import GSASIImath as G2mth
39import pytexture as ptx
40from  OpenGL.GL import *
41from OpenGL.GLU import *
42from OpenGL.GLUT import *
43from OpenGL.GLE import *
44from matplotlib.backends.backend_wx import _load_bitmap
45from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as Canvas
46from matplotlib.backends.backend_wxagg import NavigationToolbar2Wx as Toolbar
47
48# useful degree trig functions
49sind = lambda x: math.sin(x*math.pi/180.)
50cosd = lambda x: math.cos(x*math.pi/180.)
51tand = lambda x: math.tan(x*math.pi/180.)
52asind = lambda x: 180.*math.asin(x)/math.pi
53acosd = lambda x: 180.*math.acos(x)/math.pi
54atan2d = lambda x,y: 180.*math.atan2(y,x)/math.pi
55atand = lambda x: 180.*math.atan(x)/math.pi
56# numpy versions
57npsind = lambda x: np.sin(x*np.pi/180.)
58npcosd = lambda x: np.cos(x*np.pi/180.)
59nptand = lambda x: np.tan(x*np.pi/180.)
60npacosd = lambda x: 180.*np.arccos(x)/np.pi
61npasind = lambda x: 180.*np.arcsin(x)/np.pi
62npatand = lambda x: 180.*np.arctan(x)/np.pi
63npatan2d = lambda x,y: 180.*np.arctan2(x,y)/np.pi
64   
65class G2PlotMpl(wx.Panel):   
66    'needs a doc string'
67    def __init__(self,parent,id=-1,dpi=None,**kwargs):
68        wx.Panel.__init__(self,parent,id=id,**kwargs)
69        mpl.rcParams['legend.fontsize'] = 10
70        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(5,6))
71        self.canvas = Canvas(self,-1,self.figure)
72        self.toolbar = GSASIItoolbar(self.canvas)
73
74        self.toolbar.Realize()
75       
76        sizer=wx.BoxSizer(wx.VERTICAL)
77        sizer.Add(self.canvas,1,wx.EXPAND)
78        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
79        self.SetSizer(sizer)
80       
81class G2PlotOgl(wx.Panel):
82    'needs a doc string'
83    def __init__(self,parent,id=-1,dpi=None,**kwargs):
84        self.figure = wx.Panel.__init__(self,parent,id=id,**kwargs)
85        if 'win' in sys.platform:           #Windows already double buffered
86            self.canvas = wx.glcanvas.GLCanvas(self,-1,**kwargs)
87        else:                               #fix from Jim Hester for X systems
88            attribs = (wx.glcanvas.WX_GL_DOUBLEBUFFER,)         
89            self.canvas = wx.glcanvas.GLCanvas(self,-1,attribList=attribs,**kwargs)
90        self.camera = {}
91        sizer=wx.BoxSizer(wx.VERTICAL)
92        sizer.Add(self.canvas,1,wx.EXPAND)
93        self.SetSizer(sizer)
94       
95class G2Plot3D(wx.Panel):
96    'needs a doc string'
97    def __init__(self,parent,id=-1,dpi=None,**kwargs):
98        wx.Panel.__init__(self,parent,id=id,**kwargs)
99        self.figure = mpl.figure.Figure(dpi=dpi,figsize=(6,6))
100        self.canvas = Canvas(self,-1,self.figure)
101        self.toolbar = GSASIItoolbar(self.canvas)
102
103        self.toolbar.Realize()
104       
105        sizer=wx.BoxSizer(wx.VERTICAL)
106        sizer.Add(self.canvas,1,wx.EXPAND)
107        sizer.Add(self.toolbar,0,wx.LEFT|wx.EXPAND)
108        self.SetSizer(sizer)
109                             
110class G2PlotNoteBook(wx.Panel):
111    'create a tabbed window for plotting'
112    def __init__(self,parent,id=-1):
113        wx.Panel.__init__(self,parent,id=id)
114        #so one can't delete a plot page!!
115        self.nb = wx.aui.AuiNotebook(self, \
116            style=wx.aui.AUI_NB_DEFAULT_STYLE ^ wx.aui.AUI_NB_CLOSE_ON_ACTIVE_TAB)
117        sizer = wx.BoxSizer()
118        sizer.Add(self.nb,1,wx.EXPAND)
119        self.SetSizer(sizer)
120        self.status = parent.CreateStatusBar()
121        self.status.SetFieldsCount(2)
122        self.status.SetStatusWidths([150,-1])
123        self.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED, self.OnPageChanged)
124       
125        self.plotList = []
126           
127    def addMpl(self,name=""):
128        'Add a tabbed page with a matplotlib plot'
129        page = G2PlotMpl(self.nb)
130        self.nb.AddPage(page,name)
131       
132        self.plotList.append(name)
133       
134        return page.figure
135       
136    def add3D(self,name=""):
137        'Add a tabbed page with a 3D plot'
138        page = G2Plot3D(self.nb)
139        self.nb.AddPage(page,name)
140       
141        self.plotList.append(name)
142       
143        return page.figure
144       
145    def addOgl(self,name=""):
146        'Add a tabbed page with an openGL plot'
147        page = G2PlotOgl(self.nb)
148        self.nb.AddPage(page,name)
149       
150        self.plotList.append(name)
151       
152        return page.figure
153       
154    def Delete(self,name):
155        'delete a tabbed page'
156        try:
157            item = self.plotList.index(name)
158            del self.plotList[item]
159            self.nb.DeletePage(item)
160        except ValueError:          #no plot of this name - do nothing
161            return     
162               
163    def clear(self):
164        'clear all pages from plot window'
165        while self.nb.GetPageCount():
166            self.nb.DeletePage(0)
167        self.plotList = []
168        self.status.DestroyChildren()
169       
170    def Rename(self,oldName,newName):
171        'rename a tab'
172        try:
173            item = self.plotList.index(oldName)
174            self.plotList[item] = newName
175            self.nb.SetPageText(item,newName)
176        except ValueError:          #no plot of this name - do nothing
177            return     
178       
179    def OnPageChanged(self,event):
180        'respond to someone pressing a tab on the plot window'
181        if self.plotList:
182            self.status.SetStatusText('Better to select this from GSAS-II data tree',1)
183        self.status.DestroyChildren()                           #get rid of special stuff on status bar
184       
185class GSASIItoolbar(Toolbar):
186    'needs a doc string'
187    ON_MPL_HELP = wx.NewId()
188    ON_MPL_KEY = wx.NewId()
189    def __init__(self,plotCanvas):
190        Toolbar.__init__(self,plotCanvas)
191        POSITION_OF_CONFIGURE_SUBPLOTS_BTN = 6
192        self.DeleteToolByPos(POSITION_OF_CONFIGURE_SUBPLOTS_BTN)
193        parent = self.GetParent()
194        key = os.path.join(os.path.split(__file__)[0],'key.ico')
195        self.AddSimpleTool(self.ON_MPL_KEY,_load_bitmap(key),'Key press','Select key press')
196        wx.EVT_TOOL(self,self.ON_MPL_KEY,self.OnKey)
197        help = os.path.join(os.path.split(__file__)[0],'help.ico')
198        self.AddSimpleTool(self.ON_MPL_HELP,_load_bitmap(help),'Help on','Show help on')
199        wx.EVT_TOOL(self,self.ON_MPL_HELP,self.OnHelp)
200    def OnHelp(self,event):
201        'needs a doc string'
202        Page = self.GetParent().GetParent()
203        pageNo = Page.GetSelection()
204        bookmark = Page.GetPageText(pageNo)
205        bookmark = bookmark.strip(')').replace('(','_')
206        G2gd.ShowHelp(bookmark,self.TopLevelParent)
207    def OnKey(self,event):
208        'needs a doc string'
209        parent = self.GetParent()
210        if parent.Choice:
211            dlg = wx.SingleChoiceDialog(parent,'Select','Key press',list(parent.Choice))
212            if dlg.ShowModal() == wx.ID_OK:
213                sel = dlg.GetSelection()
214                event.key = parent.Choice[sel][0]
215                parent.keyPress(event)
216            dlg.Destroy()
217   
218################################################################################
219##### PlotSngl
220################################################################################
221           
222def PlotSngl(self,newPlot=False):
223    '''Single crystal structure factor plotting package - displays zone of reflections as rings proportional
224        to F, F**2, etc. as requested
225    '''
226    from matplotlib.patches import Circle
227    global HKL,HKLF
228
229    def OnSCMotion(event):
230        xpos = event.xdata
231        if xpos:
232            xpos = round(xpos)                                        #avoid out of frame mouse position
233            ypos = round(event.ydata)
234            zpos = Data['Layer']
235            if '100' in Data['Zone']:
236                HKLtxt = '(%3d,%3d,%3d)'%(zpos,xpos,ypos)
237            elif '010' in Data['Zone']:
238                HKLtxt = '(%3d,%3d,%3d)'%(xpos,zpos,ypos)
239            elif '001' in Data['Zone']:
240                HKLtxt = '(%3d,%3d,%3d)'%(xpos,ypos,zpos)
241            Page.canvas.SetToolTipString(HKLtxt)
242            self.G2plotNB.status.SetFields(['HKL = '+HKLtxt,''])
243               
244    def OnSCPick(event):
245        zpos = Data['Layer']
246        pos = event.artist.center
247        if '100' in Data['Zone']:
248            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(zpos,pos[0],pos[1]))
249            hkl = np.array([zpos,pos[0],pos[1]])
250        elif '010' in Data['Zone']:
251            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],zpos,pos[1]))
252            hkl = np.array([pos[0],zpos,pos[1]])
253        elif '001' in Data['Zone']:
254            Page.canvas.SetToolTipString('(picked:(%3d,%3d,%3d))'%(pos[0],pos[1],zpos))
255            hkl = np.array([pos[0],pos[1],zpos])
256        h,k,l = hkl
257        hklf = HKLF[np.where(np.all(HKL-hkl == [0,0,0],axis=1))]
258        if len(hklf):
259            Fosq,sig,Fcsq = hklf[0]
260            HKLtxt = '(%3d,%3d,%3d %.2f %.3f %.2f %.2f)'%(h,k,l,Fosq,sig,Fcsq,(Fosq-Fcsq)/(scale*sig))
261            self.G2plotNB.status.SetFields(['','HKL, Fosq, sig, Fcsq, delFsq/sig = '+HKLtxt])
262                                 
263    try:
264        plotNum = self.G2plotNB.plotList.index('Structure Factors')
265        Page = self.G2plotNB.nb.GetPage(plotNum)
266        if not newPlot:
267            Plot = Page.figure.gca()          #get previous powder plot & get limits
268            xylim = Plot.get_xlim(),Plot.get_ylim()
269        Page.figure.clf()
270        Page.Choice = None
271        Plot = Page.figure.gca()          #get a fresh plot after clf()
272    except ValueError:
273        Plot = self.G2plotNB.addMpl('Structure Factors').gca()
274        plotNum = self.G2plotNB.plotList.index('Structure Factors')
275        Page = self.G2plotNB.nb.GetPage(plotNum)
276        Page.canvas.mpl_connect('pick_event', OnSCPick)
277        Page.canvas.mpl_connect('motion_notify_event', OnSCMotion)
278    Page.Choice = None
279    Page.SetFocus()
280   
281    Plot.set_aspect(aspect='equal')
282    HKLref = self.PatternTree.GetItemPyData(self.Sngl)[1]
283    Data = self.PatternTree.GetItemPyData( 
284        G2gd.GetPatternTreeItemId(self,self.Sngl, 'HKL Plot Controls'))
285    Type = Data['Type']           
286    scale = Data['Scale']
287    HKLmax = Data['HKLmax']
288    HKLmin = Data['HKLmin']
289    FosqMax = Data['FoMax']
290    FoMax = math.sqrt(FosqMax)
291    xlabel = ['k, h=','h, k=','h, l=']
292    ylabel = ['l','l','k']
293    zones = ['100','010','001']
294    pzone = [[1,2],[0,2],[0,1]]
295    izone = zones.index(Data['Zone'])
296    Plot.set_title(self.PatternTree.GetItemText(self.Sngl)[5:])
297    HKL = []
298    HKLF = []
299    for refl in HKLref:
300        H = np.array(refl[:3])
301        Fosq,sig,Fcsq = refl[5:8]
302        HKL.append(H)
303        HKLF.append([Fosq,sig,Fcsq])
304        if H[izone] == Data['Layer']:
305            A = 0
306            B = 0
307            if Type == 'Fosq':
308                A = scale*Fosq/FosqMax
309                B = scale*Fcsq/FosqMax
310                C = abs(A-B)
311            elif Type == 'Fo':
312                A = scale*math.sqrt(max(0,Fosq))/FoMax
313                B = scale*math.sqrt(max(0,Fcsq))/FoMax
314                C = abs(A-B)
315            elif Type == '|DFsq|/sig':
316                if sig > 0.:
317                    A = scale*(Fosq-Fcsq)/(sig)
318                B = 0
319            elif Type == '|DFsq|>sig':
320                if sig > 0.:
321                    A = (Fosq-Fcsq)/sig
322                if abs(A) < 1.0: A = 0
323                A *= scale
324                B = 0                   
325            elif Type == '|DFsq|>3sig':
326                if sig > 0.:
327                    A = (Fosq-Fcsq)/sig
328                if abs(A) < 3.0: A = 0
329                A *= scale
330                B = 0                   
331            xy = (H[pzone[izone][0]],H[pzone[izone][1]])
332            if Type in ['|DFsq|/sig','|DFsq|>sig','|DFsq|>3sig']:
333                if A > 0.0:
334                    Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w',picker=3))
335                else:
336                    Plot.add_artist(Circle(xy,radius=-A,ec='r',fc='w',picker=3))
337            else:
338                if A > 0.0 and A > B:
339                    Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w',picker=3))
340                if B:
341                    Plot.add_artist(Circle(xy,radius=B,ec='b',fc='w'))
342                    if A < B:
343                        Plot.add_artist(Circle(xy,radius=A,ec='g',fc='w',picker=3))
344                    radius = C
345                    if radius > 0:
346                        if A > B:
347                            Plot.add_artist(Circle(xy,radius=radius,ec='g',fc='g'))
348                        else:                   
349                            Plot.add_artist(Circle(xy,radius=radius,ec='r',fc='r'))
350    HKL = np.array(HKL,dtype=np.int)
351    HKLF = np.array(HKLF)
352    Plot.set_xlabel(xlabel[izone]+str(Data['Layer']),fontsize=12)
353    Plot.set_ylabel(ylabel[izone],fontsize=12)
354    Plot.set_xlim((HKLmin[pzone[izone][0]],HKLmax[pzone[izone][0]]))
355    Plot.set_ylim((HKLmin[pzone[izone][1]],HKLmax[pzone[izone][1]]))
356    Page.canvas.draw()
357#    if not newPlot:
358#        Page.toolbar.push_current()
359#        Plot.set_xlim(xylim[0])
360#        Plot.set_ylim(xylim[1])
361#        xylim = []
362#        Page.toolbar.push_current()
363#        Page.toolbar.draw()
364#    else:
365#        Page.canvas.draw()
366       
367################################################################################
368##### PlotPatterns
369################################################################################
370           
371def PlotPatterns(G2frame,newPlot=False):
372    '''Powder pattern plotting package - displays single or multiple powder patterns as intensity vs
373    2-theta, q or TOF. Can display multiple patterns as "waterfall plots" or contour plots. Log I
374    plotting available.
375    '''
376    global HKL
377
378#    def OnKeyBox(event):
379#        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('Powder Patterns'):
380#            event.key = cb.GetValue()[0]
381#            cb.SetValue(' key press')
382#            wx.CallAfter(OnPlotKeyPress,event)
383#                       
384    def OnPlotKeyPress(event):
385        newPlot = False
386        if event.key == 'w':
387            if G2frame.Weight:
388                G2frame.Weight = False
389            else:
390                G2frame.Weight = True
391                G2frame.SinglePlot = True
392            newPlot = True
393        elif event.key == 'b':
394            if G2frame.SubBack:
395                G2frame.SubBack = False
396            else:
397                G2frame.SubBack = True
398                G2frame.SinglePlot = True               
399        elif event.key == 'n':
400            if G2frame.Contour:
401                pass
402            else:
403                if G2frame.logPlot:
404                    G2frame.logPlot = False
405                else:
406                    G2frame.Offset[0] = 0
407                    G2frame.logPlot = True
408                newPlot = True
409        elif event.key == 'u':
410            if G2frame.Contour:
411                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
412            elif G2frame.logPlot:
413                pass
414            elif G2frame.Offset[0] < 100.:
415                G2frame.Offset[0] += 1.
416        elif event.key == 'd':
417            if G2frame.Contour:
418                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
419            elif G2frame.logPlot:
420                pass
421            elif G2frame.Offset[0] > 0.:
422                G2frame.Offset[0] -= 1.
423        elif event.key == 'l':
424            G2frame.Offset[1] -= 1.
425        elif event.key == 'r':
426            G2frame.Offset[1] += 1.
427        elif event.key == 'o':
428            G2frame.Offset = [0,0]
429        elif event.key == 'c':
430            newPlot = True
431            if G2frame.Contour:
432                G2frame.Contour = False
433            else:
434                G2frame.Contour = True
435                G2frame.SinglePlot = False
436                G2frame.Offset = [0.,0.]
437        elif event.key == 'q':
438            newPlot = True
439            if G2frame.qPlot:
440                G2frame.qPlot = False
441            else:
442                G2frame.qPlot = True
443        elif event.key == 's':
444            if G2frame.Contour:
445                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
446                choice.sort()
447                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
448                if dlg.ShowModal() == wx.ID_OK:
449                    sel = dlg.GetSelection()
450                    G2frame.ContourColor = choice[sel]
451                else:
452                    G2frame.ContourColor = 'Paired'
453                dlg.Destroy()
454            else:               
455                if G2frame.SinglePlot:
456                    G2frame.SinglePlot = False
457                else:
458                    G2frame.SinglePlot = True
459            newPlot = True
460        elif event.key == '+':
461            if G2frame.PickId:
462                G2frame.PickId = False
463        elif event.key == 'i':                  #for smoothing contour plot
464            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
465               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
466               'mitchell','sinc','lanczos']
467            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
468            if dlg.ShowModal() == wx.ID_OK:
469                sel = dlg.GetSelection()
470                G2frame.Interpolate = choice[sel]
471            else:
472                G2frame.Interpolate = 'nearest'
473            dlg.Destroy()
474           
475        PlotPatterns(G2frame,newPlot=newPlot)
476       
477    def OnMotion(event):
478        xpos = event.xdata
479        if xpos:                                        #avoid out of frame mouse position
480            ypos = event.ydata
481            Page.canvas.SetCursor(wx.CROSS_CURSOR)
482            try:
483                Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
484                if 'C' in Parms['Type'][0]:
485                    wave = G2mth.getWave(Parms)
486                    if G2frame.qPlot:
487                        try:
488                            xpos = 2.0*asind(xpos*wave/(4*math.pi))
489                        except ValueError:      #avoid bad value in asin beyond upper limit
490                            pass
491                    dsp = 0.0
492                    if abs(xpos) > 0.:                  #avoid possible singularity at beam center
493                        dsp = wave/(2.*sind(abs(xpos)/2.0))
494                    if G2frame.Contour:
495                        G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f pattern ID =%5d'%(xpos,dsp,int(ypos)),1)
496                    else:
497                        G2frame.G2plotNB.status.SetStatusText('2-theta =%9.3f d =%9.5f Intensity =%9.1f'%(xpos,dsp,ypos),1)
498                else:       #TOF neutrons
499                    dsp = 0.0
500                    difC = Parms['difC'][1]
501                    dsp = xpos/difC             #rough approx.!
502                    if G2frame.Contour:
503                        G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f pattern ID =%5d'%(xpos,dsp,int(ypos)),1)
504                    else:
505                        G2frame.G2plotNB.status.SetStatusText('TOF =%9.3f d =%9.5f Intensity =%9.1f'%(xpos,dsp,ypos),1)
506                if G2frame.itemPicked:
507                    Page.canvas.SetToolTipString('%9.3f'%(xpos))
508                if G2frame.PickId:
509                    found = []
510                    if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists'] or \
511                        'PWDR' in G2frame.PatternTree.GetItemText(PickId):
512                        if len(HKL):
513                            view = Page.toolbar._views.forward()[0][:2]
514                            wid = view[1]-view[0]
515                            found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
516                        if len(found):
517                            h,k,l = found[0][:3] 
518                            Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
519                        else:
520                            Page.canvas.SetToolTipString('')
521
522            except TypeError:
523                G2frame.G2plotNB.status.SetStatusText('Select PWDR powder pattern first',1)
524                                                   
525    def OnPick(event):
526        if G2frame.itemPicked is not None: return
527        PatternId = G2frame.PatternId
528        try:
529            Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
530        except TypeError:
531            return
532        if 'C' in Parms['Type'][0]:
533            wave = G2mth.getWave(Parms)
534        else:
535            difC = Parms['difC'][1]
536        PickId = G2frame.PickId
537        pick = event.artist
538        mouse = event.mouseevent       
539        xpos = pick.get_xdata()
540        ypos = pick.get_ydata()
541        ind = event.ind
542        xy = list(zip(np.take(xpos,ind),np.take(ypos,ind))[0])
543        if G2frame.PatternTree.GetItemText(PickId) == 'Peak List':
544            if ind.all() != [0]:                                    #picked a data point
545                data = G2frame.PatternTree.GetItemPyData(G2frame.PickId)
546                XY = G2mth.setPeakparms(Parms,Parms2,xy[0],xy[1])           #what happens for a q-plot???
547                data.append(XY)
548                G2pdG.UpdatePeakGrid(G2frame,data)
549                PlotPatterns(G2frame)
550            else:                                                   #picked a peak list line
551                G2frame.itemPicked = pick
552        elif G2frame.PatternTree.GetItemText(PickId) == 'Limits':
553            if ind.all() != [0]:                                    #picked a data point
554                LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
555                data = G2frame.PatternTree.GetItemPyData(LimitId)
556                if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
557                    if G2frame.qPlot:                              #qplot - convert back to 2-theta
558                        xy[0] = 2.0*asind(xy[0]*wave/(4*math.pi))
559                if mouse.button==1:
560                    data[1][0] = min(xy[0],data[1][1])
561                if mouse.button==3:
562                    data[1][1] = max(xy[0],data[1][0])
563                G2frame.PatternTree.SetItemPyData(LimitId,data)
564                G2pdG.UpdateLimitsGrid(G2frame,data)
565                PlotPatterns(G2frame)
566            else:                                                   #picked a limit line
567                G2frame.itemPicked = pick
568        elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \
569            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
570            G2frame.itemPicked = pick
571            pick = str(pick)
572       
573    def OnRelease(event):
574        if G2frame.itemPicked is None: return
575        Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Instrument Parameters'))
576        if 'C' in Parms['Type'][0]:
577            wave = G2mth.getWave(Parms)
578        else:
579            difC = Parms['difC'][1]
580        xpos = event.xdata
581        PickId = G2frame.PickId
582        if G2frame.PatternTree.GetItemText(PickId) in ['Peak List','Limits'] and xpos:
583            lines = []
584            for line in G2frame.Lines: 
585                lines.append(line.get_xdata()[0])
586#            print G2frame.itemPicked.get_xdata()
587            lineNo = lines.index(G2frame.itemPicked.get_xdata()[0])
588            if  lineNo in [0,1]:
589                LimitId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Limits')
590                data = G2frame.PatternTree.GetItemPyData(LimitId)
591                if G2frame.qPlot:
592                    if 'C' in Parms['Type'][0]:
593                        data[1][lineNo] = 2.0*asind(wave*xpos/(4*math.pi))
594                    else:
595                        data[1][lineNo] = 2*math.pi*Parms['difC'][1]/xpos
596                else:
597                    data[1][lineNo] = xpos
598                G2frame.PatternTree.SetItemPyData(LimitId,data)
599                if G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Limits':
600                    G2pdG.UpdateLimitsGrid(G2frame,data)
601            else:
602                PeakId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Peak List')
603                data = G2frame.PatternTree.GetItemPyData(PeakId)
604#                print 'peaks',xpos
605                if event.button == 3:
606                    del data[lineNo-2]
607                else:
608                    if G2frame.qPlot:
609                        data[lineNo-2][0] = 2.0*asind(wave*xpos/(4*math.pi))
610                    else:
611                        data[lineNo-2][0] = xpos
612                G2frame.PatternTree.SetItemPyData(PeakId,data)
613                G2pdG.UpdatePeakGrid(G2frame,data)
614        elif (G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \
615            'PWDR' in G2frame.PatternTree.GetItemText(PickId)) and xpos:
616            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
617            pick = str(G2frame.itemPicked).split('(')[1].strip(')')
618            if 'line' not in pick:       #avoid data points, etc.
619                num = Phases.keys().index(pick)
620                if num:
621                    G2frame.refDelt = -(event.ydata-G2frame.refOffset)/(num*Ymax)
622                else:       #1st row of refl ticks
623                    G2frame.refOffset = event.ydata
624        PlotPatterns(G2frame)
625        G2frame.itemPicked = None   
626
627    xylim = []
628    try:
629        plotNum = G2frame.G2plotNB.plotList.index('Powder Patterns')
630        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
631        if not newPlot:
632            Plot = Page.figure.gca()          #get previous powder plot & get limits
633            xylim = Plot.get_xlim(),Plot.get_ylim()
634        Page.figure.clf()
635        Plot = Page.figure.gca()          #get a fresh plot after clf()
636    except ValueError:
637        newPlot = True
638        G2frame.Cmax = 1.0
639        Plot = G2frame.G2plotNB.addMpl('Powder Patterns').gca()
640        plotNum = G2frame.G2plotNB.plotList.index('Powder Patterns')
641        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
642        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
643        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
644        Page.canvas.mpl_connect('pick_event', OnPick)
645        Page.canvas.mpl_connect('button_release_event', OnRelease)
646    Page.SetFocus()
647    G2frame.G2plotNB.status.DestroyChildren()
648    if G2frame.Contour:
649        Page.Choice = (' key press','d: lower contour max','u: raise contour max',
650            'i: interpolation method','s: color scheme','c: contour off')
651    else:
652        if G2frame.logPlot:
653            Page.Choice = (' key press','n: log(I) off','l: offset left','r: offset right',
654                'c: contour on','q: toggle q plot','s: toggle single plot','+: no selection')
655        else:
656            Page.Choice = (' key press','l: offset left','r: offset right','d: offset down',
657                'u: offset up','o: reset offset','b: toggle subtract background','n: log(I) on','c: contour on',
658                'q: toggle q plot','s: toggle single plot','w: toggle divide by sig','+: no selection')
659    Page.keyPress = OnPlotKeyPress
660   
661    PickId = G2frame.PickId
662    PatternId = G2frame.PatternId
663    colors=['b','g','r','c','m','k']
664    Lines = []
665    if G2frame.SinglePlot:
666        Pattern = G2frame.PatternTree.GetItemPyData(PatternId)
667        Pattern.append(G2frame.PatternTree.GetItemText(PatternId))
668        PlotList = [Pattern,]
669        Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
670            G2frame.PatternId, 'Instrument Parameters'))
671        ParmList = [Parms,]
672        Title = Pattern[-1]
673    else:       
674        Title = os.path.split(G2frame.GSASprojectfile)[1]
675        PlotList = []
676        ParmList = []
677        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
678        while item:
679            if 'PWDR' in G2frame.PatternTree.GetItemText(item):
680                Pattern = G2frame.PatternTree.GetItemPyData(item)
681                if len(Pattern) < 3:                    # put name on end if needed
682                    Pattern.append(G2frame.PatternTree.GetItemText(item))
683                PlotList.append(Pattern)
684                ParmList.append(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
685                    item,'Instrument Parameters'))[0])
686            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)               
687    Ymax = 1.0
688    lenX = 0
689    if PickId:
690        if G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists']:
691            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
692            HKL = []
693            if Phases:
694                for peak in Phases[G2frame.RefList]:
695                    HKL.append(peak[:6])
696                HKL = np.array(HKL)
697        else:
698            HKL = np.array(G2frame.HKL)
699    for Pattern in PlotList:
700        xye = Pattern[1]
701        Ymax = max(Ymax,max(xye[1]))
702    offset = G2frame.Offset[0]*Ymax/100.0
703    if G2frame.logPlot:
704        Title = 'log('+Title+')'
705    Plot.set_title(Title)
706    if G2frame.qPlot:
707        Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
708    else:
709        if 'C' in ParmList[0]['Type'][0]:       
710            Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
711        else:
712            Plot.set_xlabel(r'TOF, $\mathsf{\mu}$s',fontsize=14)           
713    if G2frame.Weight:
714        Plot.set_ylabel(r'$\mathsf{I/\sigma(I)}$',fontsize=14)
715    else:
716        if 'C' in ParmList[0]['Type'][0]:
717            Plot.set_ylabel('Intensity',fontsize=14)
718        else:
719            Plot.set_ylabel('Normalized intensity',fontsize=14)
720    if G2frame.Contour:
721        ContourZ = []
722        ContourY = []
723        Nseq = 0
724    if len(PlotList) < 2:
725        G2frame.Contour = False
726    for N,Pattern in enumerate(PlotList):
727        Parms = ParmList[N]
728        if 'C' in Parms['Type'][0]:
729            wave = G2mth.getWave(Parms)
730        else:
731            difC = Parms['difC'][1]
732        ifpicked = False
733        LimitId = 0
734        xye = np.array(Pattern[1])
735        if PickId:
736            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
737            LimitId = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
738        if G2frame.qPlot:
739            Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, Pattern[2])
740            if 'C' in Parms['Type'][0]:
741                X = 4*np.pi*npsind((xye[0]-Parms['Zero'][1])/2.0)/wave
742            else:
743                X = 2*np.pi*Parms['difC'][1]/(xye[0]-Parms['Zero'][1])
744        else:
745            X = xye[0]-Parms['Zero'][1]
746        if not lenX:
747            lenX = len(X)           
748        Y = xye[1]+offset*N
749        if LimitId:
750            limits = np.array(G2frame.PatternTree.GetItemPyData(LimitId))
751            if G2frame.qPlot:
752                if 'C' in Parms['Type'][0]:
753                    limits = 4*np.pi*npsind(limits/2.0)/wave
754                else:
755                    limits = 2*np.pi*difC/limits
756            Lines.append(Plot.axvline(limits[1][0],color='g',dashes=(5,5),picker=3.))   
757            Lines.append(Plot.axvline(limits[1][1],color='r',dashes=(5,5),picker=3.))                   
758        if G2frame.Contour:
759           
760            if lenX == len(X):
761                ContourY.append(N)
762                ContourZ.append(Y)
763                ContourX = X
764                Nseq += 1
765                Plot.set_ylabel('Data sequence',fontsize=12)
766        else:
767            X += G2frame.Offset[1]*.005*N
768            if ifpicked:
769                Z = xye[3]+offset*N
770                W = xye[4]+offset*N
771                D = xye[5]-Ymax*G2frame.delOffset
772                if G2frame.logPlot:
773                    Plot.semilogy(X,Y,colors[N%6]+'+',picker=3.,clip_on=False,nonposy='mask')
774                    Plot.semilogy(X,Z,colors[(N+1)%6],picker=False,nonposy='mask')
775                    Plot.semilogy(X,W,colors[(N+2)%6],picker=False,nonposy='mask')
776                elif G2frame.Weight:
777                    DY = xye[1]*np.sqrt(xye[2])
778                    DYmax = max(DY)
779                    DZ = xye[3]*np.sqrt(xye[2])
780                    DS = xye[5]*np.sqrt(xye[2])-DYmax*G2frame.delOffset
781                    Plot.plot(X,DY,colors[N%6]+'+',picker=3.,clip_on=False)
782                    Plot.plot(X,DZ,colors[(N+1)%6],picker=False)
783                    Plot.plot(X,DS,colors[(N+3)%6],picker=False)
784                    Plot.axhline(0.,color=wx.BLACK)
785                else:
786                    if G2frame.SubBack:
787                        Plot.plot(X,Y-W,colors[N%6]+'+',picker=3.,clip_on=False)
788                        Plot.plot(X,Z-W,colors[(N+1)%6],picker=False)
789                    else:
790                        Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
791                        Plot.plot(X,Z,colors[(N+1)%6],picker=False)
792                    Plot.plot(X,W,colors[(N+2)%6],picker=False)
793                    Plot.plot(X,D,colors[(N+3)%6],picker=False)
794                    Plot.axhline(0.,color=wx.BLACK)
795                Page.canvas.SetToolTipString('')
796                if G2frame.PatternTree.GetItemText(PickId) == 'Peak List':
797                    tip = 'On data point: Pick peak - L or R MB. On line: L-move, R-delete'
798                    Page.canvas.SetToolTipString(tip)
799                    data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List'))
800                    for item in data:
801                        if G2frame.qPlot:
802                            if 'C' in Parms['Type'][0]:
803                                Lines.append(Plot.axvline(4*math.pi*sind(item[0]/2.)/wave,color=colors[N%6],picker=2.))
804                            else:
805                                Lines.append(Plot.axvline(2*math.pi*difC/item[0],color=colors[N%6],picker=2.))                               
806                        else:
807                            Lines.append(Plot.axvline(item[0],color=colors[N%6],picker=2.))
808                if G2frame.PatternTree.GetItemText(PickId) == 'Limits':
809                    tip = 'On data point: Lower limit - L MB; Upper limit - R MB. On limit: MB down to move'
810                    Page.canvas.SetToolTipString(tip)
811                    data = G2frame.LimitsTable.GetData()
812            else:
813                if G2frame.logPlot:
814                    Plot.semilogy(X,Y,colors[N%6],picker=False,nonposy='clip')
815                else:
816                    Plot.plot(X,Y,colors[N%6],picker=False)
817    if PickId and not G2frame.Contour:
818        Parms,Parms2 = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
819        if 'C' in Parms['Type'][0]:
820            wave = G2mth.getWave(Parms)
821        else:
822            difC = Parms['difC'][1]
823        if G2frame.PatternTree.GetItemText(PickId) in ['Index Peak List','Unit Cells List']:
824            peaks = np.array((G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))))
825            for peak in peaks:
826                if G2frame.qPlot:
827                    Plot.axvline(4*np.pi*sind(peak[0]/2.0)/wave,color='b')
828                else:
829                    Plot.axvline(peak[0],color='b')
830            for hkl in G2frame.HKL:
831                if G2frame.qPlot:
832                    Plot.axvline(4*np.pi*sind(hkl[5]/2.0)/wave,color='r',dashes=(5,5))
833                else:
834                    Plot.axvline(hkl[5],color='r',dashes=(5,5))
835        elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists'] or \
836            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
837            refColors=['b','r','c','g','m','k']
838            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
839            for pId,phase in enumerate(Phases):
840                peaks = Phases[phase]
841                if not peaks:
842                    continue
843                peak = np.array([[peak[4],peak[5]] for peak in peaks])
844                pos = G2frame.refOffset-pId*Ymax*G2frame.refDelt*np.ones_like(peak)
845                if G2frame.qPlot:
846                    Plot.plot(2*np.pi/peak.T[0],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
847                else:
848                    Plot.plot(peak.T[1],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
849            if len(Phases):
850                handles,legends = Plot.get_legend_handles_labels()  #got double entries in the legends for some reason
851                if handles:
852                    Plot.legend(handles[::2],legends[::2],title='Phases',loc='best')    #skip every other one
853           
854    if G2frame.Contour:
855        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
856        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
857            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
858        Page.figure.colorbar(Img)
859    else:
860        G2frame.Lines = Lines
861    if not newPlot:
862        Page.toolbar.push_current()
863        Plot.set_xlim(xylim[0])
864        Plot.set_ylim(xylim[1])
865        xylim = []
866        Page.toolbar.push_current()
867        Page.toolbar.draw()
868    else:
869        Page.canvas.draw()
870#    G2frame.Pwdr = True
871   
872################################################################################
873##### PlotDeltSig
874################################################################################
875           
876def PlotDeltSig(G2frame,kind):
877    'needs a doc string'
878    try:
879        plotNum = G2frame.G2plotNB.plotList.index('Error analysis')
880        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
881        Page.figure.clf()
882        Plot = Page.figure.gca()          #get a fresh plot after clf()
883    except ValueError:
884        newPlot = True
885        G2frame.Cmax = 1.0
886        Plot = G2frame.G2plotNB.addMpl('Error analysis').gca()
887        plotNum = G2frame.G2plotNB.plotList.index('Error analysis')
888        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
889    Page.Choice = None
890    PatternId = G2frame.PatternId
891    Pattern = G2frame.PatternTree.GetItemPyData(PatternId)
892    Pattern.append(G2frame.PatternTree.GetItemText(PatternId))
893    wtFactor = Pattern[0]['wtFactor']
894    if kind == 'PWDR':
895        limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits'))[1]
896        xye = np.array(Pattern[1])
897        xmin = np.searchsorted(xye[0],limits[0])
898        xmax = np.searchsorted(xye[0],limits[1])
899        DS = xye[5][xmin:xmax]*np.sqrt(wtFactor*xye[2][xmin:xmax])
900    elif kind == 'HKLF':
901        refl = Pattern[1]
902        DS = []
903        for ref in refl:
904            if ref[6] > 0.:
905                DS.append((ref[5]-ref[7])/ref[6])
906    Page.SetFocus()
907    G2frame.G2plotNB.status.DestroyChildren()
908    DS.sort()
909    EDS = np.zeros_like(DS)
910    DX = np.linspace(0.,1.,num=len(DS),endpoint=True)
911    oldErr = np.seterr(invalid='ignore')    #avoid problem at DX==0
912    T = np.sqrt(np.log(1.0/DX**2))
913    top = 2.515517+0.802853*T+0.010328*T**2
914    bot = 1.0+1.432788*T+0.189269*T**2+0.001308*T**3
915    EDS = np.where(DX>0,-(T-top/bot),(T-top/bot))
916    low1 = np.searchsorted(EDS,-1.)
917    hi1 = np.searchsorted(EDS,1.)
918    slp,intcp = np.polyfit(EDS[low1:hi1],DS[low1:hi1],deg=1)
919    frac = 100.*(hi1-low1)/len(DS)
920    G2frame.G2plotNB.status.SetStatusText(  \
921        'Over range -1. to 1. :'+' slope = %.3f, intercept = %.3f for %.2f%% of the fitted data'%(slp,intcp,frac),1)
922    Plot.set_title('Normal probability for '+Pattern[-1])
923    Plot.set_xlabel(r'expected $\mathsf{\Delta/\sigma}$',fontsize=14)
924    Plot.set_ylabel(r'observed $\mathsf{\Delta/\sigma}$',fontsize=14)
925    Plot.plot(EDS,DS,'r+',label='result')
926    Plot.plot([-2,2],[-2,2],'k',dashes=(5,5),label='ideal')
927    Plot.legend(loc='upper left')
928    np.seterr(invalid='warn')
929    Page.canvas.draw()
930       
931################################################################################
932##### PlotISFG
933################################################################################
934           
935def PlotISFG(G2frame,newPlot=False,type=''):
936    ''' PLotting package for PDF analysis; displays I(q), S(q), F(q) and G(r) as single
937    or multiple plots with waterfall and contour plots as options
938    '''
939    if not type:
940        type = G2frame.G2plotNB.plotList[G2frame.G2plotNB.nb.GetSelection()]
941    if type not in ['I(Q)','S(Q)','F(Q)','G(R)']:
942        return
943    superMinusOne = unichr(0xaf)+unichr(0xb9)
944   
945    def OnPlotKeyPress(event):
946        newPlot = False
947        if event.key == 'u':
948            if G2frame.Contour:
949                G2frame.Cmax = min(1.0,G2frame.Cmax*1.2)
950            elif G2frame.Offset[0] < 100.:
951                G2frame.Offset[0] += 1.
952        elif event.key == 'd':
953            if G2frame.Contour:
954                G2frame.Cmax = max(0.0,G2frame.Cmax*0.8)
955            elif G2frame.Offset[0] > 0.:
956                G2frame.Offset[0] -= 1.
957        elif event.key == 'l':
958            G2frame.Offset[1] -= 1.
959        elif event.key == 'r':
960            G2frame.Offset[1] += 1.
961        elif event.key == 'o':
962            G2frame.Offset = [0,0]
963        elif event.key == 'c':
964            newPlot = True
965            if G2frame.Contour:
966                G2frame.Contour = False
967            else:
968                G2frame.Contour = True
969                G2frame.SinglePlot = False
970                G2frame.Offset = [0.,0.]
971        elif event.key == 's':
972            if G2frame.Contour:
973                choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
974                choice.sort()
975                dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
976                if dlg.ShowModal() == wx.ID_OK:
977                    sel = dlg.GetSelection()
978                    G2frame.ContourColor = choice[sel]
979                else:
980                    G2frame.ContourColor = 'Paired'
981                dlg.Destroy()
982            else:               
983                if G2frame.SinglePlot:
984                    G2frame.SinglePlot = False
985                else:
986                    G2frame.SinglePlot = True
987        elif event.key == 'i':                  #for smoothing contour plot
988            choice = ['nearest','bilinear','bicubic','spline16','spline36','hanning',
989               'hamming','hermite','kaiser','quadric','catrom','gaussian','bessel',
990               'mitchell','sinc','lanczos']
991            dlg = wx.SingleChoiceDialog(G2frame,'Select','Interpolation',choice)
992            if dlg.ShowModal() == wx.ID_OK:
993                sel = dlg.GetSelection()
994                G2frame.Interpolate = choice[sel]
995            else:
996                G2frame.Interpolate = 'nearest'
997            dlg.Destroy()
998        elif event.key == 't' and not G2frame.Contour:
999            if G2frame.Legend:
1000                G2frame.Legend = False
1001            else:
1002                G2frame.Legend = True
1003           
1004           
1005        PlotISFG(G2frame,newPlot=newPlot,type=type)
1006       
1007    def OnKeyBox(event):
1008        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index(type):
1009            event.key = cb.GetValue()[0]
1010            cb.SetValue(' key press')
1011            wx.CallAfter(OnPlotKeyPress,event)
1012                       
1013    def OnMotion(event):
1014        xpos = event.xdata
1015        if xpos:                                        #avoid out of frame mouse position
1016            ypos = event.ydata
1017            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1018            try:
1019                if G2frame.Contour:
1020                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA pattern ID =%5d'%(xpos,int(ypos)),1)
1021                else:
1022                    G2frame.G2plotNB.status.SetStatusText('R =%.3fA %s =%.2f'%(xpos,type,ypos),1)                   
1023            except TypeError:
1024                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
1025   
1026    xylim = []
1027    try:
1028        plotNum = G2frame.G2plotNB.plotList.index(type)
1029        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1030        if not newPlot:
1031            Plot = Page.figure.gca()          #get previous plot & get limits
1032            xylim = Plot.get_xlim(),Plot.get_ylim()
1033        Page.figure.clf()
1034        Plot = Page.figure.gca()
1035    except ValueError:
1036        newPlot = True
1037        G2frame.Cmax = 1.0
1038        Plot = G2frame.G2plotNB.addMpl(type).gca()
1039        plotNum = G2frame.G2plotNB.plotList.index(type)
1040        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1041        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1042        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1043   
1044    Page.SetFocus()
1045    G2frame.G2plotNB.status.DestroyChildren()
1046    if G2frame.Contour:
1047        Page.Choice = (' key press','d: lower contour max','u: raise contour max',
1048            'i: interpolation method','s: color scheme','c: contour off')
1049    else:
1050        Page.Choice = (' key press','l: offset left','r: offset right','d: offset down','u: offset up',
1051            'o: reset offset','t: toggle legend','c: contour on','s: toggle single plot')
1052    Page.keyPress = OnPlotKeyPress
1053    PatternId = G2frame.PatternId
1054    PickId = G2frame.PickId
1055    Plot.set_title(type)
1056    if type == 'G(R)':
1057        Plot.set_xlabel(r'$R,\AA$',fontsize=14)
1058    else:
1059        Plot.set_xlabel(r'$Q,\AA$'+superMinusOne,fontsize=14)
1060    Plot.set_ylabel(r''+type,fontsize=14)
1061    colors=['b','g','r','c','m','k']
1062    name = G2frame.PatternTree.GetItemText(PatternId)[4:]
1063    Pattern = []   
1064    if G2frame.SinglePlot:
1065        name = G2frame.PatternTree.GetItemText(PatternId)
1066        name = type+name[4:]
1067        Id = G2gd.GetPatternTreeItemId(G2frame,PatternId,name)
1068        Pattern = G2frame.PatternTree.GetItemPyData(Id)
1069        if Pattern:
1070            Pattern.append(name)
1071        PlotList = [Pattern,]
1072    else:
1073        PlotList = []
1074        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
1075        while item:
1076            if 'PDF' in G2frame.PatternTree.GetItemText(item):
1077                name = type+G2frame.PatternTree.GetItemText(item)[4:]
1078                Id = G2gd.GetPatternTreeItemId(G2frame,item,name)
1079                Pattern = G2frame.PatternTree.GetItemPyData(Id)
1080                if Pattern:
1081                    Pattern.append(name)
1082                    PlotList.append(Pattern)
1083            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
1084    PDFdata = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'PDF Controls'))
1085    numbDen = G2pwd.GetNumDensity(PDFdata['ElList'],PDFdata['Form Vol'])
1086    Xb = [0.,10.]
1087    Yb = [0.,-40.*np.pi*numbDen]
1088    Ymax = 1.0
1089    lenX = 0
1090    for Pattern in PlotList:
1091        xye = Pattern[1]
1092        Ymax = max(Ymax,max(xye[1]))
1093    offset = G2frame.Offset[0]*Ymax/100.0
1094    if G2frame.Contour:
1095        ContourZ = []
1096        ContourY = []
1097        Nseq = 0
1098    for N,Pattern in enumerate(PlotList):
1099        xye = Pattern[1]
1100        if PickId:
1101            ifpicked = Pattern[2] == G2frame.PatternTree.GetItemText(PatternId)
1102        X = xye[0]
1103        if not lenX:
1104            lenX = len(X)           
1105        Y = xye[1]+offset*N
1106        if G2frame.Contour:
1107            if lenX == len(X):
1108                ContourY.append(N)
1109                ContourZ.append(Y)
1110                ContourX = X
1111                Nseq += 1
1112                Plot.set_ylabel('Data sequence',fontsize=12)
1113        else:
1114            X = xye[0]+G2frame.Offset[1]*.005*N
1115            if ifpicked:
1116                Plot.plot(X,Y,colors[N%6]+'+',picker=3.,clip_on=False)
1117                Page.canvas.SetToolTipString('')
1118            else:
1119                if G2frame.Legend:
1120                    Plot.plot(X,Y,colors[N%6],picker=False,label='Azm:'+Pattern[2].split('=')[1])
1121                else:
1122                    Plot.plot(X,Y,colors[N%6],picker=False)
1123            if type == 'G(R)':
1124                Plot.plot(Xb,Yb,color='k',dashes=(5,5))
1125            elif type == 'F(Q)':
1126                Plot.axhline(0.,color=wx.BLACK)
1127            elif type == 'S(Q)':
1128                Plot.axhline(1.,color=wx.BLACK)
1129    if G2frame.Contour:
1130        acolor = mpl.cm.get_cmap(G2frame.ContourColor)
1131        Img = Plot.imshow(ContourZ,cmap=acolor,vmin=0,vmax=Ymax*G2frame.Cmax,interpolation=G2frame.Interpolate, 
1132            extent=[ContourX[0],ContourX[-1],ContourY[0],ContourY[-1]],aspect='auto',origin='lower')
1133        Page.figure.colorbar(Img)
1134    elif G2frame.Legend:
1135        Plot.legend(loc='best')
1136    if not newPlot:
1137        Page.toolbar.push_current()
1138        Plot.set_xlim(xylim[0])
1139        Plot.set_ylim(xylim[1])
1140        xylim = []
1141        Page.toolbar.push_current()
1142        Page.toolbar.draw()
1143    else:
1144        Page.canvas.draw()
1145       
1146################################################################################
1147##### PlotXY
1148################################################################################
1149           
1150def PlotXY(G2frame,XY,newPlot=False,type=''):
1151    '''simple plot of xy data, used for diagnostic purposes
1152    '''
1153    def OnMotion(event):
1154        xpos = event.xdata
1155        if xpos:                                        #avoid out of frame mouse position
1156            ypos = event.ydata
1157            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1158            try:
1159                G2frame.G2plotNB.status.SetStatusText('X =%9.3f %s =%9.3f'%(xpos,type,ypos),1)                   
1160            except TypeError:
1161                G2frame.G2plotNB.status.SetStatusText('Select '+type+' pattern first',1)
1162
1163    try:
1164        plotNum = G2frame.G2plotNB.plotList.index(type)
1165        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1166        if not newPlot:
1167            Plot = Page.figure.gca()
1168            xylim = Plot.get_xlim(),Plot.get_ylim()
1169        Page.figure.clf()
1170        Plot = Page.figure.gca()
1171    except ValueError:
1172        newPlot = True
1173        Plot = G2frame.G2plotNB.addMpl(type).gca()
1174        plotNum = G2frame.G2plotNB.plotList.index(type)
1175        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1176        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1177   
1178    Page.Choice = None
1179    Page.SetFocus()
1180    G2frame.G2plotNB.status.DestroyChildren()
1181    Plot.set_title(type)
1182    Plot.set_xlabel(r'X',fontsize=14)
1183    Plot.set_ylabel(r''+type,fontsize=14)
1184    colors=['b','g','r','c','m','k']
1185    Ymax = 1.0
1186    lenX = 0
1187    X,Y = XY[:2]
1188    Ymax = max(Ymax,max(Y))
1189    Plot.plot(X,Y,'k',picker=False)
1190    if not newPlot:
1191        Page.toolbar.push_current()
1192        Plot.set_xlim(xylim[0])
1193        Plot.set_ylim(xylim[1])
1194        xylim = []
1195        Page.toolbar.push_current()
1196        Page.toolbar.draw()
1197    else:
1198        Page.canvas.draw()
1199
1200################################################################################
1201##### PlotPowderLines
1202################################################################################
1203           
1204def PlotPowderLines(G2frame):
1205    ''' plotting of powder lines (i.e. no powder pattern) as sticks
1206    '''
1207    global HKL
1208
1209    def OnMotion(event):
1210        xpos = event.xdata
1211        if xpos:                                        #avoid out of frame mouse position
1212            Page.canvas.SetCursor(wx.CROSS_CURSOR)
1213            G2frame.G2plotNB.status.SetFields(['','2-theta =%9.3f '%(xpos,)])
1214            if G2frame.PickId and G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List']:
1215                found = []
1216                if len(HKL):
1217                    view = Page.toolbar._views.forward()[0][:2]
1218                    wid = view[1]-view[0]
1219                    found = HKL[np.where(np.fabs(HKL.T[5]-xpos) < 0.002*wid)]
1220                if len(found):
1221                    h,k,l = found[0][:3] 
1222                    Page.canvas.SetToolTipString('%d,%d,%d'%(int(h),int(k),int(l)))
1223                else:
1224                    Page.canvas.SetToolTipString('')
1225
1226    try:
1227        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
1228        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1229        Page.figure.clf()
1230        Plot = Page.figure.gca()
1231    except ValueError:
1232        Plot = G2frame.G2plotNB.addMpl('Powder Lines').gca()
1233        plotNum = G2frame.G2plotNB.plotList.index('Powder Lines')
1234        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1235        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1236       
1237    Page.Choice = None
1238    Page.SetFocus()
1239    Plot.set_title('Powder Pattern Lines')
1240    Plot.set_xlabel(r'$\mathsf{2\theta}$',fontsize=14)
1241    PickId = G2frame.PickId
1242    PatternId = G2frame.PatternId
1243    peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))
1244    for peak in peaks:
1245        Plot.axvline(peak[0],color='b')
1246    HKL = np.array(G2frame.HKL)
1247    for hkl in G2frame.HKL:
1248        Plot.axvline(hkl[5],color='r',dashes=(5,5))
1249    xmin = peaks[0][0]
1250    xmax = peaks[-1][0]
1251    delt = xmax-xmin
1252    xlim = [max(0,xmin-delt/20.),min(180.,xmax+delt/20.)]
1253    Plot.set_xlim(xlim)
1254    Page.canvas.draw()
1255    Page.toolbar.push_current()
1256
1257################################################################################
1258##### PlotPeakWidths
1259################################################################################
1260           
1261def PlotPeakWidths(G2frame):
1262    ''' Plotting of instrument broadening terms as function of 2-theta
1263    Seen when "Instrument Parameters" chosen from powder pattern data tree
1264    '''
1265#    sig = lambda Th,U,V,W: 1.17741*math.sqrt(U*tand(Th)**2+V*tand(Th)+W)*math.pi/18000.
1266#    gam = lambda Th,X,Y: (X/cosd(Th)+Y*tand(Th))*math.pi/18000.
1267    gamFW = lambda s,g: np.exp(np.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.)
1268#    gamFW2 = lambda s,g: math.sqrt(s**2+(0.4654996*g)**2)+.5345004*g  #Ubaldo Bafile - private communication
1269    PatternId = G2frame.PatternId
1270    limitID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Limits')
1271    if limitID:
1272        limits = G2frame.PatternTree.GetItemPyData(limitID)
1273    else:
1274        return
1275    Parms,Parms2 = G2frame.PatternTree.GetItemPyData( \
1276        G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Instrument Parameters'))
1277    if 'C' in Parms['Type'][0]:
1278        lam = G2mth.getWave(Parms)
1279    else:
1280        difC = Parms['difC'][0]
1281    peakID = G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Peak List')
1282    if peakID:
1283        peaks = G2frame.PatternTree.GetItemPyData(peakID)
1284    else:
1285        peaks = []
1286   
1287    try:
1288        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
1289        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1290        Page.figure.clf()
1291        Plot = Page.figure.gca()
1292    except ValueError:
1293        Plot = G2frame.G2plotNB.addMpl('Peak Widths').gca()
1294        plotNum = G2frame.G2plotNB.plotList.index('Peak Widths')
1295        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1296    Page.Choice = None
1297    Page.SetFocus()
1298   
1299    Page.canvas.SetToolTipString('')
1300    colors=['b','g','r','c','m','k']
1301    X = []
1302    Y = []
1303    Z = []
1304    W = []
1305    if 'C' in Parms['Type'][0]:
1306        Plot.set_title('Instrument and sample peak widths')
1307        Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
1308        Plot.set_ylabel(r'$\Delta q/q, \Delta d/d$',fontsize=14)
1309        try:
1310            Xmin,Xmax = limits[1]
1311            X = np.linspace(Xmin,Xmax,num=101,endpoint=True)
1312            Q = 4.*np.pi*npsind(X/2.)/lam
1313            Z = np.ones_like(X)
1314            data = G2mth.setPeakparms(Parms,Parms2,X,Z)
1315            s = 1.17741*np.sqrt(data[4])*np.pi/18000.
1316            g = data[6]*np.pi/18000.
1317            G = gamFW(g,s)
1318            Y = s/nptand(X/2.)
1319            Z = g/nptand(X/2.)
1320            W = G/nptand(X/2.)
1321            Plot.plot(Q,Y,color='r',label='Gaussian')
1322            Plot.plot(Q,Z,color='g',label='Lorentzian')
1323            Plot.plot(Q,W,color='b',label='G+L')
1324           
1325            fit = G2mth.setPeakparms(Parms,Parms2,X,Z,useFit=True)
1326            sf = 1.17741*np.sqrt(fit[4])*np.pi/18000.
1327            gf = fit[6]*np.pi/18000.
1328            Gf = gamFW(gf,sf)
1329            Yf = sf/nptand(X/2.)
1330            Zf = gf/nptand(X/2.)
1331            Wf = Gf/nptand(X/2.)
1332            Plot.plot(Q,Yf,color='r',dashes=(5,5),label='Gaussian fit')
1333            Plot.plot(Q,Zf,color='g',dashes=(5,5),label='Lorentzian fit')
1334            Plot.plot(Q,Wf,color='b',dashes=(5,5),label='G+L fit')
1335           
1336            X = []
1337            Y = []
1338            Z = []
1339            W = []
1340            V = []
1341            for peak in peaks:
1342                X.append(4.0*math.pi*sind(peak[0]/2.0)/lam)
1343                try:
1344                    s = 1.17741*math.sqrt(peak[4])*math.pi/18000.
1345                except ValueError:
1346                    s = 0.01
1347                g = peak[6]*math.pi/18000.
1348                G = gamFW(g,s)
1349                Y.append(s/tand(peak[0]/2.))
1350                Z.append(g/tand(peak[0]/2.))
1351                W.append(G/tand(peak[0]/2.))
1352            if len(peaks):
1353                Plot.plot(X,Y,'+',color='r',label='G peak')
1354                Plot.plot(X,Z,'+',color='g',label='L peak')
1355                Plot.plot(X,W,'+',color='b',label='G+L peak')
1356            Plot.legend(loc='best')
1357            Page.canvas.draw()
1358        except ValueError:
1359            print '**** ERROR - default U,V,W profile coefficients yield sqrt of negative value at 2theta =', \
1360                '%.3f'%(2*theta)
1361            G2frame.G2plotNB.Delete('Peak Widths')
1362    else:
1363        Plot.set_title('Instrument and sample peak coefficients')
1364        Plot.set_xlabel(r'$q, \AA^{-1}$',fontsize=14)
1365        Plot.set_ylabel(r'$\alpha, \beta, \Delta q/q, \Delta d/d$',fontsize=14)
1366        Xmin,Xmax = limits[1]
1367        T = np.linspace(Xmin,Xmax,num=101,endpoint=True)
1368        Z = np.ones_like(T)
1369        data = G2mth.setPeakparms(Parms,Parms2,T,Z)
1370        ds = T/difC
1371        Q = 2.*np.pi/ds
1372        A = data[4]
1373        B = data[6]
1374        S = 1.17741*np.sqrt(data[8])/T
1375        G = data[10]/T
1376        Plot.plot(Q,A,color='r',label='Alpha')
1377        Plot.plot(Q,B,color='g',label='Beta')
1378        Plot.plot(Q,S,color='b',label='Gaussian')
1379        Plot.plot(Q,G,color='m',label='Lorentzian')
1380
1381        fit = G2mth.setPeakparms(Parms,Parms2,T,Z)
1382        ds = T/difC
1383        Q = 2.*np.pi/ds
1384        Af = fit[4]
1385        Bf = fit[6]
1386        Sf = 1.17741*np.sqrt(fit[8])/T
1387        Gf = fit[10]/T
1388        Plot.plot(Q,Af,color='r',dashes=(5,5),label='Alpha fit')
1389        Plot.plot(Q,Bf,color='g',dashes=(5,5),label='Beta fit')
1390        Plot.plot(Q,Sf,color='b',dashes=(5,5),label='Gaussian fit')
1391        Plot.plot(Q,Gf,color='m',dashes=(5,5),label='Lorentzian fit')
1392       
1393        T = []
1394        A = []
1395        B = []
1396        S = []
1397        G = []
1398        W = []
1399        Q = []
1400        V = []
1401        for peak in peaks:
1402            T.append(peak[0])
1403            A.append(peak[4])
1404            B.append(peak[6])
1405            Q.append(2.*np.pi*difC/peak[0])
1406            S.append(1.17741*np.sqrt(peak[8])/peak[0])
1407            G.append(peak[10]/peak[0])
1408           
1409       
1410        Plot.plot(Q,A,'+',color='r',label='Alpha peak')
1411        Plot.plot(Q,B,'+',color='g',label='Beta peak')
1412        Plot.plot(Q,S,'+',color='b',label='Gaussian peak')
1413        Plot.plot(Q,G,'+',color='m',label='Lorentzian peak')
1414        Plot.legend(loc='best')
1415        Page.canvas.draw()
1416
1417   
1418################################################################################
1419##### PlotSizeStrainPO
1420################################################################################
1421           
1422def PlotSizeStrainPO(G2frame,data,Start=False):
1423    '''Plot 3D mustrain/size/preferred orientation figure. In this instance data is for a phase
1424    '''
1425   
1426    PatternId = G2frame.PatternId
1427    generalData = data['General']
1428    SGData = generalData['SGData']
1429    SGLaue = SGData['SGLaue']
1430    if Start:                   #initialize the spherical harmonics qlmn arrays
1431        ptx.pyqlmninit()
1432        Start = False
1433#    MuStrKeys = G2spc.MustrainNames(SGData)
1434    cell = generalData['Cell'][1:]
1435    A,B = G2lat.cell2AB(cell[:6])
1436    Vol = cell[6]
1437    useList = data['Histograms']
1438    phase = generalData['Name']
1439    plotType = generalData['Data plot type']
1440    plotDict = {'Mustrain':'Mustrain','Size':'Size','Preferred orientation':'Pref.Ori.'}
1441    for ptype in plotDict:
1442        G2frame.G2plotNB.Delete(ptype)
1443    if plotType in ['None']:
1444        return       
1445
1446    for item in useList:
1447        if useList[item]['Show']:
1448            break
1449    else:
1450        return            #nothing to show!!
1451   
1452    numPlots = len(useList)
1453
1454    if plotType in ['Mustrain','Size']:
1455        Plot = mp3d.Axes3D(G2frame.G2plotNB.add3D(plotType))
1456    else:
1457        Plot = G2frame.G2plotNB.addMpl(plotType).gca()       
1458    plotNum = G2frame.G2plotNB.plotList.index(plotType)
1459    Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1460    Page.Choice = None
1461    Page.SetFocus()
1462    G2frame.G2plotNB.status.SetStatusText('',1)
1463    if not Page.IsShown():
1464        Page.Show()
1465   
1466    for item in useList:
1467        if useList[item]['Show']:
1468            PHI = np.linspace(0.,360.,30,True)
1469            PSI = np.linspace(0.,180.,30,True)
1470            X = np.outer(npsind(PHI),npsind(PSI))
1471            Y = np.outer(npcosd(PHI),npsind(PSI))
1472            Z = np.outer(np.ones(np.size(PHI)),npcosd(PSI))
1473            try:        #temp patch instead of 'mustrain' for old files with 'microstrain'
1474                coeff = useList[item][plotDict[plotType]]
1475            except KeyError:
1476                break
1477            if plotType in ['Mustrain','Size']:
1478                if coeff[0] == 'isotropic':
1479                    X *= coeff[1][0]
1480                    Y *= coeff[1][0]
1481                    Z *= coeff[1][0]                               
1482                elif coeff[0] == 'uniaxial':
1483                   
1484                    def uniaxCalc(xyz,iso,aniso,axes):
1485                        Z = np.array(axes)
1486                        cp = abs(np.dot(xyz,Z))
1487                        sp = np.sqrt(1.-cp**2)
1488                        R = iso*aniso/np.sqrt((iso*cp)**2+(aniso*sp)**2)
1489                        return R*xyz
1490                       
1491                    iso,aniso = coeff[1][:2]
1492                    axes = np.inner(A,np.array(coeff[3]))
1493                    axes /= nl.norm(axes)
1494                    Shkl = np.array(coeff[1])
1495                    XYZ = np.dstack((X,Y,Z))
1496                    XYZ = np.nan_to_num(np.apply_along_axis(uniaxCalc,2,XYZ,iso,aniso,axes))
1497                    X,Y,Z = np.dsplit(XYZ,3)
1498                    X = X[:,:,0]
1499                    Y = Y[:,:,0]
1500                    Z = Z[:,:,0]
1501               
1502                elif coeff[0] == 'ellipsoidal':
1503                   
1504                    def ellipseCalc(xyz,E,R):
1505                        XYZ = xyz*E.T
1506                        return np.inner(XYZ.T,R)
1507                       
1508                    S6 = coeff[4]
1509                    Sij = G2lat.U6toUij(S6)
1510                    E,R = nl.eigh(Sij)
1511                    XYZ = np.dstack((X,Y,Z))
1512                    XYZ = np.nan_to_num(np.apply_along_axis(ellipseCalc,2,XYZ,E,R))
1513                    X,Y,Z = np.dsplit(XYZ,3)
1514                    X = X[:,:,0]
1515                    Y = Y[:,:,0]
1516                    Z = Z[:,:,0]
1517                   
1518                elif coeff[0] == 'generalized':
1519                   
1520                    def genMustrain(xyz,SGData,A,Shkl):
1521                        uvw = np.inner(A.T,xyz)
1522                        Strm = np.array(G2spc.MustrainCoeff(uvw,SGData))
1523                        sum = np.sum(np.multiply(Shkl,Strm))
1524                        sum = np.where(sum > 0.01,sum,0.01)
1525                        sum = np.sqrt(sum)*math.pi/.18      #centidegrees to radians!
1526                        return sum*xyz
1527                       
1528                    Shkl = np.array(coeff[4])
1529                    if np.any(Shkl):
1530                        XYZ = np.dstack((X,Y,Z))
1531                        XYZ = np.nan_to_num(np.apply_along_axis(genMustrain,2,XYZ,SGData,A,Shkl))
1532                        X,Y,Z = np.dsplit(XYZ,3)
1533                        X = X[:,:,0]
1534                        Y = Y[:,:,0]
1535                        Z = Z[:,:,0]
1536                           
1537                if np.any(X) and np.any(Y) and np.any(Z):
1538                    errFlags = np.seterr(all='ignore')
1539                    Plot.plot_surface(X,Y,Z,rstride=1,cstride=1,color='g',linewidth=1)
1540                    np.seterr(all='ignore')
1541                    xyzlim = np.array([Plot.get_xlim3d(),Plot.get_ylim3d(),Plot.get_zlim3d()]).T
1542                    XYZlim = [min(xyzlim[0]),max(xyzlim[1])]
1543                    Plot.set_xlim3d(XYZlim)
1544                    Plot.set_ylim3d(XYZlim)
1545                    Plot.set_zlim3d(XYZlim)
1546                    Plot.set_aspect('equal')
1547                if plotType == 'Size':
1548                    Plot.set_title('Crystallite size for '+phase+'\n'+coeff[0]+' model')
1549                    Plot.set_xlabel(r'X, $\mu$m')
1550                    Plot.set_ylabel(r'Y, $\mu$m')
1551                    Plot.set_zlabel(r'Z, $\mu$m')
1552                else:   
1553                    Plot.set_title(r'$\mu$strain for '+phase+'\n'+coeff[0]+' model')
1554                    Plot.set_xlabel(r'X, $\mu$strain')
1555                    Plot.set_ylabel(r'Y, $\mu$strain')
1556                    Plot.set_zlabel(r'Z, $\mu$strain')
1557            else:
1558                h,k,l = generalData['POhkl']
1559                if coeff[0] == 'MD':
1560                    print 'March-Dollase preferred orientation plot'
1561               
1562                else:
1563                    PH = np.array(generalData['POhkl'])
1564                    phi,beta = G2lat.CrsAng(PH,cell[:6],SGData)
1565                    SHCoef = {}
1566                    for item in coeff[5]:
1567                        L,N = eval(item.strip('C'))
1568                        SHCoef['C%d,0,%d'%(L,N)] = coeff[5][item]                       
1569                    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1570                    X = np.linspace(0,90.0,26)
1571                    Y = G2lat.polfcal(ODFln,'0',X,0.0)
1572                    Plot.plot(X,Y,color='k',label=str(PH))
1573                    Plot.legend(loc='best')
1574                    Plot.set_title('Axial distribution for HKL='+str(PH)+' in '+phase)
1575                    Plot.set_xlabel(r'$\psi$',fontsize=16)
1576                    Plot.set_ylabel('MRD',fontsize=14)
1577    Page.canvas.draw()
1578   
1579################################################################################
1580##### PlotTexture
1581################################################################################
1582           
1583def PlotTexture(G2frame,data,Start=False):
1584    '''Pole figure, inverse pole figure, 3D pole distribution and 3D inverse pole distribution
1585    plotting.
1586    dict generalData contains all phase info needed which is in data
1587    '''
1588
1589    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1590    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1591    PatternId = G2frame.PatternId
1592    generalData = data['General']
1593    SGData = generalData['SGData']
1594    textureData = generalData['SH Texture']
1595    G2frame.G2plotNB.Delete('Texture')
1596    if not textureData['Order']:
1597        return                  #no plot!!
1598    SHData = generalData['SH Texture']
1599    SHCoef = SHData['SH Coeff'][1]
1600    cell = generalData['Cell'][1:7]
1601    Amat,Bmat = G2lat.cell2AB(cell)
1602    sq2 = 1.0/math.sqrt(2.0)
1603   
1604    def rp2xyz(r,p):
1605        z = npcosd(r)
1606        xy = np.sqrt(1.-z**2)
1607        return xy*npsind(p),xy*npcosd(p),z
1608           
1609    def OnMotion(event):
1610        SHData = data['General']['SH Texture']
1611        if event.xdata and event.ydata:                 #avoid out of frame errors
1612            xpos = event.xdata
1613            ypos = event.ydata
1614            if 'Inverse' in SHData['PlotType']:
1615                r = xpos**2+ypos**2
1616                if r <= 1.0:
1617                    if 'equal' in G2frame.Projection: 
1618                        r,p = 2.*npasind(np.sqrt(r)*sq2),npatan2d(ypos,xpos)
1619                    else:
1620                        r,p = 2.*npatand(np.sqrt(r)),npatan2d(ypos,xpos)
1621                    ipf = G2lat.invpolfcal(IODFln,SGData,np.array([r,]),np.array([p,]))
1622                    xyz = np.inner(Amat.T,np.array([rp2xyz(r,p)]))
1623                    y,x,z = list(xyz/np.max(np.abs(xyz)))
1624                   
1625                    G2frame.G2plotNB.status.SetFields(['',
1626                        'psi =%9.3f, beta =%9.3f, MRD =%9.3f hkl=%5.2f,%5.2f,%5.2f'%(r,p,ipf,x,y,z)])
1627                                   
1628            elif 'Axial' in SHData['PlotType']:
1629                pass
1630               
1631            else:                       #ordinary pole figure
1632                z = xpos**2+ypos**2
1633                if z <= 1.0:
1634                    z = np.sqrt(z)
1635                    if 'equal' in G2frame.Projection: 
1636                        r,p = 2.*npasind(z*sq2),npatan2d(ypos,xpos)
1637                    else:
1638                        r,p = 2.*npatand(z),npatan2d(ypos,xpos)
1639                    pf = G2lat.polfcal(ODFln,SamSym[textureData['Model']],np.array([r,]),np.array([p,]))
1640                    G2frame.G2plotNB.status.SetFields(['','phi =%9.3f, gam =%9.3f, MRD =%9.3f'%(r,p,pf)])
1641   
1642    try:
1643        plotNum = G2frame.G2plotNB.plotList.index('Texture')
1644        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1645        Page.figure.clf()
1646        Plot = Page.figure.gca()
1647        if not Page.IsShown():
1648            Page.Show()
1649    except ValueError:
1650        Plot = G2frame.G2plotNB.addMpl('Texture').gca()
1651        plotNum = G2frame.G2plotNB.plotList.index('Texture')
1652        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1653        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1654
1655    Page.Choice = None
1656    Page.SetFocus()
1657    G2frame.G2plotNB.status.SetFields(['',''])   
1658    PH = np.array(SHData['PFhkl'])
1659    phi,beta = G2lat.CrsAng(PH,cell,SGData)
1660    ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1661    PX = np.array(SHData['PFxyz'])
1662    gam = atan2d(PX[0],PX[1])
1663    xy = math.sqrt(PX[0]**2+PX[1]**2)
1664    xyz = math.sqrt(PX[0]**2+PX[1]**2+PX[2]**2)
1665    psi = asind(xy/xyz)
1666    IODFln = G2lat.Glnh(Start,SHCoef,psi,gam,SamSym[textureData['Model']])
1667    if 'Axial' in SHData['PlotType']:
1668        X = np.linspace(0,90.0,26)
1669        Y = G2lat.polfcal(ODFln,SamSym[textureData['Model']],X,0.0)
1670        Plot.plot(X,Y,color='k',label=str(SHData['PFhkl']))
1671        Plot.legend(loc='best')
1672        Plot.set_title('Axial distribution for HKL='+str(SHData['PFhkl']))
1673        Plot.set_xlabel(r'$\psi$',fontsize=16)
1674        Plot.set_ylabel('MRD',fontsize=14)
1675       
1676    else:       
1677        npts = 201
1678        if 'Inverse' in SHData['PlotType']:
1679            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
1680            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
1681            if 'equal' in G2frame.Projection:
1682                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
1683            else:
1684                R = np.where(R <= 1.,2.*npatand(R),0.0)
1685            Z = np.zeros_like(R)
1686            Z = G2lat.invpolfcal(IODFln,SGData,R,P)
1687            Z = np.reshape(Z,(npts,npts))
1688            CS = Plot.contour(Y,X,Z,aspect='equal')
1689            Plot.clabel(CS,fontsize=9,inline=1)
1690            try:
1691                Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
1692            except ValueError:
1693                pass
1694            Page.figure.colorbar(Img)
1695            Plot.set_title('Inverse pole figure for XYZ='+str(SHData['PFxyz']))
1696            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
1697                       
1698        else:
1699            X,Y = np.meshgrid(np.linspace(1.,-1.,npts),np.linspace(-1.,1.,npts))
1700            R,P = np.sqrt(X**2+Y**2).flatten(),npatan2d(X,Y).flatten()
1701            if 'equal' in G2frame.Projection:
1702                R = np.where(R <= 1.,2.*npasind(R*sq2),0.0)
1703            else:
1704                R = np.where(R <= 1.,2.*npatand(R),0.0)
1705            Z = np.zeros_like(R)
1706            Z = G2lat.polfcal(ODFln,SamSym[textureData['Model']],R,P)
1707            Z = np.reshape(Z,(npts,npts))
1708            try:
1709                CS = Plot.contour(Y,X,Z,aspect='equal')
1710                Plot.clabel(CS,fontsize=9,inline=1)
1711            except ValueError:
1712                pass
1713            Img = Plot.imshow(Z.T,aspect='equal',cmap=G2frame.ContourColor,extent=[-1,1,-1,1])
1714            Page.figure.colorbar(Img)
1715            Plot.set_title('Pole figure for HKL='+str(SHData['PFhkl']))
1716            Plot.set_xlabel(G2frame.Projection.capitalize()+' projection')
1717    Page.canvas.draw()
1718
1719################################################################################
1720##### PlotCovariance
1721################################################################################
1722           
1723def PlotCovariance(G2frame,Data):
1724    'needs a doc string'
1725    if not Data:
1726        print 'No covariance matrix available'
1727        return
1728    varyList = Data['varyList']
1729    values = Data['variables']
1730    Xmax = len(varyList)
1731    covMatrix = Data['covMatrix']
1732    sig = np.sqrt(np.diag(covMatrix))
1733    xvar = np.outer(sig,np.ones_like(sig))
1734    covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
1735    title = ' for\n'+Data['title']
1736    newAtomDict = Data['newAtomDict']
1737   
1738
1739    def OnPlotKeyPress(event):
1740        newPlot = False
1741        if event.key == 's':
1742            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
1743            choice.sort()
1744            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
1745            if dlg.ShowModal() == wx.ID_OK:
1746                sel = dlg.GetSelection()
1747                G2frame.VcovColor = choice[sel]
1748            else:
1749                G2frame.VcovColor = 'RdYlGn'
1750            dlg.Destroy()
1751        PlotCovariance(G2frame,Data)
1752
1753    def OnMotion(event):
1754        #there is a problem here - reports wrong values
1755        if event.button:
1756            ytics = imgAx.get_yticks()
1757            ytics = np.where(ytics<len(varyList),ytics,-1)
1758            ylabs = [np.where(0<=i ,varyList[int(i)],' ') for i in ytics]
1759            imgAx.set_yticklabels(ylabs)           
1760        if event.xdata and event.ydata:                 #avoid out of frame errors
1761            xpos = int(event.xdata+.5)
1762            ypos = int(event.ydata+.5)
1763            if -1 < xpos < len(varyList) and -1 < ypos < len(varyList):
1764                if xpos == ypos:
1765                    value = values[xpos]
1766                    name = varyList[xpos]
1767                    if varyList[xpos] in newAtomDict:
1768                        name,value = newAtomDict[name]                       
1769                    msg = '%s value = %.4g, esd = %.4g'%(name,value,sig[xpos])
1770                else:
1771                    msg = '%s - %s: %5.3f'%(varyList[xpos],varyList[ypos],covArray[xpos][ypos])
1772                Page.canvas.SetToolTipString(msg)
1773                G2frame.G2plotNB.status.SetFields(['',msg])
1774               
1775    try:
1776        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
1777        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1778        Page.figure.clf()
1779        Plot = Page.figure.gca()
1780        if not Page.IsShown():
1781            Page.Show()
1782    except ValueError:
1783        Plot = G2frame.G2plotNB.addMpl('Covariance').gca()
1784        plotNum = G2frame.G2plotNB.plotList.index('Covariance')
1785        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1786        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1787        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1788    Page.Choice = ['s: to change colors']
1789    Page.keyPress = OnPlotKeyPress
1790    Page.SetFocus()
1791    G2frame.G2plotNB.status.SetFields(['',''])   
1792    acolor = mpl.cm.get_cmap(G2frame.VcovColor)
1793    Img = Plot.imshow(covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower',
1794        vmin=-1.,vmax=1.)
1795    imgAx = Img.get_axes()
1796    ytics = imgAx.get_yticks()
1797    ylabs = [varyList[int(i)] for i in ytics[:-1]]
1798    imgAx.set_yticklabels(ylabs)
1799    colorBar = Page.figure.colorbar(Img)
1800    Plot.set_title('V-Cov matrix'+title)
1801    Plot.set_xlabel('Variable number')
1802    Plot.set_ylabel('Variable name')
1803    Page.canvas.draw()
1804   
1805################################################################################
1806##### PlotTorsion
1807################################################################################
1808
1809def PlotTorsion(G2frame,phaseName,Torsion,TorName,Names=[],Angles=[],Coeff=[]):
1810    'needs a doc string'
1811   
1812    global names
1813    names = Names
1814    sum = np.sum(Torsion)
1815    torsion = np.log(2*Torsion+1.)/sum
1816    tMin = np.min(torsion)
1817    tMax = np.max(torsion)
1818    torsion = 3.*(torsion-tMin)/(tMax-tMin)
1819    X = np.linspace(0.,360.,num=45)
1820   
1821    def OnMotion(event):
1822        if event.xdata and event.ydata:                 #avoid out of frame errors
1823            xpos = event.xdata
1824            ypos = event.ydata
1825            msg = 'torsion,energy: %5.3f %5.3f'%(xpos,ypos)
1826            Page.canvas.SetToolTipString(msg)
1827            G2frame.G2plotNB.status.SetFields(['',msg])
1828            if len(Angles):
1829                fX = np.where(np.fabs(Angles-xpos)<1.0)[0]
1830                if len(fX):
1831                    msg = 'atoms:'+names[fX[0]-1]
1832                    Page.canvas.SetToolTipString(msg)                   
1833
1834    try:
1835        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
1836        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1837        Page.figure.clf()
1838        Plot = Page.figure.gca()
1839        if not Page.IsShown():
1840            Page.Show()
1841    except ValueError:
1842        Plot = G2frame.G2plotNB.addMpl('Torsion').gca()
1843        plotNum = G2frame.G2plotNB.plotList.index('Torsion')
1844        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1845        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1846   
1847    Page.SetFocus()
1848    G2frame.G2plotNB.status.SetFields(['',''])
1849    Plot.plot(X,torsion,'b+')
1850    if len(Coeff):
1851        X2 = np.linspace(0.,360.,45)
1852        Y2 = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in X2])
1853        Plot.plot(X2,Y2,'r')
1854    if len(Angles):
1855        Eval = np.array([-G2mth.calcTorsionEnergy(x,Coeff)[1] for x in Angles])
1856        Plot.plot(Angles,Eval,'ro')
1857    Plot.set_xlim((0.,360.))
1858    Plot.set_title('Torsion angles for '+TorName+' in '+phaseName)
1859    Plot.set_xlabel('angle',fontsize=16)
1860    Plot.set_ylabel('Energy',fontsize=16)
1861    Page.canvas.draw()
1862   
1863################################################################################
1864##### PlotRama
1865################################################################################
1866
1867def PlotRama(G2frame,phaseName,Rama,RamaName,Names=[],PhiPsi=[],Coeff=[]):
1868    'needs a doc string'
1869
1870    global names
1871    names = Names
1872    rama = np.log(2*Rama+1.)
1873    ramaMax = np.max(rama)
1874    rama = np.reshape(rama,(45,45))
1875    global Phi,Psi
1876    Phi = []
1877    Psi = []
1878
1879    def OnPlotKeyPress(event):
1880        newPlot = False
1881        if event.key == 's':
1882            choice = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")]
1883            choice.sort()
1884            dlg = wx.SingleChoiceDialog(G2frame,'Select','Color scheme',choice)
1885            if dlg.ShowModal() == wx.ID_OK:
1886                sel = dlg.GetSelection()
1887                G2frame.RamaColor = choice[sel]
1888            else:
1889                G2frame.RamaColor = 'RdYlGn'
1890            dlg.Destroy()
1891        PlotRama(G2frame,phaseName,Rama)
1892
1893    def OnMotion(event):
1894        if event.xdata and event.ydata:                 #avoid out of frame errors
1895            xpos = event.xdata
1896            ypos = event.ydata
1897            msg = 'phi/psi: %5.3f %5.3f'%(xpos,ypos)
1898            Page.canvas.SetToolTipString(msg)
1899            G2frame.G2plotNB.status.SetFields(['',msg])
1900            if len(Phi):
1901                fPhi = np.where(np.fabs(Phi-xpos)<1.0)[0]
1902                fPsi = np.where(np.fabs(Psi-ypos)<1.0)[0]
1903                if len(fPhi) and len(fPsi) and fPhi[0] == fPsi[0]:
1904                    msg = 'atoms:'+names[fPhi[0]-1]
1905                    Page.canvas.SetToolTipString(msg)
1906           
1907    try:
1908        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
1909        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1910        Page.figure.clf()
1911        Plot = Page.figure.gca()
1912        if not Page.IsShown():
1913            Page.Show()
1914    except ValueError:
1915        Plot = G2frame.G2plotNB.addMpl('Ramachandran').gca()
1916        plotNum = G2frame.G2plotNB.plotList.index('Ramachandran')
1917        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1918        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
1919        Page.canvas.mpl_connect('key_press_event', OnPlotKeyPress)
1920
1921    Page.Choice = ['s: to change colors']
1922    Page.keyPress = OnPlotKeyPress
1923    Page.SetFocus()
1924    G2frame.G2plotNB.status.SetFields(['',''])   
1925    acolor = mpl.cm.get_cmap(G2frame.RamaColor)
1926    if RamaName == 'All' or '-1' in RamaName:
1927        if len(Coeff): 
1928            X,Y = np.meshgrid(np.linspace(-180.,180.,45),np.linspace(-180.,180.,45))
1929            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
1930            Plot.contour(X,Y,np.reshape(Z,(45,45)))
1931        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
1932            extent=[-180,180,-180,180],origin='lower')
1933        if len(PhiPsi):
1934            Phi,Psi = PhiPsi.T
1935            Phi = np.where(Phi>180.,Phi-360.,Phi)
1936            Psi = np.where(Psi>180.,Psi-360.,Psi)
1937            Plot.plot(Phi,Psi,'ro',picker=3)
1938        Plot.set_xlim((-180.,180.))
1939        Plot.set_ylim((-180.,180.))
1940    else:
1941        if len(Coeff): 
1942            X,Y = np.meshgrid(np.linspace(0.,360.,45),np.linspace(0.,360.,45))
1943            Z = np.array([-G2mth.calcRamaEnergy(x,y,Coeff)[1] for x,y in zip(X.flatten(),Y.flatten())])
1944            Plot.contour(X,Y,np.reshape(Z,(45,45)))
1945        Img = Plot.imshow(rama,aspect='equal',cmap=acolor,interpolation='nearest',
1946            extent=[0,360,0,360],origin='lower')
1947        if len(PhiPsi):
1948            Phi,Psi = PhiPsi.T
1949            Plot.plot(Phi,Psi,'ro',picker=3)
1950        Plot.set_xlim((0.,360.))
1951        Plot.set_ylim((0.,360.))
1952    Plot.set_title('Ramachandran for '+RamaName+' in '+phaseName)
1953    Plot.set_xlabel(r'$\phi$',fontsize=16)
1954    Plot.set_ylabel(r'$\psi$',fontsize=16)
1955    colorBar = Page.figure.colorbar(Img)
1956    Page.canvas.draw()
1957
1958
1959################################################################################
1960##### PlotSeq
1961################################################################################
1962           
1963def PlotSeq(G2frame,SeqData,SeqSig,SeqNames,sampleParm):
1964    'needs a doc string'
1965   
1966    def OnKeyPress(event):
1967        if event.key == 's' and sampleParm:
1968            if G2frame.xAxis:
1969                G2frame.xAxis = False
1970            else:
1971                G2frame.xAxis = True
1972            Draw(False)
1973    try:
1974        plotNum = G2frame.G2plotNB.plotList.index('Sequential refinement')
1975        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1976        Page.figure.clf()
1977        Plot = Page.figure.gca()
1978        if not Page.IsShown():
1979            Page.Show()
1980    except ValueError:
1981        Plot = G2frame.G2plotNB.addMpl('Sequential refinement').gca()
1982        plotNum = G2frame.G2plotNB.plotList.index('Sequential refinement')
1983        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
1984        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
1985        G2frame.xAxis = False
1986    Page.Choice = ['s to toggle x-axis = sample environment parameter']
1987    Page.keyPress = OnKeyPress
1988       
1989    def Draw(newPlot):
1990        Page.SetFocus()
1991        G2frame.G2plotNB.status.SetFields(['','press '])
1992        if len(SeqData):
1993            Plot.clear()
1994            if G2frame.xAxis:   
1995                xName = sampleParm.keys()[0]
1996                X = sampleParm[xName]
1997            else:
1998                X = np.arange(0,len(SeqData[0]),1)
1999                xName = 'Data sequence number'
2000            for Y,sig,name in zip(SeqData,SeqSig,SeqNames):
2001                Plot.errorbar(X,Y,yerr=sig,label=name)       
2002            Plot.legend(loc='best')
2003            Plot.set_ylabel('Parameter values')
2004            Plot.set_xlabel(xName)
2005            Page.canvas.draw()           
2006    Draw(True)
2007           
2008################################################################################
2009##### PlotExposedImage & PlotImage
2010################################################################################
2011           
2012def PlotExposedImage(G2frame,newPlot=False,event=None):
2013    '''General access module for 2D image plotting
2014    '''
2015    plotNo = G2frame.G2plotNB.nb.GetSelection()
2016    if G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Powder Image':
2017        PlotImage(G2frame,newPlot,event,newImage=True)
2018    elif G2frame.G2plotNB.nb.GetPageText(plotNo) == '2D Integration':
2019        PlotIntegration(G2frame,newPlot,event)
2020
2021def PlotImage(G2frame,newPlot=False,event=None,newImage=True):
2022    '''Plot of 2D detector images as contoured plot. Also plot calibration ellipses,
2023    masks, etc.
2024    '''
2025    from matplotlib.patches import Ellipse,Arc,Circle,Polygon
2026    import numpy.ma as ma
2027    Dsp = lambda tth,wave: wave/(2.*sind(tth/2.))
2028    global Data,Masks
2029    colors=['b','g','r','c','m','k']
2030    Data = G2frame.PatternTree.GetItemPyData(
2031        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2032    Masks = G2frame.PatternTree.GetItemPyData(
2033        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'))
2034    try:    #may be absent
2035        StrSta = G2frame.PatternTree.GetItemPyData(
2036            G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Stress/Strain'))
2037    except TypeError:   #is missing
2038        StrSta = {}
2039
2040    def OnImMotion(event):
2041        Page.canvas.SetToolTipString('')
2042        sizexy = Data['size']
2043        if event.xdata and event.ydata and len(G2frame.ImageZ):                 #avoid out of frame errors
2044            Page.canvas.SetCursor(wx.CROSS_CURSOR)
2045            item = G2frame.itemPicked
2046            pixelSize = Data['pixelSize']
2047            scalex = 1000./pixelSize[0]
2048            scaley = 1000./pixelSize[1]
2049            if item and G2frame.PatternTree.GetItemText(G2frame.PickId) == 'Image Controls':
2050                if 'Text' in str(item):
2051                    Page.canvas.SetToolTipString('%8.3f %8.3fmm'%(event.xdata,event.ydata))
2052                else:
2053                    xcent,ycent = Data['center']
2054                    xpos = event.xdata-xcent
2055                    ypos = event.ydata-ycent
2056                    tth,azm = G2img.GetTthAzm(event.xdata,event.ydata,Data)
2057                    if 'line3' in  str(item) or 'line4' in str(item) and not Data['fullIntegrate']:
2058                        Page.canvas.SetToolTipString('%6d deg'%(azm))
2059                    elif 'line1' in  str(item) or 'line2' in str(item):
2060                        Page.canvas.SetToolTipString('%8.3fdeg'%(tth))                           
2061            else:
2062                xpos = event.xdata
2063                ypos = event.ydata
2064                xpix = xpos*scalex
2065                ypix = ypos*scaley
2066                Int = 0
2067                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
2068                    Int = G2frame.ImageZ[ypix][xpix]
2069                tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
2070                Q = 2.*math.pi/dsp
2071                if G2frame.setPoly:
2072                    G2frame.G2plotNB.status.SetFields(['','Polygon mask pick - LB next point, RB close polygon'])
2073                else:
2074                    G2frame.G2plotNB.status.SetFields(\
2075                        ['','Detector 2-th =%9.3fdeg, dsp =%9.3fA, Q = %6.5fA-1, azm = %7.2fdeg, I = %6d'%(tth,dsp,Q,azm,Int)])
2076
2077    def OnImPlotKeyPress(event):
2078        try:
2079            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
2080        except TypeError:
2081            return
2082        if PickName == 'Masks':
2083            Xpos = event.xdata
2084            if not Xpos:            #got point out of frame
2085                return
2086            Ypos = event.ydata
2087            if event.key == 's':
2088                Masks['Points'].append([Xpos,Ypos,1.])
2089            elif event.key == 'r':
2090                tth = G2img.GetTth(Xpos,Ypos,Data)
2091                Masks['Rings'].append([tth,0.1])
2092            elif event.key == 'a':
2093                tth,azm = G2img.GetTthAzm(Xpos,Ypos,Data)
2094                azm = int(azm)               
2095                Masks['Arcs'].append([tth,[azm-5,azm+5],0.1])
2096            elif event.key == 'p':
2097                G2frame.setPoly = True
2098                Masks['Polygons'].append([])
2099                G2frame.G2plotNB.status.SetFields(['','Polygon mask active - LB pick next point, RB close polygon'])
2100            G2imG.UpdateMasks(G2frame,Masks)
2101        elif PickName == 'Image Controls':
2102            if event.key == 'c':
2103                Xpos = event.xdata
2104                if not Xpos:            #got point out of frame
2105                    return
2106                Ypos = event.ydata
2107                dlg = wx.MessageDialog(G2frame,'Are you sure you want to change the center?',
2108                    'Center change',style=wx.OK|wx.CANCEL)
2109                try:
2110                    if dlg.ShowModal() == wx.ID_OK:
2111                        print 'move center to: ',Xpos,Ypos
2112                        Data['center'] = [Xpos,Ypos]
2113                        G2imG.UpdateImageControls(G2frame,Data,Masks)
2114                finally:
2115                    dlg.Destroy()
2116            elif event.key == 'l':
2117                if G2frame.logPlot:
2118                    G2frame.logPlot = False
2119                else:
2120                    G2frame.logPlot = True
2121        PlotImage(G2frame,newImage=True)
2122           
2123    def OnKeyBox(event):
2124        if G2frame.G2plotNB.nb.GetSelection() == G2frame.G2plotNB.plotList.index('2D Powder Image'):
2125            event.key = cb.GetValue()[0]
2126            cb.SetValue(' key press')
2127            if event.key in 'l':
2128                wx.CallAfter(OnImPlotKeyPress,event)
2129                       
2130    def OnImPick(event):
2131        if G2frame.PatternTree.GetItemText(G2frame.PickId) not in ['Image Controls','Masks']:
2132            return
2133        if G2frame.setPoly:
2134            polygon = Masks['Polygons'][-1]
2135            xpos,ypos = event.mouseevent.xdata,event.mouseevent.ydata
2136            if xpos and ypos:                       #point inside image
2137                if len(polygon) > 2 and event.mouseevent.button == 3:
2138                    x0,y0 = polygon[0]
2139                    polygon.append([x0,y0])
2140                    G2frame.setPoly = False
2141                    G2frame.G2plotNB.status.SetFields(['','Polygon closed - RB drag a vertex to change shape'])
2142                else:
2143                    G2frame.G2plotNB.status.SetFields(['','New polygon point: %.1f,%.1f'%(xpos,ypos)])
2144                    polygon.append([xpos,ypos])
2145                G2imG.UpdateMasks(G2frame,Masks)
2146        else:
2147            if G2frame.itemPicked is not None: return
2148            G2frame.itemPicked = event.artist
2149            G2frame.mousePicked = event.mouseevent
2150       
2151    def OnImRelease(event):
2152        try:
2153            PickName = G2frame.PatternTree.GetItemText(G2frame.PickId)
2154        except TypeError:
2155            return
2156        if PickName not in ['Image Controls','Masks']:
2157            return
2158        pixelSize = Data['pixelSize']
2159        scalex = 1000./pixelSize[0]
2160        scaley = 1000./pixelSize[1]
2161        pixLimit = Data['pixLimit']
2162        if G2frame.itemPicked is None and PickName == 'Image Controls' and len(G2frame.ImageZ):
2163#            sizexy = Data['size']
2164            Xpos = event.xdata
2165            if not (Xpos and G2frame.ifGetRing):                   #got point out of frame
2166                return
2167            Ypos = event.ydata
2168            if Ypos and not Page.toolbar._active:         #make sure zoom/pan not selected
2169                if event.button == 1:
2170                    Xpix = Xpos*scalex
2171                    Ypix = Ypos*scaley
2172                    xpos,ypos,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
2173                    if I and J:
2174                        xpos += .5                              #shift to pixel center
2175                        ypos += .5
2176                        xpos /= scalex                          #convert to mm
2177                        ypos /= scaley
2178                        Data['ring'].append([xpos,ypos])
2179                elif event.button == 3:
2180                    G2frame.dataFrame.GetStatusBar().SetStatusText('Calibrating...')
2181                    if G2img.ImageCalibrate(G2frame,Data):
2182                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration successful - Show ring picks to check')
2183                        print 'Calibration successful'
2184                    else:
2185                        G2frame.dataFrame.GetStatusBar().SetStatusText('Calibration failed - Show ring picks to diagnose')
2186                        print 'Calibration failed'
2187                    G2frame.ifGetRing = False
2188                    G2imG.UpdateImageControls(G2frame,Data,Masks)
2189                    return
2190                PlotImage(G2frame,newImage=False)
2191            return
2192        else:
2193            xpos = event.xdata
2194            if xpos:                                        #avoid out of frame mouse position
2195                ypos = event.ydata
2196                if G2frame.ifGetRing:                          #delete a calibration ring pick
2197                    xypos = [xpos,ypos]
2198                    rings = Data['ring']
2199                    for ring in rings:
2200                        if np.allclose(ring,xypos,.01,0):
2201                            rings.remove(ring)                                                                       
2202                else:
2203                    tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
2204                    itemPicked = str(G2frame.itemPicked)
2205                    if 'Line2D' in itemPicked and PickName == 'Image Controls':
2206                        if 'line1' in itemPicked:
2207                            Data['IOtth'][0] = max(tth,0.001)
2208                        elif 'line2' in itemPicked:
2209                            Data['IOtth'][1] = tth
2210                        elif 'line3' in itemPicked:
2211                            Data['LRazimuth'][0] = int(azm)
2212                            if Data['fullIntegrate']:
2213                                Data['LRazimuth'][1] = Data['LRazimuth'][0]+360
2214                        elif 'line4' in itemPicked and not Data['fullIntegrate']:
2215                            Data['LRazimuth'][1] = int(azm)
2216                           
2217                        if Data['LRazimuth'][0] > Data['LRazimuth'][1]:
2218                            Data['LRazimuth'][0] -= 360
2219                           
2220                        azLim = np.array(Data['LRazimuth'])   
2221                        if np.any(azLim>360):
2222                            azLim -= 360
2223                            Data['LRazimuth'] = list(azLim)
2224                           
2225                        if  Data['IOtth'][0] > Data['IOtth'][1]:
2226                            Data['IOtth'][0],Data['IOtth'][1] = Data['IOtth'][1],Data['IOtth'][0]
2227                           
2228                        G2frame.InnerTth.SetValue("%8.2f" % (Data['IOtth'][0]))
2229                        G2frame.OuterTth.SetValue("%8.2f" % (Data['IOtth'][1]))
2230                        G2frame.Lazim.SetValue("%6d" % (Data['LRazimuth'][0]))
2231                        G2frame.Razim.SetValue("%6d" % (Data['LRazimuth'][1]))
2232                    elif 'Circle' in itemPicked and PickName == 'Masks':
2233                        spots = Masks['Points']
2234                        newPos = itemPicked.split(')')[0].split('(')[2].split(',')
2235                        newPos = np.array([float(newPos[0]),float(newPos[1])])
2236                        for spot in spots:
2237                            if np.allclose(np.array([spot[:2]]),newPos):
2238                                spot[:2] = xpos,ypos
2239                        G2imG.UpdateMasks(G2frame,Masks)
2240                    elif 'Line2D' in itemPicked and PickName == 'Masks':
2241                        Obj = G2frame.itemPicked.findobj()
2242                        rings = Masks['Rings']
2243                        arcs = Masks['Arcs']
2244                        polygons = Masks['Polygons']
2245                        for ring in G2frame.ringList:
2246                            if Obj == ring[0]:
2247                                rN = ring[1]
2248                                if ring[2] == 'o':
2249                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)-rings[rN][1]/2.
2250                                else:
2251                                    rings[rN][0] = G2img.GetTth(xpos,ypos,Data)+rings[rN][1]/2.
2252                        for arc in G2frame.arcList:
2253                            if Obj == arc[0]:
2254                                aN = arc[1]
2255                                if arc[2] == 'o':
2256                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)-arcs[aN][2]/2
2257                                elif arc[2] == 'i':
2258                                    arcs[aN][0] = G2img.GetTth(xpos,ypos,Data)+arcs[aN][2]/2
2259                                elif arc[2] == 'l':
2260                                    arcs[aN][1][0] = int(G2img.GetAzm(xpos,ypos,Data))
2261                                else:
2262                                    arcs[aN][1][1] = int(G2img.GetAzm(xpos,ypos,Data))
2263                        for poly in G2frame.polyList:
2264                            if Obj == poly[0]:
2265                                ind = G2frame.itemPicked.contains(G2frame.mousePicked)[1]['ind'][0]
2266                                oldPos = np.array([G2frame.mousePicked.xdata,G2frame.mousePicked.ydata])
2267                                pN = poly[1]
2268                                for i,xy in enumerate(polygons[pN]):
2269                                    if np.allclose(np.array([xy]),oldPos,atol=1.0):
2270                                        polygons[pN][i] = xpos,ypos
2271                        G2imG.UpdateMasks(G2frame,Masks)
2272#                    else:                  #keep for future debugging
2273#                        print str(G2frame.itemPicked),event.xdata,event.ydata,event.button
2274                PlotImage(G2frame,newImage=True)
2275            G2frame.itemPicked = None
2276           
2277    try:
2278        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
2279        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2280        if not newPlot:
2281            Plot = Page.figure.gca()          #get previous powder plot & get limits
2282            xylim = Plot.get_xlim(),Plot.get_ylim()
2283        if newImage:
2284            Page.figure.clf()
2285            Plot = Page.figure.gca()          #get a fresh plot after clf()
2286       
2287    except ValueError:
2288        Plot = G2frame.G2plotNB.addMpl('2D Powder Image').gca()
2289        plotNum = G2frame.G2plotNB.plotList.index('2D Powder Image')
2290        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2291        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
2292        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
2293        Page.canvas.mpl_connect('pick_event', OnImPick)
2294        Page.canvas.mpl_connect('button_release_event', OnImRelease)
2295        xylim = []
2296    Page.Choice = None
2297    if not event:                       #event from GUI TextCtrl - don't want focus to change to plot!!!
2298        Page.SetFocus()
2299    Title = G2frame.PatternTree.GetItemText(G2frame.Image)[4:]
2300    G2frame.G2plotNB.status.DestroyChildren()
2301    if G2frame.logPlot:
2302        Title = 'log('+Title+')'
2303    Plot.set_title(Title)
2304    try:
2305        if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Image Controls',]:
2306            Page.Choice = (' key press','l: log(I) on',)
2307            if G2frame.logPlot:
2308                Page.Choice[1] = 'l: log(I) off'
2309            Page.keyPress = OnImPlotKeyPress
2310    except TypeError:
2311        pass
2312    size,imagefile = G2frame.PatternTree.GetItemPyData(G2frame.Image)
2313    if imagefile != G2frame.oldImagefile:
2314        imagefile = G2IO.CheckImageFile(G2frame,imagefile)
2315        if not imagefile:
2316            G2frame.G2plotNB.Delete('2D Powder Image')
2317            return
2318        G2frame.PatternTree.SetItemPyData(G2frame.Image,[size,imagefile])
2319        G2frame.ImageZ = G2IO.GetImageData(G2frame,imagefile,imageOnly=True)
2320        G2frame.oldImagefile = imagefile
2321
2322    imScale = 1
2323    if len(G2frame.ImageZ) > 1024:
2324        imScale = len(G2frame.ImageZ)/1024
2325    sizexy = Data['size']
2326    pixelSize = Data['pixelSize']
2327    scalex = 1000./pixelSize[0]
2328    scaley = 1000./pixelSize[1]
2329    Xmax = sizexy[0]*pixelSize[0]/1000.
2330    Ymax = sizexy[1]*pixelSize[1]/1000.
2331    xlim = (0,Xmax)
2332    ylim = (Ymax,0)
2333    Imin,Imax = Data['range'][1]
2334    acolor = mpl.cm.get_cmap(Data['color'])
2335    xcent,ycent = Data['center']
2336    Plot.set_xlabel('Image x-axis, mm',fontsize=12)
2337    Plot.set_ylabel('Image y-axis, mm',fontsize=12)
2338    #do threshold mask - "real" mask - others are just bondaries
2339    Zlim = Masks['Thresholds'][1]
2340    wx.BeginBusyCursor()
2341    try:
2342           
2343        if newImage:                   
2344            MA = ma.masked_greater(ma.masked_less(G2frame.ImageZ,Zlim[0]),Zlim[1])
2345            MaskA = ma.getmaskarray(MA)
2346            A = G2img.ImageCompress(MA,imScale)
2347            AM = G2img.ImageCompress(MaskA,imScale)
2348            if G2frame.logPlot:
2349                A = np.where(A>0,np.log(A),0)
2350                AM = np.where(AM>0,np.log(AM),0)
2351                Imin,Imax = [np.amin(A),np.amax(A)]
2352            ImgM = Plot.imshow(AM,aspect='equal',cmap='Reds',
2353                interpolation='nearest',vmin=0,vmax=2,extent=[0,Xmax,Ymax,0])
2354            Img = Plot.imshow(A,aspect='equal',cmap=acolor,
2355                interpolation='nearest',vmin=Imin,vmax=Imax,extent=[0,Xmax,Ymax,0])
2356            if G2frame.setPoly:
2357                Img.set_picker(True)
2358   
2359        Plot.plot(xcent,ycent,'x')
2360        #G2frame.PatternTree.GetItemText(item)
2361        if Data['showLines']:
2362            LRAzim = Data['LRazimuth']                  #NB: integers
2363            Nazm = Data['outAzimuths']
2364            delAzm = float(LRAzim[1]-LRAzim[0])/Nazm
2365            AzmthOff = Data['azmthOff']
2366            IOtth = Data['IOtth']
2367            wave = Data['wavelength']
2368            dspI = wave/(2.0*sind(IOtth[0]/2.0))
2369            ellI = G2img.GetEllipse(dspI,Data)           #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)
2370            dspO = wave/(2.0*sind(IOtth[1]/2.0))
2371            ellO = G2img.GetEllipse(dspO,Data)           #Ditto & more likely for outer ellipse
2372            Azm = np.array(range(LRAzim[0],LRAzim[1]+1))-AzmthOff
2373            if ellI:
2374                xyI = []
2375                for azm in Azm:
2376                    xyI.append(G2img.GetDetectorXY(dspI,azm-90.,Data))
2377                xyI = np.array(xyI)
2378                arcxI,arcyI = xyI.T
2379                Plot.plot(arcxI,arcyI,picker=3)
2380            if ellO:
2381                xyO = []
2382                for azm in Azm:
2383                    xyO.append(G2img.GetDetectorXY(dspO,azm-90.,Data))
2384                xyO = np.array(xyO)
2385                arcxO,arcyO = xyO.T
2386                Plot.plot(arcxO,arcyO,picker=3)
2387            if ellO and ellI:
2388                Plot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3)
2389                Plot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3)
2390            for i in range(Nazm):
2391                cake = LRAzim[0]+i*delAzm-AzmthOff
2392                if Data['centerAzm']:
2393                    cake += delAzm/2.
2394                ind = np.searchsorted(Azm,cake)
2395                Plot.plot([arcxI[ind],arcxO[ind]],[arcyI[ind],arcyO[ind]],color='k',dashes=(5,5))
2396                   
2397        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Image Controls':
2398            for xring,yring in Data['ring']:
2399                Plot.plot(xring,yring,'r+',picker=3)
2400            if Data['setRings']:
2401    #            rings = np.concatenate((Data['rings']),axis=0)
2402                N = 0
2403                for ring in Data['rings']:
2404                    xring,yring = np.array(ring).T[:2]
2405                    Plot.plot(xring,yring,'+',color=colors[N%6])
2406                    N += 1           
2407            for ellipse in Data['ellipses']:
2408                cent,phi,[width,height],col = ellipse
2409                Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none'))
2410                Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
2411        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Stress/Strain':
2412            print 'plot stress/strain stuff'
2413            for ring in StrSta['d-zero']:
2414                xring,yring = ring['ImxyObs']
2415                Plot.plot(xring,yring,'r+')
2416#                for xring,yring in ring['ImxyCalc'].T:
2417#                    Plot.add_artist(Polygon(ring['ImxyCalc'].T,ec='b',fc='none'))
2418#                    Plot.plot(xring,yring)
2419        #masks - mask lines numbered after integration limit lines
2420        spots = Masks['Points']
2421        rings = Masks['Rings']
2422        arcs = Masks['Arcs']
2423        polygons = Masks['Polygons']
2424        for x,y,d in spots:
2425            Plot.add_artist(Circle((x,y),radius=d/2,fc='none',ec='r',picker=3))
2426        G2frame.ringList = []
2427        for iring,(tth,thick) in enumerate(rings):
2428            wave = Data['wavelength']
2429            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data))),2)
2430            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth-thick/2.,wave),Data))),2)
2431            G2frame.ringList.append([Plot.plot(x1,y1,'r',picker=3),iring,'o'])           
2432            G2frame.ringList.append([Plot.plot(x2,y2,'r',picker=3),iring,'i'])
2433        G2frame.arcList = []
2434        for iarc,(tth,azm,thick) in enumerate(arcs):           
2435            wave = Data['wavelength']
2436            x1,y1 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(tth+thick/2.,wave),Data),azm)),2)
2437            x2,y2 = np.hsplit(np.array(G2img.makeIdealRing(G2img.GetEllipse(Dsp(max(0.01,tth-thick/2.),wave),Data),azm)),2)
2438            G2frame.arcList.append([Plot.plot(x1,y1,'r',picker=3),iarc,'o'])           
2439            G2frame.arcList.append([Plot.plot(x2,y2,'r',picker=3),iarc,'i'])
2440            G2frame.arcList.append([Plot.plot([x1[0],x2[0]],[y1[0],y2[0]],'r',picker=3),iarc,'l'])
2441            G2frame.arcList.append([Plot.plot([x1[-1],x2[-1]],[y1[-1],y2[-1]],'r',picker=3),iarc,'u'])
2442        G2frame.polyList = []
2443        for ipoly,polygon in enumerate(polygons):
2444            x,y = np.hsplit(np.array(polygon),2)
2445            G2frame.polyList.append([Plot.plot(x,y,'r+',picker=10),ipoly])
2446            Plot.plot(x,y,'r')           
2447        if newImage:
2448            colorBar = Page.figure.colorbar(Img)
2449        Plot.set_xlim(xlim)
2450        Plot.set_ylim(ylim)
2451        Plot.invert_yaxis()
2452        if not newPlot and xylim:
2453            Page.toolbar.push_current()
2454            Plot.set_xlim(xylim[0])
2455            Plot.set_ylim(xylim[1])
2456            xylim = []
2457            Page.toolbar.push_current()
2458            Page.toolbar.draw()
2459        else:
2460            Page.canvas.draw()
2461    finally:
2462        wx.EndBusyCursor()
2463       
2464################################################################################
2465##### PlotIntegration
2466################################################################################
2467           
2468def PlotIntegration(G2frame,newPlot=False,event=None):
2469    '''Plot of 2D image after image integration with 2-theta and azimuth as coordinates
2470    '''
2471           
2472    def OnMotion(event):
2473        Page.canvas.SetToolTipString('')
2474        Page.canvas.SetCursor(wx.CROSS_CURSOR)
2475        azm = event.ydata
2476        tth = event.xdata
2477        if azm and tth:
2478            G2frame.G2plotNB.status.SetFields(\
2479                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
2480                               
2481    try:
2482        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
2483        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2484        if not newPlot:
2485            Plot = Page.figure.gca()          #get previous plot & get limits
2486            xylim = Plot.get_xlim(),Plot.get_ylim()
2487        Page.figure.clf()
2488        Plot = Page.figure.gca()          #get a fresh plot after clf()
2489       
2490    except ValueError:
2491        Plot = G2frame.G2plotNB.addMpl('2D Integration').gca()
2492        plotNum = G2frame.G2plotNB.plotList.index('2D Integration')
2493        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2494        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2495        Page.views = False
2496        view = False
2497    Page.Choice = None
2498    if not event:
2499        Page.SetFocus()
2500       
2501    Data = G2frame.PatternTree.GetItemPyData(
2502        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2503    image = G2frame.Integrate[0]
2504    xsc = G2frame.Integrate[1]
2505    ysc = G2frame.Integrate[2]
2506    Imin,Imax = Data['range'][1]
2507    acolor = mpl.cm.get_cmap(Data['color'])
2508    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
2509    Plot.set_ylabel('azimuth',fontsize=12)
2510    Plot.set_xlabel('2-theta',fontsize=12)
2511    Img = Plot.imshow(image,cmap=acolor,vmin=Imin,vmax=Imax,interpolation='nearest', \
2512        extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto')
2513    colorBar = Page.figure.colorbar(Img)
2514    if Data['ellipses']:           
2515        for ellipse in Data['ellipses']:
2516            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
2517            x,y = np.hsplit(ring,2)
2518            tth,azm = G2img.GetTthAzm(x,y,Data)
2519#            azm = np.where(azm < 0.,azm+360,azm)
2520            Plot.plot(tth,azm,'b,')
2521    if not newPlot:
2522        Page.toolbar.push_current()
2523        Plot.set_xlim(xylim[0])
2524        Plot.set_ylim(xylim[1])
2525        xylim = []
2526        Page.toolbar.push_current()
2527        Page.toolbar.draw()
2528    else:
2529        Page.canvas.draw()
2530               
2531################################################################################
2532##### PlotTRImage
2533################################################################################
2534           
2535def PlotTRImage(G2frame,tax,tay,taz,newPlot=False):
2536    '''a test plot routine - not normally used
2537    ''' 
2538           
2539    def OnMotion(event):
2540        Page.canvas.SetToolTipString('')
2541        Page.canvas.SetCursor(wx.CROSS_CURSOR)
2542        azm = event.xdata
2543        tth = event.ydata
2544        if azm and tth:
2545            G2frame.G2plotNB.status.SetFields(\
2546                ['','Detector 2-th =%9.3fdeg, azm = %7.2fdeg'%(tth,azm)])
2547                               
2548    try:
2549        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
2550        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2551        if not newPlot:
2552            Plot = Page.figure.gca()          #get previous plot & get limits
2553            xylim = Plot.get_xlim(),Plot.get_ylim()
2554        Page.figure.clf()
2555        Plot = Page.figure.gca()          #get a fresh plot after clf()
2556       
2557    except ValueError:
2558        Plot = G2frame.G2plotNB.addMpl('2D Transformed Powder Image').gca()
2559        plotNum = G2frame.G2plotNB.plotList.index('2D Transformed Powder Image')
2560        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
2561        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
2562        Page.views = False
2563        view = False
2564    Page.Choice = None
2565    Page.SetFocus()
2566       
2567    Data = G2frame.PatternTree.GetItemPyData(
2568        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
2569    Imin,Imax = Data['range'][1]
2570    step = (Imax-Imin)/5.
2571    V = np.arange(Imin,Imax,step)
2572    acolor = mpl.cm.get_cmap(Data['color'])
2573    Plot.set_title(G2frame.PatternTree.GetItemText(G2frame.Image)[4:])
2574    Plot.set_xlabel('azimuth',fontsize=12)
2575    Plot.set_ylabel('2-theta',fontsize=12)
2576    Plot.contour(tax,tay,taz,V,cmap=acolor)
2577    if Data['showLines']:
2578        IOtth = Data['IOtth']
2579        if Data['fullIntegrate']:
2580            LRAzim = [-180,180]
2581        else:
2582            LRAzim = Data['LRazimuth']                  #NB: integers
2583        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[0],IOtth[0]],picker=True)
2584        Plot.plot([LRAzim[0],LRAzim[1]],[IOtth[1],IOtth[1]],picker=True)
2585        if not Data['fullIntegrate']:
2586            Plot.plot([LRAzim[0],LRAzim[0]],[IOtth[0],IOtth[1]],picker=True)
2587            Plot.plot([LRAzim[1],LRAzim[1]],[IOtth[0],IOtth[1]],picker=True)
2588    if Data['setRings']:
2589        rings = np.concatenate((Data['rings']),axis=0)
2590        for xring,yring,dsp in rings:
2591            x,y = G2img.GetTthAzm(xring,yring,Data)
2592            Plot.plot(y,x,'r+')           
2593    if Data['ellipses']:           
2594        for ellipse in Data['ellipses']:
2595            ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color
2596            x,y = np.hsplit(ring,2)
2597            tth,azm = G2img.GetTthAzm(x,y,Data)
2598            Plot.plot(azm,tth,'b,')
2599    if not newPlot:
2600        Page.toolbar.push_current()
2601        Plot.set_xlim(xylim[0])
2602        Plot.set_ylim(xylim[1])
2603        xylim = []
2604        Page.toolbar.push_current()
2605        Page.toolbar.draw()
2606    else:
2607        Page.canvas.draw()
2608       
2609################################################################################
2610##### PlotStructure
2611################################################################################
2612           
2613def PlotStructure(G2frame,data):
2614    '''Crystal structure plotting package. Can show structures as balls, sticks, lines,
2615    thermal motion ellipsoids and polyhedra
2616    '''
2617
2618    def FindPeaksBonds(XYZ):
2619        rFact = drawingData['radiusFactor']
2620        Bonds = [[] for x in XYZ]
2621        for i,xyz in enumerate(XYZ):
2622            Dx = XYZ-xyz
2623            dist = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=1))
2624            IndB = ma.nonzero(ma.masked_greater(dist,rFact*2.2))
2625            for j in IndB[0]:
2626                Bonds[i].append(Dx[j]/2.)
2627                Bonds[j].append(-Dx[j]/2.)
2628        return Bonds
2629
2630    ForthirdPI = 4.0*math.pi/3.0
2631    generalData = data['General']
2632    cell = generalData['Cell'][1:7]
2633    Vol = generalData['Cell'][7:8][0]
2634    Amat,Bmat = G2lat.cell2AB(cell)         #Amat - crystal to cartesian, Bmat - inverse
2635    Gmat,gmat = G2lat.cell2Gmat(cell)
2636    A4mat = np.concatenate((np.concatenate((Amat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2637    B4mat = np.concatenate((np.concatenate((Bmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
2638    SGData = generalData['SGData']
2639    Mydir = generalData['Mydir']
2640    atomData = data['Atoms']
2641    mapPeaks = []
2642    drawingData = data['Drawing']   
2643    if 'Map Peaks' in data:
2644        mapPeaks = np.array(data['Map Peaks'])
2645        peakMax = 100.
2646        if len(mapPeaks):
2647            peakMax = np.max(mapPeaks.T[0])
2648    resRBData = data['RBModels'].get('Residue',[])
2649    vecRBData = data['RBModels'].get('Vector',[])
2650    rbAtmDict = {}
2651    for rbObj in resRBData+vecRBData:
2652        exclList = ['X' for i in range(len(rbObj['Ids']))]
2653        rbAtmDict.update(dict(zip(rbObj['Ids'],exclList)))
2654    testRBObj = data.get('testRBObj',{})
2655    rbObj = testRBObj.get('rbObj',{})
2656    MCSA = data.get('MCSA',{})
2657    mcsaModels = MCSA.get('Models',[])
2658    if mcsaModels:
2659            mcsaXYZ,atTypes = G2mth.UpdateMCSAxyz(Bmat,MCSA)
2660            XYZeq = []
2661            for xyz in mcsaXYZ:
2662                XYZeq += G2spc.GenAtom(xyz,SGData)[0][1:]       #skip self xyz
2663           
2664            mcsaBonds = FindPeaksBonds(mcsaXYZ)       
2665    drawAtoms = drawingData.get('Atoms',[])
2666    mapData = {}
2667    flipData = {}
2668    rhoXYZ = []
2669    showBonds = False
2670    if 'Map' in generalData:
2671        mapData = generalData['Map']
2672        showBonds = mapData.get('Show bonds',False)
2673    if 'Flip' in generalData:
2674        flipData = generalData['Flip']                       
2675        flipData['mapRoll'] = [0,0,0]
2676    Wt = np.array([255,255,255])
2677    Rd = np.array([255,0,0])
2678    Gr = np.array([0,255,0])
2679    Bl = np.array([0,0,255])
2680    Or = np.array([255,128,0])
2681    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]])
2682    uEdges = np.array([
2683        [uBox[0],uBox[1]],[uBox[0],uBox[3]],[uBox[0],uBox[4]],[uBox[1],uBox[2]], 
2684        [uBox[2],uBox[3]],[uBox[1],uBox[5]],[uBox[2],uBox[6]],[uBox[3],uBox[7]], 
2685        [uBox[4],uBox[5]],[uBox[5],uBox[6]],[uBox[6],uBox[7]],[uBox[7],uBox[4]]])
2686    mD = 0.1
2687    mV = np.array([[[-mD,0,0],[mD,0,0]],[[0,-mD,0],[0,mD,0]],[[0,0,-mD],[0,0,mD]]])
2688    mapPeakVecs = np.inner(mV,Bmat)
2689
2690    uColors = [Rd,Gr,Bl,Wt, Wt,Wt,Wt,Wt, Wt,Wt,Wt,Wt]
2691    altDown = False
2692    shiftDown = False
2693    ctrlDown = False
2694   
2695    def OnKeyBox(event):
2696        import Image
2697#        Draw()                          #make sure plot is fresh!!
2698        mode = cb.GetValue()
2699        if mode in ['jpeg','bmp','tiff',]:
2700            Fname = os.path.join(Mydir,generalData['Name']+'.'+mode)
2701            size = Page.canvas.GetSize()
2702            glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
2703            if mode in ['jpeg',]:
2704                Pix = glReadPixels(0,0,size[0],size[1],GL_RGBA, GL_UNSIGNED_BYTE)
2705                im = Image.new("RGBA", (size[0],size[1]))
2706            else:
2707                Pix = glReadPixels(0,0,size[0],size[1],GL_RGB, GL_UNSIGNED_BYTE)
2708                im = Image.new("RGB", (size[0],size[1]))
2709            im.fromstring(Pix)
2710            im.save(Fname,mode)
2711            cb.SetValue(' save as/key:')
2712            G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
2713        else:
2714            event.key = cb.GetValue()[0]
2715            cb.SetValue(' save as/key:')
2716            wx.CallAfter(OnKey,event)
2717
2718    def OnKey(event):           #on key UP!!
2719#        Draw()                          #make sure plot is fresh!!
2720        try:
2721            keyCode = event.GetKeyCode()
2722            if keyCode > 255:
2723                keyCode = 0
2724            key = chr(keyCode)
2725        except AttributeError:       #if from OnKeyBox above
2726            key = str(event.key).upper()
2727        indx = drawingData['selectedAtoms']
2728        cx,ct = drawingData['atomPtrs'][:2]
2729        if key in ['C']:
2730            drawingData['viewPoint'] = [[.5,.5,.5],[0,0]]
2731            drawingData['viewDir'] = [0,0,1]
2732            drawingData['oldxy'] = []
2733            V0 = np.array([0,0,1])
2734            V = np.inner(Amat,V0)
2735            V /= np.sqrt(np.sum(V**2))
2736            A = np.arccos(np.sum(V*V0))
2737            Q = G2mth.AV2Q(A,[0,1,0])
2738            drawingData['Quaternion'] = Q
2739            SetViewPointText(drawingData['viewPoint'][0])
2740            SetViewDirText(drawingData['viewDir'])
2741            Q = drawingData['Quaternion']
2742            G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
2743        elif key in ['N']:
2744            drawAtoms = drawingData['Atoms']
2745            pI = drawingData['viewPoint'][1]
2746            if indx:
2747                pI[0] = indx[pI[1]]
2748                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
2749                pI[1] += 1
2750                if pI[1] >= len(indx):
2751                    pI[1] = 0
2752            else:
2753                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
2754                pI[0] += 1
2755                if pI[0] >= len(drawAtoms):
2756                    pI[0] = 0
2757            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
2758            SetViewPointText(drawingData['viewPoint'][0])
2759            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
2760               
2761        elif key in ['P']:
2762            drawAtoms = drawingData['Atoms']
2763            pI = drawingData['viewPoint'][1]
2764            if indx:
2765                pI[0] = indx[pI[1]]
2766                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]
2767                pI[1] -= 1
2768                if pI[1] < 0:
2769                    pI[1] = len(indx)-1
2770            else:
2771                Tx,Ty,Tz = drawAtoms[pI[0]][cx:cx+3]               
2772                pI[0] -= 1
2773                if pI[0] < 0:
2774                    pI[0] = len(drawAtoms)-1
2775            drawingData['viewPoint'] = [[Tx,Ty,Tz],pI]
2776            SetViewPointText(drawingData['viewPoint'][0])           
2777            G2frame.G2plotNB.status.SetStatusText('View point at atom '+drawAtoms[pI[0]][ct-1]+str(pI),1)
2778        elif key in ['U','D','L','R'] and mapData['Flip'] == True:
2779            dirDict = {'U':[0,1],'D':[0,-1],'L':[-1,0],'R':[1,0]}
2780            SetMapRoll(dirDict[key])
2781            SetPeakRoll(dirDict[key])
2782            SetMapPeaksText(mapPeaks)
2783        Draw('key')
2784           
2785    def GetTruePosition(xy,Add=False):
2786        View = glGetIntegerv(GL_VIEWPORT)
2787        Proj = glGetDoublev(GL_PROJECTION_MATRIX)
2788        Model = glGetDoublev(GL_MODELVIEW_MATRIX)
2789        Zmax = 1.
2790        if Add:
2791            Indx = GetSelectedAtoms()
2792        if G2frame.dataDisplay.GetPageText(getSelection()) == 'Map peaks':
2793            for i,peak in enumerate(mapPeaks):
2794                x,y,z = peak[1:4]
2795                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
2796                XY = [int(X),int(View[3]-Y)]
2797                if np.allclose(xy,XY,atol=10) and Z < Zmax:
2798                    Zmax = Z
2799                    try:
2800                        Indx.remove(i)
2801                        ClearSelectedAtoms()
2802                        for id in Indx:
2803                            SetSelectedAtoms(id,Add)
2804                    except:
2805                        SetSelectedAtoms(i,Add)
2806        else:
2807            cx = drawingData['atomPtrs'][0]
2808            for i,atom in enumerate(drawAtoms):
2809                x,y,z = atom[cx:cx+3]
2810                X,Y,Z = gluProject(x,y,z,Model,Proj,View)
2811                XY = [int(X),int(View[3]-Y)]
2812                if np.allclose(xy,XY,atol=10) and Z < Zmax:
2813                    Zmax = Z
2814                    try:
2815                        Indx.remove(i)
2816                        ClearSelectedAtoms()
2817                        for id in Indx:
2818                            SetSelectedAtoms(id,Add)
2819                    except:
2820                        SetSelectedAtoms(i,Add)
2821                                       
2822    def OnMouseDown(event):
2823        xy = event.GetPosition()
2824        if event.ShiftDown():
2825            if event.LeftIsDown():
2826                GetTruePosition(xy)
2827            elif event.RightIsDown():
2828                GetTruePosition(xy,True)
2829        else:
2830            drawingData['oldxy'] = list(xy)
2831       
2832    def OnMouseMove(event):
2833        if event.ShiftDown():           #don't want any inadvertant moves when picking
2834            return
2835        newxy = event.GetPosition()
2836                               
2837        if event.Dragging():
2838            if event.AltDown() and rbObj:
2839                if event.LeftIsDown():
2840                    SetRBRotation(newxy)
2841                    Q = rbObj['Orient'][0]
2842                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
2843                elif event.RightIsDown():
2844                    SetRBTranslation(newxy)
2845                    Tx,Ty,Tz = rbObj['Orig'][0]
2846                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
2847                elif event.MiddleIsDown():
2848                    SetRBRotationZ(newxy)
2849                    Q = rbObj['Orient'][0]
2850                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
2851                Draw('move')
2852            elif not event.ControlDown():
2853                if event.LeftIsDown():
2854                    SetRotation(newxy)
2855                    Q = drawingData['Quaternion']
2856                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
2857                elif event.RightIsDown():
2858                    SetTranslation(newxy)
2859                    Tx,Ty,Tz = drawingData['viewPoint'][0]
2860                    G2frame.G2plotNB.status.SetStatusText('New view point: %.4f, %.4f, %.4f'%(Tx,Ty,Tz),1)
2861                elif event.MiddleIsDown():
2862                    SetRotationZ(newxy)
2863                    Q = drawingData['Quaternion']
2864                    G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
2865                Draw('move')
2866           
2867       
2868    def OnMouseWheel(event):
2869        if event.ShiftDown():
2870            return
2871        drawingData['cameraPos'] += event.GetWheelRotation()/24.
2872        drawingData['cameraPos'] = max(10,min(500,drawingData['cameraPos']))
2873        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(drawingData['cameraPos']),1)
2874        page = getSelection()
2875        if page:
2876            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2877                G2frame.dataDisplay.cameraPosTxt.SetLabel('Camera Position: '+'%.2f'%(drawingData['cameraPos']))
2878                G2frame.dataDisplay.cameraSlider.SetValue(drawingData['cameraPos'])
2879            Draw('wheel')
2880       
2881    def getSelection():
2882        try:
2883            return G2frame.dataDisplay.GetSelection()
2884        except AttributeError:
2885            G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!')
2886            return 0
2887           
2888    def SetViewPointText(VP):
2889        page = getSelection()
2890        if page:
2891            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2892                G2frame.dataDisplay.viewPoint.SetValue('%.3f %.3f %.3f'%(VP[0],VP[1],VP[2]))
2893               
2894    def SetRBOrigText():
2895        page = getSelection()
2896        if page:
2897            if G2frame.dataDisplay.GetPageText(page) == 'RB Models':
2898                for i,sizer in enumerate(testRBObj['Sizers']['Xsizers']):
2899                    sizer.SetValue('%8.5f'%(testRBObj['rbObj']['Orig'][0][i]))
2900                   
2901    def SetRBOrienText():
2902        page = getSelection()
2903        if page:
2904            if G2frame.dataDisplay.GetPageText(page) == 'RB Models':
2905                for i,sizer in enumerate(testRBObj['Sizers']['Osizers']):
2906                    sizer.SetValue('%8.5f'%(testRBObj['rbObj']['Orient'][0][i]))
2907               
2908    def SetViewDirText(VD):
2909        page = getSelection()
2910        if page:
2911            if G2frame.dataDisplay.GetPageText(page) == 'Draw Options':
2912                G2frame.dataDisplay.viewDir.SetValue('%.3f %.3f %.3f'%(VD[0],VD[1],VD[2]))
2913               
2914    def SetMapPeaksText(mapPeaks):
2915        page = getSelection()
2916        if page:
2917            if G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
2918                G2frame.MapPeaksTable.SetData(mapPeaks)
2919                panel = G2frame.dataDisplay.GetPage(page).GetChildren()
2920                names = [child.GetName() for child in panel]
2921                panel[names.index('grid window')].Refresh()
2922           
2923    def ClearSelectedAtoms():
2924        page = getSelection()
2925        if page:
2926            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2927                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Draw Atoms
2928            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
2929                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
2930            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2931                G2frame.dataDisplay.GetPage(page).ClearSelection()      #this is the Atoms grid in Atoms
2932               
2933                   
2934    def SetSelectedAtoms(ind,Add=False):
2935        page = getSelection()
2936        if page:
2937            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2938                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)      #this is the Atoms grid in Draw Atoms
2939            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
2940                G2frame.dataDisplay.GetPage(page).SelectRow(ind,Add)                 
2941            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2942                Id = drawAtoms[ind][-3]
2943                for i,atom in enumerate(atomData):
2944                    if atom[-1] == Id:
2945                        G2frame.dataDisplay.GetPage(page).SelectRow(i)      #this is the Atoms grid in Atoms
2946                 
2947    def GetSelectedAtoms():
2948        page = getSelection()
2949        Ind = []
2950        if page:
2951            if G2frame.dataDisplay.GetPageText(page) == 'Draw Atoms':
2952                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Draw Atoms
2953            elif G2frame.dataDisplay.GetPageText(page) == 'Map peaks':
2954                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()
2955            elif G2frame.dataDisplay.GetPageText(page) == 'Atoms':
2956                Ind = G2frame.dataDisplay.GetPage(page).GetSelectedRows()      #this is the Atoms grid in Atoms
2957        return Ind
2958                                       
2959    def SetBackground():
2960        R,G,B,A = Page.camera['backColor']
2961        glClearColor(R,G,B,A)
2962        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
2963       
2964    def SetLights():
2965        glEnable(GL_DEPTH_TEST)
2966        glShadeModel(GL_SMOOTH)
2967        glEnable(GL_LIGHTING)
2968        glEnable(GL_LIGHT0)
2969        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
2970        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
2971        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
2972       
2973    def GetRoll(newxy,rho):
2974        Q = drawingData['Quaternion']
2975        dxy = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,newxy+[0,]))
2976        dxy = np.array(dxy*rho.shape)       
2977        roll = np.where(dxy>0.5,1,np.where(dxy<-.5,-1,0))
2978        return roll
2979               
2980    def SetMapRoll(newxy):
2981        rho = mapData['rho']
2982        roll = GetRoll(newxy,rho)
2983        mapData['rho'] = np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2984        drawingData['oldxy'] = list(newxy)
2985       
2986    def SetPeakRoll(newxy):
2987        rho = mapData['rho']
2988        roll = GetRoll(newxy,rho)
2989        steps = 1./np.array(rho.shape)
2990        dxy = roll*steps
2991        for peak in mapPeaks:
2992            peak[1:4] += dxy
2993            peak[1:4] %= 1.
2994            peak[4] = np.sqrt(np.sum(np.inner(Amat,peak[1:4])**2))
2995               
2996    def SetTranslation(newxy):
2997#first get translation vector in screen coords.       
2998        oldxy = drawingData['oldxy']
2999        if not len(oldxy): oldxy = list(newxy)
3000        dxy = newxy-oldxy
3001        drawingData['oldxy'] = list(newxy)
3002        V = np.array([-dxy[0],dxy[1],0.])
3003#then transform to rotated crystal coordinates & apply to view point       
3004        Q = drawingData['Quaternion']
3005        V = np.inner(Bmat,G2mth.prodQVQ(G2mth.invQ(Q),V))
3006        Tx,Ty,Tz = drawingData['viewPoint'][0]
3007        Tx += V[0]*0.01
3008        Ty += V[1]*0.01
3009        Tz += V[2]*0.01
3010        drawingData['viewPoint'][0] =  Tx,Ty,Tz
3011        SetViewPointText([Tx,Ty,Tz])
3012       
3013    def SetRBTranslation(newxy):
3014#first get translation vector in screen coords.       
3015        oldxy = drawingData['oldxy']
3016        if not len(oldxy): oldxy = list(newxy)
3017        dxy = newxy-oldxy
3018        drawingData['oldxy'] = list(newxy)
3019        V = np.array([-dxy[0],dxy[1],0.])
3020#then transform to rotated crystal coordinates & apply to RB origin       
3021        Q = drawingData['Quaternion']
3022        V = np.inner(Bmat,G2mth.prodQVQ(G2mth.invQ(Q),V))
3023        Tx,Ty,Tz = rbObj['Orig'][0]
3024        Tx -= V[0]*0.01
3025        Ty -= V[1]*0.01
3026        Tz -= V[2]*0.01
3027        rbObj['Orig'][0] =  Tx,Ty,Tz
3028        SetRBOrigText()
3029       
3030    def SetRotation(newxy):
3031#first get rotation vector in screen coords. & angle increment       
3032        oldxy = drawingData['oldxy']
3033        if not len(oldxy): oldxy = list(newxy)
3034        dxy = newxy-oldxy
3035        drawingData['oldxy'] = list(newxy)
3036        V = np.array([dxy[1],dxy[0],0.])
3037        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
3038# next transform vector back to xtal coordinates via inverse quaternion
3039# & make new quaternion
3040        Q = drawingData['Quaternion']
3041        V = G2mth.prodQVQ(G2mth.invQ(Q),np.inner(Bmat,V))
3042        DQ = G2mth.AVdeg2Q(A,V)
3043        Q = G2mth.prodQQ(Q,DQ)
3044        drawingData['Quaternion'] = Q
3045# finally get new view vector - last row of rotation matrix
3046        VD = np.inner(Bmat,G2mth.Q2Mat(Q)[2])
3047        VD /= np.sqrt(np.sum(VD**2))
3048        drawingData['viewDir'] = VD
3049        SetViewDirText(VD)
3050       
3051    def SetRotationZ(newxy):                       
3052#first get rotation vector (= view vector) in screen coords. & angle increment       
3053        View = glGetIntegerv(GL_VIEWPORT)
3054        cent = [View[2]/2,View[3]/2]
3055        oldxy = drawingData['oldxy']
3056        if not len(oldxy): oldxy = list(newxy)
3057        dxy = newxy-oldxy
3058        drawingData['oldxy'] = list(newxy)
3059        V = drawingData['viewDir']
3060        A = [0,0]
3061        A[0] = dxy[1]*.25
3062        A[1] = dxy[0]*.25
3063        if newxy[0] > cent[0]:
3064            A[0] *= -1
3065        if newxy[1] < cent[1]:
3066            A[1] *= -1       
3067# next transform vector back to xtal coordinates & make new quaternion
3068        Q = drawingData['Quaternion']
3069        V = np.inner(Amat,V)
3070        Qx = G2mth.AVdeg2Q(A[0],V)
3071        Qy = G2mth.AVdeg2Q(A[1],V)
3072        Q = G2mth.prodQQ(Q,Qx)
3073        Q = G2mth.prodQQ(Q,Qy)
3074        drawingData['Quaternion'] = Q
3075
3076    def SetRBRotation(newxy):
3077#first get rotation vector in screen coords. & angle increment       
3078        oldxy = drawingData['oldxy']
3079        if not len(oldxy): oldxy = list(newxy)
3080        dxy = newxy-oldxy
3081        drawingData['oldxy'] = list(newxy)
3082        V = np.array([dxy[1],dxy[0],0.])
3083        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
3084# next transform vector back to xtal coordinates via inverse quaternion
3085# & make new quaternion
3086        Q = rbObj['Orient'][0]              #rotate RB to Cart
3087        QC = drawingData['Quaternion']      #rotate Cart to drawing
3088        V = G2mth.prodQVQ(G2mth.invQ(QC),V)
3089        V = G2mth.prodQVQ(G2mth.invQ(Q),V)
3090        DQ = G2mth.AVdeg2Q(A,V)
3091        Q = G2mth.prodQQ(Q,DQ)
3092        rbObj['Orient'][0] = Q
3093        SetRBOrienText()
3094       
3095    def SetRBRotationZ(newxy):                       
3096#first get rotation vector (= view vector) in screen coords. & angle increment       
3097        View = glGetIntegerv(GL_VIEWPORT)
3098        cent = [View[2]/2,View[3]/2]
3099        oldxy = drawingData['oldxy']
3100        if not len(oldxy): oldxy = list(newxy)
3101        dxy = newxy-oldxy
3102        drawingData['oldxy'] = list(newxy)
3103        V = drawingData['viewDir']
3104        A = [0,0]
3105        A[0] = dxy[1]*.25
3106        A[1] = dxy[0]*.25
3107        if newxy[0] < cent[0]:
3108            A[0] *= -1
3109        if newxy[1] > cent[1]:
3110            A[1] *= -1       
3111# next transform vector back to RB coordinates & make new quaternion
3112        Q = rbObj['Orient'][0]              #rotate RB to cart
3113        V = np.inner(Amat,V)
3114        V = -G2mth.prodQVQ(G2mth.invQ(Q),V)
3115        Qx = G2mth.AVdeg2Q(A[0],V)
3116        Qy = G2mth.AVdeg2Q(A[1],V)
3117        Q = G2mth.prodQQ(Q,Qx)
3118        Q = G2mth.prodQQ(Q,Qy)
3119        rbObj['Orient'][0] = Q
3120        SetRBOrienText()
3121
3122    def RenderBox():
3123        glEnable(GL_COLOR_MATERIAL)
3124        glLineWidth(2)
3125        glEnable(GL_BLEND)
3126        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
3127        glEnable(GL_LINE_SMOOTH)
3128        glBegin(GL_LINES)
3129        for line,color in zip(uEdges,uColors):
3130            glColor3ubv(color)
3131            glVertex3fv(line[0])
3132            glVertex3fv(line[1])
3133        glEnd()
3134        glColor4ubv([0,0,0,0])
3135        glDisable(GL_LINE_SMOOTH)
3136        glDisable(GL_BLEND)
3137        glDisable(GL_COLOR_MATERIAL)
3138       
3139    def RenderUnitVectors(x,y,z):
3140        xyz = np.array([x,y,z])
3141        glEnable(GL_COLOR_MATERIAL)
3142        glLineWidth(1)
3143        glPushMatrix()
3144        glTranslate(x,y,z)
3145        glScalef(1/cell[0],1/cell[1],1/cell[2])
3146        glBegin(GL_LINES)
3147        for line,color in zip(uEdges,uColors)[:3]:
3148            glColor3ubv(color)
3149            glVertex3fv(-line[1]/2.)
3150            glVertex3fv(line[1]/2.)
3151        glEnd()
3152        glPopMatrix()
3153        glColor4ubv([0,0,0,0])
3154        glDisable(GL_COLOR_MATERIAL)
3155               
3156    def RenderSphere(x,y,z,radius,color):
3157        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3158        glPushMatrix()
3159        glTranslate(x,y,z)
3160        glMultMatrixf(B4mat.T)
3161        q = gluNewQuadric()
3162        gluSphere(q,radius,20,10)
3163        glPopMatrix()
3164       
3165    def RenderDots(XYZ,RC):
3166        glEnable(GL_COLOR_MATERIAL)
3167        XYZ = np.array(XYZ)
3168        glPushMatrix()
3169        for xyz,rc in zip(XYZ,RC):
3170            x,y,z = xyz
3171            r,c = rc
3172            glColor3ubv(c)
3173            glPointSize(r*50)
3174            glBegin(GL_POINTS)
3175            glVertex3fv(xyz)
3176            glEnd()
3177        glPopMatrix()
3178        glColor4ubv([0,0,0,0])
3179        glDisable(GL_COLOR_MATERIAL)
3180       
3181    def RenderSmallSphere(x,y,z,radius,color):
3182        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3183        glPushMatrix()
3184        glTranslate(x,y,z)
3185        glMultMatrixf(B4mat.T)
3186        q = gluNewQuadric()
3187        gluSphere(q,radius,4,2)
3188        glPopMatrix()
3189               
3190    def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color):
3191        s1,s2,s3 = E
3192        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3193        glPushMatrix()
3194        glTranslate(x,y,z)
3195        glMultMatrixf(B4mat.T)
3196        glMultMatrixf(R4.T)
3197        glEnable(GL_NORMALIZE)
3198        glScale(s1,s2,s3)
3199        q = gluNewQuadric()
3200        gluSphere(q,ellipseProb,20,10)
3201        glDisable(GL_NORMALIZE)
3202        glPopMatrix()
3203       
3204    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
3205        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3206        glPushMatrix()
3207        glTranslate(x,y,z)
3208        glMultMatrixf(B4mat.T)
3209        for bond in Bonds:
3210            glPushMatrix()
3211            Dx = np.inner(Amat,bond)
3212            Z = np.sqrt(np.sum(Dx**2))
3213            if Z:
3214                azm = atan2d(-Dx[1],-Dx[0])
3215                phi = acosd(Dx[2]/Z)
3216                glRotate(-azm,0,0,1)
3217                glRotate(phi,1,0,0)
3218                q = gluNewQuadric()
3219                gluCylinder(q,radius,radius,Z,slice,2)
3220            glPopMatrix()           
3221        glPopMatrix()
3222               
3223    def RenderLines(x,y,z,Bonds,color):
3224        glShadeModel(GL_FLAT)
3225        xyz = np.array([x,y,z])
3226        glEnable(GL_COLOR_MATERIAL)
3227        glLineWidth(1)
3228        glColor3fv(color)
3229        glPushMatrix()
3230        glBegin(GL_LINES)
3231        for bond in Bonds:
3232            glVertex3fv(xyz)
3233            glVertex3fv(xyz+bond)
3234        glEnd()
3235        glColor4ubv([0,0,0,0])
3236        glPopMatrix()
3237        glDisable(GL_COLOR_MATERIAL)
3238        glShadeModel(GL_SMOOTH)
3239       
3240    def RenderPolyhedra(x,y,z,Faces,color):
3241        glPushMatrix()
3242        glTranslate(x,y,z)
3243        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3244        glShadeModel(GL_SMOOTH)
3245        glMultMatrixf(B4mat.T)
3246        for face,norm in Faces:
3247            glPolygonMode(GL_FRONT_AND_BACK,GL_FILL)
3248            glFrontFace(GL_CW)
3249            glNormal3fv(norm)
3250            glBegin(GL_TRIANGLES)
3251            for vert in face:
3252                glVertex3fv(vert)
3253            glEnd()
3254        glPopMatrix()
3255
3256    def RenderMapPeak(x,y,z,color,den):
3257        glShadeModel(GL_FLAT)
3258        xyz = np.array([x,y,z])
3259        glEnable(GL_COLOR_MATERIAL)
3260        glLineWidth(3)
3261        glColor3fv(color*den/255)
3262        glPushMatrix()
3263        glBegin(GL_LINES)
3264        for vec in mapPeakVecs:
3265            glVertex3fv(vec[0]+xyz)
3266            glVertex3fv(vec[1]+xyz)
3267        glEnd()
3268        glColor4ubv([0,0,0,0])
3269        glPopMatrix()
3270        glDisable(GL_COLOR_MATERIAL)
3271        glShadeModel(GL_SMOOTH)
3272       
3273    def RenderBackbone(Backbone,BackboneColor,radius):
3274        glPushMatrix()
3275        glMultMatrixf(B4mat.T)
3276        glEnable(GL_COLOR_MATERIAL)
3277        glShadeModel(GL_SMOOTH)
3278        gleSetJoinStyle(TUBE_NORM_EDGE | TUBE_JN_ANGLE | TUBE_JN_CAP)
3279        glePolyCylinder(Backbone,BackboneColor,radius)
3280        glPopMatrix()       
3281        glDisable(GL_COLOR_MATERIAL)
3282       
3283    def RenderLabel(x,y,z,label,r,color):       
3284        glPushMatrix()
3285        glTranslate(x,y,z)
3286        glMultMatrixf(B4mat.T)
3287        glDisable(GL_LIGHTING)
3288        glColor3fv(color)
3289        glRasterPos3f(0,0,0)
3290        for c in list(label):
3291            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
3292        glEnable(GL_LIGHTING)
3293        glPopMatrix()
3294       
3295    def RenderMap(rho,rhoXYZ,indx,Rok):
3296        glShadeModel(GL_FLAT)
3297        cLevel = drawingData['contourLevel']
3298        XYZ = []
3299        RC = []
3300        for i,xyz in enumerate(rhoXYZ):
3301            if not Rok[i]:
3302                x,y,z = xyz
3303                I,J,K = indx[i]
3304                alpha = 1.0
3305                if cLevel < 1.:
3306                    alpha = (abs(rho[I,J,K])/mapData['rhoMax']-cLevel)/(1.-cLevel)
3307                if rho[I,J,K] < 0.:
3308                    XYZ.append(xyz)
3309                    RC.append([0.1*alpha,Rd])
3310                else:
3311                    XYZ.append(xyz)
3312                    RC.append([0.1*alpha,Gr])
3313        RenderDots(XYZ,RC)
3314        glShadeModel(GL_SMOOTH)
3315                           
3316    def Draw(caller=''):
3317#useful debug?       
3318#        if caller:
3319#            print caller
3320# end of useful debug
3321        mapData = generalData['Map']
3322        pageName = ''
3323        page = getSelection()
3324        if page:
3325            pageName = G2frame.dataDisplay.GetPageText(page)
3326        rhoXYZ = []
3327        if len(mapData['rho']):
3328            VP = np.array(drawingData['viewPoint'][0])-np.array([.5,.5,.5])
3329            contLevel = drawingData['contourLevel']*mapData['rhoMax']
3330            if 'delt-F' in mapData['MapType']:
3331                rho = ma.array(mapData['rho'],mask=(np.abs(mapData['rho'])<contLevel))
3332            else:
3333                rho = ma.array(mapData['rho'],mask=(mapData['rho']<contLevel))
3334            steps = 1./np.array(rho.shape)
3335            incre = np.where(VP>=0,VP%steps,VP%steps-steps)
3336            Vsteps = -np.array(VP/steps,dtype='i')
3337            rho = np.roll(np.roll(np.roll(rho,Vsteps[0],axis=0),Vsteps[1],axis=1),Vsteps[2],axis=2)
3338            indx = np.array(ma.nonzero(rho)).T
3339            rhoXYZ = indx*steps+VP-incre
3340            Nc = len(rhoXYZ)
3341            rcube = 2000.*Vol/(ForthirdPI*Nc)
3342            rmax = math.exp(math.log(rcube)/3.)**2
3343            radius = min(drawingData['mapSize']**2,rmax)
3344            view = np.array(drawingData['viewPoint'][0])
3345            Rok = np.sum(np.inner(Amat,rhoXYZ-view).T**2,axis=1)>radius
3346        Ind = GetSelectedAtoms()
3347        VS = np.array(Page.canvas.GetSize())
3348        aspect = float(VS[0])/float(VS[1])
3349        cPos = drawingData['cameraPos']
3350        Zclip = drawingData['Zclip']*cPos/200.
3351        Q = drawingData['Quaternion']
3352        Tx,Ty,Tz = drawingData['viewPoint'][0]
3353        cx,ct,cs,ci = drawingData['atomPtrs']
3354        bondR = drawingData['bondRadius']
3355        G,g = G2lat.cell2Gmat(cell)
3356        GS = G
3357        GS[0][1] = GS[1][0] = math.sqrt(GS[0][0]*GS[1][1])
3358        GS[0][2] = GS[2][0] = math.sqrt(GS[0][0]*GS[2][2])
3359        GS[1][2] = GS[2][1] = math.sqrt(GS[1][1]*GS[2][2])
3360        ellipseProb = G2lat.criticalEllipse(drawingData['ellipseProb']/100.)
3361       
3362        SetBackground()
3363        glInitNames()
3364        glPushName(0)
3365       
3366        glMatrixMode(GL_PROJECTION)
3367        glLoadIdentity()
3368        glViewport(0,0,VS[0],VS[1])
3369        gluPerspective(20.,aspect,cPos-Zclip,cPos+Zclip)
3370        gluLookAt(0,0,cPos,0,0,0,0,1,0)
3371        SetLights()           
3372           
3373        glMatrixMode(GL_MODELVIEW)
3374        glLoadIdentity()
3375        matRot = G2mth.Q2Mat(Q)
3376        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3377        glMultMatrixf(matRot.T)
3378        glMultMatrixf(A4mat.T)
3379        glTranslate(-Tx,-Ty,-Tz)
3380        if drawingData['unitCellBox']:
3381            RenderBox()
3382        if drawingData['showABC']:
3383            x,y,z = drawingData['viewPoint'][0]
3384            RenderUnitVectors(x,y,z)
3385        Backbones = {}
3386        BackboneColor = []
3387        time0 = time.time()
3388#        glEnable(GL_BLEND)
3389#        glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
3390        for iat,atom in enumerate(drawingData['Atoms']):
3391            x,y,z = atom[cx:cx+3]
3392            Bonds = atom[-2]
3393            Faces = atom[-1]
3394            try:
3395                atNum = generalData['AtomTypes'].index(atom[ct])
3396            except ValueError:
3397                atNum = -1
3398            CL = atom[cs+2]
3399            atColor = np.array(CL)/255.
3400            if drawingData['showRigidBodies'] and atom[ci] in rbAtmDict:
3401                bndColor = Or
3402            else:
3403                bndColor = atColor
3404            if iat in Ind and G2frame.dataDisplay.GetPageText(getSelection()) != 'Map peaks':
3405                atColor = np.array(Gr)/255.
3406#            color += [.25,]
3407            radius = 0.5
3408            if atom[cs] != '':
3409                try:
3410                    glLoadName(atom[-3])
3411                except: #problem with old files - missing code
3412                    pass                   
3413            if 'balls' in atom[cs]:
3414                vdwScale = drawingData['vdwScale']
3415                ballScale = drawingData['ballScale']
3416                if atNum < 0:
3417                    radius = 0.2
3418                elif 'H' == atom[ct]:
3419                    if drawingData['showHydrogen']:
3420                        if 'vdW' in atom[cs] and atNum >= 0:
3421                            radius = vdwScale*generalData['vdWRadii'][atNum]
3422                        else:
3423                            radius = ballScale*drawingData['sizeH']
3424                    else:
3425                        radius = 0.0
3426                else:
3427                    if 'vdW' in atom[cs]:
3428                        radius = vdwScale*generalData['vdWRadii'][atNum]
3429                    else:
3430                        radius = ballScale*generalData['BondRadii'][atNum]
3431                RenderSphere(x,y,z,radius,atColor)
3432                if 'sticks' in atom[cs]:
3433                    RenderBonds(x,y,z,Bonds,bondR,bndColor)
3434            elif 'ellipsoids' in atom[cs]:
3435                RenderBonds(x,y,z,Bonds,bondR,bndColor)
3436                if atom[cs+3] == 'A':                   
3437                    Uij = atom[cs+5:cs+11]
3438                    U = np.multiply(G2spc.Uij2U(Uij),GS)
3439                    U = np.inner(Amat,np.inner(U,Amat).T)
3440                    E,R = nl.eigh(U)
3441                    R4 = np.concatenate((np.concatenate((R,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3442                    E = np.sqrt(E)
3443                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
3444                        pass
3445                    else:
3446                        RenderEllipsoid(x,y,z,ellipseProb,E,R4,atColor)                   
3447                else:
3448                    if atom[ct] == 'H' and not drawingData['showHydrogen']:
3449                        pass
3450                    else:
3451                        radius = ellipseProb*math.sqrt(abs(atom[cs+4]))
3452                        RenderSphere(x,y,z,radius,atColor)
3453            elif 'lines' in atom[cs]:
3454                radius = 0.1
3455                RenderLines(x,y,z,Bonds,bndColor)
3456#                RenderBonds(x,y,z,Bonds,0.05,color,6)
3457            elif atom[cs] == 'sticks':
3458                radius = 0.1
3459                RenderBonds(x,y,z,Bonds,bondR,bndColor)
3460            elif atom[cs] == 'polyhedra':
3461                RenderPolyhedra(x,y,z,Faces,atColor)
3462            elif atom[cs] == 'backbone':
3463                if atom[ct-1].split()[0] in ['C','N']:
3464                    if atom[2] not in Backbones:
3465                        Backbones[atom[2]] = []
3466                    Backbones[atom[2]].append(list(np.inner(Amat,np.array([x,y,z]))))
3467                    BackboneColor.append(list(atColor))
3468                   
3469            if atom[cs+1] == 'type':
3470                RenderLabel(x,y,z,'  '+atom[ct],radius,Gr)
3471            elif atom[cs+1] == 'name':
3472                RenderLabel(x,y,z,'  '+atom[ct-1],radius,Gr)
3473            elif atom[cs+1] == 'number':
3474                RenderLabel(x,y,z,'  '+str(iat),radius,Gr)
3475            elif atom[cs+1] == 'residue' and atom[ct-1] == 'CA':
3476                RenderLabel(x,y,z,'  '+atom[ct-4],radius,Gr)
3477            elif atom[cs+1] == '1-letter' and atom[ct-1] == 'CA':
3478                RenderLabel(x,y,z,'  '+atom[ct-3],radius,Gr)
3479            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
3480                RenderLabel(x,y,z,'  '+atom[ct-2],radius,Gr)
3481#        glDisable(GL_BLEND)
3482        if len(rhoXYZ):
3483            RenderMap(rho,rhoXYZ,indx,Rok)
3484        if len(mapPeaks):
3485            XYZ = mapPeaks.T[1:4].T
3486            mapBonds = FindPeaksBonds(XYZ)
3487            for ind,[mag,x,y,z,d] in enumerate(mapPeaks):
3488                if ind in Ind and pageName == 'Map peaks':
3489                    RenderMapPeak(x,y,z,Gr,1.0)
3490                else:
3491                    RenderMapPeak(x,y,z,Wt,mag/peakMax)
3492                if showBonds:
3493                    RenderLines(x,y,z,mapBonds[ind],Wt)
3494        if len(testRBObj) and pageName == 'RB Models':
3495            XYZ = G2mth.UpdateRBXYZ(Bmat,testRBObj['rbObj'],testRBObj['rbData'],testRBObj['rbType'])[0]
3496            rbBonds = FindPeaksBonds(XYZ)
3497            for ind,[x,y,z] in enumerate(XYZ):
3498                aType = testRBObj['rbAtTypes'][ind]
3499                name = '  '+aType+str(ind)
3500                color = np.array(testRBObj['AtInfo'][aType][1])
3501                RenderSphere(x,y,z,0.2,color/255.)
3502#                RenderMapPeak(x,y,z,color,1.0)
3503                RenderBonds(x,y,z,rbBonds[ind],0.03,Gr)
3504                RenderLabel(x,y,z,name,0.2,Or)
3505        if len(mcsaModels) > 1 and pageName == 'MC/SA':             #skip the default MD entry
3506            for ind,[x,y,z] in enumerate(mcsaXYZ):
3507                aType = atTypes[ind]
3508                name = '  '+aType+str(ind)
3509                color = np.array(MCSA['AtInfo'][aType][1])
3510                RenderSphere(x,y,z,0.2,color/255.)
3511                RenderBonds(x,y,z,mcsaBonds[ind],0.03,Gr)
3512                RenderLabel(x,y,z,name,0.2,Or)
3513        if Backbones:
3514            for chain in Backbones:
3515                Backbone = Backbones[chain]
3516                RenderBackbone(Backbone,BackboneColor,bondR)
3517#        print time.time()-time0
3518        Page.canvas.SwapBuffers()
3519       
3520    def OnSize(event):
3521        Draw('size')
3522       
3523    def OnFocus(event):         #not needed?? Bind commented out below
3524        Draw('focus')
3525       
3526    try:
3527        plotNum = G2frame.G2plotNB.plotList.index(generalData['Name'])
3528        Page = G2frame.G2plotNB.nb.GetPage(plotNum)       
3529    except ValueError:
3530        Plot = G2frame.G2plotNB.addOgl(generalData['Name'])
3531        plotNum = G2frame.G2plotNB.plotList.index(generalData['Name'])
3532        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3533        Page.views = False
3534        view = False
3535        altDown = False
3536    Page.SetFocus()
3537    Page.Choice = None
3538    if mapData['Flip']:
3539        choice = [' save as/key:','jpeg','tiff','bmp','u: roll up','d: roll down','l: roll left','r: roll right']
3540    else:
3541        choice = [' save as/key:','jpeg','tiff','bmp','c: center on 1/2,1/2,1/2','n: next','p: previous']
3542    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice)
3543    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
3544    cb.SetValue(' save as/key:')
3545    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
3546    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
3547    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
3548    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
3549    Page.canvas.Bind(wx.EVT_KEY_UP, OnKey)
3550    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
3551    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
3552#    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
3553    Page.camera['position'] = drawingData['cameraPos']
3554    Page.camera['viewPoint'] = np.inner(Amat,drawingData['viewPoint'][0])
3555    Page.camera['backColor'] = np.array(list(drawingData['backColor'])+[0,])/255.
3556    try:
3557        Page.canvas.SetCurrent()
3558    except:
3559        pass
3560    Draw('main')
3561       
3562################################################################################
3563#### Plot Rigid Body
3564################################################################################
3565
3566def PlotRigidBody(G2frame,rbType,AtInfo,rbData,defaults):
3567    '''RB plotting package. Can show rigid body structures as balls & sticks
3568    '''
3569
3570    def FindBonds(XYZ):
3571        rbTypes = rbData['rbTypes']
3572        Radii = []
3573        for Atype in rbTypes:
3574            Radii.append(AtInfo[Atype][0])
3575            if Atype == 'H':
3576                Radii[-1] = 0.5
3577        Radii = np.array(Radii)
3578        Bonds = [[] for i in range(len(Radii))]
3579        for i,xyz in enumerate(XYZ):
3580            Dx = XYZ-xyz
3581            dist = np.sqrt(np.sum(Dx**2,axis=1))
3582            sumR = Radii[i]+Radii
3583            IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
3584            for j in IndB[0]:
3585                Bonds[i].append(Dx[j]*Radii[i]/sumR[j])
3586                Bonds[j].append(-Dx[j]*Radii[j]/sumR[j])
3587        return Bonds
3588                       
3589    Wt = np.array([255,255,255])
3590    Rd = np.array([255,0,0])
3591    Gr = np.array([0,255,0])
3592    Bl = np.array([0,0,255])
3593    uBox = np.array([[0,0,0],[1,0,0],[0,1,0],[0,0,1]])
3594    uEdges = np.array([[uBox[0],uBox[1]],[uBox[0],uBox[2]],[uBox[0],uBox[3]]])
3595    uColors = [Rd,Gr,Bl]
3596    if rbType == 'Vector':
3597        atNames = [str(i)+':'+Ty for i,Ty in enumerate(rbData['rbTypes'])]
3598        XYZ = np.array([[0.,0.,0.] for Ty in rbData['rbTypes']])
3599        for imag,mag in enumerate(rbData['VectMag']):
3600            XYZ += mag*rbData['rbVect'][imag]
3601        Bonds = FindBonds(XYZ)
3602    elif rbType == 'Residue':
3603#        atNames = [str(i)+':'+Ty for i,Ty in enumerate(rbData['atNames'])]
3604        atNames = rbData['atNames']
3605        XYZ = np.copy(rbData['rbXYZ'])      #don't mess with original!
3606        Seq = rbData['rbSeq']
3607        for ia,ib,ang,mv in Seq:
3608            va = XYZ[ia]-XYZ[ib]
3609            Q = G2mth.AVdeg2Q(ang,va)
3610            for im in mv:
3611                vb = XYZ[im]-XYZ[ib]
3612                vb = G2mth.prodQVQ(Q,vb)
3613                XYZ[im] = XYZ[ib]+vb
3614        Bonds = FindBonds(XYZ)
3615    elif rbType == 'Z-matrix':
3616        pass
3617
3618#    def SetRBOrigin():
3619#        page = getSelection()
3620#        if page:
3621#            if G2frame.dataDisplay.GetPageText(page) == 'Rigid bodies':
3622#                G2frame.MapPeaksTable.SetData(mapPeaks)
3623#                panel = G2frame.dataDisplay.GetPage(page).GetChildren()
3624#                names = [child.GetName() for child in panel]
3625#                panel[names.index('grid window')].Refresh()
3626           
3627    def OnMouseDown(event):
3628        xy = event.GetPosition()
3629        defaults['oldxy'] = list(xy)
3630
3631    def OnMouseMove(event):
3632        newxy = event.GetPosition()
3633                               
3634        if event.Dragging():
3635            if event.LeftIsDown():
3636                SetRotation(newxy)
3637                Q = defaults['Quaternion']
3638                G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
3639#            elif event.RightIsDown():
3640#                SetRBOrigin(newxy)
3641            elif event.MiddleIsDown():
3642                SetRotationZ(newxy)
3643                Q = defaults['Quaternion']
3644                G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
3645            Draw('move')
3646       
3647    def OnMouseWheel(event):
3648        defaults['cameraPos'] += event.GetWheelRotation()/24
3649        defaults['cameraPos'] = max(10,min(500,defaults['cameraPos']))
3650        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(defaults['cameraPos']),1)
3651        Draw('wheel')
3652       
3653    def SetBackground():
3654        R,G,B,A = Page.camera['backColor']
3655        glClearColor(R,G,B,A)
3656        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
3657       
3658    def SetLights():
3659        glEnable(GL_DEPTH_TEST)
3660        glShadeModel(GL_FLAT)
3661        glEnable(GL_LIGHTING)
3662        glEnable(GL_LIGHT0)
3663        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,0)
3664        glLightfv(GL_LIGHT0,GL_AMBIENT,[1,1,1,.8])
3665        glLightfv(GL_LIGHT0,GL_DIFFUSE,[1,1,1,1])
3666       
3667#    def SetRBOrigin(newxy):
3668##first get translation vector in screen coords.       
3669#        oldxy = defaults['oldxy']
3670#        if not len(oldxy): oldxy = list(newxy)
3671#        dxy = newxy-oldxy
3672#        defaults['oldxy'] = list(newxy)
3673#        V = np.array([dxy[0],-dxy[1],0.])/100.
3674#        Q = defaults['Quaternion']
3675#        V = G2mth.prodQVQ(G2mth.invQ(Q),V)
3676#        rbData['rbXYZ'] += V
3677#        PlotRigidBody(G2frame,rbType,AtInfo,rbData,defaults)
3678#               
3679    def SetRotation(newxy):
3680#first get rotation vector in screen coords. & angle increment       
3681        oldxy = defaults['oldxy']
3682        if not len(oldxy): oldxy = list(newxy)
3683        dxy = newxy-oldxy
3684        defaults['oldxy'] = list(newxy)
3685        V = np.array([dxy[1],dxy[0],0.])
3686        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
3687# next transform vector back to xtal coordinates via inverse quaternion
3688# & make new quaternion
3689        Q = defaults['Quaternion']
3690        V = G2mth.prodQVQ(G2mth.invQ(Q),V)
3691        DQ = G2mth.AVdeg2Q(A,V)
3692        Q = G2mth.prodQQ(Q,DQ)
3693        defaults['Quaternion'] = Q
3694# finally get new view vector - last row of rotation matrix
3695        VD = G2mth.Q2Mat(Q)[2]
3696        VD /= np.sqrt(np.sum(VD**2))
3697        defaults['viewDir'] = VD
3698       
3699    def SetRotationZ(newxy):                       
3700#first get rotation vector (= view vector) in screen coords. & angle increment       
3701        View = glGetIntegerv(GL_VIEWPORT)
3702        cent = [View[2]/2,View[3]/2]
3703        oldxy = defaults['oldxy']
3704        if not len(oldxy): oldxy = list(newxy)
3705        dxy = newxy-oldxy
3706        defaults['oldxy'] = list(newxy)
3707        V = defaults['viewDir']
3708        A = [0,0]
3709        A[0] = dxy[1]*.25
3710        A[1] = dxy[0]*.25
3711        if newxy[0] > cent[0]:
3712            A[0] *= -1
3713        if newxy[1] < cent[1]:
3714            A[1] *= -1       
3715# next transform vector back to xtal coordinates & make new quaternion
3716        Q = defaults['Quaternion']
3717        Qx = G2mth.AVdeg2Q(A[0],V)
3718        Qy = G2mth.AVdeg2Q(A[1],V)
3719        Q = G2mth.prodQQ(Q,Qx)
3720        Q = G2mth.prodQQ(Q,Qy)
3721        defaults['Quaternion'] = Q
3722
3723    def RenderUnitVectors(x,y,z):
3724        xyz = np.array([x,y,z])
3725        glEnable(GL_COLOR_MATERIAL)
3726        glLineWidth(1)
3727        glPushMatrix()
3728        glTranslate(x,y,z)
3729        glBegin(GL_LINES)
3730        for line,color in zip(uEdges,uColors):
3731            glColor3ubv(color)
3732            glVertex3fv(-line[1])
3733            glVertex3fv(line[1])
3734        glEnd()
3735        glPopMatrix()
3736        glColor4ubv([0,0,0,0])
3737        glDisable(GL_COLOR_MATERIAL)
3738               
3739    def RenderSphere(x,y,z,radius,color):
3740        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3741        glPushMatrix()
3742        glTranslate(x,y,z)
3743        q = gluNewQuadric()
3744        gluSphere(q,radius,20,10)
3745        glPopMatrix()
3746       
3747    def RenderBonds(x,y,z,Bonds,radius,color,slice=20):
3748        glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color)
3749        glPushMatrix()
3750        glTranslate(x,y,z)
3751        for Dx in Bonds:
3752            glPushMatrix()
3753            Z = np.sqrt(np.sum(Dx**2))
3754            if Z:
3755                azm = atan2d(-Dx[1],-Dx[0])
3756                phi = acosd(Dx[2]/Z)
3757                glRotate(-azm,0,0,1)
3758                glRotate(phi,1,0,0)
3759                q = gluNewQuadric()
3760                gluCylinder(q,radius,radius,Z,slice,2)
3761            glPopMatrix()           
3762        glPopMatrix()
3763               
3764    def RenderLabel(x,y,z,label):       
3765        glPushMatrix()
3766        glTranslate(x,y,z)
3767        glDisable(GL_LIGHTING)
3768        glColor3f(1.0,1.0,1.0)
3769        glRasterPos3f(0,0,0)
3770        for c in list(label):
3771            glutBitmapCharacter(GLUT_BITMAP_8_BY_13,ord(c))
3772        glEnable(GL_LIGHTING)
3773        glPopMatrix()
3774       
3775    def Draw(caller=''):
3776#useful debug?       
3777#        if caller:
3778#            print caller
3779# end of useful debug
3780        cPos = defaults['cameraPos']
3781        VS = np.array(Page.canvas.GetSize())
3782        aspect = float(VS[0])/float(VS[1])
3783        Zclip = 500.0
3784        Q = defaults['Quaternion']
3785        SetBackground()
3786        glInitNames()
3787        glPushName(0)
3788       
3789        glMatrixMode(GL_PROJECTION)
3790        glLoadIdentity()
3791        glViewport(0,0,VS[0],VS[1])
3792        gluPerspective(20.,aspect,0.1,500.)
3793        gluLookAt(0,0,cPos,0,0,0,0,1,0)
3794        SetLights()           
3795           
3796        glMatrixMode(GL_MODELVIEW)
3797        glLoadIdentity()
3798        matRot = G2mth.Q2Mat(Q)
3799        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
3800        glMultMatrixf(matRot.T)
3801        RenderUnitVectors(0.,0.,0.)
3802        radius = 0.2
3803        for iat,atom in enumerate(XYZ):
3804            x,y,z = atom
3805            CL = AtInfo[rbData['rbTypes'][iat]][1]
3806            color = np.array(CL)/255.
3807            RenderSphere(x,y,z,radius,color)
3808            RenderBonds(x,y,z,Bonds[iat],0.05,color)
3809            RenderLabel(x,y,z,'  '+atNames[iat])
3810        Page.canvas.SwapBuffers()
3811
3812    def OnSize(event):
3813        Draw('size')
3814       
3815    try:
3816        plotNum = G2frame.G2plotNB.plotList.index('Rigid body')
3817        Page = G2frame.G2plotNB.nb.GetPage(plotNum)       
3818    except ValueError:
3819        Plot = G2frame.G2plotNB.addOgl('Rigid body')
3820        plotNum = G2frame.G2plotNB.plotList.index('Rigid body')
3821        Page = G2frame.G2plotNB.nb.GetPage(plotNum)
3822        Page.views = False
3823        view = False
3824        altDown = False
3825    Page.SetFocus()
3826    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
3827    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
3828    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
3829    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
3830    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
3831    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
3832    Page.camera['position'] = defaults['cameraPos']
3833    Page.camera['backColor'] = np.array([0,0,0,0])
3834    Page.canvas.SetCurrent()
3835    Draw('main')
Note: See TracBrowser for help on using the repository browser.