Context Navigation

source: trunk/GSASIIIO.py @ 163

Last change on this file since 163 was 163, checked in by vondreele, 15 years ago
new versions
File size: 38.0 KB

Line
1	"""GSASIIIO: functions for IO of data
2	Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3	"""
4
5	import wx
6	import math
7	import numpy as np
8	import cPickle
9	import sys
10	import GSASIIpath
11	import GSASIIgrid as G2gd
12	import os.path as ospath
13
14	def sfloat(S):
15	if S.strip():
16	return float(S)
17	else:
18	return 0.0
19
20	def sint(S):
21	if S.strip():
22	return int(S)
23	else:
24	return 0
25
26	def SelectPowderData(self, filename):
27	"""Selects banks of data from a filename of any GSAS powder data format
28	Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
29	Returns - a list of banks to be read; each entry in list is a tuple containing:
30	filename: same as input filename
31	Pos: position for start of data; record just after BANK record
32	Bank: the BANK record
33	"""
34	File = open(filename,'Ur')
35	Title = '''
36	First line of this file:
37	'''+File.readline()
38	dlg = wx.MessageDialog(self, Title, 'Is this the file you want?',
39	wx.YES_NO \| wx.ICON_QUESTION)
40	try:
41	result = dlg.ShowModal()
42	finally:
43	dlg.Destroy()
44	if result == wx.ID_NO: return (0,0)
45
46	self.IparmName = GetInstrumentFile(self,filename)
47	if self.IparmName:
48	Iparm = GetInstrumentData(self.IparmName)
49	else:
50	Iparm = {} #Assume CuKa lab data if no iparm file
51	Iparm['INS HTYPE '] = 'PXC '
52	Iparm['INS 1 ICONS'] = ' 1.540500 1.544300 0.0 0 0.7 0 0.5 '
53	Iparm['INS 1PRCF1 '] = ' 3 8 0.01 '
54	Iparm['INS 1PRCF11'] = ' 2.000000E+00 -2.000000E+00 5.000000E+00 0.000000E+00 '
55	Iparm['INS 1PRCF12'] = ' 0.000000E+00 0.000000E+00 0.150000E-01 0.150000E-01 '
56	S = 1
57	Banks = []
58	Pos = []
59	FoundData = []
60	Comments = []
61	wx.BeginBusyCursor()
62	try:
63	while S:
64	S = File.readline()
65	if S[:1] != '#':
66	if S[:4] == 'BANK':
67	Banks.append(S)
68	Pos.append(File.tell())
69	else:
70	Comments.append(S[:-1])
71	File.close()
72	finally:
73	wx.EndBusyCursor()
74	if Comments:
75	print 'Comments on file:'
76	for Comment in Comments: print Comment
77	if Banks:
78	result = [0]
79	if len(Banks) >= 2:
80	dlg = wx.MultiChoiceDialog(self, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
81	try:
82	if dlg.ShowModal() == wx.ID_OK:
83	result = dlg.GetSelections()
84	else:
85	result = []
86	finally:
87	dlg.Destroy()
88	for i in result:
89	FoundData.append((filename,Pos[i],Banks[i]))
90	else:
91	dlg = wx.MessageDialog(self, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK \| wx.ICON_ERROR)
92	try:
93	result = dlg.ShowModal()
94	finally:
95	dlg.Destroy()
96	return FoundData,Iparm,Comments
97
98	def GetInstrumentFile(self,filename):
99	import os.path as op
100	dlg = wx.FileDialog(self,'Choose an instrument file','.', '', 'GSAS iparm file (.prm)\|.prm\|All files(.)\|.', wx.OPEN)
101	if self.dirname:
102	dlg.SetDirectory(self.dirname)
103	Tname = filename[:filename.index('.')]+'.prm'
104	if op.exists(Tname):
105	self.IparmName = Tname
106	if self.IparmName: dlg.SetFilename(self.IparmName)
107	filename = ''
108	try:
109	if dlg.ShowModal() == wx.ID_OK:
110	filename = dlg.GetPath()
111	finally:
112	dlg.Destroy()
113	return filename
114
115	def GetInstrumentData(IparmName):
116	file = open(IparmName, 'Ur')
117	S = 1
118	Iparm = {}
119	while S:
120	S = file.readline()
121	Iparm[S[:12]] = S[12:-1]
122	return Iparm
123
124	def GetPowderPeaks(fileName):
125	sind = lambda x: math.sin(x*math.pi/180.)
126	asind = lambda x: 180.*math.asin(x)/math.pi
127	Cuka = 1.54052
128	File = open(fileName,'Ur')
129	Comments = []
130	peaks = []
131	S = File.readline()
132	while S:
133	if S[:1] == '#':
134	Comments.append(S[:-1])
135	else:
136	item = S.split()
137	if len(item) == 1:
138	peaks.append([float(item[0]),1.0])
139	elif len(item) > 1:
140	peaks.append([float(item[0]),float(item[0])])
141	S = File.readline()
142	File.close()
143	if Comments:
144	print 'Comments on file:'
145	for Comment in Comments: print Comment
146	Peaks = []
147	if peaks[0][0] > peaks[-1][0]: # d-spacings - assume CuKa
148	for peak in peaks:
149	dsp = peak[0]
150	sth = Cuka/(2.0*dsp)
151	if sth < 1.0:
152	tth = 2.0*asind(sth)
153	else:
154	break
155	Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
156	else: #2-thetas - assume Cuka (for now)
157	for peak in peaks:
158	tth = peak[0]
159	dsp = Cuka/(2.0*sind(tth/2.0))
160	Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
161	return Comments,Peaks
162
163	def GetHKLData(filename):
164	print 'Reading: '+filename
165	File = open(filename,'Ur')
166	HKLref = []
167	HKLmin = [1000,1000,1000]
168	HKLmax = [0,0,0]
169	FoMax = 0
170	ifFc = False
171	S = File.readline()
172	while '#' in S[0]: #get past comments if any
173	S = File.readline()
174	if '_' in S: #cif style .hkl file
175	while 'loop_' not in S: #skip preliminaries if any - can't have 'loop_' in them!
176	S = File.readline()
177	S = File.readline() #get past 'loop_' line
178	pos = 0
179	hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
180	while S:
181	if '_' in S:
182	if 'index_h' in S:
183	hpos = pos
184	elif 'index_k' in S:
185	kpos = pos
186	elif 'index_l' in S:
187	lpos = pos
188	elif 'F_squared_meas' in S:
189	Fosqpos = pos
190	elif 'F_squared_calc' in S:
191	Fcsqpos = pos
192	elif 'F_squared_sigma' in S:
193	sigpos = pos
194	pos += 1
195	else:
196	data = S.split()
197	if data: #avoid blank lines
198	HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
199	h,k,l = HKL
200	Fosq = float(data[Fosqpos])
201	if sigpos != -1:
202	sigFosq = float(data[sigpos])
203	else:
204	sigFosq = 1.
205	if Fcsqpos != -1:
206	Fcsq = float(data[Fcsqpos])
207	if Fcsq:
208	ifFc = True
209	else:
210	Fcsq = 0.
211
212	HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
213	HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
214	FoMax = max(FoMax,Fosq)
215	HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0]) #room for Fc, Fcp, Fcpp & phase
216	S = File.readline()
217	else: #dumb h,k,l,Fo,sigFo .hkl file
218	while S:
219	h,k,l,Fo,sigFo = S.split()
220	HKL = np.array([int(h),int(k),int(l)])
221	h,k,l = HKL
222	Fo = float(Fo)
223	sigFo = float(sigFo)
224	HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
225	HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
226	FoMax = max(FoMax,Fo)
227	HKLref.append([HKL,Fo*2,2.Fo*sigFo,0,0,0,0]) #room for Fc, Fcp, Fcpp & phase
228	S = File.readline()
229	File.close()
230	return HKLref,HKLmin,HKLmax,FoMax,ifFc
231
232	def GetPowderData(filename,Pos,Bank,DataType):
233	'''Reads one BANK of data from GSAS raw powder data file
234	input:
235	filename: GSAS raw powder file dataname
236	Pos: start of data in file just after BANK record
237	Bank: the BANK record
238	DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
239	returns: list [x,y,e,yc,yb]
240	x: np.array of x-axis values
241	y: np.array of powder pattern intensities
242	w: np.array of w=sig(intensity)^2 values
243	yc: np.array of calc. intensities (zero)
244	yb: np.array of calc. background (zero)
245	yd: np.array of obs-calc profiles
246	'''
247	print 'Reading: '+filename
248	print 'Bank: '+Bank[:-1]
249	if 'FXYE' in Bank:
250	return GetFXYEdata(filename,Pos,Bank,DataType)
251	elif 'FXY' in Bank:
252	return GetFXYdata(filename,Pos,Bank,DataType)
253	elif 'ESD' in Bank:
254	return GetESDdata(filename,Pos,Bank,DataType)
255	elif 'STD' in Bank:
256	return GetSTDdata(filename,Pos,Bank,DataType)
257	else:
258	return GetSTDdata(filename,Pos,Bank,DataType)
259	return []
260
261	def GetFXYEdata(filename,Pos,Bank,DataType):
262	File = open(filename,'Ur')
263	File.seek(Pos)
264	x = []
265	y = []
266	w = []
267	S = File.readline()
268	while S and S[:4] != 'BANK':
269	vals = S.split()
270	if DataType[2] == 'C':
271	x.append(float(vals[0])/100.) #CW: from centidegrees to degrees
272	elif DataType[2] == 'T':
273	x.append(float(vals[0])/1000.0) #TOF: from musec to millisec
274	f = float(vals[1])
275	if f <= 0.0:
276	y.append(0.0)
277	w.append(1.0)
278	else:
279	y.append(float(vals[1]))
280	w.append(1.0/float(vals[2])**2)
281	S = File.readline()
282	File.close()
283	N = len(x)
284	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
285
286	def GetFXYdata(filename,Pos,Bank,DataType):
287	File = open(filename,'Ur')
288	File.seek(Pos)
289	x = []
290	y = []
291	w = []
292	S = File.readline()
293	while S and S[:4] != 'BANK':
294	vals = S.split()
295	if DataType[2] == 'C':
296	x.append(float(vals[0])/100.) #CW: from centidegrees to degrees
297	elif DataType[2] == 'T':
298	x.append(float(vals[0])/1000.0) #TOF: from musec to millisec
299	f = float(vals[1])
300	if f > 0.0:
301	y.append(f)
302	w.append(1.0/f)
303	else:
304	y.append(0.0)
305	w.append(1.0)
306	S = File.readline()
307	File.close()
308	N = len(x)
309	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
310
311	def GetESDdata(filename,Pos,Bank,DataType):
312	File = open(filename,'Ur')
313	cons = Bank.split()
314	if DataType[2] == 'C':
315	start = float(cons[5])/100.0 #CW: from centidegrees to degrees
316	step = float(cons[6])/100.0
317	elif DataType[2] == 'T':
318	start = float(cons[5])/1000.0 #TOF: from musec to millisec
319	step = float(cons[6])/1000.0
320	File.seek(Pos)
321	x = []
322	y = []
323	w = []
324	S = File.readline()
325	j = 0
326	while S and S[:4] != 'BANK':
327	for i in range(0,80,16):
328	xi = start+step*j
329	yi = sfloat(S[i:i+8])
330	ei = sfloat(S[i+8:i+16])
331	x.append(xi)
332	if yi > 0.0:
333	y.append(yi)
334	w.append(1.0/ei**2)
335	else:
336	y.append(0.0)
337	w.append(1.0)
338	j += 1
339	S = File.readline()
340	File.close()
341	N = len(x)
342	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
343
344	def GetSTDdata(filename,Pos,Bank,DataType):
345	File = open(filename,'Ur')
346	cons = Bank.split()
347	Nch = cons[2]
348	if DataType[2] == 'C':
349	start = float(cons[5])/100.0 #CW: from centidegrees to degrees
350	step = float(cons[6])/100.0
351	elif DataType[2] == 'T':
352	start = float(cons[5])/1000.0 #TOF: from musec to millisec - not likely!
353	step = float(cons[6])/1000.0
354	File.seek(Pos)
355	x = []
356	y = []
357	w = []
358	S = File.readline()
359	j = 0
360	while S and S[:4] != 'BANK':
361	for i in range(0,80,8):
362	xi = start+step*j
363	ni = max(sint(S[i:i+2]),1)
364	yi = max(sfloat(S[i+2:i+8]),0.0)
365	if yi:
366	ei = math.sqrt(yi*ni)
367	else:
368	yi = 0.0
369	ei = 1.0
370	j += 1
371	if j < Nch:
372	x.append(xi)
373	y.append(yi)
374	w.append(1.0/ei**2)
375	S = File.readline()
376	File.close()
377	N = len(x)
378	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
379
380	def GetImageData(imagefile,imageOnly=False):
381	ext = ospath.splitext(imagefile)[1]
382	Comments = []
383	if ext == '.tif':
384	Comments,Data,Size,Image = GetTifData(imagefile)
385	elif ext == '.img':
386	Comments,Data,Size,Image = GetImgData(imagefile)
387	Image[0][0] = 0
388	elif ext == '.mar3450' or ext == '.mar2300':
389	Comments,Data,Size,Image = GetMAR345Data(imagefile)
390	elif ext in ['.sum','.avg']:
391	Comments,Data,Size,Image = GetGEsumData(imagefile)
392	elif ext == '.G2img':
393	return GetG2Image(imagefile)
394	if imageOnly:
395	return Image
396	else:
397	return Comments,Data,Size,Image
398
399	def PutG2Image(filename,image):
400	File = open(filename,'wb')
401	cPickle.dump(image,File,1)
402	File.close()
403	return
404
405	def GetG2Image(filename):
406	File = open(filename,'rb')
407	image = cPickle.load(File)
408	File.close()
409	return image
410
411	def GetGEsumData(filename,imageOnly=False):
412	import array as ar
413	if not imageOnly:
414	print 'Read GE sum file: ',filename
415	File = open(filename,'rb')
416	size = 2048
417	if '.sum' in filename:
418	head = ['GE detector sum data from APS 1-ID',]
419	if '.avg' in filename:
420	head = ['GE detector avg data from APS 1-ID',]
421	image = np.zeros(shape=(size,size),dtype=np.int32)
422	row = 0
423	pos = 0
424	while row < size:
425	File.seek(pos)
426	if '.sum' in filename:
427	line = ar.array('f',File.read(4*size))
428	pos += 4*size
429	elif '.avg' in filename:
430	line = ar.array('H',File.read(2*size))
431	pos += 2*size
432	image[row] = np.asarray(line)
433	row += 1
434	data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8]}
435	File.close()
436	if imageOnly:
437	return image
438	else:
439	return head,data,size,image
440
441	def GetImgData(filename,imageOnly=False):
442	import struct as st
443	import array as ar
444	if not imageOnly:
445	print 'Read ADSC img file: ',filename
446	File = open(filename,'rb')
447	head = File.read(511)
448	lines = head.split('\n')
449	head = []
450	center = [0,0]
451	for line in lines[1:-2]:
452	line = line.strip()[:-1]
453	if line:
454	if 'SIZE1' in line:
455	size = int(line.split('=')[1])
456	elif 'WAVELENGTH' in line:
457	wave = float(line.split('=')[1])
458	elif 'BIN' in line:
459	if line.split('=')[1] == '2x2':
460	pixel=(102,102)
461	else:
462	pixel = (51,51)
463	elif 'DISTANCE' in line:
464	distance = float(line.split('=')[1])
465	elif 'CENTER_X' in line:
466	center[0] = float(line.split('=')[1])
467	elif 'CENTER_Y' in line:
468	center[1] = float(line.split('=')[1])
469	head.append(line)
470	data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
471	image = []
472	row = 0
473	pos = 512
474	image = np.zeros(shape=(size,size),dtype=np.int32)
475	while row < size:
476	File.seek(pos)
477	line = ar.array('H',File.read(2*size))
478	image[row] = np.asarray(line)
479	row += 1
480	pos += 2*size
481	File.close()
482	if imageOnly:
483	return image
484	else:
485	return lines[1:-2],data,size,image
486
487	def GetMAR345Data(filename,imageOnly=False):
488	import array as ar
489	import struct as st
490	try:
491	import pack_f as pf
492	except:
493	msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
494	caption="Import Error",
495	style=wx.ICON_ERROR \| wx.OK \| wx.STAY_ON_TOP)
496	msg.ShowModal()
497	return None,None,None,None
498
499	if not imageOnly:
500	print 'Read Mar345 file: ',filename
501	File = open(filename,'rb')
502	head = File.read(4095)
503	numbers = st.unpack('<iiiiiiiiii',head[:40])
504	lines = head[128:].split('\n')
505	head = []
506	for line in lines:
507	line = line.strip()
508	if 'PIXEL' in line:
509	values = line.split()
510	pixel = (int(values[2]),int(values[4])) #in microns
511	elif 'WAVELENGTH' in line:
512	wave = float(line.split()[1])
513	elif 'DISTANCE' in line:
514	distance = float(line.split()[1]) #in mm
515	elif 'CENTER' in line:
516	values = line.split()
517	center = [float(values[2])/10.,float(values[4])/10.] #make in mm from pixels
518	if line:
519	head.append(line)
520	data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
521	for line in head:
522	if 'FORMAT' in line[0:6]:
523	items = line.split()
524	size = int(items[1])
525	pos = 4096
526	File.seek(pos)
527	line = File.read(8)
528	while 'CCP4' not in line: #get past overflow list for now
529	line = File.read(8)
530	pos += 8
531	pos += 37
532	File.seek(pos)
533	raw = File.read()
534	File.close()
535	image = np.zeros(shape=(size,size),dtype=np.int32)
536	image = pf.pack_f(len(raw),raw,size,image)
537	if imageOnly:
538	return image.T #transpose to get it right way around
539	else:
540	return head,data,size,image.T
541
542	def GetTifData(filename,imageOnly=False):
543	# only works for APS Perkin-Elmer detector data files in "TIFF" format that are readable by Fit2D
544	import struct as st
545	import array as ar
546	File = open(filename,'rb')
547	dataType = 5
548	try:
549	Meta = open(filename+'.metadata','Ur')
550	head = Meta.readlines()
551	for line in head:
552	line = line.strip()
553	if 'dataType' in line:
554	dataType = int(line.split('=')[1])
555	Meta.close()
556	except IOError:
557	print 'no metadata file found - will try to read file anyway'
558	head = 'no metadata file found'
559	tag = File.read(3)
560	if tag != 'II*':
561	lines = ['not a detector tiff file',]
562	return lines,0,0
563	size,Ityp = st.unpack('<ii',File.read(8))
564	if Ityp == 0:
565	tifType = 'Pilatus'
566	pixy = (172,172)
567	pos = 4096
568	if not imageOnly:
569	print 'Read Pilatus tiff file: ',filename
570	elif Ityp == 1:
571	tifType = 'PE'
572	pixy = (200,200)
573	pos = 8
574	if not imageOnly:
575	print 'Read APS PE-detector tiff file: ',filename
576	elif Ityp == 3328:
577	tifType = 'MAR'
578	pixy = (79,79)
579	pos = 4096
580	if not imageOnly:
581	print 'Read MAR CCD tiff file: ',filename
582	else:
583	lines = 'unknown tif type'
584	return lines,0,0
585	image = np.zeros(shape=(size,size),dtype=np.int32)
586	row = 0
587	while row < size:
588	File.seek(pos)
589	if 'PE' in tifType:
590	if dataType == 5:
591	line = ar.array('f',File.read(4*size))
592	else:
593	line = ar.array('l',File.read(4*size))
594	pos += 4*size
595	elif 'MAR' in tifType:
596	line = ar.array('H',File.read(2*size))
597	pos += 2*size
598	image[row] = np.asarray(line)
599	row += 1
600	data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':[204.8,204.8]}
601	File.close()
602	if imageOnly:
603	return image
604	else:
605	return head,data,size,image
606
607	def ProjFileOpen(self):
608	file = open(self.GSASprojectfile,'rb')
609	print 'load from file: ',self.GSASprojectfile
610	wx.BeginBusyCursor()
611	try:
612	while True:
613	try:
614	data = cPickle.load(file)
615	except EOFError:
616	break
617	datum = data[0]
618	print 'load: ',datum[0]
619
620	#temporary fixes to old project files
621	#fix to convert old style list arrays to numpy arrays
622	if 'PWDR' in datum[0] and 'list' in str(type(datum[1][1][0])):
623	X = datum[1][1]
624	X = [np.array(X[0]),np.array(X[1]),np.array(X[2]),np.array(X[3]),np.array(X[4]),np.array(X[5])]
625	datum[1] = [datum[1][0],X]
626	print 'powder data converted to numpy arrays'
627	#temporary fix to insert 'PWDR' in front of powder names
628	if 'PKS' not in datum[0] and 'IMG' not in datum[0] and 'SNGL' not in datum[0]:
629	if datum[0] not in ['Notebook','Controls','Phases'] and 'PWDR' not in datum[0]:
630	datum[0] = 'PWDR '+datum[0]
631	print 'add PWDR to powder names'
632	#end of temporary fixes
633
634	Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
635	self.PatternTree.SetItemPyData(Id,datum[1])
636	for datus in data[1:]:
637	print ' load: ',datus[0]
638
639	#temporary fix to add azimuthal angle to instrument parameters
640	if 'PWDR' in datum[0] and 'Instrument Parameters' in datus[0]:
641	if len(datus[1][0]) == 10 or len(datus[1][0]) == 12:
642	datus[1][0] += (0.0,) #add missing azimuthal angle
643	datus[1][1].append(0.0)
644	datus[1][2].append(0.0)
645	datus[1][3].append('Azimuth')
646	print 'add azimuth to instrument parameters'
647	#end of temporary fix
648
649	sub = self.PatternTree.AppendItem(Id,datus[0])
650	self.PatternTree.SetItemPyData(sub,datus[1])
651
652	#temporary fix to add Comments to powder data sets
653	if 'PWDR' in datum[0] and not G2gd.GetPatternTreeItemId(self,Id, 'Comments'):
654	print 'no comments - add comments'
655	sub = self.PatternTree.AppendItem(Id,'Comments')
656	self.PatternTree.SetItemPyData(sub,['no comments'])
657	#end of temporary fix
658
659	if 'IMG' in datum[0]: #retreive image default flag & data if set
660	Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
661	if Data['setDefault']:
662	self.imageDefault = Data
663	#temporary fix to add masks
664	if not G2gd.GetPatternTreeItemId(self,Id, 'Masks'):
665	Imin,Imax = Data['range'][0]
666	Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]}
667	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Masks'),Masks)
668	#end of temporary fix
669
670	file.close()
671
672	#temporary fix to add Notebook & Controls to project
673	if not G2gd.GetPatternTreeItemId(self,self.root,'Notebook'):
674	sub = self.PatternTree.AppendItem(parent=self.root,text='Notebook')
675	self.PatternTree.SetItemPyData(sub,[''])
676	sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
677	self.PatternTree.SetItemPyData(sub,[0])
678	print 'add Notebook and Controls to project'
679
680	finally:
681	wx.EndBusyCursor()
682	print 'project load successful'
683	self.NewPlot = True
684
685	def ProjFileSave(self):
686	if not self.PatternTree.IsEmpty():
687	file = open(self.GSASprojectfile,'wb')
688	print 'save to file: ',self.GSASprojectfile
689	wx.BeginBusyCursor()
690	try:
691	item, cookie = self.PatternTree.GetFirstChild(self.root)
692	while item:
693	data = []
694	name = self.PatternTree.GetItemText(item)
695	print 'save: ',name
696	data.append([name,self.PatternTree.GetItemPyData(item)])
697	item2, cookie2 = self.PatternTree.GetFirstChild(item)
698	while item2:
699	name = self.PatternTree.GetItemText(item2)
700	print ' save: ',name
701	data.append([name,self.PatternTree.GetItemPyData(item2)])
702	item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)
703	item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
704	cPickle.dump(data,file,1)
705	file.close()
706	finally:
707	wx.EndBusyCursor()
708	print 'project save successful'
709
710	def SaveIntegration(self,PickId,data):
711	azms = self.Integrate[1]
712	X = self.Integrate[2][:-1]
713	Xminmax = [X[0],X[-1]]
714	N = len(X)
715	Id = self.PatternTree.GetItemParent(PickId)
716	name = self.PatternTree.GetItemText(Id)
717	name = name.replace('IMG ','PWDR ')
718	Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
719	names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth']
720	codes = [0 for i in range(11)]
721	parms = ['PXC',data['wavelength'],0.0,0.0,1.0,-1.0,0.3,0.0,1.0,0.0,0.0]
722	Azms = [(azms[i+1]+azms[i])/2. for i in range(len(azms)-1)]
723	for i,azm in enumerate(Azms):
724	item, cookie = self.PatternTree.GetFirstChild(self.root)
725	Id = 0
726	while item:
727	Name = self.PatternTree.GetItemText(item)
728	if name == Name:
729	Id = item
730	item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
731	parms[10] = azm
732	Y = self.Integrate[0][i]
733	W = np.sqrt(Y)
734	if Id:
735	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)
736	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
737	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
738	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms,codes,names])
739	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
740	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
741	self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])
742	else:
743	Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(azm))
744	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)
745	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
746	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
747	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms,codes,names])
748	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
749	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
750	self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
751	self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
752	self.PatternTree.SelectItem(Id)
753	self.PatternTree.Expand(Id)
754	self.PatternId = Id
755
756	def powderFxyeSave(self,powderfile):
757	file = open(powderfile,'w')
758	prm = open(powderfile.strip('fxye')+'prm','w') #old style GSAS parm file
759	print 'save powder pattern to file: ',powderfile
760	wx.BeginBusyCursor()
761	Inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
762	self.PickId, 'Instrument Parameters'))[1]
763	if len(Inst) == 11: #single wavelength
764	lam1 = Inst[1]
765	lam2 = 0.0
766	GU,GV,GW = Inst[4:7]
767	LX,LY = Inst[7:9]
768	SL = HL = Inst[9]/2.0
769	else: #Ka1 & Ka2
770	lam1 = Inst[1]
771	lam2 = Inst[2]
772	GU,GV,GW = Inst[6:9]
773	LX,LY = Inst[9:11]
774	SL = HL = Inst[11]/2.0
775	prm.write( ' 123456789012345678901234567890123456789012345678901234567890 '+'\n')
776	prm.write( 'INS BANK 1 '+'\n')
777	prm.write( 'INS HTYPE PXCR '+'\n')
778	prm.write(('INS 1 ICONS%10.7f%10.7f 0.0000 0.990 0 0.500 '+'\n')%(lam1,lam2))
779	prm.write( 'INS 1 IRAD 0 '+'\n')
780	prm.write( 'INS 1I HEAD '+'\n')
781	prm.write( 'INS 1I ITYP 0 0.0000 180.0000 1 '+'\n')
782	prm.write( 'INS 1PRCF1 3 8 0.00100 '+'\n')
783	prm.write(('INS 1PRCF11 %15.6g%15.6g%15.6g%15.6g '+'\n')%(GU,GV,GW,0.0))
784	prm.write(('INS 1PRCF12 %15.6g%15.6g%15.6g%15.6g '+'\n')%(LX,LY,SL,HL))
785	prm.close()
786	try:
787	x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
788	file.write(powderfile+'\n')
789	file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
790	100.x[0],100.(x[1]-x[0])))
791	XYW = zip(x,y,w)
792	for X,Y,W in XYW:
793	file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,W))
794	file.close()
795	finally:
796	wx.EndBusyCursor()
797	print 'powder pattern file written'
798
799	def powderXyeSave(self,powderfile):
800	file = open(powderfile,'w')
801	print 'save powder pattern to file: ',powderfile
802	wx.BeginBusyCursor()
803	try:
804	x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
805	XYW = zip(x,y,w)
806	for X,Y,W in XYW:
807	file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
808	file.close()
809	finally:
810	wx.EndBusyCursor()
811	print 'powder pattern file written'
812
813	def PeakListSave(self,file,peaks):
814	print 'save peak list to file: ',self.peaklistfile
815	if not peaks:
816	dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
817	try:
818	result = dlg.ShowModal()
819	finally:
820	dlg.Destroy()
821	return
822	for peak in peaks:
823	file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
824	(peak[0],peak[2],peak[4],peak[6]))
825	print 'peak list saved'
826
827	def IndexPeakListSave(self,peaks):
828	file = open(self.peaklistfile,'wa')
829	print 'save index peak list to file: ',self.peaklistfile
830	wx.BeginBusyCursor()
831	try:
832	if not peaks:
833	dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
834	try:
835	result = dlg.ShowModal()
836	finally:
837	dlg.Destroy()
838	return
839	for peak in peaks:
840	file.write("%12.6f\n" % (peak[7]))
841	file.close()
842	finally:
843	wx.EndBusyCursor()
844	print 'index peak list saved'
845
846	def ReadEXPPhase(self,filename):
847	import GSASIIspc as G2spc
848	import GSASIIlattice as G2lat
849	import GSASIIElem as G2el
850	file = open(filename, 'Ur')
851	Phase = {}
852	S = 1
853	Expr = [{},{},{},{},{},{},{},{},{}]
854	while S:
855	S = file.readline()
856	if 'EXPR NPHAS' in S[:12]:
857	Num = S[12:-1].count('0')
858	NPhas = S[12:-1].split()
859	if 'CRS' in S[:3]:
860	N = int(S[3:4])-1
861	Expr[N][S[:12]] = S[12:-1]
862	file.close()
863	PNames = []
864	for n,N in enumerate(NPhas):
865	if N != '0':
866	result = n
867	key = 'CRS'+str(n+1)+' PNAM'
868	PNames.append(Expr[n][key])
869	if Num < 8:
870	dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
871	try:
872	if dlg.ShowModal() == wx.ID_OK:
873	result = dlg.GetSelection()
874	finally:
875	dlg.Destroy()
876	EXPphase = Expr[result]
877	keyList = EXPphase.keys()
878	keyList.sort()
879	SGData = {}
880	if NPhas[result] == '1':
881	Ptype = 'nuclear'
882	elif NPhas[result] in ['2','3']:
883	Ptype = 'magnetic'
884	elif NPhas[result] == '4':
885	Ptype = 'macromolecular'
886	elif NPhas[result] == '10':
887	Ptype = 'Pawley'
888	for key in keyList:
889	if 'PNAM' in key:
890	PhaseName = EXPphase[key].strip()
891	elif 'ABC ' in key:
892	abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]
893	elif 'ANGLES' in key:
894	angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]
895	elif 'SG SYM' in key:
896	SpGrp = EXPphase[key][:15].strip()
897	E,SGData = G2spc.SpcGroup(SpGrp)
898	Atoms = []
899	if Ptype == 'nuclear':
900	for key in keyList:
901	if 'AT' in key:
902	if key[11:] == 'A':
903	S = EXPphase[key]
904	elif key[11:] == 'B':
905	S += EXPphase[key]
906	Atom = [S[50:58].strip(),S[:10].strip(),'',
907	float(S[10:20]),float(S[20:30]),float(S[30:40]),
908	float(S[40:50]),'',int(S[60:62]),S[130:131]]
909	if Atom[9] == 'I':
910	Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
911	elif Atom[9] == 'A':
912	Atom += [0.0,float(S[68:78]),float(S[78:88]),
913	float(S[88:98]),float(S[98:108]),
914	float(S[108:118]),float(S[118:128])]
915	XYZ = Atom[3:6]
916	Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
917	Atoms.append(Atom)
918	elif Ptype == 'macromolecular':
919	for key in keyList:
920	if 'AT' in key[6:8]:
921	S = EXPphase[key]
922	Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
923	S[46:51].strip(),S[:8].strip(),'',
924	float(S[16:24]),float(S[24:32]),float(S[32:40]),
925	float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
926	XYZ = Atom[6:9]
927	Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
928	Atoms.append(Atom)
929	Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
930	Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
931	'Cell':[False,]+abc+angles+[Volume,]}
932	Phase['Atoms'] = Atoms
933	Phase['Drawing'] = {}
934	Phase['Histograms'] = {}
935	return Phase
936
937	def ReadPDBPhase(filename):
938	import GSASIIspc as G2spc
939	import GSASIIlattice as G2lat
940	import GSASIIElem as G2el
941	EightPiSq = 8.math.pi*2
942	file = open(filename, 'Ur')
943	Phase = {}
944	Title = ''
945	Compnd = ''
946	Atoms = []
947	A = np.zeros(shape=(3,3))
948	S = file.readline()
949	while S:
950	Atom = []
951	if 'TITLE' in S[:5]:
952	Title = S[10:72].strip()
953	S = file.readline()
954	elif 'COMPND ' in S[:10]:
955	Compnd = S[10:72].strip()
956	S = file.readline()
957	elif 'CRYST' in S[:5]:
958	abc = S[7:34].split()
959	angles = S[34:55].split()
960	cell=[float(abc[0]),float(abc[1]),float(abc[2]),
961	float(angles[0]),float(angles[1]),float(angles[2])]
962	Volume = G2lat.calc_V(G2lat.cell2A(cell))
963	AA,AB = G2lat.cell2AB(cell)
964	SpGrp = S[55:65]
965	E,SGData = G2spc.SpcGroup(SpGrp)
966	S = file.readline()
967	elif 'SCALE' in S[:5]:
968	V = (S[10:41].split())
969	A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
970	S = file.readline()
971	elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
972	XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
973	XYZ = np.inner(AB,XYZ)
974	SytSym,Mult = G2spc.SytSym(XYZ,SGData)
975	Uiso = float(S[61:67])/EightPiSq
976	Type = S[12:14].upper()
977	if Type[0] in '123456789':
978	Type = Type[1:]
979	Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
980	S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
981	float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
982	S = file.readline()
983	if 'ANISOU' in S[:6]:
984	Uij = S[30:72].split()
985	Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
986	float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
987	Atom = Atom[:14]+Uij
988	Atom[12] = 'A'
989	S = file.readline()
990	Atoms.append(Atom)
991	else:
992	S = file.readline()
993	file.close()
994	if Title:
995	PhaseName = Title
996	elif Compnd:
997	PhaseName = Compnd
998	else:
999	PhaseName = 'None'
1000	Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1001	'Cell':[False,]+cell+[Volume,]}
1002	Phase['Atoms'] = Atoms
1003	Phase['Drawing'] = {}
1004	Phase['Histograms'] = {}
1005
1006	return Phase
1007
1008	def ReadCIFAtoms(self,data):
1009	print data

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

Download in other formats: