Context Navigation

source: trunk/GSASIIIO.py @ 55

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

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

Download in other formats: