source: trunk/GSASIIstrIO.py @ 3814

Last change on this file since 3814 was 3814, checked in by toby, 3 years ago

refactor to move some IO-only routines; add initial image support to scriptable

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