Context Navigation

source: trunk/GSASIIIO.py @ 615

Last change on this file since 615 was 615, checked in by toby, 12 years ago
add powder CIF importer, minor other updates to readers
Property svn:keywords set to `Date Author Revision URL Id`
File size: 64.9 KB

Line
1	# -- coding: utf-8 --
2	"""GSASIIIO: functions for IO of data
3	Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
4	"""
5	########### SVN repository information ###################
6	# $Date: 2012-05-16 03:54:23 +0000 (Wed, 16 May 2012) $
7	# $Author: toby $
8	# $Revision: 615 $
9	# $URL: trunk/GSASIIIO.py $
10	# $Id: GSASIIIO.py 615 2012-05-16 03:54:23Z toby $
11	########### SVN repository information ###################
12	import wx
13	import math
14	import numpy as np
15	import cPickle
16	import sys
17	import random as ran
18	import GSASIIpath
19	import GSASIIgrid as G2gd
20	import GSASIIspc as G2spc
21	import GSASIIlattice as G2lat
22	import GSASIIpwdGUI as G2pdG
23	import GSASIIElem as G2el
24	import os
25	import os.path as ospath
26
27	def sfloat(S):
28	if S.strip():
29	return float(S)
30	else:
31	return 0.0
32
33	def sint(S):
34	if S.strip():
35	return int(S)
36	else:
37	return 0
38
39	def FileDlgFixExt(dlg,file):
40	#this is needed to fix a problem in linux wx.FileDialog
41	ext = dlg.GetWildcard().split('\|')[2dlg.GetFilterIndex()+1].strip('')
42	if ext not in file:
43	file += ext
44	return file
45
46	# to be removed
47	def SelectPowderData(G2frame, filename):
48	"""Selects banks of data from a filename of any GSAS powder data format
49	Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
50	Returns - a list of banks to be read; each entry in list is a tuple containing:
51	filename: same as input filename
52	Pos: position for start of data; record just after BANK record
53	Bank: the BANK record
54	"""
55	File = open(filename,'Ur')
56	Title = '''
57	First line of this file:
58	'''+File.readline()
59	dlg = wx.MessageDialog(G2frame, Title, 'Is this the file you want?',
60	wx.YES_NO \| wx.ICON_QUESTION)
61	try:
62	result = dlg.ShowModal()
63	finally:
64	dlg.Destroy()
65	if result == wx.ID_NO: return (0,0)
66	Temperature = 300
67
68	if '.xye' in filename: #Topas style xye file (e.g. 2-th, I, sig) - no iparm file/no BANK record
69	dlg = wx.MessageDialog(G2frame,'''Is this laboratory Cu Ka1/Ka2 data?
70	(No = 0.6A wavelength synchrotron data)
71	Change wavelength in Instrument Parameters if needed''','Data type?',
72	wx.YES_NO \| wx.ICON_QUESTION)
73	try:
74	result = dlg.ShowModal()
75	finally:
76	dlg.Destroy()
77	print result
78	if result == wx.ID_YES:
79	Iparm = {} #Assume CuKa lab data
80	Iparm['INS HTYPE '] = 'PXC '
81	Iparm['INS 1 ICONS'] = ' 1.540500 1.544300 0.0 0 0.7 0 0.5 '
82	Iparm['INS 1PRCF1 '] = ' 3 8 0.01 '
83	Iparm['INS 1PRCF11'] = ' 2.000000E+00 -2.000000E+00 5.000000E+00 0.000000E+00 '
84	Iparm['INS 1PRCF12'] = ' 0.000000E+00 0.000000E+00 0.100000E-01 0.100000E-01 '
85	else:
86	Iparm = {} #Assume 0.6A synchrotron data
87	Iparm['INS HTYPE '] = 'PXC '
88	Iparm['INS 1 ICONS'] = ' 0.600000 0.000000 0.0 0 0.99 0 0.5 '
89	Iparm['INS 1PRCF1 '] = ' 3 8 0.01 '
90	Iparm['INS 1PRCF11'] = ' 1.000000E+00 -1.000000E+00 0.300000E+00 0.000000E+00 '
91	Iparm['INS 1PRCF12'] = ' 0.000000E+00 0.000000E+00 0.100000E-02 0.100000E-02 '
92
93
94	else: #GSAS style fxye or fxy file (e.g. 100*2-th, I, sig)
95	G2frame.IparmName = GetInstrumentFile(G2frame,filename)
96	if G2frame.IparmName:
97	Iparm = GetInstrumentData(G2frame.IparmName)
98	else:
99	Iparm = {} #Assume CuKa lab data if no iparm file
100	Iparm['INS HTYPE '] = 'PXC '
101	Iparm['INS 1 ICONS'] = ' 1.540500 1.544300 0.0 0 0.7 0 0.5 '
102	Iparm['INS 1PRCF1 '] = ' 3 8 0.01 '
103	Iparm['INS 1PRCF11'] = ' 2.000000E+00 -2.000000E+00 5.000000E+00 0.000000E+00 '
104	Iparm['INS 1PRCF12'] = ' 0.000000E+00 0.000000E+00 0.100000E-01 0.100000E-01 '
105	S = 1
106	Banks = []
107	Pos = []
108	FoundData = []
109	Comments = []
110	wx.BeginBusyCursor()
111	try:
112	while S:
113	S = File.readline()
114	if S[:1] != '#':
115	if S[:4] == 'BANK':
116	Banks.append(S)
117	Pos.append(File.tell())
118	elif '.xye' in filename: #No BANK in a xye file
119	Banks.append('BANK 1 XYE')
120	Pos.append(File.tell())
121	break
122	else:
123	Comments.append(S[:-1])
124	if 'Temp' in S.split('=')[0]:
125	Temperature = float(S.split('=')[1])
126	File.close()
127	finally:
128	wx.EndBusyCursor()
129	if Comments:
130	print 'Comments on file:'
131	for Comment in Comments: print Comment
132	if Banks:
133	result = [0]
134	if len(Banks) >= 2:
135	dlg = wx.MultiChoiceDialog(G2frame, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
136	try:
137	if dlg.ShowModal() == wx.ID_OK:
138	result = dlg.GetSelections()
139	else:
140	result = []
141	finally:
142	dlg.Destroy()
143	for i in result:
144	FoundData.append((filename,Pos[i],Banks[i]))
145	else:
146	dlg = wx.MessageDialog(G2frame, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK \| wx.ICON_ERROR)
147	try:
148	result = dlg.ShowModal()
149	finally:
150	dlg.Destroy()
151	return FoundData,Iparm,Comments,Temperature
152
153	# to be removed
154	def GetInstrumentFile(G2frame,filename):
155	import os.path as op
156	dlg = wx.FileDialog(G2frame,'Choose an instrument file','.', '', 'GSAS iparm file (.prm)\|.prm\|All files(.)\|.',
157	wx.OPEN\|wx.CHANGE_DIR)
158	Tname = filename[:filename.index('.')]+'.prm'
159	if op.exists(Tname):
160	G2frame.IparmName = Tname
161	if G2frame.IparmName: dlg.SetFilename(G2frame.IparmName)
162	filename = ''
163	try:
164	if dlg.ShowModal() == wx.ID_OK:
165	filename = dlg.GetPath()
166	finally:
167	dlg.Destroy()
168	return filename
169
170	# to be removed
171	def GetInstrumentData(IparmName):
172	file = open(IparmName, 'Ur')
173	S = 1
174	Iparm = {}
175	while S:
176	S = file.readline()
177	Iparm[S[:12]] = S[12:-1]
178	return Iparm
179
180	def GetPowderPeaks(fileName):
181	sind = lambda x: math.sin(x*math.pi/180.)
182	asind = lambda x: 180.*math.asin(x)/math.pi
183	Cuka = 1.54052
184	File = open(fileName,'Ur')
185	Comments = []
186	peaks = []
187	S = File.readline()
188	while S:
189	if S[:1] == '#':
190	Comments.append(S[:-1])
191	else:
192	item = S.split()
193	if len(item) == 1:
194	peaks.append([float(item[0]),1.0])
195	elif len(item) > 1:
196	peaks.append([float(item[0]),float(item[0])])
197	S = File.readline()
198	File.close()
199	if Comments:
200	print 'Comments on file:'
201	for Comment in Comments: print Comment
202	Peaks = []
203	if peaks[0][0] > peaks[-1][0]: # d-spacings - assume CuKa
204	for peak in peaks:
205	dsp = peak[0]
206	sth = Cuka/(2.0*dsp)
207	if sth < 1.0:
208	tth = 2.0*asind(sth)
209	else:
210	break
211	Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
212	else: #2-thetas - assume Cuka (for now)
213	for peak in peaks:
214	tth = peak[0]
215	dsp = Cuka/(2.0*sind(tth/2.0))
216	Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
217	return Comments,Peaks
218
219	def GetPawleyPeaks(filename):
220	rt2ln2x2 = 2.35482
221	File = open(filename,'Ur')
222	PawleyPeaks = []
223	S = File.readline() #skip header
224	S = File.readline()
225	item = S.split()
226	while S:
227	h,k,l = int(item[0]),int(item[1]),int(item[2])
228	mult = int(item[3])
229	tth = float(item[5])
230	sig = float(item[6])/rt2ln2x2
231	Iobs = float(item[7])*mult
232	PawleyPeaks.append([h,k,l,mult,tth,False,Iobs,0.0])
233	S = File.readline()
234	item = S.split()
235	if item[3] == '-100.0000': #find trailer
236	break
237	File.close()
238	return PawleyPeaks
239
240	# this will be removed eventually
241	def GetHKLData(filename):
242	print 'Reading: '+filename
243	File = open(filename,'Ur')
244	HKLref = []
245	HKLmin = [1000,1000,1000]
246	HKLmax = [0,0,0]
247	FoMax = 0
248	ifFc = False
249	S = File.readline()
250	while '#' in S[0]: #get past comments if any
251	S = File.readline()
252	if '_' in S: #cif style .hkl file
253	while 'loop_' not in S: #skip preliminaries if any - can't have 'loop_' in them!
254	S = File.readline()
255	S = File.readline() #get past 'loop_' line
256	pos = 0
257	hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
258	while S:
259	if '_' in S:
260	if 'index_h' in S:
261	hpos = pos
262	elif 'index_k' in S:
263	kpos = pos
264	elif 'index_l' in S:
265	lpos = pos
266	elif 'F_squared_meas' in S:
267	Fosqpos = pos
268	elif 'F_squared_calc' in S:
269	Fcsqpos = pos
270	elif 'F_squared_sigma' in S:
271	sigpos = pos
272	pos += 1
273	else:
274	data = S.split()
275	if data: #avoid blank lines
276	HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
277	h,k,l = HKL
278	Fosq = float(data[Fosqpos])
279	if sigpos != -1:
280	sigFosq = float(data[sigpos])
281	else:
282	sigFosq = 1.
283	if Fcsqpos != -1:
284	Fcsq = float(data[Fcsqpos])
285	if Fcsq:
286	ifFc = True
287	else:
288	Fcsq = 0.
289
290	HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
291	HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
292	FoMax = max(FoMax,Fosq)
293	HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0]) #room for Fcp, Fcpp & phase
294	S = File.readline()
295	else: #dumb h,k,l,Fo,sigFo .hkl file
296	while S:
297	h,k,l,Fo,sigFo = S.split()
298	HKL = np.array([int(h),int(k),int(l)])
299	h,k,l = HKL
300	Fo = float(Fo)
301	sigFo = float(sigFo)
302	HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
303	HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
304	FoMax = max(FoMax,Fo)
305	HKLref.append([HKL,Fo*2,2.Fo*sigFo,0,0,0,0]) #room for Fc, Fcp, Fcpp & phase
306	S = File.readline()
307	File.close()
308	return HKLref,HKLmin,HKLmax,FoMax,ifFc
309
310	# to be removed
311	def GetPowderData(filename,Pos,Bank,DataType):
312	'''Reads one BANK of data from GSAS raw powder data file
313	input:
314	filename: GSAS raw powder file dataname
315	Pos: start of data in file just after BANK record
316	Bank: the BANK record
317	DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
318	returns: list [x,y,e,yc,yb]
319	x: np.array of x-axis values
320	y: np.array of powder pattern intensities
321	w: np.array of w=sig(intensity)^2 values
322	yc: np.array of calc. intensities (zero)
323	yb: np.array of calc. background (zero)
324	yd: np.array of obs-calc profiles
325	'''
326	print 'Reading: '+filename
327	print 'Bank: '+Bank[:-1]
328	if 'FXYE' in Bank:
329	return GetFXYEdata(filename,Pos,Bank,DataType)
330	elif ' XYE' in Bank:
331	return GetXYEdata(filename,Pos,Bank,DataType)
332	elif 'FXY' in Bank:
333	return GetFXYdata(filename,Pos,Bank,DataType)
334	elif 'ESD' in Bank:
335	return GetESDdata(filename,Pos,Bank,DataType)
336	elif 'STD' in Bank:
337	return GetSTDdata(filename,Pos,Bank,DataType)
338	else:
339	return GetSTDdata(filename,Pos,Bank,DataType)
340	return []
341
342	# to be removed
343	def GetFXYEdata(filename,Pos,Bank,DataType):
344	File = open(filename,'Ur')
345	File.seek(Pos)
346	x = []
347	y = []
348	w = []
349	S = File.readline()
350	while S and S[:4] != 'BANK':
351	vals = S.split()
352	if DataType[2] == 'C':
353	x.append(float(vals[0])/100.) #CW: from centidegrees to degrees
354	elif DataType[2] == 'T':
355	x.append(float(vals[0])/1000.0) #TOF: from musec to millisec
356	f = float(vals[1])
357	if f <= 0.0:
358	y.append(0.0)
359	w.append(1.0)
360	else:
361	y.append(float(vals[1]))
362	w.append(1.0/float(vals[2])**2)
363	S = File.readline()
364	File.close()
365	N = len(x)
366	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
367
368	# to be removed
369	def GetXYEdata(filename,Pos,Bank,DataType):
370	File = open(filename,'Ur')
371	File.seek(Pos)
372	x = []
373	y = []
374	w = []
375	S = File.readline()
376	while S:
377	vals = S.split()
378	try:
379	x.append(float(vals[0]))
380	f = float(vals[1])
381	if f <= 0.0:
382	y.append(0.0)
383	w.append(1.0)
384	else:
385	y.append(float(vals[1]))
386	w.append(1.0/float(vals[2])**2)
387	S = File.readline()
388	except ValueError:
389	break
390	File.close()
391	N = len(x)
392	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
393
394
395	# to be removed
396	def GetFXYdata(filename,Pos,Bank,DataType):
397	File = open(filename,'Ur')
398	File.seek(Pos)
399	x = []
400	y = []
401	w = []
402	S = File.readline()
403	while S and S[:4] != 'BANK':
404	vals = S.split()
405	if DataType[2] == 'C':
406	x.append(float(vals[0])/100.) #CW: from centidegrees to degrees
407	elif DataType[2] == 'T':
408	x.append(float(vals[0])/1000.0) #TOF: from musec to millisec
409	f = float(vals[1])
410	if f > 0.0:
411	y.append(f)
412	w.append(1.0/f)
413	else:
414	y.append(0.0)
415	w.append(1.0)
416	S = File.readline()
417	File.close()
418	N = len(x)
419	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
420
421	# to be removed
422	def GetESDdata(filename,Pos,Bank,DataType):
423	File = open(filename,'Ur')
424	cons = Bank.split()
425	if DataType[2] == 'C':
426	start = float(cons[5])/100.0 #CW: from centidegrees to degrees
427	step = float(cons[6])/100.0
428	elif DataType[2] == 'T':
429	start = float(cons[5])/1000.0 #TOF: from musec to millisec
430	step = float(cons[6])/1000.0
431	File.seek(Pos)
432	x = []
433	y = []
434	w = []
435	S = File.readline()
436	j = 0
437	while S and S[:4] != 'BANK':
438	for i in range(0,80,16):
439	xi = start+step*j
440	yi = sfloat(S[i:i+8])
441	ei = sfloat(S[i+8:i+16])
442	x.append(xi)
443	if yi > 0.0:
444	y.append(yi)
445	w.append(1.0/ei**2)
446	else:
447	y.append(0.0)
448	w.append(1.0)
449	j += 1
450	S = File.readline()
451	File.close()
452	N = len(x)
453	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
454
455	# to be removed
456	def GetSTDdata(filename,Pos,Bank,DataType):
457	File = open(filename,'Ur')
458	cons = Bank.split()
459	Nch = int(cons[2])
460	if DataType[2] == 'C':
461	start = float(cons[5])/100.0 #CW: from centidegrees to degrees
462	step = float(cons[6])/100.0
463	elif DataType[2] == 'T':
464	start = float(cons[5])/1000.0 #TOF: from musec to millisec - not likely!
465	step = float(cons[6])/1000.0
466	File.seek(Pos)
467	x = []
468	y = []
469	w = []
470	S = File.readline()
471	j = 0
472	while S and S[:4] != 'BANK':
473	for i in range(0,80,8):
474	xi = start+step*j
475	ni = max(sint(S[i:i+2]),1)
476	yi = max(sfloat(S[i+2:i+8]),0.0)
477	if yi:
478	vi = yi/ni
479	else:
480	yi = 0.0
481	vi = 1.0
482	j += 1
483	if j < Nch:
484	x.append(xi)
485	y.append(yi)
486	w.append(1.0/vi)
487	S = File.readline()
488	File.close()
489	N = len(x)
490	return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
491
492	def CheckImageFile(G2frame,imagefile):
493	if not ospath.exists(imagefile):
494	dlg = wx.FileDialog(G2frame, 'Bad image file name; choose name', '.', '',\
495	'Any image file (.tif;.tiff;.mar;.avg;.sum;*.img)\
496	\|.tif;.tiff;.mar;.avg;.sum;*.img\|\
497	Any detector tif (.tif;.tiff)\|.tif;.tiff\|\
498	MAR file (.mar)\|.mar\|\
499	GE Image (.avg;.sum)\|.avg;.sum\|\
500	ADSC Image (.img)\|.img\|\
501	All files (.)\|.',wx.OPEN\|wx.CHANGE_DIR)
502	try:
503	dlg.SetFilename(''+ospath.split(imagefile)[1])
504	if dlg.ShowModal() == wx.ID_OK:
505	imagefile = dlg.GetPath()
506	else:
507	imagefile = False
508	finally:
509	dlg.Destroy()
510	return imagefile
511
512	def GetImageData(G2frame,imagefile,imageOnly=False):
513	ext = ospath.splitext(imagefile)[1]
514	Comments = []
515	if ext == '.tif' or ext == '.tiff':
516	Comments,Data,Npix,Image = GetTifData(imagefile)
517	elif ext == '.img':
518	Comments,Data,Npix,Image = GetImgData(imagefile)
519	Image[0][0] = 0
520	elif ext == '.mar3450' or ext == '.mar2300':
521	Comments,Data,Npix,Image = GetMAR345Data(imagefile)
522	elif ext in ['.sum','.avg','']:
523	Comments,Data,Npix,Image = GetGEsumData(imagefile)
524	elif ext == '.G2img':
525	Comments,Data,Npix,Image = GetG2Image(imagefile)
526	if imageOnly:
527	return Image
528	else:
529	return Comments,Data,Npix,Image
530
531	def PutG2Image(filename,Comments,Data,Npix,image):
532	File = open(filename,'wb')
533	cPickle.dump([Comments,Data,Npix,image],File,1)
534	File.close()
535	return
536
537	def GetG2Image(filename):
538	File = open(filename,'rb')
539	Comments,Data,Npix,image = cPickle.load(File)
540	File.close()
541	return Comments,Data,Npix,image
542
543	def GetGEsumData(filename,imageOnly=False):
544	import struct as st
545	import array as ar
546	if not imageOnly:
547	print 'Read GE sum file: ',filename
548	File = open(filename,'rb')
549	if '.sum' in filename:
550	head = ['GE detector sum data from APS 1-ID',]
551	sizexy = [2048,2048]
552	elif '.avg' in filename:
553	head = ['GE detector avg data from APS 1-ID',]
554	sizexy = [2048,2048]
555	else:
556	head = ['GE detector raw data from APS 1-ID',]
557	File.seek(18)
558	size,nframes = st.unpack('<ih',File.read(6))
559	sizexy = [2048,2048]
560	pos = 8192
561	File.seek(pos)
562	Npix = sizexy[0]*sizexy[1]
563	if '.sum' in filename:
564	image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
565	elif '.avg' in filename:
566	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
567	else:
568	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
569	while nframes > 1:
570	image += np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
571	nframes -= 1
572	image = np.reshape(image,(sizexy[1],sizexy[0]))
573	data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':sizexy}
574	File.close()
575	if imageOnly:
576	return image
577	else:
578	return head,data,Npix,image
579
580	def GetImgData(filename,imageOnly=False):
581	import struct as st
582	import array as ar
583	if not imageOnly:
584	print 'Read ADSC img file: ',filename
585	File = open(filename,'rb')
586	head = File.read(511)
587	lines = head.split('\n')
588	head = []
589	center = [0,0]
590	for line in lines[1:-2]:
591	line = line.strip()[:-1]
592	if line:
593	if 'SIZE1' in line:
594	size = int(line.split('=')[1])
595	Npix = size*size
596	elif 'WAVELENGTH' in line:
597	wave = float(line.split('=')[1])
598	elif 'BIN' in line:
599	if line.split('=')[1] == '2x2':
600	pixel=(102,102)
601	else:
602	pixel = (51,51)
603	elif 'DISTANCE' in line:
604	distance = float(line.split('=')[1])
605	elif 'CENTER_X' in line:
606	center[0] = float(line.split('=')[1])
607	elif 'CENTER_Y' in line:
608	center[1] = float(line.split('=')[1])
609	head.append(line)
610	data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
611	image = []
612	row = 0
613	pos = 512
614	File.seek(pos)
615	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
616	image = np.reshape(image,(sizexy[1],sizexy[0]))
617	# image = np.zeros(shape=(size,size),dtype=np.int32)
618	# while row < size:
619	# File.seek(pos)
620	# line = ar.array('H',File.read(2*size))
621	# image[row] = np.asarray(line)
622	# row += 1
623	# pos += 2*size
624	File.close()
625	if imageOnly:
626	return image
627	else:
628	return lines[1:-2],data,Npix,image
629
630	def GetMAR345Data(filename,imageOnly=False):
631	import array as ar
632	import struct as st
633	try:
634	import pack_f as pf
635	except:
636	msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
637	caption="Import Error",
638	style=wx.ICON_ERROR \| wx.OK \| wx.STAY_ON_TOP)
639	msg.ShowModal()
640	return None,None,None,None
641
642	if not imageOnly:
643	print 'Read Mar345 file: ',filename
644	File = open(filename,'rb')
645	head = File.read(4095)
646	numbers = st.unpack('<iiiiiiiiii',head[:40])
647	lines = head[128:].split('\n')
648	head = []
649	for line in lines:
650	line = line.strip()
651	if 'PIXEL' in line:
652	values = line.split()
653	pixel = (int(values[2]),int(values[4])) #in microns
654	elif 'WAVELENGTH' in line:
655	wave = float(line.split()[1])
656	elif 'DISTANCE' in line:
657	distance = float(line.split()[1]) #in mm
658	elif 'CENTER' in line:
659	values = line.split()
660	center = [float(values[2])/10.,float(values[4])/10.] #make in mm from pixels
661	if line:
662	head.append(line)
663	data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
664	for line in head:
665	if 'FORMAT' in line[0:6]:
666	items = line.split()
667	size = int(items[1])
668	Npix = size*size
669	pos = 4096
670	data['size'] = [size,size]
671	File.seek(pos)
672	line = File.read(8)
673	while 'CCP4' not in line: #get past overflow list for now
674	line = File.read(8)
675	pos += 8
676	pos += 37
677	File.seek(pos)
678	raw = File.read()
679	File.close()
680	image = np.zeros(shape=(size,size),dtype=np.int32)
681	image = pf.pack_f(len(raw),raw,size,image)
682	if imageOnly:
683	return image.T #transpose to get it right way around
684	else:
685	return head,data,Npix,image.T
686
687	def GetTifData(filename,imageOnly=False):
688	import struct as st
689	import array as ar
690	File = open(filename,'rb')
691	dataType = 5
692	try:
693	Meta = open(filename+'.metadata','Ur')
694	head = Meta.readlines()
695	for line in head:
696	line = line.strip()
697	if 'dataType=' in line:
698	dataType = int(line.split('=')[1])
699	Meta.close()
700	except IOError:
701	print 'no metadata file found - will try to read file anyway'
702	head = ['no metadata file found',]
703
704	tag = File.read(2)
705	byteOrd = '<'
706	if tag == 'II' and int(st.unpack('<h',File.read(2))[0]) == 42: #little endian
707	IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
708	elif tag == 'MM' and int(st.unpack('>h',File.read(2))[0]) == 42: #big endian
709	byteOrd = '>'
710	IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
711	else:
712	lines = ['not a detector tiff file',]
713	return lines,0,0,0
714	File.seek(IFD) #get number of directory entries
715	NED = int(st.unpack(byteOrd+'h',File.read(2))[0])
716	IFD = {}
717	for ied in range(NED):
718	Tag,Type = st.unpack(byteOrd+'Hh',File.read(4))
719	nVal = st.unpack(byteOrd+'i',File.read(4))[0]
720	if Type == 1:
721	Value = st.unpack(byteOrd+nVal*'b',File.read(nVal))
722	elif Type == 2:
723	Value = st.unpack(byteOrd+'i',File.read(4))
724	elif Type == 3:
725	Value = st.unpack(byteOrd+nVal'h',File.read(nVal2))
726	x = st.unpack(byteOrd+nVal'h',File.read(nVal2))
727	elif Type == 4:
728	Value = st.unpack(byteOrd+nVal'i',File.read(nVal4))
729	elif Type == 5:
730	Value = st.unpack(byteOrd+nVal'i',File.read(nVal4))
731	elif Type == 11:
732	Value = st.unpack(byteOrd+nVal'f',File.read(nVal4))
733	IFD[Tag] = [Type,nVal,Value]
734	# for key in IFD:
735	# print key,IFD[key]
736	sizexy = [IFD[256][2][0],IFD[257][2][0]]
737	[nx,ny] = sizexy
738	Npix = nx*ny
739	if 272 in IFD:
740	ifd = IFD[272]
741	File.seek(ifd[2][0])
742	S = File.read(ifd[1])
743	if 'PILATUS' in S:
744	tifType = 'Pilatus'
745	dataType = 0
746	pixy = (172,172)
747	File.seek(4096)
748	if not imageOnly:
749	print 'Read Pilatus tiff file: ',filename
750	image = ar.array('L',File.read(4*Npix))
751	image = np.array(np.asarray(image),dtype=np.int32)
752	elif 262 in IFD and IFD[262][2][0] > 4:
753	tifType = 'DND'
754	pixy = (158,158)
755	File.seek(512)
756	if not imageOnly:
757	print 'Read DND SAX/WAX-detector tiff file: ',filename
758	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
759	elif sizexy == [1536,1536]:
760	tifType = 'APS Gold'
761	pixy = (150,150)
762	File.seek(64)
763	if not imageOnly:
764	print 'Read Gold tiff file:',filename
765	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
766	elif sizexy == [2048,2048] or sizexy == [1024,1024]:
767	if IFD[273][2][0] == 8:
768	if IFD[258][2][0] == 32:
769	tifType = 'PE'
770	pixy = (200,200)
771	File.seek(8)
772	if not imageOnly:
773	print 'Read APS PE-detector tiff file: ',filename
774	if dataType == 5:
775	image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
776	else:
777	image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
778	elif IFD[273][2][0] == 4096:
779	tifType = 'MAR'
780	pixy = (158,158)
781	File.seek(4096)
782	if not imageOnly:
783	print 'Read MAR CCD tiff file: ',filename
784	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
785	elif sizexy == [4096,4096]:
786	if IFD[273][2][0] == 8:
787	if IFD[258][2][0] == 16:
788	tifType = 'scanCCD'
789	pixy = (9,9)
790	File.seek(8)
791	if not imageOnly:
792	print 'Read APS scanCCD tiff file: ',filename
793	image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
794	# elif sizexy == [960,960]:
795	# tiftype = 'PE-BE'
796	# pixy = (200,200)
797	# File.seek(8)
798	# if not imageOnly:
799	# print 'Read Gold tiff file:',filename
800	# image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
801
802	else:
803	lines = ['not a known detector tiff file',]
804	return lines,0,0,0
805
806	image = np.reshape(image,(sizexy[1],sizexy[0]))
807	center = [pixy[0]sizexy[0]/2000,pixy[1]sizexy[1]/2000]
808	data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
809	File.close()
810	if imageOnly:
811	return image
812	else:
813	return head,data,Npix,image
814
815	def ProjFileOpen(G2frame):
816	file = open(G2frame.GSASprojectfile,'rb')
817	print 'load from file: ',G2frame.GSASprojectfile
818	wx.BeginBusyCursor()
819	try:
820	while True:
821	try:
822	data = cPickle.load(file)
823	except EOFError:
824	break
825	datum = data[0]
826
827	Id = G2frame.PatternTree.AppendItem(parent=G2frame.root,text=datum[0])
828	if 'PWDR' in datum[0]:
829	G2frame.PatternTree.SetItemPyData(Id,datum[1][:3]) #temp. trim off junk
830	else:
831	G2frame.PatternTree.SetItemPyData(Id,datum[1])
832	for datus in data[1:]:
833	sub = G2frame.PatternTree.AppendItem(Id,datus[0])
834	G2frame.PatternTree.SetItemPyData(sub,datus[1])
835	if 'IMG' in datum[0]: #retreive image default flag & data if set
836	Data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Image Controls'))
837	if Data['setDefault']:
838	G2frame.imageDefault = Data
839	file.close()
840
841	finally:
842	wx.EndBusyCursor()
843	print 'project load successful'
844	G2frame.NewPlot = True
845
846	def ProjFileSave(G2frame):
847	if not G2frame.PatternTree.IsEmpty():
848	file = open(G2frame.GSASprojectfile,'wb')
849	print 'save to file: ',G2frame.GSASprojectfile
850	wx.BeginBusyCursor()
851	try:
852	item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
853	while item:
854	data = []
855	name = G2frame.PatternTree.GetItemText(item)
856	data.append([name,G2frame.PatternTree.GetItemPyData(item)])
857	item2, cookie2 = G2frame.PatternTree.GetFirstChild(item)
858	while item2:
859	name = G2frame.PatternTree.GetItemText(item2)
860	data.append([name,G2frame.PatternTree.GetItemPyData(item2)])
861	item2, cookie2 = G2frame.PatternTree.GetNextChild(item, cookie2)
862	item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
863	cPickle.dump(data,file,1)
864	file.close()
865	finally:
866	wx.EndBusyCursor()
867	print 'project save successful'
868
869	def SaveIntegration(G2frame,PickId,data):
870	azms = G2frame.Integrate[1]
871	X = G2frame.Integrate[2][:-1]
872	Xminmax = [X[0],X[-1]]
873	N = len(X)
874	Id = G2frame.PatternTree.GetItemParent(PickId)
875	name = G2frame.PatternTree.GetItemText(Id)
876	name = name.replace('IMG ','PWDR ')
877	Comments = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id, 'Comments'))
878	names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth']
879	codes = [0 for i in range(11)]
880	LRazm = data['LRazimuth']
881	Azms = []
882	if data['fullIntegrate'] and data['outAzimuths'] == 1:
883	Azms = [45.0,] #a poor man's average?
884	else:
885	for i,azm in enumerate(azms[:-1]):
886	Azms.append((azms[i+1]+azm)/2.)
887	for i,azm in enumerate(azms[:-1]):
888	item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
889	Id = 0
890	while item:
891	Name = G2frame.PatternTree.GetItemText(item)
892	if name == Name:
893	Id = item
894	item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
895	parms = ['PXC',data['wavelength'],0.0,0.99,1.0,-0.10,0.4,0.30,1.0,0.0001,Azms[i]] #set polarization for synchrotron radiation!
896	Y = G2frame.Integrate[0][i]
897	W = 1./Y #probably not true
898	Sample = G2pdG.SetDefaultSample()
899	if Id:
900	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id, 'Comments'),Comments)
901	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Limits'),[tuple(Xminmax),Xminmax])
902	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0],
903	{'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
904	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
905	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Peak List'),[])
906	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Index Peak List'),[])
907	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Unit Cells List'),[])
908	G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Reflection Lists'),{})
909	else:
910	Id = G2frame.PatternTree.AppendItem(parent=G2frame.root,text=name+" Azm= %.2f"%(Azms[i]))
911	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Comments'),Comments)
912	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
913	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0],
914	{'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
915	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
916	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
917	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Peak List'),[])
918	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Index Peak List'),[])
919	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
920	G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Reflection Lists'),{})
921	G2frame.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
922	G2frame.PatternTree.SelectItem(Id)
923	G2frame.PatternTree.Expand(Id)
924	G2frame.PatternId = Id
925
926	def powderFxyeSave(G2frame,exports,powderfile):
927	head,tail = ospath.split(powderfile)
928	name,ext = tail.split('.')
929	for i,export in enumerate(exports):
930	filename = ospath.join(head,name+'-%03d.'%(i)+ext)
931	prmname = filename.strip(ext)+'prm'
932	prm = open(prmname,'w') #old style GSAS parm file
933	PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
934	Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, \
935	PickId, 'Instrument Parameters'))[1::2] #get values & names
936	Inst = dict(zip(Names,Values))
937	print Inst['Type']
938	prm.write( ' 123456789012345678901234567890123456789012345678901234567890 '+'\n')
939	prm.write( 'INS BANK 1 '+'\n')
940	prm.write(('INS HTYPE %sR '+'\n')%(Inst['Type']))
941	if 'Lam1' in Inst: #Ka1 & Ka2
942	prm.write(('INS 1 ICONS%10.7f%10.7f 0.0000 0.990 0 0.500 '+'\n')%(Inst['Lam1'],Inst['Lam2']))
943	elif 'Lam' in Inst: #single wavelength
944	prm.write(('INS 1 ICONS%10.7f%10.7f 0.0000 0.990 0 0.500 '+'\n')%(Inst['Lam'],0.0))
945	prm.write( 'INS 1 IRAD 0 '+'\n')
946	prm.write( 'INS 1I HEAD '+'\n')
947	prm.write( 'INS 1I ITYP 0 0.0000 180.0000 1 '+'\n')
948	prm.write(('INS 1DETAZM%10.3f '+'\n')%(Inst['Azimuth']))
949	prm.write( 'INS 1PRCF1 3 8 0.00100 '+'\n')
950	prm.write(('INS 1PRCF11 %15.6g%15.6g%15.6g%15.6g '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
951	prm.write(('INS 1PRCF12 %15.6g%15.6g%15.6g%15.6g '+'\n')%(Inst['X'],Inst['Y'],Inst['SH/L']/2.,Inst['SH/L']/2.))
952	prm.close()
953	file = open(filename,'w')
954	print 'save powder pattern to file: ',filename
955	x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
956	file.write(powderfile+'\n')
957	file.write('Instrument parameter file:'+ospath.split(prmname)[1]+'\n')
958	file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
959	100.x[0],100.(x[1]-x[0])))
960	s = list(np.sqrt(1./np.array(w)))
961	XYW = zip(x,y,s)
962	for X,Y,S in XYW:
963	file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
964	file.close()
965	print 'powder pattern file '+filename+' written'
966
967	def powderXyeSave(G2frame,exports,powderfile):
968	head,tail = ospath.split(powderfile)
969	name,ext = tail.split('.')
970	for i,export in enumerate(exports):
971	filename = ospath.join(head,name+'-%03d.'%(i)+ext)
972	PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
973	file = open(filename,'w')
974	file.write('#%s\n'%(export))
975	print 'save powder pattern to file: ',filename
976	x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
977	s = list(np.sqrt(1./np.array(w)))
978	XYW = zip(x,y,s)
979	for X,Y,W in XYW:
980	file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
981	file.close()
982	print 'powder pattern file '+filename+' written'
983
984	def PDFSave(G2frame,exports):
985	for export in exports:
986	PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
987	SQname = 'S(Q)'+export[4:]
988	GRname = 'G(R)'+export[4:]
989	sqfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.sq')
990	grfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.gr')
991	sqId = G2gd.GetPatternTreeItemId(G2frame, PickId, SQname)
992	grId = G2gd.GetPatternTreeItemId(G2frame, PickId, GRname)
993	sqdata = np.array(G2frame.PatternTree.GetItemPyData(sqId)[1][:2]).T
994	grdata = np.array(G2frame.PatternTree.GetItemPyData(grId)[1][:2]).T
995	sqfile = open(sqfilename,'w')
996	grfile = open(grfilename,'w')
997	sqfile.write('#T S(Q) %s\n'%(export))
998	grfile.write('#T G(R) %s\n'%(export))
999	sqfile.write('#L Q S(Q)\n')
1000	grfile.write('#L R G(R)\n')
1001	for q,sq in sqdata:
1002	sqfile.write("%15.6g %15.6g\n" % (q,sq))
1003	sqfile.close()
1004	for r,gr in grdata:
1005	grfile.write("%15.6g %15.6g\n" % (r,gr))
1006	grfile.close()
1007
1008	def PeakListSave(G2frame,file,peaks):
1009	print 'save peak list to file: ',G2frame.peaklistfile
1010	if not peaks:
1011	dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1012	try:
1013	result = dlg.ShowModal()
1014	finally:
1015	dlg.Destroy()
1016	return
1017	for peak in peaks:
1018	file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1019	(peak[0],peak[2],peak[4],peak[6]))
1020	print 'peak list saved'
1021
1022	def IndexPeakListSave(G2frame,peaks):
1023	file = open(G2frame.peaklistfile,'wa')
1024	print 'save index peak list to file: ',G2frame.peaklistfile
1025	wx.BeginBusyCursor()
1026	try:
1027	if not peaks:
1028	dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1029	try:
1030	result = dlg.ShowModal()
1031	finally:
1032	dlg.Destroy()
1033	return
1034	for peak in peaks:
1035	file.write("%12.6f\n" % (peak[7]))
1036	file.close()
1037	finally:
1038	wx.EndBusyCursor()
1039	print 'index peak list saved'
1040
1041	def SetNewPhase(Name='New Phase',SGData=G2spc.SpcGroup('P 1')[1],cell=[1.0,1.0,1.0,90.,90,90.,1.]):
1042	phaseData = {
1043	'General':{
1044	'Name':Name,
1045	'Type':'nuclear',
1046	'SGData':SGData,
1047	'Cell':[False,]+cell,
1048	'Pawley dmin':1.0,
1049	'Data plot type':'None',
1050	'SH Texture':{
1051	'Order':0,
1052	'Model':'cylindrical',
1053	'Sample omega':[False,0.0],
1054	'Sample chi':[False,0.0],
1055	'Sample phi':[False,0.0],
1056	'SH Coeff':[False,{}],
1057	'SHShow':False,
1058	'PFhkl':[0,0,1],
1059	'PFxyz':[0,0,1],
1060	'PlotType':'Pole figure'}},
1061	'Atoms':[],
1062	'Drawing':{},
1063	'Histograms':{},
1064	'Pawley ref':[],
1065	'Models':{},
1066	}
1067	return phaseData
1068
1069	def ReadEXPPhase(G2frame,filename):
1070	shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1071	textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1072	'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1073	'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1074	shNcof = 0
1075	file = open(filename, 'Ur')
1076	S = 1
1077	Expr = [{},{},{},{},{},{},{},{},{}]
1078	while S:
1079	S = file.readline()
1080	if 'EXPR NPHAS' in S[:12]:
1081	Num = S[12:-1].count('0')
1082	NPhas = S[12:-1].split()
1083	if 'CRS' in S[:3]:
1084	N = int(S[3:4])-1
1085	Expr[N][S[:12]] = S[12:-1]
1086	file.close()
1087	PNames = []
1088	for n,N in enumerate(NPhas):
1089	if N != '0':
1090	result = n
1091	key = 'CRS'+str(n+1)+' PNAM'
1092	PNames.append(Expr[n][key])
1093	if Num < 8:
1094	dlg = wx.SingleChoiceDialog(G2frame, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1095	try:
1096	if dlg.ShowModal() == wx.ID_OK:
1097	result = dlg.GetSelection()
1098	finally:
1099	dlg.Destroy()
1100	EXPphase = Expr[result]
1101	keyList = EXPphase.keys()
1102	keyList.sort()
1103	SGData = {}
1104	if NPhas[result] == '1':
1105	Ptype = 'nuclear'
1106	elif NPhas[result] in ['2','3']:
1107	Ptype = 'magnetic'
1108	elif NPhas[result] == '4':
1109	Ptype = 'macromolecular'
1110	elif NPhas[result] == '10':
1111	Ptype = 'Pawley'
1112	for key in keyList:
1113	if 'PNAM' in key:
1114	PhaseName = EXPphase[key].strip()
1115	elif 'ABC ' in key:
1116	abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]
1117	elif 'ANGLES' in key:
1118	angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]
1119	elif 'SG SYM' in key:
1120	SpGrp = EXPphase[key][:15].strip()
1121	E,SGData = G2spc.SpcGroup(SpGrp)
1122	elif 'OD ' in key:
1123	SHdata = EXPphase[key].split() # may not have all 9 values
1124	SHvals = 9*[0]
1125	for i in range(9):
1126	try:
1127	float(SHdata[i])
1128	SHvals[i] = SHdata[i]
1129	except:
1130	pass
1131	textureData['Order'] = int(SHvals[0])
1132	textureData['Model'] = shModels[int(SHvals[2])]
1133	textureData['Sample omega'] = [False,float(SHvals[6])]
1134	textureData['Sample chi'] = [False,float(SHvals[7])]
1135	textureData['Sample phi'] = [False,float(SHvals[8])]
1136	shNcof = int(SHvals[1])
1137	Atoms = []
1138	if Ptype == 'nuclear':
1139	for key in keyList:
1140	if 'AT' in key:
1141	if key[11:] == 'A':
1142	S = EXPphase[key]
1143	elif key[11:] == 'B':
1144	S += EXPphase[key]
1145	Atom = [S[50:58].strip(),S[:10].strip(),'',
1146	float(S[10:20]),float(S[20:30]),float(S[30:40]),
1147	float(S[40:50]),'',int(S[60:62]),S[130:131]]
1148	if Atom[9] == 'I':
1149	Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1150	elif Atom[9] == 'A':
1151	Atom += [0.0,float(S[68:78]),float(S[78:88]),
1152	float(S[88:98]),float(S[98:108]),
1153	float(S[108:118]),float(S[118:128])]
1154	XYZ = Atom[3:6]
1155	Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1156	Atom.append(ran.randint(0,sys.maxint))
1157	Atoms.append(Atom)
1158	elif Ptype == 'macromolecular':
1159	for key in keyList:
1160	if 'AT' in key[6:8]:
1161	S = EXPphase[key]
1162	Atom = [S[56:60],S[50:54].strip().upper(),S[54:56],
1163	S[46:51].strip(),S[:8].strip(),'',
1164	float(S[16:24]),float(S[24:32]),float(S[32:40]),
1165	float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1166	XYZ = Atom[6:9]
1167	Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1168	Atom.append(ran.randint(0,sys.maxint))
1169	Atoms.append(Atom)
1170	Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1171	if shNcof:
1172	shCoef = {}
1173	nRec = [i+1 for i in range((shNcof-1)/6+1)]
1174	for irec in nRec:
1175	ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1176	indx = EXPphase[ODkey].split()
1177	ODkey = ODkey[:-1]+'B'
1178	vals = EXPphase[ODkey].split()
1179	for i,val in enumerate(vals):
1180	key = 'C(%s,%s,%s)'%(indx[3i],indx[3i+1],indx[3*i+2])
1181	shCoef[key] = float(val)
1182	textureData['SH Coeff'] = [False,shCoef]
1183
1184	Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=abc+angles+[Volume,])
1185	general = Phase['General']
1186	general['Type'] = Ptype
1187	general['SH Texture'] = textureData
1188	Phase['Atoms'] = Atoms
1189	return Phase
1190
1191	def ReadPDBPhase(filename):
1192	EightPiSq = 8.math.pi*2
1193	file = open(filename, 'Ur')
1194	Phase = {}
1195	Title = ''
1196	Compnd = ''
1197	Atoms = []
1198	A = np.zeros(shape=(3,3))
1199	S = file.readline()
1200	while S:
1201	Atom = []
1202	if 'TITLE' in S[:5]:
1203	Title = S[10:72].strip()
1204	S = file.readline()
1205	elif 'COMPND ' in S[:10]:
1206	Compnd = S[10:72].strip()
1207	S = file.readline()
1208	elif 'CRYST' in S[:5]:
1209	abc = S[7:34].split()
1210	angles = S[34:55].split()
1211	cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1212	float(angles[0]),float(angles[1]),float(angles[2])]
1213	Volume = G2lat.calc_V(G2lat.cell2A(cell))
1214	AA,AB = G2lat.cell2AB(cell)
1215	SpGrp = S[55:65]
1216	E,SGData = G2spc.SpcGroup(SpGrp)
1217	while E:
1218	print G2spc.SGErrors(E)
1219	dlg = wx.TextEntryDialog(None,
1220	SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
1221	'ERROR in space group symbol','',style=wx.OK)
1222	if dlg.ShowModal() == wx.ID_OK:
1223	SpGrp = dlg.GetValue()
1224	E,SGData = G2spc.SpcGroup(SpGrp)
1225	else:
1226	return None
1227	dlg.Destroy()
1228	SGlines = G2spc.SGPrint(SGData)
1229	for line in SGlines: print line
1230	S = file.readline()
1231	elif 'SCALE' in S[:5]:
1232	V = (S[10:41].split())
1233	A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1234	S = file.readline()
1235	elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1236	XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1237	XYZ = np.inner(AB,XYZ)
1238	XYZ = np.where(abs(XYZ)<0.00001,0,XYZ)
1239	SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1240	Uiso = float(S[61:67])/EightPiSq
1241	Type = S[12:14].upper()
1242	if Type[0] in '123456789':
1243	Type = Type[1:]
1244	Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1245	S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1246	float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1247	S = file.readline()
1248	if 'ANISOU' in S[:6]:
1249	Uij = S[30:72].split()
1250	Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1251	float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1252	Atom = Atom[:14]+Uij
1253	Atom[12] = 'A'
1254	S = file.readline()
1255	Atom.append(ran.randint(0,sys.maxint))
1256	Atoms.append(Atom)
1257	else:
1258	S = file.readline()
1259	file.close()
1260	if Title:
1261	PhaseName = Title
1262	elif Compnd:
1263	PhaseName = Compnd
1264	else:
1265	PhaseName = 'None'
1266	Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=cell+[Volume,])
1267	Phase['General']['Type'] = 'macromolecular'
1268	Phase['Atoms'] = Atoms
1269
1270	return Phase
1271
1272	class MultipleChoicesDialog(wx.Dialog):
1273	'''A dialog that offers a series of choices, each with a title and a wx.Choice
1274	widget. Intended to be used Modally.
1275	typical input:
1276	choicelist=[ ('a','b','c'), ('test1','test2'),('no choice',)]
1277	headinglist = [ 'select a, b or c', 'select 1 of 2', 'No option here']
1278	selections are placed in self.chosen when OK is pressed
1279	'''
1280	def __init__(self,choicelist,headinglist,
1281	head='Select options',
1282	title='Please select from options below',
1283	parent=None):
1284	self.chosen = []
1285	wx.Dialog.__init__(
1286	self,parent,wx.ID_ANY,head,
1287	pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
1288	panel = wx.Panel(self)
1289	mainSizer = wx.BoxSizer(wx.VERTICAL)
1290	mainSizer.Add((10,10),1)
1291	topLabl = wx.StaticText(panel,wx.ID_ANY,title)
1292	mainSizer.Add(topLabl,0,wx.ALIGN_CENTER_VERTICAL\|wx.CENTER,10)
1293	self.ChItems = []
1294	for choice,lbl in zip(choicelist,headinglist):
1295	mainSizer.Add((10,10),1)
1296	self.chosen.append(0)
1297	topLabl = wx.StaticText(panel,wx.ID_ANY,' '+lbl)
1298	mainSizer.Add(topLabl,0,wx.ALIGN_LEFT,10)
1299	self.ChItems.append(wx.Choice(self, wx.ID_ANY, (100, 50), choices = choice))
1300	mainSizer.Add(self.ChItems[-1],0,wx.ALIGN_CENTER,10)
1301
1302	OkBtn = wx.Button(panel,-1,"Ok")
1303	OkBtn.Bind(wx.EVT_BUTTON, self.OnOk)
1304	cancelBtn = wx.Button(panel,-1,"Cancel")
1305	cancelBtn.Bind(wx.EVT_BUTTON, self.OnCancel)
1306	btnSizer = wx.BoxSizer(wx.HORIZONTAL)
1307	btnSizer.Add((20,20),1)
1308	btnSizer.Add(OkBtn)
1309	btnSizer.Add((20,20),1)
1310	btnSizer.Add(cancelBtn)
1311	btnSizer.Add((20,20),1)
1312	mainSizer.Add(btnSizer,0,wx.EXPAND\|wx.BOTTOM\|wx.TOP, 10)
1313	panel.SetSizer(mainSizer)
1314	panel.Fit()
1315	self.Fit()
1316
1317	def OnOk(self,event):
1318	parent = self.GetParent()
1319	if parent is not None: parent.Raise()
1320	# save the results from the choice widgets
1321	self.chosen = []
1322	for w in self.ChItems:
1323	self.chosen.append(w.GetSelection())
1324	self.EndModal(wx.ID_OK)
1325
1326	def OnCancel(self,event):
1327	parent = self.GetParent()
1328	if parent is not None: parent.Raise()
1329	self.chosen = []
1330	self.EndModal(wx.ID_CANCEL)
1331
1332	######################################################################
1333	# base classes for reading various types of data files
1334	# not used directly, only by subclassing
1335	######################################################################
1336	E,SGData = G2spc.SpcGroup('P 1') # data structure for default space group
1337	class ImportBaseclass(object):
1338	'''Defines a base class for the importing of data files (diffraction
1339	data, coordinates,...
1340	'''
1341	def __init__(self,
1342	formatName,
1343	longFormatName=None,
1344	extensionlist=[],
1345	strictExtension=False,
1346	):
1347	self.formatName = formatName # short string naming file type
1348	if longFormatName: # longer string naming file type
1349	self.longFormatName = longFormatName
1350	else:
1351	self.longFormatName = formatName
1352	# define extensions that are allowed for the file type
1353	# for windows, remove any extensions that are duplicate, as case is ignored
1354	if sys.platform == 'windows' and extensionlist:
1355	extensionlist = list(set([s.lower() for s in extensionlist]))
1356	self.extensionlist = extensionlist
1357	# If strictExtension is True, the file will not be read, unless
1358	# the extension matches one in the extensionlist
1359	self.strictExtension = strictExtension
1360	self.warnings = ''
1361	self.errors = ''
1362	# used for readers that will use multiple passes to read
1363	# more than one data block
1364	self.repeat = False
1365	self.repeatcount = 0
1366	#print 'created',self.__class__
1367
1368	def BlockSelector(self, ChoiceList, ParentFrame=None,
1369	title='Select a block',
1370	size=None, header='Block Selector'):
1371	''' Provide a wx dialog to select a block if the file contains more
1372	than one set of data and one must be selected
1373	'''
1374	dlg = wx.SingleChoiceDialog(
1375	ParentFrame,
1376	title, header,
1377	ChoiceList,
1378	)
1379	if size: dlg.SetSize(size)
1380	if dlg.ShowModal() == wx.ID_OK:
1381	sel = dlg.GetSelection()
1382	dlg.Destroy()
1383	return sel
1384	else:
1385	dlg.Destroy()
1386	return None
1387
1388	def MultipleBlockSelector(self, ChoiceList, ParentFrame=None,
1389	title='Select a block',
1390	size=None, header='Block Selector'):
1391	''' Provide a wx dialog to select a block of data if the file contains more
1392	than one set of data and one must be selected.
1393	Returns a list of the selected blocks
1394	'''
1395	dlg = wx.MultiChoiceDialog(
1396	ParentFrame,
1397	title, header,
1398	ChoiceList+['Select all'],
1399	wx.CHOICEDLG_STYLE
1400	)
1401	if size: dlg.SetSize(size)
1402	if dlg.ShowModal() == wx.ID_OK:
1403	sel = dlg.GetSelections()
1404	dlg.Destroy()
1405	else:
1406	dlg.Destroy()
1407	return []
1408	selected = []
1409	if len(ChoiceList) in sel:
1410	return range(len(ChoiceList))
1411	else:
1412	return sel
1413	return selected
1414
1415	def MultipleChoicesDialog(self, choicelist, headinglist, ParentFrame=None, **kwargs):
1416	'''A modal dialog that offers a series of choices, each with a title and a wx.Choice
1417	widget.
1418	typical input:
1419	choicelist=[ ('a','b','c'), ('test1','test2'),('no choice',)]
1420	headinglist = [ 'select a, b or c', 'select 1 of 2', 'No option here']
1421	optional keyword parameters are: head (window title) and title
1422	returns a list of selected indicies for each choice (or None)
1423	'''
1424	result = None
1425	dlg = MultipleChoicesDialog(choicelist,headinglist,
1426	parent=ParentFrame, **kwargs)
1427	if dlg.ShowModal() == wx.ID_OK:
1428	result = dlg.chosen
1429	dlg.Destroy()
1430	return result
1431
1432	def ShowBusy(self):
1433	wx.BeginBusyCursor()
1434
1435	def DoneBusy(self):
1436	wx.EndBusyCursor()
1437
1438	# def Reader(self, filename, filepointer, ParentFrame=None, **unused):
1439	# '''This method must be supplied in the child class
1440	# it will read the file
1441	# '''
1442	# return True # if read OK
1443	# return False # if an error occurs
1444
1445	def ExtensionValidator(self, filename):
1446	'''This methods checks if the file has the correct extension
1447	Return False if this filename will not be supported by this reader
1448	Return True if the extension matches the list supplied by the reader
1449	Return None if the reader allows un-registered extensions
1450	'''
1451	if filename:
1452	ext = os.path.splitext(filename)[1]
1453	if sys.platform == 'windows': ext = ext.lower()
1454	if ext in self.extensionlist: return True
1455	if self.strictExtension: return False
1456	return None
1457
1458	def ContentsValidator(self, filepointer):
1459	'''This routine will attempt to determine if the file can be read
1460	with the current format.
1461	This will typically be overridden with a method that
1462	takes a quick scan of [some of]
1463	the file contents to do a "sanity" check if the file
1464	appears to match the selected format.
1465	Expected to be called via self.Validator()
1466	'''
1467	#filepointer.seek(0) # rewind the file pointer
1468	return True
1469
1470	class ImportPhase(ImportBaseclass):
1471	'''Defines a base class for the reading of files with coordinates
1472	'''
1473	def __init__(self,
1474	formatName,
1475	longFormatName=None,
1476	extensionlist=[],
1477	strictExtension=False,
1478	):
1479	# call parent __init__
1480	ImportBaseclass.__init__(self,formatName,
1481	longFormatName,
1482	extensionlist,
1483	strictExtension)
1484	# define a default Phase structure
1485	self.Phase = SetNewPhase(Name='new phase',SGData=SGData)
1486
1487	def PhaseSelector(self, ChoiceList, ParentFrame=None,
1488	title='Select a phase', size=None,
1489	header='Phase Selector'):
1490	''' Provide a wx dialog to select a phase if the file contains more
1491	than one phase
1492	'''
1493	return self.BlockSelector(ChoiceList,
1494	ParentFrame,
1495	title,
1496	size,
1497	header)
1498
1499	######################################################################
1500	class ImportStructFactor(ImportBaseclass):
1501	'''Defines a base class for the reading of files with tables
1502	of structure factors
1503	'''
1504	def __init__(self,
1505	formatName,
1506	longFormatName=None,
1507	extensionlist=[],
1508	strictExtension=False,
1509	):
1510	ImportBaseclass.__init__(self,formatName,
1511	longFormatName,
1512	extensionlist,
1513	strictExtension)
1514
1515	# define contents of Structure Factor entry
1516	self.Controls = { # dictionary with plotting controls
1517	'Type' : 'Fosq',
1518	'ifFc' : False, #
1519	'HKLmax' : [None,None,None],
1520	'HKLmin' : [None,None,None],
1521	'FoMax' : None, # maximum observed structure factor as Fo
1522	'Zone' : '001',
1523	'Layer' : 0,
1524	'Scale' : 1.0,
1525	'log-lin' : 'lin',
1526	}
1527	self.Parameters = [ # list with data collection parameters
1528	('SXC',1.5428),
1529	['SXC',1.5428],
1530	['Type','Lam']
1531	]
1532	self.RefList = []
1533
1534	def UpdateParameters(self,Type=None,Wave=None):
1535	HistType = self.Parameters[0][0]
1536	HistWave = self.Parameters[0][1]
1537	if Type is not None:
1538	HistType = Type
1539	if Wave is not None:
1540	HistWave = Wave
1541	self.Parameters = [ # overwrite entire list
1542	(HistType,HistWave),
1543	[HistType,HistWave],
1544	['Type','Lam']
1545	]
1546
1547	def UpdateControls(self,Type='Fosq',FcalcPresent=False):
1548	'''Scan through the reflections to update the Controls dictionary
1549	'''
1550	self.Controls['Type'] = Type
1551	self.Controls['iffc'] = FcalcPresent
1552	HKLmax = [None,None,None]
1553	HKLmin = [None,None,None]
1554	Fo2max = None
1555	for HKL,Fo2,SFo2,Fc,Fcp,Fcpp,phase in self.RefList:
1556	if Fo2max is None:
1557	Fo2max = Fo2
1558	else:
1559	Fo2max = max(Fo2max,Fo2)
1560	for i,hkl in enumerate(HKL):
1561	if HKLmax[i] is None:
1562	HKLmax[i] = hkl
1563	HKLmin[i] = hkl
1564	else:
1565	HKLmax[i] = max(HKLmax[i],hkl)
1566	HKLmin[i] = min(HKLmin[i],hkl)
1567	self.Controls['HKLmax'] = HKLmax
1568	self.Controls['HKLmin'] = HKLmin
1569	if Type == 'Fosq':
1570	self.Controls['FoMax'] = np.sqrt(Fo2max)
1571	elif Type == 'Fo':
1572	self.Controls['FoMax'] = Fo2max
1573	else:
1574	print "Unsupported Stract Fact type in ImportStructFactor.UpdateControls"
1575	raise Exception,"Unsupported Stract Fact type in ImportStructFactor.UpdateControls"
1576
1577	######################################################################
1578	class ImportPowderData(ImportBaseclass):
1579	'''Defines a base class for the reading of files with powder data
1580	'''
1581	# define some default instrument parameter files
1582	# just like GSAS, sigh
1583	Iparm_CuKa12 = { # Default Inst. parms for CuKa lab data
1584	'INS HTYPE ':'PXC ',
1585	'INS 1 ICONS':' 1.540500 1.544300 0.0 0 0.7 0 0.5 ',
1586	'INS 1PRCF1 ':' 3 8 0.01 ',
1587	'INS 1PRCF11':' 2.000000E+00 -2.000000E+00 5.000000E+00 0.000000E+00 ',
1588	'INS 1PRCF12':' 0.000000E+00 0.000000E+00 0.150000E-01 0.150000E-01 ',
1589	}
1590	Iparm_Sync06 = { # Default Inst. parms for 0.6A synchrotron data
1591	'INS HTYPE ':'PXC ',
1592	'INS 1 ICONS':' 0.600000 0.000000 0.0 0 0.99 0 0.5 ',
1593	'INS 1PRCF1 ':' 3 8 0.01 ',
1594	'INS 1PRCF11':' 1.000000E+00 -1.000000E+00 0.300000E+00 0.000000E+00 ',
1595	'INS 1PRCF12':' 0.000000E+00 0.000000E+00 0.100000E-01 0.100000E-01 ',
1596	}
1597	def __init__(self,
1598	formatName,
1599	longFormatName=None,
1600	extensionlist=[],
1601	strictExtension=False,
1602	):
1603	ImportBaseclass.__init__(self,formatName,
1604	longFormatName,
1605	extensionlist,
1606	strictExtension)
1607	self.powderentry = ['',None,None] # (filename,Pos,Bank)
1608	self.powderdata = [] # Powder dataset
1609	'''A powder data set is a list with items [x,y,w,yc,yb,yd]:
1610	np.array(x), # x-axis values
1611	np.array(y), # powder pattern intensities
1612	np.array(w), # 1/sig(intensity)^2 values (weights)
1613	np.array(yc), # calc. intensities (zero)
1614	np.array(yb), # calc. background (zero)
1615	np.array(yd), # obs-calc profiles
1616	'''
1617	self.comments = []
1618	self.idstring = ''
1619	self.Sample = G2pdG.SetDefaultSample()
1620	self.instparm = None # name hint
1621	self.instfile = '' # full path name to instrument parameter file
1622	self.instbank = '' # inst parm bank number
1623	self.instmsg = '' # a label that gets printed to show
1624	# where instrument parameters are from
1625	self.numbanks = 1
1626	self.instdict = {} # place items here that will be transferred to the instrument parameters
1627
1628	if __name__ == '__main__':
1629	app = wx.PySimpleApp()
1630	frm = wx.Frame(None) # create a frame
1631	choicelist=[ ('a','b','c'),
1632	('test1','test2'),('no choice',)]
1633	titles = [ 'a, b or c', 'tests', 'No option here']
1634	dlg = MultipleChoicesDialog(
1635	choicelist,titles,
1636	parent=frm)
1637	if dlg.ShowModal() == wx.ID_OK:
1638	print 'Got OK'

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

Download in other formats: