source: trunk/GSASIIstrIO.py @ 3813

Last change on this file since 3813 was 3813, checked in by vondreele, 3 years ago

fix problem with update of modulation vector after refinement
change start Fc value to 1.0 from 100.0; gives better stability in starting LeBail? refinements.

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