Context Navigation

source: MPbranch/GSASIIIO.py @ 568

Last change on this file since 568 was 485, checked in by vondreele, 13 years ago

Property svn:keywords set to `Date Author Revision URL Id`
File size: 53.1 KB

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

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

Download in other formats: