source: trunk/GSASIIstrIO.py @ 1143

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

rework constraints to handle names and refine flag for new var (input linear constraints); redo powder 'Sample Parameters'

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