source: trunk/GSASIIstrIO.py @ 1616

Last change on this file since 1616 was 1616, checked in by vondreele, 7 years ago

oops forgot to remove a raise Exception

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