Context Navigation

source: trunk/GSASIIplot.py @ 1204

Last change on this file since 1204 was 1204, checked in by vondreele, 9 years ago
SASD now in q only (not 2-theta); fix minor bugs
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 172.6 KB

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

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

Download in other formats: