Context Navigation

source: trunk/GSASIIIO.py @ 223

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

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

Download in other formats: