source: trunk/GSASIIstrIO.py @ 989

Last change on this file since 989 was 989, checked in by toby, 9 years ago

add constraint derivs for single xtal; allow cellFill to work without esds; update docs; fix controls copy bug

File size: 110.2 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrIO: structure I/O routines*
4-------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2013-05-17 12:13:15 -0500 (Fri, 17 May 2013) $
9# $Author: vondreele $
10# $Revision: 920 $
11# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIstrIO.py $
12# $Id: GSASIIstrIO.py 920 2013-05-17 17:13:15Z vondreele $
13########### SVN repository information ###################
14import sys
15import os
16import os.path as ospath
17import time
18import math
19import copy
20import random
21import cPickle
22import numpy as np
23import numpy.ma as ma
24import numpy.linalg as nl
25import scipy.optimize as so
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 920 $")
28import GSASIIElem as G2el
29import GSASIIlattice as G2lat
30import GSASIIspc as G2spc
31import GSASIImapvars as G2mv
32import GSASIImath as G2mth
33
34sind = lambda x: np.sin(x*np.pi/180.)
35cosd = lambda x: np.cos(x*np.pi/180.)
36tand = lambda x: np.tan(x*np.pi/180.)
37asind = lambda x: 180.*np.arcsin(x)/np.pi
38acosd = lambda x: 180.*np.arccos(x)/np.pi
39atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
40   
41ateln2 = 8.0*math.log(2.0)
42
43def GetControls(GPXfile):
44    ''' Returns dictionary of control items found in GSASII gpx file
45
46    :param str GPXfile: full .gpx file name
47    :return: dictionary of control items
48    '''
49    Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
50        'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
51    fl = open(GPXfile,'rb')
52    while True:
53        try:
54            data = cPickle.load(fl)
55        except EOFError:
56            break
57        datum = data[0]
58        if datum[0] == 'Controls':
59            Controls.update(datum[1])
60    fl.close()
61    return Controls
62   
63def GetConstraints(GPXfile):
64    '''Read the constraints from the GPX file and interpret them
65    '''
66    constList = []
67    fl = open(GPXfile,'rb')
68    while True:
69        try:
70            data = cPickle.load(fl)
71        except EOFError:
72            break
73        datum = data[0]
74        if datum[0] == 'Constraints':
75            constDict = datum[1]
76            for item in constDict:
77                constList += constDict[item]
78    fl.close()
79    constDict,fixedList,ignored = ProcessConstraints(constList)
80    if ignored:
81        print ignored,'old-style Constraints were rejected'
82    return constDict,fixedList
83   
84def ProcessConstraints(constList):
85    """Interpret the constraints in the constList input into a dictionary, etc.
86   
87    :param list constList: a list of lists where each item in the outer list
88      specifies a constraint of some form. The last item in each inner list
89      determines which of the four constraints types has been input:
90
91        * h (hold): a single variable that will not be varied. It
92          will be removed from the varyList later.
93        * c (constraint): specifies a linear relationship that
94          can be varied as a new grouped variable
95          a fixed value.
96        * f (fixed): specifies a linear relationship that is assigned
97          a fixed value.
98        * e (equivalence): specifies a series of variables where the
99          first variable in the last can be used to generate the
100          values for all the remaining variables.
101
102    :returns:  a tuple of (constDict,fixedList,ignored) where:
103     
104      * constDict (list) contains the constraint relationships
105      * fixedList (list) contains the fixed values for type
106        of constraint.
107      * ignored (int) counts the number of invalid constraint items
108        (should always be zero!)
109
110    """
111    constDict = []
112    fixedList = []
113    ignored = 0
114    for item in constList:
115        if item[-1] == 'h':
116            # process a hold
117            fixedList.append('0')
118            constDict.append({item[0][1]:0.0})
119        elif item[-1] == 'f':
120            # process a new variable
121            fixedList.append(None)
122            constDict.append({})
123            for term in item[:-3]:
124                constDict[-1][term[1]] = term[0]
125            #constFlag[-1] = ['Vary']
126        elif item[-1] == 'c': 
127            # process a contraint relationship
128            fixedList.append(str(item[-3]))
129            constDict.append({})
130            for term in item[:-3]:
131                constDict[-1][term[1]] = term[0]
132            #constFlag[-1] = ['VaryFree']
133        elif item[-1] == 'e':
134            # process an equivalence
135            firstmult = None
136            eqlist = []
137            for term in item[:-3]:
138                if term[0] == 0: term[0] = 1.0
139                if firstmult is None:
140                    firstmult,firstvar = term
141                else:
142                    eqlist.append([term[1],firstmult/term[0]])
143            G2mv.StoreEquivalence(firstvar,eqlist)
144        else:
145            ignored += 1
146    return constDict,fixedList,ignored
147
148def CheckConstraints(GPXfile):
149    '''Load constraints and related info and return any error or warning messages'''
150    # init constraints
151    G2mv.InitVars()   
152    # get variables
153    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
154    if not Phases:
155        return 'Error: No Phases!',''
156    if not Histograms:
157        return 'Error: no diffraction data',''
158    rigidbodyDict = GetRigidBodies(GPXfile)
159    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
160    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
161    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
162    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
163    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
164    varyList = rbVary+phaseVary+hapVary+histVary
165    constrDict,fixedList = GetConstraints(GPXfile)
166    return G2mv.CheckConstraints(varyList,constrDict,fixedList)
167   
168def GetRestraints(GPXfile):
169    '''Read the restraints from the GPX file.
170    Throws an exception if not found in the .GPX file
171    '''
172    fl = open(GPXfile,'rb')
173    while True:
174        try:
175            data = cPickle.load(fl)
176        except EOFError:
177            break
178        datum = data[0]
179        if datum[0] == 'Restraints':
180            restraintDict = datum[1]
181    fl.close()
182    return restraintDict
183   
184def GetRigidBodies(GPXfile):
185    '''Read the rigid body models from the GPX file
186    '''
187    fl = open(GPXfile,'rb')
188    while True:
189        try:
190            data = cPickle.load(fl)
191        except EOFError:
192            break
193        datum = data[0]
194        if datum[0] == 'Rigid bodies':
195            rigidbodyDict = datum[1]
196    fl.close()
197    return rigidbodyDict
198       
199def GetFprime(controlDict,Histograms):
200    'Needs a doc string'
201    FFtables = controlDict['FFtables']
202    if not FFtables:
203        return
204    histoList = Histograms.keys()
205    histoList.sort()
206    for histogram in histoList:
207        if histogram[:4] in ['PWDR','HKLF']:
208            Histogram = Histograms[histogram]
209            hId = Histogram['hId']
210            hfx = ':%d:'%(hId)
211            keV = controlDict[hfx+'keV']
212            for El in FFtables:
213                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
214                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
215                FFtables[El][hfx+'FP'] = FP
216                FFtables[El][hfx+'FPP'] = FPP               
217           
218def GetPhaseNames(GPXfile):
219    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
220
221    :param str GPXfile: full .gpx file name
222    :return: list of phase names
223    '''
224    fl = open(GPXfile,'rb')
225    PhaseNames = []
226    while True:
227        try:
228            data = cPickle.load(fl)
229        except EOFError:
230            break
231        datum = data[0]
232        if 'Phases' == datum[0]:
233            for datus in data[1:]:
234                PhaseNames.append(datus[0])
235    fl.close()
236    return PhaseNames
237
238def GetAllPhaseData(GPXfile,PhaseName):
239    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
240
241    :param str GPXfile: full .gpx file name
242    :param str PhaseName: phase name
243    :return: phase dictionary
244    '''       
245    fl = open(GPXfile,'rb')
246    General = {}
247    Atoms = []
248    while True:
249        try:
250            data = cPickle.load(fl)
251        except EOFError:
252            break
253        datum = data[0]
254        if 'Phases' == datum[0]:
255            for datus in data[1:]:
256                if datus[0] == PhaseName:
257                    break
258    fl.close()
259    return datus[1]
260   
261def GetHistograms(GPXfile,hNames):
262    """ Returns a dictionary of histograms found in GSASII gpx file
263
264    :param str GPXfile: full .gpx file name
265    :param str hNames: list of histogram names
266    :return: dictionary of histograms (types = PWDR & HKLF)
267
268    """
269    fl = open(GPXfile,'rb')
270    Histograms = {}
271    while True:
272        try:
273            data = cPickle.load(fl)
274        except EOFError:
275            break
276        datum = data[0]
277        hist = datum[0]
278        if hist in hNames:
279            if 'PWDR' in hist[:4]:
280                PWDRdata = {}
281                try:
282                    PWDRdata.update(datum[1][0])        #weight factor
283                except ValueError:
284                    PWDRdata['wtFactor'] = 1.0          #patch
285                PWDRdata['Data'] = datum[1][1]          #powder data arrays
286                PWDRdata[data[2][0]] = data[2][1]       #Limits
287                PWDRdata[data[3][0]] = data[3][1]       #Background
288                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
289                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
290                try:
291                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
292                except IndexError:
293                    PWDRdata['Reflection Lists'] = {}
294   
295                Histograms[hist] = PWDRdata
296            elif 'HKLF' in hist[:4]:
297                HKLFdata = {}
298                try:
299                    HKLFdata.update(datum[1][0])        #weight factor
300                except ValueError:
301                    HKLFdata['wtFactor'] = 1.0          #patch
302                HKLFdata['Data'] = datum[1][1]
303                HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
304                HKLFdata['Reflection Lists'] = None
305                Histograms[hist] = HKLFdata           
306    fl.close()
307    return Histograms
308   
309def GetHistogramNames(GPXfile,hType):
310    """ Returns a list of histogram names found in GSASII gpx file
311
312    :param str GPXfile: full .gpx file name
313    :param str hNames: list of histogram names
314    :return: list of histogram names (types = PWDR & HKLF)
315
316    """
317    fl = open(GPXfile,'rb')
318    HistogramNames = []
319    while True:
320        try:
321            data = cPickle.load(fl)
322        except EOFError:
323            break
324        datum = data[0]
325        if datum[0][:4] in hType:
326            HistogramNames.append(datum[0])
327    fl.close()
328    return HistogramNames
329   
330def GetUsedHistogramsAndPhases(GPXfile):
331    ''' Returns all histograms that are found in any phase
332    and any phase that uses a histogram
333
334    :param str GPXfile: full .gpx file name
335    :return: (Histograms,Phases)
336
337     * Histograms = dictionary of histograms as {name:data,...}
338     * Phases = dictionary of phases that use histograms
339
340    '''
341    phaseNames = GetPhaseNames(GPXfile)
342    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
343    allHistograms = GetHistograms(GPXfile,histoList)
344    phaseData = {}
345    for name in phaseNames: 
346        phaseData[name] =  GetAllPhaseData(GPXfile,name)
347    Histograms = {}
348    Phases = {}
349    for phase in phaseData:
350        Phase = phaseData[phase]
351        if Phase['Histograms']:
352            if phase not in Phases:
353                pId = phaseNames.index(phase)
354                Phase['pId'] = pId
355                Phases[phase] = Phase
356            for hist in Phase['Histograms']:
357                if 'Use' not in Phase['Histograms'][hist]:      #patch
358                    Phase['Histograms'][hist]['Use'] = True         
359                if hist not in Histograms and Phase['Histograms'][hist]['Use']:
360                    Histograms[hist] = allHistograms[hist]
361                    hId = histoList.index(hist)
362                    Histograms[hist]['hId'] = hId
363    return Histograms,Phases
364   
365def getBackupName(GPXfile,makeBack):
366    '''
367    Get the name for the backup .gpx file name
368   
369    :param str GPXfile: full .gpx file name
370    :param bool makeBack: if True the name of a new file is returned, if
371      False the name of the last file that exists is returned
372    :returns: the name of a backup file
373   
374    '''
375    GPXpath,GPXname = ospath.split(GPXfile)
376    if GPXpath == '': GPXpath = '.'
377    Name = ospath.splitext(GPXname)[0]
378    files = os.listdir(GPXpath)
379    last = 0
380    for name in files:
381        name = name.split('.')
382        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
383            if makeBack:
384                last = max(last,int(name[1].strip('bak'))+1)
385            else:
386                last = max(last,int(name[1].strip('bak')))
387    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
388    return GPXback   
389       
390def GPXBackup(GPXfile,makeBack=True):
391    '''
392    makes a backup of the current .gpx file (?)
393   
394    :param str GPXfile: full .gpx file name
395    :param bool makeBack: if True (default), the backup is written to
396      a new file; if False, the last backup is overwritten
397    :returns: the name of the backup file that was written
398    '''
399    import distutils.file_util as dfu
400    GPXback = getBackupName(GPXfile,makeBack)
401    dfu.copy_file(GPXfile,GPXback)
402    return GPXback
403
404def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
405    ''' Updates gpxfile from all histograms that are found in any phase
406    and any phase that used a histogram. Also updates rigid body definitions.
407
408
409    :param str GPXfile: full .gpx file name
410    :param dict Histograms: dictionary of histograms as {name:data,...}
411    :param dict Phases: dictionary of phases that use histograms
412    :param dict RigidBodies: dictionary of rigid bodies
413    :param dict CovData: dictionary of refined variables, varyList, & covariance matrix
414    :param bool makeBack: True if new backup of .gpx file is to be made; else use the last one made
415
416    '''
417                       
418    GPXback = GPXBackup(GPXfile,makeBack)
419    print 'Read from file:',GPXback
420    print 'Save to file  :',GPXfile
421    infile = open(GPXback,'rb')
422    outfile = open(GPXfile,'wb')
423    while True:
424        try:
425            data = cPickle.load(infile)
426        except EOFError:
427            break
428        datum = data[0]
429#        print 'read: ',datum[0]
430        if datum[0] == 'Phases':
431            for iphase in range(len(data)):
432                if data[iphase][0] in Phases:
433                    phaseName = data[iphase][0]
434                    data[iphase][1].update(Phases[phaseName])
435        elif datum[0] == 'Covariance':
436            data[0][1] = CovData
437        elif datum[0] == 'Rigid bodies':
438            data[0][1] = RigidBodies
439        try:
440            histogram = Histograms[datum[0]]
441#            print 'found ',datum[0]
442            data[0][1][1] = histogram['Data']
443            for datus in data[1:]:
444#                print '    read: ',datus[0]
445                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
446                    datus[1] = histogram[datus[0]]
447        except KeyError:
448            pass
449                               
450        cPickle.dump(data,outfile,1)
451    infile.close()
452    outfile.close()
453    print 'GPX file save successful'
454   
455def SetSeqResult(GPXfile,Histograms,SeqResult):
456    '''
457    Needs doc string
458   
459    :param str GPXfile: full .gpx file name
460    '''
461    GPXback = GPXBackup(GPXfile)
462    print 'Read from file:',GPXback
463    print 'Save to file  :',GPXfile
464    infile = open(GPXback,'rb')
465    outfile = open(GPXfile,'wb')
466    while True:
467        try:
468            data = cPickle.load(infile)
469        except EOFError:
470            break
471        datum = data[0]
472        if datum[0] == 'Sequential results':
473            data[0][1] = SeqResult
474        try:
475            histogram = Histograms[datum[0]]
476            data[0][1][1] = histogram['Data']
477            for datus in data[1:]:
478                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
479                    datus[1] = histogram[datus[0]]
480        except KeyError:
481            pass
482                               
483        cPickle.dump(data,outfile,1)
484    infile.close()
485    outfile.close()
486    print 'GPX file save successful'
487                       
488def ShowBanner(pFile=None):
489    'Print authorship, copyright and citation notice'
490    print >>pFile,80*'*'
491    print >>pFile,'   General Structure Analysis System-II Crystal Structure Refinement'
492    print >>pFile,'              by Robert B. Von Dreele & Brian H. Toby'
493    print >>pFile,'                Argonne National Laboratory(C), 2010'
494    print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,' 
495    print >>pFile,'            as Operator of Argonne National Laboratory.'
496    print >>pFile,'                          Please cite:'
497    print >>pFile,'   B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
498
499    print >>pFile,80*'*','\n'
500
501def ShowControls(Controls,pFile=None):
502    'Print controls information'
503    print >>pFile,' Least squares controls:'
504    print >>pFile,' Refinement type: ',Controls['deriv type']
505    if 'Hessian' in Controls['deriv type']:
506        print >>pFile,' Maximum number of cycles:',Controls['max cyc']
507    else:
508        print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
509    print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
510   
511def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
512    'needs a doc string'
513#    if SGLaue in ['-1','2/m','mmm']:
514#        return                      #no Pawley symmetry required constraints
515    eqvDict = {}
516    for i,varyI in enumerate(pawleyVary):
517        eqvDict[varyI] = []
518        refI = int(varyI.split(':')[-1])
519        ih,ik,il = PawleyRef[refI][:3]
520        dspI = PawleyRef[refI][4]
521        for varyJ in pawleyVary[i+1:]:
522            refJ = int(varyJ.split(':')[-1])
523            jh,jk,jl = PawleyRef[refJ][:3]
524            dspJ = PawleyRef[refJ][4]
525            if SGLaue in ['4/m','4/mmm']:
526                isum = ih**2+ik**2
527                jsum = jh**2+jk**2
528                if abs(il) == abs(jl) and isum == jsum:
529                    eqvDict[varyI].append(varyJ) 
530            elif SGLaue in ['3R','3mR']:
531                isum = ih**2+ik**2+il**2
532                jsum = jh**2+jk**2*jl**2
533                isum2 = ih*ik+ih*il+ik*il
534                jsum2 = jh*jk+jh*jl+jk*jl
535                if isum == jsum and isum2 == jsum2:
536                    eqvDict[varyI].append(varyJ) 
537            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
538                isum = ih**2+ik**2+ih*ik
539                jsum = jh**2+jk**2+jh*jk
540                if abs(il) == abs(jl) and isum == jsum:
541                    eqvDict[varyI].append(varyJ) 
542            elif SGLaue in ['m3','m3m']:
543                isum = ih**2+ik**2+il**2
544                jsum = jh**2+jk**2+jl**2
545                if isum == jsum:
546                    eqvDict[varyI].append(varyJ)
547            elif abs(dspI-dspJ)/dspI < 1.e-4:
548                eqvDict[varyI].append(varyJ) 
549    for item in pawleyVary:
550        if eqvDict[item]:
551            for item2 in pawleyVary:
552                if item2 in eqvDict[item]:
553                    eqvDict[item2] = []
554            G2mv.StoreEquivalence(item,eqvDict[item])
555                   
556def cellVary(pfx,SGData): 
557    'needs a doc string'
558    if SGData['SGLaue'] in ['-1',]:
559        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
560    elif SGData['SGLaue'] in ['2/m',]:
561        if SGData['SGUniq'] == 'a':
562            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
563        elif SGData['SGUniq'] == 'b':
564            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
565        else:
566            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
567    elif SGData['SGLaue'] in ['mmm',]:
568        return [pfx+'A0',pfx+'A1',pfx+'A2']
569    elif SGData['SGLaue'] in ['4/m','4/mmm']:
570        return [pfx+'A0',pfx+'A2']
571    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
572        return [pfx+'A0',pfx+'A2']
573    elif SGData['SGLaue'] in ['3R', '3mR']:
574        return [pfx+'A0',pfx+'A3']                       
575    elif SGData['SGLaue'] in ['m3m','m3']:
576        return [pfx+'A0',]
577       
578################################################################################
579##### Rigid Body Models and not General.get('doPawley')
580################################################################################
581       
582def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
583    'needs a doc string'
584   
585    def PrintResRBModel(RBModel):
586        atNames = RBModel['atNames']
587        rbRef = RBModel['rbRef']
588        rbSeq = RBModel['rbSeq']
589        print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
590            'No. times used: ',RBModel['useCount']
591        print >>pFile,'    At name       x          y          z'
592        for name,xyz in zip(atNames,RBModel['rbXYZ']):
593            print >>pFile,%8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
594        print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
595            ' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
596        if rbSeq:
597            for i,rbseq in enumerate(rbSeq):
598                print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
599                    atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
600       
601    def PrintVecRBModel(RBModel):
602        rbRef = RBModel['rbRef']
603        atTypes = RBModel['rbTypes']
604        print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
605            'No. times used: ',RBModel['useCount']
606        for i in range(len(RBModel['VectMag'])):
607            print >>pFile,'Vector no.: ',i,' Magnitude: ', \
608                '%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
609            print >>pFile,'  No. Type     vx         vy         vz'
610            for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
611                print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
612        print >>pFile,'  No. Type      x          y          z'
613        for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
614            print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
615        print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
616            ' & atom ',rbRef[0],' -> atom ',rbRef[2]
617    rbVary = []
618    rbDict = {}
619    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
620    if len(rbIds['Vector']):
621        for irb,item in enumerate(rbIds['Vector']):
622            if rigidbodyDict['Vector'][item]['useCount']:
623                RBmags = rigidbodyDict['Vector'][item]['VectMag']
624                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
625                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
626                    pid = '::RBV;'+str(i)+':'+str(irb)
627                    rbDict[pid] = mag
628                    if ref:
629                        rbVary.append(pid)
630                if Print:
631                    print >>pFile,'\nVector rigid body model:'
632                    PrintVecRBModel(rigidbodyDict['Vector'][item])
633    if len(rbIds['Residue']):
634        for item in rbIds['Residue']:
635            if rigidbodyDict['Residue'][item]['useCount']:
636                if Print:
637                    print >>pFile,'\nResidue rigid body model:'
638                    PrintResRBModel(rigidbodyDict['Residue'][item])
639    return rbVary,rbDict
640   
641def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
642    'needs a doc string'
643   
644    def PrintRBVectandSig(VectRB,VectSig):
645        print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
646        namstr = '  names :'
647        valstr = '  values:'
648        sigstr = '  esds  :'
649        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
650            namstr += '%12s'%('Vect '+str(i))
651            valstr += '%12.4f'%(val)
652            if sig:
653                sigstr += '%12.4f'%(sig)
654            else:
655                sigstr += 12*' '
656        print >>pFile,namstr
657        print >>pFile,valstr
658        print >>pFile,sigstr       
659       
660    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
661    if not RBIds['Vector']:
662        return
663    for irb,item in enumerate(RBIds['Vector']):
664        if rigidbodyDict['Vector'][item]['useCount']:
665            VectSig = []
666            RBmags = rigidbodyDict['Vector'][item]['VectMag']
667            for i,mag in enumerate(RBmags):
668                name = '::RBV;'+str(i)+':'+str(irb)
669                mag = parmDict[name]
670                if name in sigDict:
671                    VectSig.append(sigDict[name])
672            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
673       
674################################################################################
675##### Phase data
676################################################################################       
677                   
678def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
679    'needs a doc string'
680           
681    def PrintFFtable(FFtable):
682        print >>pFile,'\n X-ray scattering factors:'
683        print >>pFile,'   Symbol     fa                                      fb                                      fc'
684        print >>pFile,99*'-'
685        for Ename in FFtable:
686            ffdata = FFtable[Ename]
687            fa = ffdata['fa']
688            fb = ffdata['fb']
689            print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
690                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
691               
692    def PrintBLtable(BLtable):
693        print >>pFile,'\n Neutron scattering factors:'
694        print >>pFile,'   Symbol   isotope       mass       b       resonant terms'
695        print >>pFile,99*'-'
696        for Ename in BLtable:
697            bldata = BLtable[Ename]
698            isotope = bldata[0]
699            mass = bldata[1][0]
700            blen = bldata[1][1]
701            bres = []
702            if len(bldata[1]) > 2:
703                bres = bldata[1][2:]
704            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
705            for item in bres:
706                line += '%10.5g'%(item)
707            print >>pFile,line
708           
709    def PrintRBObjects(resRBData,vecRBData):
710       
711        def PrintRBThermals():
712            tlstr = ['11','22','33','12','13','23']
713            sstr = ['12','13','21','23','31','32','AA','BB']
714            TLS = RB['ThermalMotion'][1]
715            TLSvar = RB['ThermalMotion'][2]
716            if 'T' in RB['ThermalMotion'][0]:
717                print >>pFile,'TLS data'
718                text = ''
719                for i in range(6):
720                    text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
721                print >>pFile,text
722                if 'L' in RB['ThermalMotion'][0]: 
723                    text = ''
724                    for i in range(6,12):
725                        text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
726                    print >>pFile,text
727                if 'S' in RB['ThermalMotion'][0]:
728                    text = ''
729                    for i in range(12,20):
730                        text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
731                    print >>pFile,text
732            if 'U' in RB['ThermalMotion'][0]:
733                print >>pFile,'Uiso data'
734                text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])           
735           
736        if len(resRBData):
737            for RB in resRBData:
738                Oxyz = RB['Orig'][0]
739                Qrijk = RB['Orient'][0]
740                Angle = 2.0*acosd(Qrijk[0])
741                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
742                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
743                print >>pFile,'Orientation angle,vector:',      \
744                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
745                Torsions = RB['Torsions']
746                if len(Torsions):
747                    text = 'Torsions: '
748                    for torsion in Torsions:
749                        text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
750                    print >>pFile,text
751                PrintRBThermals()
752        if len(vecRBData):
753            for RB in vecRBData:
754                Oxyz = RB['Orig'][0]
755                Qrijk = RB['Orient'][0]
756                Angle = 2.0*acosd(Qrijk[0])
757                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
758                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]           
759                print >>pFile,'Orientation angle,vector:',      \
760                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
761                PrintRBThermals()
762               
763    def PrintAtoms(General,Atoms):
764        cx,ct,cs,cia = General['AtomPtrs']
765        print >>pFile,'\n Atoms:'
766        line = '   name    type  refine?   x         y         z    '+ \
767            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
768        if General['Type'] == 'magnetic':
769            line += '   Mx     My     Mz'
770        elif General['Type'] == 'macromolecular':
771            line = ' res no residue chain'+line
772        print >>pFile,line
773        if General['Type'] == 'nuclear':
774            print >>pFile,135*'-'
775            for i,at in enumerate(Atoms):
776                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
777                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
778                if at[cia] == 'I':
779                    line += '%8.4f'%(at[cia+1])+48*' '
780                else:
781                    line += 8*' '
782                    for j in range(6):
783                        line += '%8.4f'%(at[cia+1+j])
784                print >>pFile,line
785        elif General['Type'] == 'macromolecular':
786            print >>pFile,135*'-'           
787            for i,at in enumerate(Atoms):
788                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
789                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
790                if at[cia] == 'I':
791                    line += '%8.4f'%(at[cia+1])+48*' '
792                else:
793                    line += 8*' '
794                    for j in range(6):
795                        line += '%8.4f'%(at[cia+1+j])
796                print >>pFile,line
797       
798    def PrintTexture(textureData):
799        topstr = '\n Spherical harmonics texture: Order:' + \
800            str(textureData['Order'])
801        if textureData['Order']:
802            print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
803        else:
804            print >>pFile,topstr
805            return
806        names = ['omega','chi','phi']
807        line = '\n'
808        for name in names:
809            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
810        print >>pFile,line
811        print >>pFile,'\n Texture coefficients:'
812        ptlbls = ' names :'
813        ptstr =  ' values:'
814        SHcoeff = textureData['SH Coeff'][1]
815        for item in SHcoeff:
816            ptlbls += '%12s'%(item)
817            ptstr += '%12.4f'%(SHcoeff[item]) 
818        print >>pFile,ptlbls
819        print >>pFile,ptstr
820       
821    def MakeRBParms(rbKey,phaseVary,phaseDict):
822        rbid = str(rbids.index(RB['RBId']))
823        pfxRB = pfx+'RB'+rbKey+'P'
824        pstr = ['x','y','z']
825        ostr = ['a','i','j','k']
826        for i in range(3):
827            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
828            phaseDict[name] = RB['Orig'][0][i]
829            if RB['Orig'][1]:
830                phaseVary += [name,]
831        pfxRB = pfx+'RB'+rbKey+'O'
832        for i in range(4):
833            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
834            phaseDict[name] = RB['Orient'][0][i]
835            if RB['Orient'][1] == 'AV' and i:
836                phaseVary += [name,]
837            elif RB['Orient'][1] == 'A' and not i:
838                phaseVary += [name,]
839           
840    def MakeRBThermals(rbKey,phaseVary,phaseDict):
841        rbid = str(rbids.index(RB['RBId']))
842        tlstr = ['11','22','33','12','13','23']
843        sstr = ['12','13','21','23','31','32','AA','BB']
844        if 'T' in RB['ThermalMotion'][0]:
845            pfxRB = pfx+'RB'+rbKey+'T'
846            for i in range(6):
847                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
848                phaseDict[name] = RB['ThermalMotion'][1][i]
849                if RB['ThermalMotion'][2][i]:
850                    phaseVary += [name,]
851        if 'L' in RB['ThermalMotion'][0]:
852            pfxRB = pfx+'RB'+rbKey+'L'
853            for i in range(6):
854                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
855                phaseDict[name] = RB['ThermalMotion'][1][i+6]
856                if RB['ThermalMotion'][2][i+6]:
857                    phaseVary += [name,]
858        if 'S' in RB['ThermalMotion'][0]:
859            pfxRB = pfx+'RB'+rbKey+'S'
860            for i in range(8):
861                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
862                phaseDict[name] = RB['ThermalMotion'][1][i+12]
863                if RB['ThermalMotion'][2][i+12]:
864                    phaseVary += [name,]
865        if 'U' in RB['ThermalMotion'][0]:
866            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
867            phaseDict[name] = RB['ThermalMotion'][1][0]
868            if RB['ThermalMotion'][2][0]:
869                phaseVary += [name,]
870               
871    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
872        rbid = str(rbids.index(RB['RBId']))
873        pfxRB = pfx+'RB'+rbKey+'Tr;'
874        for i,torsion in enumerate(RB['Torsions']):
875            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
876            phaseDict[name] = torsion[0]
877            if torsion[1]:
878                phaseVary += [name,]
879                   
880    if Print:
881        print  >>pFile,'\n Phases:'
882    phaseVary = []
883    phaseDict = {}
884    phaseConstr = {}
885    pawleyLookup = {}
886    FFtables = {}                   #scattering factors - xrays
887    BLtables = {}                   # neutrons
888    Natoms = {}
889    AtMults = {}
890    AtIA = {}
891    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
892    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
893    atomIndx = {}
894    for name in PhaseData:
895        General = PhaseData[name]['General']
896        pId = PhaseData[name]['pId']
897        pfx = str(pId)+'::'
898        FFtable = G2el.GetFFtable(General['AtomTypes'])
899        BLtable = G2el.GetBLtable(General)
900        FFtables.update(FFtable)
901        BLtables.update(BLtable)
902        Atoms = PhaseData[name]['Atoms']
903        AtLookup = G2mth.FillAtomLookUp(Atoms)
904        PawleyRef = PhaseData[name].get('Pawley ref',[])
905        SGData = General['SGData']
906        SGtext = G2spc.SGPrint(SGData)
907        cell = General['Cell']
908        A = G2lat.cell2A(cell[1:7])
909        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
910            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
911        if cell[0]:
912            phaseVary += cellVary(pfx,SGData)
913        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
914        if resRBData:
915            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
916            for iRB,RB in enumerate(resRBData):
917                MakeRBParms('R',phaseVary,phaseDict)
918                MakeRBThermals('R',phaseVary,phaseDict)
919                MakeRBTorsions('R',phaseVary,phaseDict)
920       
921        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
922        if vecRBData:
923            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
924            for iRB,RB in enumerate(vecRBData):
925                MakeRBParms('V',phaseVary,phaseDict)
926                MakeRBThermals('V',phaseVary,phaseDict)
927                   
928        Natoms[pfx] = 0
929        if Atoms and not General.get('doPawley'):
930            cx,ct,cs,cia = General['AtomPtrs']
931            if General['Type'] in ['nuclear','macromolecular']:
932                Natoms[pfx] = len(Atoms)
933                for i,at in enumerate(Atoms):
934                    atomIndx[at[-1]] = [pfx,i]      #lookup table for restraints
935                    phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
936                        pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
937                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
938                        pfx+'AI/A:'+str(i):at[cia],})
939                    if at[cia] == 'I':
940                        phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
941                    else:
942                        phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
943                            pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
944                    if 'F' in at[ct+1]:
945                        phaseVary.append(pfx+'Afrac:'+str(i))
946                    if 'X' in at[ct+1]:
947                        xId,xCoef = G2spc.GetCSxinel(at[cs])
948                        names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
949                        equivs = [[],[],[]]
950                        for j in range(3):
951                            if xId[j] > 0:                               
952                                phaseVary.append(names[j])
953                                equivs[xId[j]-1].append([names[j],xCoef[j]])
954                        for equiv in equivs:
955                            if len(equiv) > 1:
956                                name = equiv[0][0]
957                                for eqv in equiv[1:]:
958                                    G2mv.StoreEquivalence(name,(eqv,))
959                    if 'U' in at[ct+1]:
960                        if at[9] == 'I':
961                            phaseVary.append(pfx+'AUiso:'+str(i))
962                        else:
963                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
964                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
965                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
966                            equivs = [[],[],[],[],[],[]]
967                            for j in range(6):
968                                if uId[j] > 0:                               
969                                    phaseVary.append(names[j])
970                                    equivs[uId[j]-1].append([names[j],uCoef[j]])
971                            for equiv in equivs:
972                                if len(equiv) > 1:
973                                    name = equiv[0][0]
974                                    for eqv in equiv[1:]:
975                                        G2mv.StoreEquivalence(name,(eqv,))
976#            elif General['Type'] == 'magnetic':
977#            elif General['Type'] == 'macromolecular':
978            textureData = General['SH Texture']
979            if textureData['Order']:
980                phaseDict[pfx+'SHorder'] = textureData['Order']
981                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
982                for item in ['omega','chi','phi']:
983                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
984                    if textureData['Sample '+item][0]:
985                        phaseVary.append(pfx+'SH '+item)
986                for item in textureData['SH Coeff'][1]:
987                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
988                    if textureData['SH Coeff'][0]:
989                        phaseVary.append(pfx+item)
990               
991            if Print:
992                print >>pFile,'\n Phase name: ',General['Name']
993                print >>pFile,135*'-'
994                PrintFFtable(FFtable)
995                PrintBLtable(BLtable)
996                print >>pFile,''
997                for line in SGtext: print >>pFile,line
998                PrintRBObjects(resRBData,vecRBData)
999                PrintAtoms(General,Atoms)
1000                print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1001                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1002                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1003                PrintTexture(textureData)
1004                if name in RestraintDict:
1005                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1006                        textureData,RestraintDict[name],pFile)
1007                   
1008        elif PawleyRef:
1009            pawleyVary = []
1010            for i,refl in enumerate(PawleyRef):
1011                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1012                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1013                if refl[5]:
1014                    pawleyVary.append(pfx+'PWLref:'+str(i))
1015            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
1016            phaseVary += pawleyVary
1017               
1018    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1019   
1020def cellFill(pfx,SGData,parmDict,sigDict): 
1021    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1022    from the parameter and sig dictionaries.
1023
1024    :param str pfx: parameter prefix ("n::", where n is a phase number)
1025    :param dict SGdata: a symmetry object
1026    :param dict parmDict: a dictionary of parameters
1027    :param dict sigDict:  a dictionary of uncertainties on parameters
1028
1029    :returns: A,sigA where each is a list of six terms with the A terms
1030    '''
1031    if SGData['SGLaue'] in ['-1',]:
1032        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1033            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1034    elif SGData['SGLaue'] in ['2/m',]:
1035        if SGData['SGUniq'] == 'a':
1036            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1037                parmDict[pfx+'A3'],0,0]
1038        elif SGData['SGUniq'] == 'b':
1039            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1040                0,parmDict[pfx+'A4'],0]
1041        else:
1042            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1043                0,0,parmDict[pfx+'A5']]
1044    elif SGData['SGLaue'] in ['mmm',]:
1045        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1046    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1047        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1048    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1049        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1050            parmDict[pfx+'A0'],0,0]
1051    elif SGData['SGLaue'] in ['3R', '3mR']:
1052        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1053            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1054    elif SGData['SGLaue'] in ['m3m','m3']:
1055        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1056
1057    try:
1058        if SGData['SGLaue'] in ['-1',]:
1059            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1060                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1061        elif SGData['SGLaue'] in ['2/m',]:
1062            if SGData['SGUniq'] == 'a':
1063                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1064                    sigDict[pfx+'A3'],0,0]
1065            elif SGData['SGUniq'] == 'b':
1066                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1067                    0,sigDict[pfx+'A4'],0]
1068            else:
1069                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1070                    0,0,sigDict[pfx+'A5']]
1071        elif SGData['SGLaue'] in ['mmm',]:
1072            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1073        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1074            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1075        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1076            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1077        elif SGData['SGLaue'] in ['3R', '3mR']:
1078            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1079        elif SGData['SGLaue'] in ['m3m','m3']:
1080            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1081    except KeyError:
1082        sigA = [0,0,0,0,0,0]
1083
1084    return A,sigA
1085       
1086def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1087    'needs a doc string'
1088    if phaseRest:
1089        Amat = G2lat.cell2AB(cell)[0]
1090        cx,ct,cs = AtPtrs[:3]
1091        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1092            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1093            ['ChemComp','Sites'],['Texture','HKLs']]
1094        for name,rest in names:
1095            itemRest = phaseRest[name]
1096            if itemRest[rest] and itemRest['Use']:
1097                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1098                if name in ['Bond','Angle','Plane','Chiral']:
1099                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1100                    for indx,ops,obs,esd in itemRest[rest]:
1101                        try:
1102                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1103                            AtName = ''
1104                            for i,Aname in enumerate(AtNames):
1105                                AtName += Aname
1106                                if ops[i] == '1':
1107                                    AtName += '-'
1108                                else:
1109                                    AtName += '+('+ops[i]+')-'
1110                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1111                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1112                            if name == 'Bond':
1113                                calc = G2mth.getRestDist(XYZ,Amat)
1114                            elif name == 'Angle':
1115                                calc = G2mth.getRestAngle(XYZ,Amat)
1116                            elif name == 'Plane':
1117                                calc = G2mth.getRestPlane(XYZ,Amat)
1118                            elif name == 'Chiral':
1119                                calc = G2mth.getRestChiral(XYZ,Amat)
1120                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1121                        except KeyError:
1122                            del itemRest[rest]
1123                elif name in ['Torsion','Rama']:
1124                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1125                    coeffDict = itemRest['Coeff']
1126                    for indx,ops,cofName,esd in itemRest[rest]:
1127                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1128                        AtName = ''
1129                        for i,Aname in enumerate(AtNames):
1130                            AtName += Aname+'+('+ops[i]+')-'
1131                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1132                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1133                        if name == 'Torsion':
1134                            tor = G2mth.getRestTorsion(XYZ,Amat)
1135                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1136                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1137                        else:
1138                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1139                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1140                            print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])
1141                elif name == 'ChemComp':
1142                    print >>pFile,'     atoms   mul*frac  factor     prod'
1143                    for indx,factors,obs,esd in itemRest[rest]:
1144                        try:
1145                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1146                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1147                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1148                            mulfrac = mul*frac
1149                            calcs = mul*frac*factors
1150                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1151                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1152                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1153                        except KeyError:
1154                            del itemRest[rest]
1155                elif name == 'Texture' and textureData['Order']:
1156                    Start = False
1157                    SHCoef = textureData['SH Coeff'][1]
1158                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1159                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1160                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1161                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1162                        PH = np.array(hkl)
1163                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1164                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1165                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1166                        Z = ma.masked_greater(Z,0.0)
1167                        num = ma.count(Z)
1168                        sum = 0
1169                        if num:
1170                            sum = np.sum(Z)
1171                        print '   %d %d %d  %d %8.3f %8.3f %8d   %s    %8.3f'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2)
1172       
1173def getCellEsd(pfx,SGData,A,covData):
1174    'needs a doc string'
1175    dpr = 180./np.pi
1176    rVsq = G2lat.calc_rVsq(A)
1177    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1178    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1179    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1180    scos = cosd(cellst[3:6])
1181    ssin = sind(cellst[3:6])
1182    scot = scos/ssin
1183    rcos = cosd(cell[3:6])
1184    rsin = sind(cell[3:6])
1185    rcot = rcos/rsin
1186    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1187    varyList = covData['varyList']
1188    covMatrix = covData['covMatrix']
1189    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1190    Ax = np.array(A)
1191    Ax[3:] /= 2.
1192    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
1193        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1194    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1195    Vol = 1/np.sqrt(rVsq)
1196    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
1197    R123 = Ax[0]*Ax[1]*Ax[2]
1198    dsasdg = np.zeros((3,6))
1199    dadg = np.zeros((6,6))
1200    for i0 in range(3):         #0  1   2
1201        i1 = (i0+1)%3           #1  2   0
1202        i2 = (i1+1)%3           #2  0   1
1203        i3 = 5-i2               #3  5   4
1204        i4 = 5-i1               #4  3   5
1205        i5 = 5-i0               #5  4   3
1206        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1207        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1208        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1209        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1210        denom = np.sqrt(denmsq)
1211        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1212        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1213        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1214        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1215        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1216        dadg[i5][i5] = -Ax[i0]/denom
1217    for i0 in range(3):
1218        i1 = (i0+1)%3
1219        i2 = (i1+1)%3
1220        i3 = 5-i2
1221        for ij in range(6):
1222            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1223            if ij == i0:
1224                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1225            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1226    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1227    var = np.diag(sigMat)
1228    CS = np.where(var>0.,np.sqrt(var),0.)
1229    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
1230    return cellSig           
1231   
1232def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1233    'needs a doc string'
1234   
1235    def PrintAtomsAndSig(General,Atoms,atomsSig):
1236        print >>pFile,'\n Atoms:'
1237        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1238        if General['Type'] == 'magnetic':
1239            line += '   Mx     My     Mz'
1240        elif General['Type'] == 'macromolecular':
1241            line = ' res no  residue  chain '+line
1242        print >>pFile,line
1243        if General['Type'] == 'nuclear':
1244            print >>pFile,135*'-'
1245            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1246                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1247            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
1248            for atyp in General['AtomTypes']:       #zero composition
1249                General['NoAtoms'][atyp] = 0.
1250            for i,at in enumerate(Atoms):
1251                General['NoAtoms'][at[1]] += at[6]*float(at[8])     #new composition
1252                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1253                valstr = ' values:'
1254                sigstr = ' sig   :'
1255                for ind in [3,4,5,6]:
1256                    sigind = str(i)+':'+str(ind)
1257                    valstr += fmt[ind]%(at[ind])                   
1258                    if sigind in atomsSig:
1259                        sigstr += fmt[ind]%(atomsSig[sigind])
1260                    else:
1261                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1262                if at[9] == 'I':
1263                    valstr += fmt[10]%(at[10])
1264                    if str(i)+':10' in atomsSig:
1265                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1266                    else:
1267                        sigstr += 8*' '
1268                else:
1269                    valstr += 8*' '
1270                    sigstr += 8*' '
1271                    for ind in [11,12,13,14,15,16]:
1272                        sigind = str(i)+':'+str(ind)
1273                        valstr += fmt[ind]%(at[ind])
1274                        if sigind in atomsSig:                       
1275                            sigstr += fmt[ind]%(atomsSig[sigind])
1276                        else:
1277                            sigstr += 8*' '
1278                print >>pFile,name
1279                print >>pFile,valstr
1280                print >>pFile,sigstr
1281               
1282    def PrintRBObjPOAndSig(rbfx,rbsx):
1283        namstr = '  names :'
1284        valstr = '  values:'
1285        sigstr = '  esds  :'
1286        for i,px in enumerate(['Px:','Py:','Pz:']):
1287            name = pfx+rbfx+px+rbsx
1288            namstr += '%12s'%('Pos '+px[1])
1289            valstr += '%12.5f'%(parmDict[name])
1290            if name in sigDict:
1291                sigstr += '%12.5f'%(sigDict[name])
1292            else:
1293                sigstr += 12*' '
1294        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1295            name = pfx+rbfx+po+rbsx
1296            namstr += '%12s'%('Ori '+po[1])
1297            valstr += '%12.5f'%(parmDict[name])
1298            if name in sigDict:
1299                sigstr += '%12.5f'%(sigDict[name])
1300            else:
1301                sigstr += 12*' '
1302        print >>pFile,namstr
1303        print >>pFile,valstr
1304        print >>pFile,sigstr
1305       
1306    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1307        namstr = '  names :'
1308        valstr = '  values:'
1309        sigstr = '  esds  :'
1310        if 'N' not in TLS:
1311            print >>pFile,' Thermal motion:'
1312        if 'T' in TLS:
1313            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1314                name = pfx+rbfx+pt+rbsx
1315                namstr += '%12s'%(pt[:3])
1316                valstr += '%12.5f'%(parmDict[name])
1317                if name in sigDict:
1318                    sigstr += '%12.5f'%(sigDict[name])
1319                else:
1320                    sigstr += 12*' '
1321            print >>pFile,namstr
1322            print >>pFile,valstr
1323            print >>pFile,sigstr
1324        if 'L' in TLS:
1325            namstr = '  names :'
1326            valstr = '  values:'
1327            sigstr = '  esds  :'
1328            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1329                name = pfx+rbfx+pt+rbsx
1330                namstr += '%12s'%(pt[:3])
1331                valstr += '%12.3f'%(parmDict[name])
1332                if name in sigDict:
1333                    sigstr += '%12.3f'%(sigDict[name])
1334                else:
1335                    sigstr += 12*' '
1336            print >>pFile,namstr
1337            print >>pFile,valstr
1338            print >>pFile,sigstr
1339        if 'S' in TLS:
1340            namstr = '  names :'
1341            valstr = '  values:'
1342            sigstr = '  esds  :'
1343            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1344                name = pfx+rbfx+pt+rbsx
1345                namstr += '%12s'%(pt[:3])
1346                valstr += '%12.4f'%(parmDict[name])
1347                if name in sigDict:
1348                    sigstr += '%12.4f'%(sigDict[name])
1349                else:
1350                    sigstr += 12*' '
1351            print >>pFile,namstr
1352            print >>pFile,valstr
1353            print >>pFile,sigstr
1354        if 'U' in TLS:
1355            name = pfx+rbfx+'U:'+rbsx
1356            namstr = '  names :'+'%12s'%('Uiso')
1357            valstr = '  values:'+'%12.5f'%(parmDict[name])
1358            if name in sigDict:
1359                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1360            else:
1361                sigstr = '  esds  :'+12*' '
1362            print >>pFile,namstr
1363            print >>pFile,valstr
1364            print >>pFile,sigstr
1365       
1366    def PrintRBObjTorAndSig(rbsx):
1367        namstr = '  names :'
1368        valstr = '  values:'
1369        sigstr = '  esds  :'
1370        nTors = len(RBObj['Torsions'])
1371        if nTors:
1372            print >>pFile,' Torsions:'
1373            for it in range(nTors):
1374                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1375                namstr += '%12s'%('Tor'+str(it))
1376                valstr += '%12.4f'%(parmDict[name])
1377                if name in sigDict:
1378                    sigstr += '%12.4f'%(sigDict[name])
1379            print >>pFile,namstr
1380            print >>pFile,valstr
1381            print >>pFile,sigstr
1382               
1383    def PrintSHtextureAndSig(textureData,SHtextureSig):
1384        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1385        names = ['omega','chi','phi']
1386        namstr = '  names :'
1387        ptstr =  '  values:'
1388        sigstr = '  esds  :'
1389        for name in names:
1390            namstr += '%12s'%(name)
1391            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1392            if 'Sample '+name in SHtextureSig:
1393                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1394            else:
1395                sigstr += 12*' '
1396        print >>pFile,namstr
1397        print >>pFile,ptstr
1398        print >>pFile,sigstr
1399        print >>pFile,'\n Texture coefficients:'
1400        namstr = '  names :'
1401        ptstr =  '  values:'
1402        sigstr = '  esds  :'
1403        SHcoeff = textureData['SH Coeff'][1]
1404        for name in SHcoeff:
1405            namstr += '%12s'%(name)
1406            ptstr += '%12.3f'%(SHcoeff[name])
1407            if name in SHtextureSig:
1408                sigstr += '%12.3f'%(SHtextureSig[name])
1409            else:
1410                sigstr += 12*' '
1411        print >>pFile,namstr
1412        print >>pFile,ptstr
1413        print >>pFile,sigstr
1414           
1415    print >>pFile,'\n Phases:'
1416    for phase in Phases:
1417        print >>pFile,' Result for phase: ',phase
1418        Phase = Phases[phase]
1419        General = Phase['General']
1420        SGData = General['SGData']
1421        Atoms = Phase['Atoms']
1422        AtLookup = G2mth.FillAtomLookUp(Atoms)
1423        cell = General['Cell']
1424        pId = Phase['pId']
1425        pfx = str(pId)+'::'
1426        if cell[0]:
1427            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1428            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1429            print >>pFile,' Reciprocal metric tensor: '
1430            ptfmt = "%15.9f"
1431            names = ['A11','A22','A33','A12','A13','A23']
1432            namstr = '  names :'
1433            ptstr =  '  values:'
1434            sigstr = '  esds  :'
1435            for name,a,siga in zip(names,A,sigA):
1436                namstr += '%15s'%(name)
1437                ptstr += ptfmt%(a)
1438                if siga:
1439                    sigstr += ptfmt%(siga)
1440                else:
1441                    sigstr += 15*' '
1442            print >>pFile,namstr
1443            print >>pFile,ptstr
1444            print >>pFile,sigstr
1445            cell[1:7] = G2lat.A2cell(A)
1446            cell[7] = G2lat.calc_V(A)
1447            print >>pFile,' New unit cell:'
1448            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1449            names = ['a','b','c','alpha','beta','gamma','Volume']
1450            namstr = '  names :'
1451            ptstr =  '  values:'
1452            sigstr = '  esds  :'
1453            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1454                namstr += '%12s'%(name)
1455                ptstr += fmt%(a)
1456                if siga:
1457                    sigstr += fmt%(siga)
1458                else:
1459                    sigstr += 12*' '
1460            print >>pFile,namstr
1461            print >>pFile,ptstr
1462            print >>pFile,sigstr
1463           
1464        General['Mass'] = 0.
1465        if Phase['General'].get('doPawley'):
1466            pawleyRef = Phase['Pawley ref']
1467            for i,refl in enumerate(pawleyRef):
1468                key = pfx+'PWLref:'+str(i)
1469                refl[6] = parmDict[key]
1470                if key in sigDict:
1471                    refl[7] = sigDict[key]
1472                else:
1473                    refl[7] = 0
1474        else:
1475            VRBIds = RBIds['Vector']
1476            RRBIds = RBIds['Residue']
1477            RBModels = Phase['RBModels']
1478            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1479                jrb = VRBIds.index(RBObj['RBId'])
1480                rbsx = str(irb)+':'+str(jrb)
1481                print >>pFile,' Vector rigid body parameters:'
1482                PrintRBObjPOAndSig('RBV',rbsx)
1483                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1484            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1485                jrb = RRBIds.index(RBObj['RBId'])
1486                rbsx = str(irb)+':'+str(jrb)
1487                print >>pFile,' Residue rigid body parameters:'
1488                PrintRBObjPOAndSig('RBR',rbsx)
1489                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1490                PrintRBObjTorAndSig(rbsx)
1491            atomsSig = {}
1492            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
1493                for i,at in enumerate(Atoms):
1494                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1495                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1496                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1497                    for ind in [3,4,5,6]:
1498                        at[ind] = parmDict[names[ind]]
1499                        if ind in [3,4,5]:
1500                            name = names[ind].replace('A','dA')
1501                        else:
1502                            name = names[ind]
1503                        if name in sigDict:
1504                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1505                    if at[9] == 'I':
1506                        at[10] = parmDict[names[10]]
1507                        if names[10] in sigDict:
1508                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1509                    else:
1510                        for ind in [11,12,13,14,15,16]:
1511                            at[ind] = parmDict[names[ind]]
1512                            if names[ind] in sigDict:
1513                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1514                    ind = General['AtomTypes'].index(at[1])
1515                    General['Mass'] += General['AtomMass'][ind]*at[6]*at[8]
1516            PrintAtomsAndSig(General,Atoms,atomsSig)
1517       
1518        textureData = General['SH Texture']   
1519        if textureData['Order']:
1520            SHtextureSig = {}
1521            for name in ['omega','chi','phi']:
1522                aname = pfx+'SH '+name
1523                textureData['Sample '+name][1] = parmDict[aname]
1524                if aname in sigDict:
1525                    SHtextureSig['Sample '+name] = sigDict[aname]
1526            for name in textureData['SH Coeff'][1]:
1527                aname = pfx+name
1528                textureData['SH Coeff'][1][name] = parmDict[aname]
1529                if aname in sigDict:
1530                    SHtextureSig[name] = sigDict[aname]
1531            PrintSHtextureAndSig(textureData,SHtextureSig)
1532        if phase in RestraintDict:
1533            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1534                textureData,RestraintDict[phase],pFile)
1535                   
1536################################################################################
1537##### Histogram & Phase data
1538################################################################################       
1539                   
1540def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1541    '''Loads the HAP histogram/phase information into dicts
1542
1543    :param dict Phases: phase information
1544    :param dict Histograms: Histogram information
1545    :param bool Print: prints information as it is read
1546    :param file pFile: file object to print to (the default, None causes printing to the console)
1547    :param bool resetRefList: Should the contents of the Reflection List be initialized
1548      on loading. The default, True, initializes the Reflection List as it is loaded.
1549
1550    :returns: (hapVary,hapDict,controlDict)
1551      * hapVary: list of refined variables
1552      * hapDict: dict with refined variables and their values
1553      * controlDict: dict with computation controls (?)
1554    '''
1555   
1556    def PrintSize(hapData):
1557        if hapData[0] in ['isotropic','uniaxial']:
1558            line = '\n Size model    : %9s'%(hapData[0])
1559            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1560            if hapData[0] == 'uniaxial':
1561                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1562            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1563            print >>pFile,line
1564        else:
1565            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
1566                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1567            Snames = ['S11','S22','S33','S12','S13','S23']
1568            ptlbls = ' names :'
1569            ptstr =  ' values:'
1570            varstr = ' refine:'
1571            for i,name in enumerate(Snames):
1572                ptlbls += '%12s' % (name)
1573                ptstr += '%12.6f' % (hapData[4][i])
1574                varstr += '%12s' % (str(hapData[5][i]))
1575            print >>pFile,ptlbls
1576            print >>pFile,ptstr
1577            print >>pFile,varstr
1578       
1579    def PrintMuStrain(hapData,SGData):
1580        if hapData[0] in ['isotropic','uniaxial']:
1581            line = '\n Mustrain model: %9s'%(hapData[0])
1582            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1583            if hapData[0] == 'uniaxial':
1584                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1585            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1586            print >>pFile,line
1587        else:
1588            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1589                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1590            Snames = G2spc.MustrainNames(SGData)
1591            ptlbls = ' names :'
1592            ptstr =  ' values:'
1593            varstr = ' refine:'
1594            for i,name in enumerate(Snames):
1595                ptlbls += '%12s' % (name)
1596                ptstr += '%12.6f' % (hapData[4][i])
1597                varstr += '%12s' % (str(hapData[5][i]))
1598            print >>pFile,ptlbls
1599            print >>pFile,ptstr
1600            print >>pFile,varstr
1601
1602    def PrintHStrain(hapData,SGData):
1603        print >>pFile,'\n Hydrostatic/elastic strain: '
1604        Hsnames = G2spc.HStrainNames(SGData)
1605        ptlbls = ' names :'
1606        ptstr =  ' values:'
1607        varstr = ' refine:'
1608        for i,name in enumerate(Hsnames):
1609            ptlbls += '%12s' % (name)
1610            ptstr += '%12.6f' % (hapData[0][i])
1611            varstr += '%12s' % (str(hapData[1][i]))
1612        print >>pFile,ptlbls
1613        print >>pFile,ptstr
1614        print >>pFile,varstr
1615
1616    def PrintSHPO(hapData):
1617        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1618            str(hapData[4])+' Refine? '+str(hapData[2])
1619        ptlbls = ' names :'
1620        ptstr =  ' values:'
1621        for item in hapData[5]:
1622            ptlbls += '%12s'%(item)
1623            ptstr += '%12.3f'%(hapData[5][item]) 
1624        print >>pFile,ptlbls
1625        print >>pFile,ptstr
1626   
1627    def PrintBabinet(hapData):
1628        print >>pFile,'\n Babinet form factor modification: '
1629        ptlbls = ' names :'
1630        ptstr =  ' values:'
1631        varstr = ' refine:'
1632        for name in ['BabA','BabU']:
1633            ptlbls += '%12s' % (name)
1634            ptstr += '%12.6f' % (hapData[name][0])
1635            varstr += '%12s' % (str(hapData[name][1]))
1636        print >>pFile,ptlbls
1637        print >>pFile,ptstr
1638        print >>pFile,varstr
1639       
1640   
1641    hapDict = {}
1642    hapVary = []
1643    controlDict = {}
1644    poType = {}
1645    poAxes = {}
1646    spAxes = {}
1647    spType = {}
1648   
1649    for phase in Phases:
1650        HistoPhase = Phases[phase]['Histograms']
1651        SGData = Phases[phase]['General']['SGData']
1652        cell = Phases[phase]['General']['Cell'][1:7]
1653        A = G2lat.cell2A(cell)
1654        pId = Phases[phase]['pId']
1655        histoList = HistoPhase.keys()
1656        histoList.sort()
1657        for histogram in histoList:
1658            try:
1659                Histogram = Histograms[histogram]
1660            except KeyError:                       
1661                #skip if histogram not included e.g. in a sequential refinement
1662                continue
1663            hapData = HistoPhase[histogram]
1664            hId = Histogram['hId']
1665            if 'PWDR' in histogram:
1666                limits = Histogram['Limits'][1]
1667                inst = Histogram['Instrument Parameters'][0]
1668                Zero = inst['Zero'][1]
1669                if 'C' in inst['Type'][1]:
1670                    try:
1671                        wave = inst['Lam'][1]
1672                    except KeyError:
1673                        wave = inst['Lam1'][1]
1674                    dmin = wave/(2.0*sind(limits[1]/2.0))
1675                pfx = str(pId)+':'+str(hId)+':'
1676                for item in ['Scale','Extinction']:
1677                    hapDict[pfx+item] = hapData[item][0]
1678                    if hapData[item][1]:
1679                        hapVary.append(pfx+item)
1680                names = G2spc.HStrainNames(SGData)
1681                for i,name in enumerate(names):
1682                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1683                    if hapData['HStrain'][1][i]:
1684                        hapVary.append(pfx+name)
1685                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1686                if hapData['Pref.Ori.'][0] == 'MD':
1687                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1688                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1689                    if hapData['Pref.Ori.'][2]:
1690                        hapVary.append(pfx+'MD')
1691                else:                           #'SH' spherical harmonics
1692                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1693                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1694                    for item in hapData['Pref.Ori.'][5]:
1695                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1696                        if hapData['Pref.Ori.'][2]:
1697                            hapVary.append(pfx+item)
1698                for item in ['Mustrain','Size']:
1699                    controlDict[pfx+item+'Type'] = hapData[item][0]
1700                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
1701                    if hapData[item][2][2]:
1702                        hapVary.append(pfx+item+';mx')
1703                    if hapData[item][0] in ['isotropic','uniaxial']:
1704                        hapDict[pfx+item+';i'] = hapData[item][1][0]
1705                        if hapData[item][2][0]:
1706                            hapVary.append(pfx+item+';i')
1707                        if hapData[item][0] == 'uniaxial':
1708                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1709                            hapDict[pfx+item+';a'] = hapData[item][1][1]
1710                            if hapData[item][2][1]:
1711                                hapVary.append(pfx+item+';a')
1712                    else:       #generalized for mustrain or ellipsoidal for size
1713                        Nterms = len(hapData[item][4])
1714                        if item == 'Mustrain':
1715                            names = G2spc.MustrainNames(SGData)
1716                            pwrs = []
1717                            for name in names:
1718                                h,k,l = name[1:]
1719                                pwrs.append([int(h),int(k),int(l)])
1720                            controlDict[pfx+'MuPwrs'] = pwrs
1721                        for i in range(Nterms):
1722                            sfx = ':'+str(i)
1723                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1724                            if hapData[item][5][i]:
1725                                hapVary.append(pfx+item+sfx)
1726                for bab in ['BabA','BabU']:
1727                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1728                    if hapData['Babinet'][bab][1]:
1729                        hapVary.append(pfx+bab)
1730                               
1731                if Print: 
1732                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1733                    print >>pFile,135*'-'
1734                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1735                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1736                    if hapData['Pref.Ori.'][0] == 'MD':
1737                        Ax = hapData['Pref.Ori.'][3]
1738                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1739                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1740                    else: #'SH' for spherical harmonics
1741                        PrintSHPO(hapData['Pref.Ori.'])
1742                    PrintSize(hapData['Size'])
1743                    PrintMuStrain(hapData['Mustrain'],SGData)
1744                    PrintHStrain(hapData['HStrain'],SGData)
1745                    if hapData['Babinet']['BabA'][0]:
1746                        PrintBabinet(hapData['Babinet'])
1747                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1748                refList = []
1749                for h,k,l,d in HKLd:
1750                    ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1751                    mul *= 2      # for powder overlap of Friedel pairs
1752                    if ext:
1753                        continue
1754                    if 'C' in inst['Type'][0]:
1755                        pos = 2.0*asind(wave/(2.0*d))+Zero
1756                        if limits[0] < pos < limits[1]:
1757                            refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1758                            #last item should contain form factor values by atom type
1759                    else:
1760                        raise ValueError 
1761                if resetRefList: Histogram['Reflection Lists'][phase] = refList
1762            elif 'HKLF' in histogram:
1763                inst = Histogram['Instrument Parameters'][0]
1764                hId = Histogram['hId']
1765                hfx = ':%d:'%(hId)
1766                for item in inst:
1767                    hapDict[hfx+item] = inst[item][1]
1768                pfx = str(pId)+':'+str(hId)+':'
1769                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1770                if hapData['Scale'][1]:
1771                    hapVary.append(pfx+'Scale')
1772                               
1773                extApprox,extType,extParms = hapData['Extinction']
1774                controlDict[pfx+'EType'] = extType
1775                controlDict[pfx+'EApprox'] = extApprox
1776                controlDict[pfx+'Tbar'] = extParms['Tbar']
1777                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1778                if 'Primary' in extType:
1779                    Ekey = ['Ep',]
1780                elif 'I & II' in extType:
1781                    Ekey = ['Eg','Es']
1782                elif 'Secondary Type II' == extType:
1783                    Ekey = ['Es',]
1784                elif 'Secondary Type I' == extType:
1785                    Ekey = ['Eg',]
1786                else:   #'None'
1787                    Ekey = []
1788                for eKey in Ekey:
1789                    hapDict[pfx+eKey] = extParms[eKey][0]
1790                    if extParms[eKey][1]:
1791                        hapVary.append(pfx+eKey)
1792                for bab in ['BabA','BabU']:
1793                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1794                    if hapData['Babinet'][bab][1]:
1795                        hapVary.append(pfx+bab)
1796                if Print: 
1797                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1798                    print >>pFile,135*'-'
1799                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1800                    if extType != 'None':
1801                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1802                        text = ' Parameters       :'
1803                        for eKey in Ekey:
1804                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1805                        print >>pFile,text
1806                    if hapData['Babinet']['BabA'][0]:
1807                        PrintBabinet(hapData['Babinet'])
1808                Histogram['Reflection Lists'] = phase       
1809               
1810    return hapVary,hapDict,controlDict
1811   
1812def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1813    'needs a doc string'
1814   
1815    def PrintSizeAndSig(hapData,sizeSig):
1816        line = '\n Size model:     %9s'%(hapData[0])
1817        refine = False
1818        if hapData[0] in ['isotropic','uniaxial']:
1819            line += ' equatorial:%12.4f'%(hapData[1][0])
1820            if sizeSig[0][0]:
1821                line += ', sig:%8.4f'%(sizeSig[0][0])
1822                refine = True
1823            if hapData[0] == 'uniaxial':
1824                line += ' axial:%12.4f'%(hapData[1][1])
1825                if sizeSig[0][1]:
1826                    refine = True
1827                    line += ', sig:%8.4f'%(sizeSig[0][1])
1828            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1829            if sizeSig[0][2]:
1830                refine = True
1831                line += ', sig:%8.4f'%(sizeSig[0][2])
1832            if refine:
1833                print >>pFile,line
1834        else:
1835            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1836            if sizeSig[0][2]:
1837                refine = True
1838                line += ', sig:%8.4f'%(sizeSig[0][2])
1839            Snames = ['S11','S22','S33','S12','S13','S23']
1840            ptlbls = ' name  :'
1841            ptstr =  ' value :'
1842            sigstr = ' sig   :'
1843            for i,name in enumerate(Snames):
1844                ptlbls += '%12s' % (name)
1845                ptstr += '%12.6f' % (hapData[4][i])
1846                if sizeSig[1][i]:
1847                    refine = True
1848                    sigstr += '%12.6f' % (sizeSig[1][i])
1849                else:
1850                    sigstr += 12*' '
1851            if refine:
1852                print >>pFile,line
1853                print >>pFile,ptlbls
1854                print >>pFile,ptstr
1855                print >>pFile,sigstr
1856       
1857    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1858        line = '\n Mustrain model: %9s'%(hapData[0])
1859        refine = False
1860        if hapData[0] in ['isotropic','uniaxial']:
1861            line += ' equatorial:%12.1f'%(hapData[1][0])
1862            if mustrainSig[0][0]:
1863                line += ', sig:%8.1f'%(mustrainSig[0][0])
1864                refine = True
1865            if hapData[0] == 'uniaxial':
1866                line += ' axial:%12.1f'%(hapData[1][1])
1867                if mustrainSig[0][1]:
1868                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1869            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1870            if mustrainSig[0][2]:
1871                refine = True
1872                line += ', sig:%8.3f'%(mustrainSig[0][2])
1873            if refine:
1874                print >>pFile,line
1875        else:
1876            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1877            if mustrainSig[0][2]:
1878                refine = True
1879                line += ', sig:%8.3f'%(mustrainSig[0][2])
1880            Snames = G2spc.MustrainNames(SGData)
1881            ptlbls = ' name  :'
1882            ptstr =  ' value :'
1883            sigstr = ' sig   :'
1884            for i,name in enumerate(Snames):
1885                ptlbls += '%12s' % (name)
1886                ptstr += '%12.6f' % (hapData[4][i])
1887                if mustrainSig[1][i]:
1888                    refine = True
1889                    sigstr += '%12.6f' % (mustrainSig[1][i])
1890                else:
1891                    sigstr += 12*' '
1892            if refine:
1893                print >>pFile,line
1894                print >>pFile,ptlbls
1895                print >>pFile,ptstr
1896                print >>pFile,sigstr
1897           
1898    def PrintHStrainAndSig(hapData,strainSig,SGData):
1899        Hsnames = G2spc.HStrainNames(SGData)
1900        ptlbls = ' name  :'
1901        ptstr =  ' value :'
1902        sigstr = ' sig   :'
1903        refine = False
1904        for i,name in enumerate(Hsnames):
1905            ptlbls += '%12s' % (name)
1906            ptstr += '%12.6g' % (hapData[0][i])
1907            if name in strainSig:
1908                refine = True
1909                sigstr += '%12.6g' % (strainSig[name])
1910            else:
1911                sigstr += 12*' '
1912        if refine:
1913            print >>pFile,'\n Hydrostatic/elastic strain: '
1914            print >>pFile,ptlbls
1915            print >>pFile,ptstr
1916            print >>pFile,sigstr
1917       
1918    def PrintSHPOAndSig(pfx,hapData,POsig):
1919        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1920        ptlbls = ' names :'
1921        ptstr =  ' values:'
1922        sigstr = ' sig   :'
1923        for item in hapData[5]:
1924            ptlbls += '%12s'%(item)
1925            ptstr += '%12.3f'%(hapData[5][item])
1926            if pfx+item in POsig:
1927                sigstr += '%12.3f'%(POsig[pfx+item])
1928            else:
1929                sigstr += 12*' ' 
1930        print >>pFile,ptlbls
1931        print >>pFile,ptstr
1932        print >>pFile,sigstr
1933       
1934    def PrintExtAndSig(pfx,hapData,ScalExtSig):
1935        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1936        text = ''
1937        for item in hapData[2]:
1938            if pfx+item in ScalExtSig:
1939                text += '       %s: '%(item)
1940                text += '%12.2e'%(hapData[2][item][0])
1941                if pfx+item in ScalExtSig:
1942                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1943        print >>pFile,text       
1944       
1945    def PrintBabinetAndSig(pfx,hapData,BabSig):
1946        print >>pFile,'\n Babinet form factor modification: '
1947        ptlbls = ' names :'
1948        ptstr =  ' values:'
1949        sigstr = ' sig   :'
1950        for item in hapData:
1951            ptlbls += '%12s'%(item)
1952            ptstr += '%12.3f'%(hapData[item][0])
1953            if pfx+item in BabSig:
1954                sigstr += '%12.3f'%(BabSig[pfx+item])
1955            else:
1956                sigstr += 12*' ' 
1957        print >>pFile,ptlbls
1958        print >>pFile,ptstr
1959        print >>pFile,sigstr
1960   
1961    PhFrExtPOSig = {}
1962    SizeMuStrSig = {}
1963    ScalExtSig = {}
1964    BabSig = {}
1965    wtFrSum = {}
1966    for phase in Phases:
1967        HistoPhase = Phases[phase]['Histograms']
1968        General = Phases[phase]['General']
1969        SGData = General['SGData']
1970        pId = Phases[phase]['pId']
1971        histoList = HistoPhase.keys()
1972        histoList.sort()
1973        for histogram in histoList:
1974            try:
1975                Histogram = Histograms[histogram]
1976            except KeyError:                       
1977                #skip if histogram not included e.g. in a sequential refinement
1978                continue
1979            hapData = HistoPhase[histogram]
1980            hId = Histogram['hId']
1981            pfx = str(pId)+':'+str(hId)+':'
1982            if hId not in wtFrSum:
1983                wtFrSum[hId] = 0.
1984            if 'PWDR' in histogram:
1985                for item in ['Scale','Extinction']:
1986                    hapData[item][0] = parmDict[pfx+item]
1987                    if pfx+item in sigDict:
1988                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1989                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
1990                if hapData['Pref.Ori.'][0] == 'MD':
1991                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1992                    if pfx+'MD' in sigDict:
1993                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
1994                else:                           #'SH' spherical harmonics
1995                    for item in hapData['Pref.Ori.'][5]:
1996                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1997                        if pfx+item in sigDict:
1998                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1999                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2000                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2001                    pfx+'HStrain':{}})                 
2002                for item in ['Mustrain','Size']:
2003                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2004                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2005                    if pfx+item+';mx' in sigDict:
2006                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2007                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2008                        hapData[item][1][0] = parmDict[pfx+item+';i']
2009                        if item == 'Size':
2010                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2011                        if pfx+item+';i' in sigDict: 
2012                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2013                        if hapData[item][0] == 'uniaxial':
2014                            hapData[item][1][1] = parmDict[pfx+item+';a']
2015                            if item == 'Size':
2016                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2017                            if pfx+item+';a' in sigDict:
2018                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2019                    else:       #generalized for mustrain or ellipsoidal for size
2020                        Nterms = len(hapData[item][4])
2021                        for i in range(Nterms):
2022                            sfx = ':'+str(i)
2023                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2024                            if pfx+item+sfx in sigDict:
2025                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2026                names = G2spc.HStrainNames(SGData)
2027                for i,name in enumerate(names):
2028                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2029                    if pfx+name in sigDict:
2030                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2031                for name in ['BabA','BabU']:
2032                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2033                    if pfx+name in sigDict:
2034                        BabSig[pfx+name] = sigDict[pfx+name]               
2035               
2036            elif 'HKLF' in histogram:
2037                for item in ['Scale',]:
2038                    if parmDict.get(pfx+item):
2039                        hapData[item][0] = parmDict[pfx+item]
2040                        if pfx+item in sigDict:
2041                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2042                for item in ['Ep','Eg','Es']:
2043                    if parmDict.get(pfx+item):
2044                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2045                        if pfx+item in sigDict:
2046                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2047                for name in ['BabA','BabU']:
2048                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2049                    if pfx+name in sigDict:
2050                        BabSig[pfx+name] = sigDict[pfx+name]               
2051
2052    if Print:
2053        for phase in Phases:
2054            HistoPhase = Phases[phase]['Histograms']
2055            General = Phases[phase]['General']
2056            SGData = General['SGData']
2057            pId = Phases[phase]['pId']
2058            histoList = HistoPhase.keys()
2059            histoList.sort()
2060            for histogram in histoList:
2061                try:
2062                    Histogram = Histograms[histogram]
2063                except KeyError:                       
2064                    #skip if histogram not included e.g. in a sequential refinement
2065                    continue
2066                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2067                print >>pFile,130*'-'
2068                hapData = HistoPhase[histogram]
2069                hId = Histogram['hId']
2070                pfx = str(pId)+':'+str(hId)+':'
2071                if 'PWDR' in histogram:
2072                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2073                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2074               
2075                    if pfx+'Scale' in PhFrExtPOSig:
2076                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2077                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2078                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2079                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2080                    if pfx+'Extinction' in PhFrExtPOSig:
2081                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2082                    if hapData['Pref.Ori.'][0] == 'MD':
2083                        if pfx+'MD' in PhFrExtPOSig:
2084                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2085                    else:
2086                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2087                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2088                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2089                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2090                    if len(BabSig):
2091                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2092                   
2093                elif 'HKLF' in histogram:
2094                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2095                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2096                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2097                    if pfx+'Scale' in ScalExtSig:
2098                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2099                    if hapData['Extinction'][0] != 'None':
2100                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2101                    if len(BabSig):
2102                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2103
2104################################################################################
2105##### Histogram data
2106################################################################################       
2107                   
2108def GetHistogramData(Histograms,Print=True,pFile=None):
2109    'needs a doc string'
2110   
2111    def GetBackgroundParms(hId,Background):
2112        Back = Background[0]
2113        DebyePeaks = Background[1]
2114        bakType,bakFlag = Back[:2]
2115        backVals = Back[3:]
2116        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2117        backDict = dict(zip(backNames,backVals))
2118        backVary = []
2119        if bakFlag:
2120            backVary = backNames
2121        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2122        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2123        debyeDict = {}
2124        debyeList = []
2125        for i in range(DebyePeaks['nDebye']):
2126            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2127            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2128            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2129        debyeVary = []
2130        for item in debyeList:
2131            if item[1]:
2132                debyeVary.append(item[0])
2133        backDict.update(debyeDict)
2134        backVary += debyeVary
2135        peakDict = {}
2136        peakList = []
2137        for i in range(DebyePeaks['nPeaks']):
2138            peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
2139                ':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
2140            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2141            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2142        peakVary = []
2143        for item in peakList:
2144            if item[1]:
2145                peakVary.append(item[0])
2146        backDict.update(peakDict)
2147        backVary += peakVary
2148        return bakType,backDict,backVary           
2149       
2150    def GetInstParms(hId,Inst):     
2151        dataType = Inst['Type'][0]
2152        instDict = {}
2153        insVary = []
2154        pfx = ':'+str(hId)+':'
2155        insKeys = Inst.keys()
2156        insKeys.sort()
2157        for item in insKeys:
2158            insName = pfx+item
2159            instDict[insName] = Inst[item][1]
2160            if Inst[item][2]:
2161                insVary.append(insName)
2162#        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2163#        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2164        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2165        return dataType,instDict,insVary
2166       
2167    def GetSampleParms(hId,Sample):
2168        sampVary = []
2169        hfx = ':'+str(hId)+':'       
2170        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2171            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2172        Type = Sample['Type']
2173        if 'Bragg' in Type:             #Bragg-Brentano
2174            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
2175                sampDict[hfx+item] = Sample[item][0]
2176                if Sample[item][1]:
2177                    sampVary.append(hfx+item)
2178        elif 'Debye' in Type:        #Debye-Scherrer
2179            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2180                sampDict[hfx+item] = Sample[item][0]
2181                if Sample[item][1]:
2182                    sampVary.append(hfx+item)
2183        return Type,sampDict,sampVary
2184       
2185    def PrintBackground(Background):
2186        Back = Background[0]
2187        DebyePeaks = Background[1]
2188        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2189        line = ' Coefficients: '
2190        for i,back in enumerate(Back[3:]):
2191            line += '%10.3f'%(back)
2192            if i and not i%10:
2193                line += '\n'+15*' '
2194        print >>pFile,line
2195        if DebyePeaks['nDebye']:
2196            print >>pFile,'\n Debye diffuse scattering coefficients'
2197            parms = ['DebyeA','DebyeR','DebyeU']
2198            line = ' names :  '
2199            for parm in parms:
2200                line += '%8s refine?'%(parm)
2201            print >>pFile,line
2202            for j,term in enumerate(DebyePeaks['debyeTerms']):
2203                line = ' term'+'%2d'%(j)+':'
2204                for i in range(3):
2205                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2206                print >>pFile,line
2207        if DebyePeaks['nPeaks']:
2208            print >>pFile,'\n Single peak coefficients'
2209            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2210            line = ' names :  '
2211            for parm in parms:
2212                line += '%8s refine?'%(parm)
2213            print >>pFile,line
2214            for j,term in enumerate(DebyePeaks['peaksList']):
2215                line = ' peak'+'%2d'%(j)+':'
2216                for i in range(4):
2217                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2218                print >>pFile,line
2219       
2220    def PrintInstParms(Inst):
2221        print >>pFile,'\n Instrument Parameters:'
2222        ptlbls = ' name  :'
2223        ptstr =  ' value :'
2224        varstr = ' refine:'
2225        insKeys = Inst.keys()
2226        insKeys.sort()
2227        for item in insKeys:
2228            if item != 'Type':
2229                ptlbls += '%12s' % (item)
2230                ptstr += '%12.6f' % (Inst[item][1])
2231                if item in ['Lam1','Lam2','Azimuth']:
2232                    varstr += 12*' '
2233                else:
2234                    varstr += '%12s' % (str(bool(Inst[item][2])))
2235        print >>pFile,ptlbls
2236        print >>pFile,ptstr
2237        print >>pFile,varstr
2238       
2239    def PrintSampleParms(Sample):
2240        print >>pFile,'\n Sample Parameters:'
2241        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2242            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2243        ptlbls = ' name  :'
2244        ptstr =  ' value :'
2245        varstr = ' refine:'
2246        if 'Bragg' in Sample['Type']:
2247            for item in ['Scale','Shift','Transparency']:
2248                ptlbls += '%14s'%(item)
2249                ptstr += '%14.4f'%(Sample[item][0])
2250                varstr += '%14s'%(str(bool(Sample[item][1])))
2251           
2252        elif 'Debye' in Type:        #Debye-Scherrer
2253            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2254                ptlbls += '%14s'%(item)
2255                ptstr += '%14.4f'%(Sample[item][0])
2256                varstr += '%14s'%(str(bool(Sample[item][1])))
2257
2258        print >>pFile,ptlbls
2259        print >>pFile,ptstr
2260        print >>pFile,varstr
2261       
2262    histDict = {}
2263    histVary = []
2264    controlDict = {}
2265    histoList = Histograms.keys()
2266    histoList.sort()
2267    for histogram in histoList:
2268        if 'PWDR' in histogram:
2269            Histogram = Histograms[histogram]
2270            hId = Histogram['hId']
2271            pfx = ':'+str(hId)+':'
2272            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2273            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2274           
2275            Background = Histogram['Background']
2276            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2277            controlDict[pfx+'bakType'] = Type
2278            histDict.update(bakDict)
2279            histVary += bakVary
2280           
2281            Inst = Histogram['Instrument Parameters'][0]
2282            Type,instDict,insVary = GetInstParms(hId,Inst)
2283            controlDict[pfx+'histType'] = Type
2284            if pfx+'Lam1' in instDict:
2285                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2286            else:
2287                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2288            histDict.update(instDict)
2289            histVary += insVary
2290           
2291            Sample = Histogram['Sample Parameters']
2292            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2293            controlDict[pfx+'instType'] = Type
2294            histDict.update(sampDict)
2295            histVary += sampVary
2296           
2297   
2298            if Print: 
2299                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2300                print >>pFile,135*'-'
2301                Units = {'C':' deg','T':' msec'}
2302                units = Units[controlDict[pfx+'histType'][2]]
2303                Limits = controlDict[pfx+'Limits']
2304                print >>pFile,' Instrument type: ',Sample['Type']
2305                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
2306                PrintSampleParms(Sample)
2307                PrintInstParms(Inst)
2308                PrintBackground(Background)
2309        elif 'HKLF' in histogram:
2310            Histogram = Histograms[histogram]
2311            hId = Histogram['hId']
2312            pfx = ':'+str(hId)+':'
2313            controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
2314            Inst = Histogram['Instrument Parameters'][0]
2315            controlDict[pfx+'histType'] = Inst['Type'][0]
2316            histDict[pfx+'Lam'] = Inst['Lam'][1]
2317            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2318    return histVary,histDict,controlDict
2319   
2320def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2321    'needs a doc string'
2322   
2323    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2324        Back = Background[0]
2325        DebyePeaks = Background[1]
2326        lenBack = len(Back[3:])
2327        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2328        for i in range(lenBack):
2329            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2330            if pfx+'Back:'+str(i) in sigDict:
2331                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2332        if DebyePeaks['nDebye']:
2333            for i in range(DebyePeaks['nDebye']):
2334                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2335                for j,name in enumerate(names):
2336                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2337                    if name in sigDict:
2338                        backSig[lenBack+3*i+j] = sigDict[name]           
2339        if DebyePeaks['nPeaks']:
2340            for i in range(DebyePeaks['nPeaks']):
2341                names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
2342                    pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
2343                for j,name in enumerate(names):
2344                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2345                    if name in sigDict:
2346                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2347        return backSig
2348       
2349    def SetInstParms(pfx,Inst,parmDict,sigDict):
2350        instSig = {}
2351        insKeys = Inst.keys()
2352        insKeys.sort()
2353        for item in insKeys:
2354            insName = pfx+item
2355            Inst[item][1] = parmDict[insName]
2356            if insName in sigDict:
2357                instSig[item] = sigDict[insName]
2358            else:
2359                instSig[item] = 0
2360        return instSig
2361       
2362    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2363        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2364            sampSig = [0 for i in range(3)]
2365            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
2366                Sample[item][0] = parmDict[pfx+item]
2367                if pfx+item in sigDict:
2368                    sampSig[i] = sigDict[pfx+item]
2369        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2370            sampSig = [0 for i in range(4)]
2371            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2372                Sample[item][0] = parmDict[pfx+item]
2373                if pfx+item in sigDict:
2374                    sampSig[i] = sigDict[pfx+item]
2375        return sampSig
2376       
2377    def PrintBackgroundSig(Background,backSig):
2378        Back = Background[0]
2379        DebyePeaks = Background[1]
2380        lenBack = len(Back[3:])
2381        valstr = ' value : '
2382        sigstr = ' sig   : '
2383        refine = False
2384        for i,back in enumerate(Back[3:]):
2385            valstr += '%10.4g'%(back)
2386            if Back[1]:
2387                refine = True
2388                sigstr += '%10.4g'%(backSig[i])
2389            else:
2390                sigstr += 10*' '
2391        if refine:
2392            print >>pFile,'\n Background function: ',Back[0]
2393            print >>pFile,valstr
2394            print >>pFile,sigstr
2395        if DebyePeaks['nDebye']:
2396            ifAny = False
2397            ptfmt = "%12.3f"
2398            names =  ' names :'
2399            ptstr =  ' values:'
2400            sigstr = ' esds  :'
2401            for item in sigDict:
2402                if 'Debye' in item:
2403                    ifAny = True
2404                    names += '%12s'%(item)
2405                    ptstr += ptfmt%(parmDict[item])
2406                    sigstr += ptfmt%(sigDict[item])
2407            if ifAny:
2408                print >>pFile,'\n Debye diffuse scattering coefficients'
2409                print >>pFile,names
2410                print >>pFile,ptstr
2411                print >>pFile,sigstr
2412        if DebyePeaks['nPeaks']:
2413            ifAny = False
2414            ptfmt = "%14.3f"
2415            names =  ' names :'
2416            ptstr =  ' values:'
2417            sigstr = ' esds  :'
2418            for item in sigDict:
2419                if 'BkPk' in item:
2420                    ifAny = True
2421                    names += '%14s'%(item)
2422                    ptstr += ptfmt%(parmDict[item])
2423                    sigstr += ptfmt%(sigDict[item])
2424            if ifAny:
2425                print >>pFile,'\n Single peak coefficients'
2426                print >>pFile,names
2427                print >>pFile,ptstr
2428                print >>pFile,sigstr
2429       
2430    def PrintInstParmsSig(Inst,instSig):
2431        ptlbls = ' names :'
2432        ptstr =  ' value :'
2433        sigstr = ' sig   :'
2434        refine = False
2435        insKeys = instSig.keys()
2436        insKeys.sort()
2437        for name in insKeys:
2438            if name not in  ['Type','Lam1','Lam2','Azimuth']:
2439                ptlbls += '%12s' % (name)
2440                ptstr += '%12.6f' % (Inst[name][1])
2441                if instSig[name]:
2442                    refine = True
2443                    sigstr += '%12.6f' % (instSig[name])
2444                else:
2445                    sigstr += 12*' '
2446        if refine:
2447            print >>pFile,'\n Instrument Parameters:'
2448            print >>pFile,ptlbls
2449            print >>pFile,ptstr
2450            print >>pFile,sigstr
2451       
2452    def PrintSampleParmsSig(Sample,sampleSig):
2453        ptlbls = ' names :'
2454        ptstr =  ' values:'
2455        sigstr = ' sig   :'
2456        refine = False
2457        if 'Bragg' in Sample['Type']:
2458            for i,item in enumerate(['Scale','Shift','Transparency']):
2459                ptlbls += '%14s'%(item)
2460                ptstr += '%14.4f'%(Sample[item][0])
2461                if sampleSig[i]:
2462                    refine = True
2463                    sigstr += '%14.4f'%(sampleSig[i])
2464                else:
2465                    sigstr += 14*' '
2466           
2467        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2468            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2469                ptlbls += '%14s'%(item)
2470                ptstr += '%14.4f'%(Sample[item][0])
2471                if sampleSig[i]:
2472                    refine = True
2473                    sigstr += '%14.4f'%(sampleSig[i])
2474                else:
2475                    sigstr += 14*' '
2476
2477        if refine:
2478            print >>pFile,'\n Sample Parameters:'
2479            print >>pFile,ptlbls
2480            print >>pFile,ptstr
2481            print >>pFile,sigstr
2482       
2483    histoList = Histograms.keys()
2484    histoList.sort()
2485    for histogram in histoList:
2486        if 'PWDR' in histogram:
2487            Histogram = Histograms[histogram]
2488            hId = Histogram['hId']
2489            pfx = ':'+str(hId)+':'
2490            Background = Histogram['Background']
2491            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2492           
2493            Inst = Histogram['Instrument Parameters'][0]
2494            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2495       
2496            Sample = Histogram['Sample Parameters']
2497            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2498
2499            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2500            print >>pFile,135*'-'
2501            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2502            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2503            if Print:
2504                print >>pFile,' Instrument type: ',Sample['Type']
2505                PrintSampleParmsSig(Sample,sampSig)
2506                PrintInstParmsSig(Inst,instSig)
2507                PrintBackgroundSig(Background,backSig)
2508               
Note: See TracBrowser for help on using the repository browser.