Context Navigation

source: trunk/GSASIIIO.py @ 614

Last change on this file since 614 was 614, checked in by toby, 11 years ago
import GSAS powder data tested; import routines moved to imports
Property svn:keywords set to `Date Author Revision URL Id`
File size: 61.1 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-11 17:03:59 +0000 (Fri, 11 May 2012) $
7	# $Author: toby $
8	# $Revision: 614 $
9	# $URL: trunk/GSASIIIO.py $
10	# $Id: GSASIIIO.py 614 2012-05-11 17:03:59Z 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	######################################################################
1273	# base classes for reading various types of data files
1274	# not used directly, only by subclassing
1275	######################################################################
1276	E,SGData = G2spc.SpcGroup('P 1') # data structure for default space group
1277	class ImportBaseclass(object):
1278	'''Defines a base class for the importing of data files (diffraction
1279	data, coordinates,...
1280	'''
1281	def __init__(self,
1282	formatName,
1283	longFormatName=None,
1284	extensionlist=[],
1285	strictExtension=False,
1286	):
1287	self.formatName = formatName # short string naming file type
1288	if longFormatName: # longer string naming file type
1289	self.longFormatName = longFormatName
1290	else:
1291	self.longFormatName = formatName
1292	# define extensions that are allowed for the file type
1293	# for windows, remove any extensions that are duplicate, as case is ignored
1294	if sys.platform == 'windows' and extensionlist:
1295	extensionlist = list(set([s.lower() for s in extensionlist]))
1296	self.extensionlist = extensionlist
1297	# If strictExtension is True, the file will not be read, unless
1298	# the extension matches one in the extensionlist
1299	self.strictExtension = strictExtension
1300	self.warnings = ''
1301	self.errors = ''
1302	# used for readers that will use multiple passes to read
1303	# more than one data block
1304	self.repeat = False
1305	self.repeatcount = 0
1306	#print 'created',self.__class__
1307
1308	def BlockSelector(self, ChoiceList, ParentFrame=None,
1309	title='Select a block',
1310	size=None, header='Block Selector'):
1311	''' Provide a wx dialog to select a block if the file contains more
1312	than one set of data and one must be selected
1313	'''
1314	dlg = wx.SingleChoiceDialog(
1315	ParentFrame,
1316	title, header,
1317	ChoiceList,
1318	)
1319	if size: dlg.SetSize(size)
1320	if dlg.ShowModal() == wx.ID_OK:
1321	sel = dlg.GetSelection()
1322	dlg.Destroy()
1323	return sel
1324	else:
1325	dlg.Destroy()
1326	return None
1327
1328	def MultipleBlockSelector(self, ChoiceList, ParentFrame=None,
1329	title='Select a block',
1330	size=None, header='Block Selector'):
1331	''' Provide a wx dialog to select a block of data if the file contains more
1332	than one set of data and one must be selected.
1333	Returns a list of the selected blocks
1334	'''
1335	dlg = wx.MultiChoiceDialog(
1336	ParentFrame,
1337	title, header,
1338	ChoiceList+['Select all'],
1339	wx.CHOICEDLG_STYLE
1340	)
1341	if size: dlg.SetSize(size)
1342	if dlg.ShowModal() == wx.ID_OK:
1343	sel = dlg.GetSelections()
1344	dlg.Destroy()
1345	else:
1346	dlg.Destroy()
1347	return []
1348	selected = []
1349	if len(ChoiceList) in sel:
1350	return range(len(ChoiceList))
1351	else:
1352	return sel
1353	return selected
1354
1355	def ShowBusy(self):
1356	wx.BeginBusyCursor()
1357
1358	def DoneBusy(self):
1359	wx.EndBusyCursor()
1360
1361	# def Reader(self, filename, filepointer, ParentFrame=None, **unused):
1362	# '''This method must be supplied in the child class
1363	# it will read the file
1364	# '''
1365	# return True # if read OK
1366	# return False # if an error occurs
1367
1368	def ExtensionValidator(self, filename):
1369	'''This methods checks if the file has the correct extension
1370	Return False if this filename will not be supported by this reader
1371	Return True if the extension matches the list supplied by the reader
1372	Return None if the reader allows un-registered extensions
1373	'''
1374	if filename:
1375	ext = os.path.splitext(filename)[1]
1376	if sys.platform == 'windows': ext = ext.lower()
1377	if ext in self.extensionlist: return True
1378	if self.strictExtension: return False
1379	return None
1380
1381	def ContentsValidator(self, filepointer):
1382	'''This routine will attempt to determine if the file can be read
1383	with the current format.
1384	This will typically be overridden with a method that
1385	takes a quick scan of [some of]
1386	the file contents to do a "sanity" check if the file
1387	appears to match the selected format.
1388	Expected to be called via self.Validator()
1389	'''
1390	#filepointer.seek(0) # rewind the file pointer
1391	return True
1392
1393	class ImportPhase(ImportBaseclass):
1394	'''Defines a base class for the reading of files with coordinates
1395	'''
1396	def __init__(self,
1397	formatName,
1398	longFormatName=None,
1399	extensionlist=[],
1400	strictExtension=False,
1401	):
1402	# call parent __init__
1403	ImportBaseclass.__init__(self,formatName,
1404	longFormatName,
1405	extensionlist,
1406	strictExtension)
1407	# define a default Phase structure
1408	self.Phase = SetNewPhase(Name='new phase',SGData=SGData)
1409
1410	def PhaseSelector(self, ChoiceList, ParentFrame=None,
1411	title='Select a phase', size=None,
1412	header='Phase Selector'):
1413	''' Provide a wx dialog to select a phase if the file contains more
1414	than one phase
1415	'''
1416	return self.BlockSelector(ChoiceList,
1417	ParentFrame,
1418	title,
1419	size,
1420	header)
1421
1422	######################################################################
1423	class ImportStructFactor(ImportBaseclass):
1424	'''Defines a base class for the reading of files with tables
1425	of structure factors
1426	'''
1427	def __init__(self,
1428	formatName,
1429	longFormatName=None,
1430	extensionlist=[],
1431	strictExtension=False,
1432	):
1433	ImportBaseclass.__init__(self,formatName,
1434	longFormatName,
1435	extensionlist,
1436	strictExtension)
1437
1438	# define contents of Structure Factor entry
1439	self.Controls = { # dictionary with plotting controls
1440	'Type' : 'Fosq',
1441	'ifFc' : False, #
1442	'HKLmax' : [None,None,None],
1443	'HKLmin' : [None,None,None],
1444	'FoMax' : None, # maximum observed structure factor as Fo
1445	'Zone' : '001',
1446	'Layer' : 0,
1447	'Scale' : 1.0,
1448	'log-lin' : 'lin',
1449	}
1450	self.Parameters = [ # list with data collection parameters
1451	('SXC',1.5428),
1452	['SXC',1.5428],
1453	['Type','Lam']
1454	]
1455	self.RefList = []
1456
1457	def UpdateParameters(self,Type=None,Wave=None):
1458	HistType = self.Parameters[0][0]
1459	HistWave = self.Parameters[0][1]
1460	if Type is not None:
1461	HistType = Type
1462	if Wave is not None:
1463	HistWave = Wave
1464	self.Parameters = [ # overwrite entire list
1465	(HistType,HistWave),
1466	[HistType,HistWave],
1467	['Type','Lam']
1468	]
1469
1470	def UpdateControls(self,Type='Fosq',FcalcPresent=False):
1471	'''Scan through the reflections to update the Controls dictionary
1472	'''
1473	self.Controls['Type'] = Type
1474	self.Controls['iffc'] = FcalcPresent
1475	HKLmax = [None,None,None]
1476	HKLmin = [None,None,None]
1477	Fo2max = None
1478	for HKL,Fo2,SFo2,Fc,Fcp,Fcpp,phase in self.RefList:
1479	if Fo2max is None:
1480	Fo2max = Fo2
1481	else:
1482	Fo2max = max(Fo2max,Fo2)
1483	for i,hkl in enumerate(HKL):
1484	if HKLmax[i] is None:
1485	HKLmax[i] = hkl
1486	HKLmin[i] = hkl
1487	else:
1488	HKLmax[i] = max(HKLmax[i],hkl)
1489	HKLmin[i] = min(HKLmin[i],hkl)
1490	self.Controls['HKLmax'] = HKLmax
1491	self.Controls['HKLmin'] = HKLmin
1492	if Type == 'Fosq':
1493	self.Controls['FoMax'] = np.sqrt(Fo2max)
1494	elif Type == 'Fo':
1495	self.Controls['FoMax'] = Fo2max
1496	else:
1497	print "Unsupported Stract Fact type in ImportStructFactor.UpdateControls"
1498	raise Exception,"Unsupported Stract Fact type in ImportStructFactor.UpdateControls"
1499
1500	######################################################################
1501	class ImportPowderData(ImportBaseclass):
1502	'''Defines a base class for the reading of files with powder data
1503	'''
1504	# define some default instrument parameter files
1505	# just like GSAS, sigh
1506	Iparm_CuKa12 = { # Default Inst. parms for CuKa lab data
1507	'INS HTYPE ':'PXC ',
1508	'INS 1 ICONS':' 1.540500 1.544300 0.0 0 0.7 0 0.5 ',
1509	'INS 1PRCF1 ':' 3 8 0.01 ',
1510	'INS 1PRCF11':' 2.000000E+00 -2.000000E+00 5.000000E+00 0.000000E+00 ',
1511	'INS 1PRCF12':' 0.000000E+00 0.000000E+00 0.150000E-01 0.150000E-01 ',
1512	}
1513	Iparm_Sync06 = { # Default Inst. parms for 0.6A synchrotron data
1514	'INS HTYPE ':'PXC ',
1515	'INS 1 ICONS':' 0.600000 0.000000 0.0 0 0.99 0 0.5 ',
1516	'INS 1PRCF1 ':' 3 8 0.01 ',
1517	'INS 1PRCF11':' 1.000000E+00 -1.000000E+00 0.300000E+00 0.000000E+00 ',
1518	'INS 1PRCF12':' 0.000000E+00 0.000000E+00 0.100000E-01 0.100000E-01 ',
1519	}
1520	def __init__(self,
1521	formatName,
1522	longFormatName=None,
1523	extensionlist=[],
1524	strictExtension=False,
1525	):
1526	ImportBaseclass.__init__(self,formatName,
1527	longFormatName,
1528	extensionlist,
1529	strictExtension)
1530	self.powderentry = ['',None,None] # (filename,Pos,Bank)
1531	self.powderdata = [] # Powder dataset
1532	'''A powder data set is a list with items [x,y,w,yc,yb,yd]:
1533	np.array(x), # x-axis values
1534	np.array(y), # powder pattern intensities
1535	np.array(w), # 1/sig(intensity)^2 values (weights)
1536	np.array(yc), # calc. intensities (zero)
1537	np.array(yb), # calc. background (zero)
1538	np.array(yd), # obs-calc profiles
1539	'''
1540	self.comments = []
1541	self.idstring = ''
1542	self.Sample = G2pdG.SetDefaultSample()
1543	self.instparm = None # name hint
1544	self.instfile = '' # full path name to instrument parameter file
1545	self.instbank = '' # inst parm bank number
1546	self.instmsg = '' # a label that gets printed to show
1547	# where instrument parameters are from
1548	self.numbanks = 1

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

Download in other formats: