Context Navigation

source: trunk/GSASIIIO.py @ 607

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

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

Download in other formats: