Context Navigation

source: trunk/GSASIIplot.py @ 1032

Last change on this file since 1032 was 1032, checked in by toby, 10 years ago
Use Q not q; only read .gpx files
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 161.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: