Context Navigation

source: trunk/GSASIIIO.py @ 59

Last change on this file since 59 was 56, checked in by vondreel, 15 years ago
major modification for images - not stored in project file. Now read as needed in plotting routine. Reduces .gpx file size from 10s of MB to 10s of kb.
File size: 31.1 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 GetImageData(imagefile,imageOnly=False):
391	ext = ospath.splitext(imagefile)[1]
392	Comments = []
393	if ext == '.tif':
394	Comments,Data,Size,Image = GetTifData(imagefile)
395	elif ext == '.img':
396	Comments,Data,Size,Image = GetImgData(imagefile)
397	Image[0][0] = 0
398	elif ext == '.mar3450' or ext == '.mar2300':
399	Comments,Data,Size,Image = GetMAR345Data(imagefile)
400	elif ext in ['.sum','.avg']:
401	Comments,Data,Size,Image = GetGEsumData(imagefile)
402	if imageOnly:
403	return Image
404	else:
405	return Comments,Data,Size,Image
406
407
408	def GetGEsumData(filename,imageOnly=False):
409	import array as ar
410	print 'Read GE sum file: ',filename
411	File = open(filename,'rb')
412	size = 2048
413	if '.sum' in filename:
414	head = ['GE detector sum data from APS 1-ID',]
415	if '.avg' in filename:
416	head = ['GE detector avg data from APS 1-ID',]
417	image = np.zeros(shape=(size,size),dtype=np.int32)
418	row = 0
419	pos = 0
420	while row < size:
421	File.seek(pos)
422	if '.sum' in filename:
423	line = ar.array('f',File.read(4*size))
424	pos += 4*size
425	elif '.avg' in filename:
426	line = ar.array('H',File.read(2*size))
427	pos += 2*size
428	image[row] = np.asarray(line)
429	row += 1
430	data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8]}
431	File.close()
432	if imageOnly:
433	return image
434	else:
435	return head,data,size,image
436
437	def GetImgData(filename,imageOnly=False):
438	import struct as st
439	import array as ar
440	print 'Read ADSC img file: ',filename
441	File = open(filename,'rb')
442	head = File.read(511)
443	lines = head.split('\n')
444	head = []
445	center = [0,0]
446	for line in lines[1:-2]:
447	line = line.strip()[:-1]
448	if line:
449	if 'SIZE1' in line:
450	size = int(line.split('=')[1])
451	elif 'WAVELENGTH' in line:
452	wave = float(line.split('=')[1])
453	elif 'BIN' in line:
454	if line.split('=')[1] == '2x2':
455	pixel=(102,102)
456	else:
457	pixel = (51,51)
458	elif 'DISTANCE' in line:
459	distance = float(line.split('=')[1])
460	elif 'CENTER_X' in line:
461	center[0] = float(line.split('=')[1])
462	elif 'CENTER_Y' in line:
463	center[1] = float(line.split('=')[1])
464	head.append(line)
465	data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
466	image = []
467	row = 0
468	pos = 512
469	image = np.zeros(shape=(size,size),dtype=np.int32)
470	while row < size:
471	File.seek(pos)
472	line = ar.array('H',File.read(2*size))
473	image[row] = np.asarray(line)
474	row += 1
475	pos += 2*size
476	File.close()
477	if imageOnly:
478	return image
479	else:
480	return lines[1:-2],data,size,image
481
482	def GetMAR345Data(filename,imageOnly=False):
483	import array as ar
484	import struct as st
485	try:
486	import pack_f as pf
487	except:
488	msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
489	caption="Import Error",
490	style=wx.ICON_ERROR \| wx.OK \| wx.STAY_ON_TOP)
491	msg.ShowModal()
492	return None,None,None,None
493
494	print 'Read Mar345 file: ',filename
495	File = open(filename,'rb')
496	head = File.read(4095)
497	numbers = st.unpack('<iiiiiiiiii',head[:40])
498	lines = head[128:].split('\n')
499	head = []
500	for line in lines:
501	line = line.strip()
502	if 'PIXEL' in line:
503	values = line.split()
504	pixel = (int(values[2]),int(values[4])) #in microns
505	elif 'WAVELENGTH' in line:
506	wave = float(line.split()[1])
507	elif 'DISTANCE' in line:
508	distance = float(line.split()[1]) #in mm
509	elif 'CENTER' in line:
510	values = line.split()
511	center = [float(values[2])/10.,float(values[4])/10.] #make in mm from pixels
512	if line:
513	head.append(line)
514	data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
515	for line in head:
516	if 'FORMAT' in line[0:6]:
517	items = line.split()
518	size = int(items[1])
519	pos = 4096
520	File.seek(pos)
521	line = File.read(8)
522	while 'CCP4' not in line: #get past overflow list for now
523	line = File.read(8)
524	pos += 8
525	pos += 37
526	File.seek(pos)
527	raw = File.read()
528	File.close()
529	image = np.zeros(shape=(size,size),dtype=np.int32)
530	image = pf.pack_f(len(raw),raw,size,image)
531	if imageOnly:
532	return image.T
533	else:
534	return head,data,size,image.T
535
536	def GetTifData(filename,imageOnly=False):
537	# only works for APS Perkin-Elmer detector data files in "TIFF" format that are readable by Fit2D
538	import struct as st
539	import array as ar
540	print 'Read APS PE-detector tiff file: ',filename
541	File = open(filename,'Ur')
542	dataType = 5
543	try:
544	Meta = open(filename+'.metadata','Ur')
545	head = Meta.readlines()
546	for line in head:
547	line = line.strip()
548	if 'dataType' in line:
549	dataType = int(line.split('=')[1])
550	Meta.close()
551	except IOError:
552	print 'no metadata file found - will try to read file anyway'
553	head = 'no metadata file found'
554	tag = File.read(3)
555	if tag != 'II*':
556	lines = ['not a APS PE-detector tiff file',]
557	return lines,0,0
558	size = st.unpack('<i',File.read(4))[0]
559	image = np.zeros(shape=(size,size),dtype=np.int32)
560	row = 0
561	pos = 8
562	while row < size:
563	File.seek(pos)
564	if dataType == 5:
565	line = ar.array('f',File.read(4*size))
566	else:
567	line = ar.array('l',File.read(4*size))
568	image[row] = np.asarray(line)
569	row += 1
570	pos += 4*size
571	data = {'pixelSize':(200,200),'wavelength':0.10,'distance':100.0,'center':[204.8,204.8]}
572	File.close()
573	if imageOnly:
574	return image
575	else:
576	return head,data,size,image
577
578
579	def ProjFileOpen(self):
580	file = open(self.GSASprojectfile,'rb')
581	print 'load from file: ',self.GSASprojectfile
582	wx.BeginBusyCursor()
583	try:
584	while True:
585	try:
586	data = cPickle.load(file)
587	except EOFError:
588	break
589	datum = data[0]
590	print 'load: ',datum[0]
591	if 'PWDR' in datum[0] and 'list' in str(type(datum[1][1][0])): #fix to convert old style list arrays to numpy arrays
592	X = datum[1][1]
593	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])]
594	datum[1] = [datum[1][0],X]
595	print 'powder data converted to numpy arrays'
596	if 'PKS' not in datum[0] and 'IMG' not in datum[0] and 'SNGL' not in datum[0]:
597	if datum[0] not in ['Notebook','Controls','Phases'] and 'PWDR' not in datum[0]: #temporary fix
598	datum[0] = 'PWDR '+datum[0]
599	Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
600	self.PatternTree.SetItemPyData(Id,datum[1])
601	for datus in data[1:]:
602	print ' load: ',datus[0]
603	if 'PWDR' in datum[0] and 'Instrument Parameters' in datus[0]:
604	if len(datus[1][0]) == 10 or len(datus[1][0]) == 12:
605	datus[1][0] += (0.0,) #add missing azimuthal angle
606	datus[1][1].append(0.0)
607	datus[1][2].append(0.0)
608	datus[1][3].append('Azimuth')
609	sub = self.PatternTree.AppendItem(Id,datus[0])
610	self.PatternTree.SetItemPyData(sub,datus[1])
611	if 'PWDR' in datum[0] and not G2gd.GetPatternTreeItemId(self,Id, 'Comments'):
612	print 'no comments - add comments'
613	sub = self.PatternTree.AppendItem(Id,'Comments')
614	self.PatternTree.SetItemPyData(sub,['no comments'])
615
616	file.close()
617	if not G2gd.GetPatternTreeItemId(self,self.root,'Notebook'):
618	sub = self.PatternTree.AppendItem(parent=self.root,text='Notebook')
619	self.PatternTree.SetItemPyData(sub,[''])
620	sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
621	self.PatternTree.SetItemPyData(sub,[0])
622
623	finally:
624	wx.EndBusyCursor()
625	print 'project load successful'
626	self.NewPlot = True
627
628	def ProjFileSave(self):
629	if not self.PatternTree.IsEmpty():
630	file = open(self.GSASprojectfile,'wb')
631	print 'save to file: ',self.GSASprojectfile
632	wx.BeginBusyCursor()
633	try:
634	item, cookie = self.PatternTree.GetFirstChild(self.root)
635	while item:
636	data = []
637	name = self.PatternTree.GetItemText(item)
638	print 'save: ',name
639	data.append([name,self.PatternTree.GetItemPyData(item)])
640	item2, cookie2 = self.PatternTree.GetFirstChild(item)
641	while item2:
642	name = self.PatternTree.GetItemText(item2)
643	print ' save: ',name
644	data.append([name,self.PatternTree.GetItemPyData(item2)])
645	item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)
646	item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
647	cPickle.dump(data,file,1)
648	file.close()
649	finally:
650	wx.EndBusyCursor()
651	print 'project save successful'
652
653	def powderFxyeSave(self,powderfile):
654	file = open(powderfile,'w')
655	print 'save powder pattern to file: ',powderfile
656	wx.BeginBusyCursor()
657	try:
658	x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
659	x = x*100.
660	XYW = zip(x,y,w)
661	for X,Y,W in XYW:
662	file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
663	file.close()
664	finally:
665	wx.EndBusyCursor()
666	print 'powder pattern file written'
667
668	def powderXyeSave(self,powderfile):
669	file = open(powderfile,'w')
670	print 'save powder pattern to file: ',powderfile
671	wx.BeginBusyCursor()
672	try:
673	x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
674	XYW = zip(x,y,w)
675	for X,Y,W in XYW:
676	file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
677	file.close()
678	finally:
679	wx.EndBusyCursor()
680	print 'powder pattern file written'
681
682	def PeakListSave(self,file,peaks):
683	print 'save peak list to file: ',self.peaklistfile
684	if not peaks:
685	dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
686	try:
687	result = dlg.ShowModal()
688	finally:
689	dlg.Destroy()
690	return
691	for peak in peaks:
692	file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
693	(peak[0],peak[2],peak[4],peak[6]))
694	print 'peak list saved'
695
696	def IndexPeakListSave(self,peaks):
697	file = open(self.peaklistfile,'wa')
698	print 'save index peak list to file: ',self.peaklistfile
699	wx.BeginBusyCursor()
700	try:
701	if not peaks:
702	dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
703	try:
704	result = dlg.ShowModal()
705	finally:
706	dlg.Destroy()
707	return
708	for peak in peaks:
709	file.write("%12.6f\n" % (peak[7]))
710	file.close()
711	finally:
712	wx.EndBusyCursor()
713	print 'index peak list saved'
714
715	def ReadEXPPhase(filename):
716	import GSASIIspc as G2spc
717	import GSASIIcomp as G2cmp
718	import GSASIIElem as G2el
719	file = open(filename, 'Ur')
720	Phase = {}
721	S = 1
722	Expr = [{},{},{},{},{},{},{},{},{}]
723	while S:
724	S = file.readline()
725	if 'EXPR NPHAS' in S[:12]:
726	Num = S[12:-1].count('0')
727	NPhas = S[12:-1].split()
728	if 'CRS' in S[:3]:
729	N = int(S[3:4])-1
730	Expr[N][S[:12]] = S[12:-1]
731	file.close()
732	PNames = []
733	for n,N in enumerate(NPhas):
734	if N != '0':
735	result = n
736	key = 'CRS'+str(n+1)+' PNAM'
737	PNames.append(Expr[n][key])
738	if Num < 8:
739	dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
740	try:
741	if dlg.ShowModal() == wx.ID_OK:
742	result = dlg.GetSelection()
743	finally:
744	dlg.Destroy()
745	EXPphase = Expr[result]
746	keyList = EXPphase.keys()
747	keyList.sort()
748	SGData = {}
749	if NPhas[result] == '1':
750	Ptype = 'nuclear'
751	elif NPhas[result] in ['2','3']:
752	Ptype = 'magnetic'
753	elif NPhas[result] == '4':
754	Ptype = 'macromolecular'
755	elif NPhas[result] == '10':
756	Ptype = 'Pawley'
757	for key in keyList:
758	if 'PNAM' in key:
759	PhaseName = EXPphase[key].strip()
760	elif 'ABC ' in key:
761	abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]
762	elif 'ANGLES' in key:
763	angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]
764	elif 'SG SYM' in key:
765	SpGrp = EXPphase[key][:15].strip()
766	E,SGData = G2spc.SpcGroup(SpGrp)
767	Atoms = []
768	if Ptype == 'nuclear':
769	for key in keyList:
770	if 'AT' in key:
771	if key[11:] == 'A':
772	S = EXPphase[key]
773	elif key[11:] == 'B':
774	S += EXPphase[key]
775	Atom = [S[50:58].strip(),S[:10].strip(),'',
776	float(S[10:20]),float(S[20:30]),float(S[30:40]),
777	float(S[40:50]),'',int(S[60:62]),S[130:131]]
778	if Atom[9] == 'I':
779	Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
780	elif Atom[9] == 'A':
781	Atom += [0.0,float(S[68:78]),float(S[78:88]),
782	float(S[88:98]),float(S[98:108]),
783	float(S[108:118]),float(S[118:128])]
784	XYZ = Atom[3:6]
785	Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
786	Atoms.append(Atom)
787	elif Ptype == 'macromolecular':
788	for key in keyList:
789	if 'AT' in key[6:8]:
790	S = EXPphase[key]
791	Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
792	S[46:51].strip(),S[:8].strip(),'',
793	float(S[16:24]),float(S[24:32]),float(S[32:40]),
794	float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
795	XYZ = Atom[6:9]
796	Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
797	Atoms.append(Atom)
798	Volume = G2cmp.calc_V(G2cmp.cell2A(abc+angles))
799	Phase['General'] = [PhaseName,Ptype,SGData,[False,]+abc+angles+[Volume,],[False,1.0],
800	0,0,0,0,0]
801	Phase['Atoms'] = Atoms
802	return Phase
803
804	def ReadPDBPhase(filename):
805	import GSASIIspc as G2spc
806	import GSASIIcomp as G2cmp
807	import GSASIIElem as G2el
808	EightPiSq = 8.math.pi*2
809	file = open(filename, 'Ur')
810	Phase = {}
811	Title = ''
812	Compnd = ''
813	Atoms = []
814	A = np.zeros(shape=(3,3))
815	S = file.readline()
816	while S:
817	Atom = []
818	if 'TITLE' in S[:5]:
819	Title = S[10:72].strip()
820	S = file.readline()
821	elif 'COMPND ' in S[:10]:
822	Compnd = S[10:72].strip()
823	S = file.readline()
824	elif 'CRYST' in S[:5]:
825	abc = S[7:34].split()
826	angles = S[34:55].split()
827	cell=[float(abc[0]),float(abc[1]),float(abc[2]),
828	float(angles[0]),float(angles[1]),float(angles[2])]
829	Volume = G2cmp.calc_V(G2cmp.cell2A(cell))
830	AA,AB = G2cmp.cell2AB(cell)
831	SpGrp = S[55:65]
832	E,SGData = G2spc.SpcGroup(SpGrp)
833	S = file.readline()
834	elif 'SCALE' in S[:5]:
835	V = (S[10:41].split())
836	A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
837	S = file.readline()
838	elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
839	XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
840	XYZ = np.sum(AA*XYZ,axis=1)
841	SytSym,Mult = G2spc.SytSym(XYZ,SGData)
842	Uiso = float(S[61:67])/EightPiSq
843	Type = S[12:14].upper()
844	if Type[0] in '123456789':
845	Type = Type[1:]
846	Atom = [int(S[22:27]),S[17:20].upper(),S[20:22],
847	S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
848	float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
849	S = file.readline()
850	if 'ANISOU' in S[:6]:
851	Uij = S[30:72].split()
852	Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
853	float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
854	Atom = Atom[:14]+Uij
855	Atom[12] = 'A'
856	S = file.readline()
857	Atoms.append(Atom)
858	else:
859	S = file.readline()
860	file.close()
861	if Title:
862	PhaseName = Title
863	elif Compnd:
864	PhaseName = Compnd
865	else:
866	PhaseName = 'None'
867	Phase['General'] = [PhaseName,'macromolecular',SGData,[False,]+cell+[Volume,],[False,1.0],
868	0,0,0,0,0]
869	Phase['Atoms'] = Atoms
870
871	return Phase
872
873	def ReadCIFAtoms(self,data):
874	print data

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

Download in other formats: