source: trunk/GSASIIstrIO.py @ 1137

Last change on this file since 1137 was 1137, checked in by vondreele, 8 years ago

allow a 'None' for image calibrant to facilitate fitting single rings
fix an atom hard reference to use atm pointer

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