Context Navigation

source: trunk/GSASIIplot.py @ 998

Last change on this file since 998 was 998, checked in by vondreele, 10 years ago
make widths correct in Export HKLs from PWD data
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 159.3 KB

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

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

Download in other formats: