source: trunk/GSASIIstrIO.py @ 3811

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

scriptable: fix reading py2 gpx in py3 & use of multi-set instprm files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 152.0 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2019-02-04 21:35:43 +0000 (Mon, 04 Feb 2019) $
4# $Author: toby $
5# $Revision: 3811 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 3811 2019-02-04 21:35:43Z toby $
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: 3811 $")
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 var in ['mV0','mV1','mV2']:
2019                    namstr += '%12s'%(pfx+var)
2020                    ptstr += '%12.6f'%(parmDict[pfx+var])
2021                    if pfx+var in sigDict:
2022                        sigstr += '%12.6f'%(sigDict[pfx+var])
2023                    else:
2024                        sigstr += 12*' '
2025                pFile.write(namstr+'\n')
2026                pFile.write(ptstr+'\n')
2027                pFile.write(sigstr+'\n')
2028           
2029        General['Mass'] = 0.
2030        if Phase['General'].get('doPawley'):
2031            pawleyRef = Phase['Pawley ref']
2032            for i,refl in enumerate(pawleyRef):
2033                key = pfx+'PWLref:'+str(i)
2034                refl[ik] = parmDict[key]
2035                if key in sigDict:
2036                    refl[ik+1] = sigDict[key]
2037                else:
2038                    refl[ik+1] = 0
2039        else:
2040            VRBIds = RBIds['Vector']
2041            RRBIds = RBIds['Residue']
2042            RBModels = Phase['RBModels']
2043            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2044                jrb = VRBIds.index(RBObj['RBId'])
2045                rbsx = str(irb)+':'+str(jrb)
2046                pFile.write(' Vector rigid body parameters:\n')
2047                PrintRBObjPOAndSig('RBV',rbsx)
2048                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2049            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2050                jrb = RRBIds.index(RBObj['RBId'])
2051                rbsx = str(irb)+':'+str(jrb)
2052                pFile.write(' Residue rigid body parameters:\n')
2053                PrintRBObjPOAndSig('RBR',rbsx)
2054                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2055                PrintRBObjTorAndSig(rbsx)
2056            atomsSig = {}
2057            wavesSig = {}
2058            cx,ct,cs,cia = General['AtomPtrs']
2059            for i,at in enumerate(Atoms):
2060                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2061                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2062                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2063                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2064                for ind in range(cx,cx+4):
2065                    at[ind] = parmDict[names[ind]]
2066                    if ind in range(cx,cx+3):
2067                        name = names[ind].replace('A','dA')
2068                    else:
2069                        name = names[ind]
2070                    if name in sigDict:
2071                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2072                if at[cia] == 'I':
2073                    at[cia+1] = parmDict[names[cia+1]]
2074                    if names[cia+1] in sigDict:
2075                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2076                else:
2077                    for ind in range(cia+2,cia+8):
2078                        at[ind] = parmDict[names[ind]]
2079                        if names[ind] in sigDict:
2080                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2081                if General['Type'] == 'magnetic':
2082                    for ind in range(cx+4,cx+7):
2083                        at[ind] = parmDict[names[ind]]
2084                        if names[ind] in sigDict:
2085                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2086                ind = General['AtomTypes'].index(at[ct])
2087                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2088                if General.get('Modulated',False):
2089                    AtomSS = at[-1]['SS1']
2090                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2091                        Waves = AtomSS[Stype]
2092                        if len(Waves):
2093                            waveType = Waves[0]
2094                        else:
2095                            continue
2096                        for iw,wave in enumerate(Waves[1:]):
2097                            stiw = str(i)+':'+str(iw)
2098                            if Stype == 'Spos':
2099                                if waveType in ['ZigZag','Block',] and not iw:
2100                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2101                                else:
2102                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2103                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2104                            elif Stype == 'Sadp':
2105                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2106                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2107                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2108                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2109                            elif Stype == 'Sfrac':
2110                                if 'Crenel' in waveType and not iw:
2111                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2112                                else:
2113                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2114                            elif Stype == 'Smag':
2115                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2116                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2117                            for iname,name in enumerate(names):
2118                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2119                                if pfx+name in sigDict:
2120                                    wavesSig[name] = sigDict[pfx+name]
2121                   
2122            PrintAtomsAndSig(General,Atoms,atomsSig)
2123            if General['Type'] == 'magnetic':
2124                PrintMomentsAndSig(General,Atoms,atomsSig)
2125            if General.get('Modulated',False):
2126                PrintWavesAndSig(General,Atoms,wavesSig)
2127               
2128            density = G2mth.getDensity(General)[0]
2129            pFile.write('\n Density: %f.4 g/cm**3\n'%density)
2130           
2131       
2132        textureData = General['SH Texture']   
2133        if textureData['Order']:
2134            SHtextureSig = {}
2135            for name in ['omega','chi','phi']:
2136                aname = pfx+'SH '+name
2137                textureData['Sample '+name][1] = parmDict[aname]
2138                if aname in sigDict:
2139                    SHtextureSig['Sample '+name] = sigDict[aname]
2140            for name in textureData['SH Coeff'][1]:
2141                aname = pfx+name
2142                textureData['SH Coeff'][1][name] = parmDict[aname]
2143                if aname in sigDict:
2144                    SHtextureSig[name] = sigDict[aname]
2145            PrintSHtextureAndSig(textureData,SHtextureSig)
2146        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2147            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2148                textureData,RestraintDict[phase],pFile)
2149                   
2150################################################################################
2151##### Histogram & Phase data
2152################################################################################       
2153                   
2154def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2155    '''Loads the HAP histogram/phase information into dicts
2156
2157    :param dict Phases: phase information
2158    :param dict Histograms: Histogram information
2159    :param bool Print: prints information as it is read
2160    :param file pFile: file object to print to (the default, None causes printing to the console)
2161    :param bool resetRefList: Should the contents of the Reflection List be initialized
2162      on loading. The default, True, initializes the Reflection List as it is loaded.
2163
2164    :returns: (hapVary,hapDict,controlDict)
2165      * hapVary: list of refined variables
2166      * hapDict: dict with refined variables and their values
2167      * controlDict: dict with fixed parameters
2168    '''
2169   
2170    def PrintSize(hapData):
2171        if hapData[0] in ['isotropic','uniaxial']:
2172            line = '\n Size model    : %9s'%(hapData[0])
2173            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2174            if hapData[0] == 'uniaxial':
2175                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2176            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2177            pFile.write(line+'\n')
2178        else:
2179            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2180                (hapData[0],hapData[1][2],hapData[2][2]))
2181            Snames = ['S11','S22','S33','S12','S13','S23']
2182            ptlbls = ' names :'
2183            ptstr =  ' values:'
2184            varstr = ' refine:'
2185            for i,name in enumerate(Snames):
2186                ptlbls += '%12s' % (name)
2187                ptstr += '%12.3f' % (hapData[4][i])
2188                varstr += '%12s' % (str(hapData[5][i]))
2189            pFile.write(ptlbls+'\n')
2190            pFile.write(ptstr+'\n')
2191            pFile.write(varstr+'\n')
2192       
2193    def PrintMuStrain(hapData,SGData):
2194        if hapData[0] in ['isotropic','uniaxial']:
2195            line = '\n Mustrain model: %9s'%(hapData[0])
2196            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2197            if hapData[0] == 'uniaxial':
2198                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2199            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2200            pFile.write(line+'\n')
2201        else:
2202            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2203                (hapData[0],hapData[1][2],hapData[2][2]))
2204            Snames = G2spc.MustrainNames(SGData)
2205            ptlbls = ' names :'
2206            ptstr =  ' values:'
2207            varstr = ' refine:'
2208            for i,name in enumerate(Snames):
2209                ptlbls += '%12s' % (name)
2210                ptstr += '%12.1f' % (hapData[4][i])
2211                varstr += '%12s' % (str(hapData[5][i]))
2212            pFile.write(ptlbls+'\n')
2213            pFile.write(ptstr+'\n')
2214            pFile.write(varstr+'\n')
2215
2216    def PrintHStrain(hapData,SGData):
2217        pFile.write('\n Hydrostatic/elastic strain:\n')
2218        Hsnames = G2spc.HStrainNames(SGData)
2219        ptlbls = ' names :'
2220        ptstr =  ' values:'
2221        varstr = ' refine:'
2222        for i,name in enumerate(Hsnames):
2223            ptlbls += '%12s' % (name)
2224            ptstr += '%12.4g' % (hapData[0][i])
2225            varstr += '%12s' % (str(hapData[1][i]))
2226        pFile.write(ptlbls+'\n')
2227        pFile.write(ptstr+'\n')
2228        pFile.write(varstr+'\n')
2229
2230    def PrintSHPO(hapData):
2231        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2232        ptlbls = ' names :'
2233        ptstr =  ' values:'
2234        for item in hapData[5]:
2235            ptlbls += '%12s'%(item)
2236            ptstr += '%12.3f'%(hapData[5][item]) 
2237        pFile.write(ptlbls+'\n')
2238        pFile.write(ptstr+'\n')
2239   
2240    def PrintBabinet(hapData):
2241        pFile.write('\n Babinet form factor modification:\n')
2242        ptlbls = ' names :'
2243        ptstr =  ' values:'
2244        varstr = ' refine:'
2245        for name in ['BabA','BabU']:
2246            ptlbls += '%12s' % (name)
2247            ptstr += '%12.6f' % (hapData[name][0])
2248            varstr += '%12s' % (str(hapData[name][1]))
2249        pFile.write(ptlbls+'\n')
2250        pFile.write(ptstr+'\n')
2251        pFile.write(varstr+'\n')
2252       
2253    hapDict = {}
2254    hapVary = []
2255    controlDict = {}
2256   
2257    for phase in Phases:
2258        HistoPhase = Phases[phase]['Histograms']
2259        SGData = Phases[phase]['General']['SGData']
2260        cell = Phases[phase]['General']['Cell'][1:7]
2261        A = G2lat.cell2A(cell)
2262        if Phases[phase]['General'].get('Modulated',False):
2263            SSGData = Phases[phase]['General']['SSGData']
2264            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2265        pId = Phases[phase]['pId']
2266        histoList = list(HistoPhase.keys())
2267        histoList.sort()
2268        for histogram in histoList:
2269            try:
2270                Histogram = Histograms[histogram]
2271            except KeyError:                       
2272                #skip if histogram not included e.g. in a sequential refinement
2273                continue
2274            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2275                if phase in Histograms[histogram]['Reflection Lists']:
2276                    del Histograms[histogram]['Reflection Lists'][phase]
2277                continue
2278            hapData = HistoPhase[histogram]
2279            hId = Histogram['hId']
2280            if 'PWDR' in histogram:
2281                limits = Histogram['Limits'][1]
2282                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2283                if 'C' in inst['Type'][1]:
2284                    try:
2285                        wave = inst['Lam'][1]
2286                    except KeyError:
2287                        wave = inst['Lam1'][1]
2288                    dmin = wave/(2.0*sind(limits[1]/2.0))
2289                elif 'T' in inst['Type'][0]:
2290                    dmin = limits[0]/inst['difC'][1]
2291                pfx = str(pId)+':'+str(hId)+':'
2292                if Phases[phase]['General']['doPawley']:
2293                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2294                    hapDict[pfx+'newLeBail'] = True
2295                    Tmin = G2lat.Dsp2pos(inst,dmin)
2296                    if 'C' in inst['Type'][1]:
2297                        limits[1] = min(limits[1],Tmin)
2298                    else:
2299                        limits[0] = max(limits[0],Tmin)
2300                else:
2301                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2302                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2303                if Phases[phase]['General']['Type'] == 'magnetic':
2304                    dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2305                for item in ['Scale','Extinction']:
2306                    hapDict[pfx+item] = hapData[item][0]
2307                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2308                        hapVary.append(pfx+item)
2309                names = G2spc.HStrainNames(SGData)
2310                HSvals = []
2311                for i,name in enumerate(names):
2312                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2313                    HSvals.append(hapDict[pfx+name])
2314                    if hapData['HStrain'][1][i]:
2315#                    if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
2316                        hapVary.append(pfx+name)
2317                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2318                if hapData['Pref.Ori.'][0] == 'MD':
2319                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2320                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2321                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2322                        hapVary.append(pfx+'MD')
2323                else:                           #'SH' spherical harmonics
2324                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2325                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2326                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2327                    controlDict[pfx+'SHhkl'] = []
2328                    try: #patch for old Pref.Ori. items
2329                        controlDict[pfx+'SHtoler'] = 0.1
2330                        if hapData['Pref.Ori.'][6][0] != '':
2331                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2332                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2333                    except IndexError:
2334                        pass
2335                    for item in hapData['Pref.Ori.'][5]:
2336                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2337                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2338                            hapVary.append(pfx+item)
2339                for item in ['Mustrain','Size']:
2340                    controlDict[pfx+item+'Type'] = hapData[item][0]
2341                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2342                    if hapData[item][2][2]:
2343                        hapVary.append(pfx+item+';mx')
2344                    if hapData[item][0] in ['isotropic','uniaxial']:
2345                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2346                        if hapData[item][2][0]:
2347                            hapVary.append(pfx+item+';i')
2348                        if hapData[item][0] == 'uniaxial':
2349                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2350                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2351                            if hapData[item][2][1]:
2352                                hapVary.append(pfx+item+';a')
2353                    else:       #generalized for mustrain or ellipsoidal for size
2354                        Nterms = len(hapData[item][4])
2355                        if item == 'Mustrain':
2356                            names = G2spc.MustrainNames(SGData)
2357                            pwrs = []
2358                            for name in names:
2359                                h,k,l = name[1:]
2360                                pwrs.append([int(h),int(k),int(l)])
2361                            controlDict[pfx+'MuPwrs'] = pwrs
2362                        for i in range(Nterms):
2363                            sfx = ';'+str(i)
2364                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2365                            if hapData[item][5][i]:
2366                                hapVary.append(pfx+item+sfx)
2367                if Phases[phase]['General']['Type'] != 'magnetic':
2368                    for bab in ['BabA','BabU']:
2369                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2370                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2371                            hapVary.append(pfx+bab)
2372                               
2373                if Print: 
2374                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2375                    pFile.write(135*'='+'\n')
2376                    if hapDict[pfx+'LeBail']:
2377                        pFile.write(' Perform LeBail extraction\n')                     
2378                    else:
2379                        pFile.write(' Phase fraction  : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2380                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2381                        if hapData['Pref.Ori.'][0] == 'MD':
2382                            Ax = hapData['Pref.Ori.'][3]
2383                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2384                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2385                        else: #'SH' for spherical harmonics
2386                            PrintSHPO(hapData['Pref.Ori.'])
2387                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2388                    PrintSize(hapData['Size'])
2389                    PrintMuStrain(hapData['Mustrain'],SGData)
2390                    PrintHStrain(hapData['HStrain'],SGData)
2391                    if Phases[phase]['General']['Type'] != 'magnetic':
2392                        if hapData['Babinet']['BabA'][0]:
2393                            PrintBabinet(hapData['Babinet'])
2394                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2395                    hapData['newLeBail'] = True
2396                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2397                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2398                        hapData['newLeBail'] = False
2399                    refList = []
2400#                    Uniq = []
2401#                    Phi = []
2402                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2403                    if Phases[phase]['General'].get('Modulated',False):
2404                        ifSuper = True
2405                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2406                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2407                        for h,k,l,m,d in HKLd:
2408                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2409                            mul *= 2      # for powder overlap of Friedel pairs
2410                            if m or not ext or useExt:
2411                                if 'C' in inst['Type'][0]:
2412                                    pos = G2lat.Dsp2pos(inst,d)
2413                                    if limits[0] < pos < limits[1]:
2414                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,100., 0.0,0.0,1.0,1.0,1.0])
2415                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2416#                                        Uniq.append(uniq)
2417#                                        Phi.append(phi)
2418                                elif 'T' in inst['Type'][0]:
2419                                    pos = G2lat.Dsp2pos(inst,d)
2420                                    if limits[0] < pos < limits[1]:
2421                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2422                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,100., 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
2423                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2424                                        #TODO - if tabulated put alp & bet in here
2425#                                        Uniq.append(uniq)
2426#                                        Phi.append(phi)
2427                    else:
2428                        ifSuper = False
2429                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2430                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2431                        for h,k,l,d in HKLd:
2432                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2433                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2434                                ext = G2spc.checkMagextc([h,k,l],SGData)
2435                            mul *= 2      # for powder overlap of Friedel pairs
2436                            if ext and not useExt:
2437                                continue
2438                            if 'C' in inst['Type'][0]:
2439                                pos = G2lat.Dsp2pos(inst,d)
2440                                if limits[0] < pos < limits[1]:
2441                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,100., 0.0,0.0,1.0,1.0,1.0])
2442                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2443#                                    Uniq.append(uniq)
2444#                                    Phi.append(phi)
2445                            elif 'T' in inst['Type'][0]:
2446                                pos = G2lat.Dsp2pos(inst,d)
2447                                if limits[0] < pos < limits[1]:
2448                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2449                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,100., 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
2450                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2451#                                    Uniq.append(uniq)
2452#                                    Phi.append(phi)
2453                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2454            elif 'HKLF' in histogram:
2455                inst = Histogram['Instrument Parameters'][0]
2456                hId = Histogram['hId']
2457                hfx = ':%d:'%(hId)
2458                for item in inst:
2459                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2460                    hapDict[hfx+item] = inst[item][1]
2461                pfx = str(pId)+':'+str(hId)+':'
2462                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2463                if hapData['Scale'][1]:
2464                    hapVary.append(pfx+'Scale')
2465                               
2466                extApprox,extType,extParms = hapData['Extinction']
2467                controlDict[pfx+'EType'] = extType
2468                controlDict[pfx+'EApprox'] = extApprox
2469                if 'C' in inst['Type'][0]:
2470                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2471                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2472                if 'Primary' in extType:
2473                    Ekey = ['Ep',]
2474                elif 'I & II' in extType:
2475                    Ekey = ['Eg','Es']
2476                elif 'Secondary Type II' == extType:
2477                    Ekey = ['Es',]
2478                elif 'Secondary Type I' == extType:
2479                    Ekey = ['Eg',]
2480                else:   #'None'
2481                    Ekey = []
2482                for eKey in Ekey:
2483                    hapDict[pfx+eKey] = extParms[eKey][0]
2484                    if extParms[eKey][1]:
2485                        hapVary.append(pfx+eKey)
2486                for bab in ['BabA','BabU']:
2487                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2488                    if hapData['Babinet'][bab][1]:
2489                        hapVary.append(pfx+bab)
2490                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2491                if len(Twins) == 1:
2492                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2493                    if hapData.get('Flack',[0,False])[1]:
2494                        hapVary.append(pfx+'Flack')
2495                sumTwFr = 0.
2496                controlDict[pfx+'TwinLaw'] = []
2497                controlDict[pfx+'TwinInv'] = []
2498                NTL = 0           
2499                for it,twin in enumerate(Twins):
2500                    if 'bool' in str(type(twin[0])):
2501                        controlDict[pfx+'TwinInv'].append(twin[0])
2502                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2503                    else:
2504                        NTL += 1
2505                        controlDict[pfx+'TwinInv'].append(False)
2506                        controlDict[pfx+'TwinLaw'].append(twin[0])
2507                    if it:
2508                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2509                        sumTwFr += twin[1]
2510                    else:
2511                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2512                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2513                    if Twins[0][1][1]:
2514                        hapVary.append(pfx+'TwinFr:'+str(it))
2515                controlDict[pfx+'NTL'] = NTL
2516                #need to make constraint on TwinFr
2517                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2518                if len(Twins) > 1:    #force sum to unity
2519                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2520                if Print: 
2521                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2522                    pFile.write(135*'='+'\n')
2523                    pFile.write(' Scale factor     : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2524                    if extType != 'None':
2525                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2526                        text = ' Parameters       :'
2527                        for eKey in Ekey:
2528                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2529                        pFile.write(text+'\n')
2530                    if hapData['Babinet']['BabA'][0]:
2531                        PrintBabinet(hapData['Babinet'])
2532                    if not SGData['SGInv'] and len(Twins) == 1:
2533                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2534                    if len(Twins) > 1:
2535                        for it,twin in enumerate(Twins):
2536                            if 'bool' in str(type(twin[0])):
2537                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2538                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2539                            else:
2540                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2541                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2542                       
2543                Histogram['Reflection Lists'] = phase       
2544               
2545    return hapVary,hapDict,controlDict
2546   
2547def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2548    'needs a doc string'
2549   
2550    def PrintSizeAndSig(hapData,sizeSig):
2551        line = '\n Size model:     %9s'%(hapData[0])
2552        refine = False
2553        if hapData[0] in ['isotropic','uniaxial']:
2554            line += ' equatorial:%12.5f'%(hapData[1][0])
2555            if sizeSig[0][0]:
2556                line += ', sig:%8.4f'%(sizeSig[0][0])
2557                refine = True
2558            if hapData[0] == 'uniaxial':
2559                line += ' axial:%12.4f'%(hapData[1][1])
2560                if sizeSig[0][1]:
2561                    refine = True
2562                    line += ', sig:%8.4f'%(sizeSig[0][1])
2563            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2564            if sizeSig[0][2]:
2565                refine = True
2566                line += ', sig:%8.4f'%(sizeSig[0][2])
2567            if refine:
2568                pFile.write(line+'\n')
2569        else:
2570            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2571            if sizeSig[0][2]:
2572                refine = True
2573                line += ', sig:%8.4f'%(sizeSig[0][2])
2574            Snames = ['S11','S22','S33','S12','S13','S23']
2575            ptlbls = ' name  :'
2576            ptstr =  ' value :'
2577            sigstr = ' sig   :'
2578            for i,name in enumerate(Snames):
2579                ptlbls += '%12s' % (name)
2580                ptstr += '%12.6f' % (hapData[4][i])
2581                if sizeSig[1][i]:
2582                    refine = True
2583                    sigstr += '%12.6f' % (sizeSig[1][i])
2584                else:
2585                    sigstr += 12*' '
2586            if refine:
2587                pFile.write(line+'\n')
2588                pFile.write(ptlbls+'\n')
2589                pFile.write(ptstr+'\n')
2590                pFile.write(sigstr+'\n')
2591       
2592    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2593        line = '\n Mustrain model: %9s\n'%(hapData[0])
2594        refine = False
2595        if hapData[0] in ['isotropic','uniaxial']:
2596            line += ' equatorial:%12.1f'%(hapData[1][0])
2597            if mustrainSig[0][0]:
2598                line += ', sig:%8.1f'%(mustrainSig[0][0])
2599                refine = True
2600            if hapData[0] == 'uniaxial':
2601                line += ' axial:%12.1f'%(hapData[1][1])
2602                if mustrainSig[0][1]:
2603                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2604            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2605            if mustrainSig[0][2]:
2606                refine = True
2607                line += ', sig:%8.3f'%(mustrainSig[0][2])
2608            if refine:
2609                pFile.write(line+'\n')
2610        else:
2611            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2612            if mustrainSig[0][2]:
2613                refine = True
2614                line += ', sig:%8.3f'%(mustrainSig[0][2])
2615            Snames = G2spc.MustrainNames(SGData)
2616            ptlbls = ' name  :'
2617            ptstr =  ' value :'
2618            sigstr = ' sig   :'
2619            for i,name in enumerate(Snames):
2620                ptlbls += '%12s' % (name)
2621                ptstr += '%12.1f' % (hapData[4][i])
2622                if mustrainSig[1][i]:
2623                    refine = True
2624                    sigstr += '%12.1f' % (mustrainSig[1][i])
2625                else:
2626                    sigstr += 12*' '
2627            if refine:
2628                pFile.write(line+'\n')
2629                pFile.write(ptlbls+'\n')
2630                pFile.write(ptstr+'\n')
2631                pFile.write(sigstr+'\n')
2632           
2633    def PrintHStrainAndSig(hapData,strainSig,SGData):
2634        Hsnames = G2spc.HStrainNames(SGData)
2635        ptlbls = ' name  :'
2636        ptstr =  ' value :'
2637        sigstr = ' sig   :'
2638        refine = False
2639        for i,name in enumerate(Hsnames):
2640            ptlbls += '%12s' % (name)
2641            ptstr += '%12.4g' % (hapData[0][i])
2642            if name in strainSig:
2643                refine = True
2644                sigstr += '%12.4g' % (strainSig[name])
2645            else:
2646                sigstr += 12*' '
2647        if refine:
2648            pFile.write('\n Hydrostatic/elastic strain:\n')
2649            pFile.write(ptlbls+'\n')
2650            pFile.write(ptstr+'\n')
2651            pFile.write(sigstr+'\n')
2652       
2653    def PrintSHPOAndSig(pfx,hapData,POsig):
2654        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2655        ptlbls = ' names :'
2656        ptstr =  ' values:'
2657        sigstr = ' sig   :'
2658        for item in hapData[5]:
2659            ptlbls += '%12s'%(item)
2660            ptstr += '%12.3f'%(hapData[5][item])
2661            if pfx+item in POsig:
2662                sigstr += '%12.3f'%(POsig[pfx+item])
2663            else:
2664                sigstr += 12*' ' 
2665        pFile.write(ptlbls+'\n')
2666        pFile.write(ptstr+'\n')
2667        pFile.write(sigstr+'\n')
2668       
2669    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2670        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2671        text = ''
2672        for item in hapData[2]:
2673            if pfx+item in ScalExtSig:
2674                text += '       %s: '%(item)
2675                text += '%12.2e'%(hapData[2][item][0])
2676                if pfx+item in ScalExtSig:
2677                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2678        pFile.write(text+'\n')   
2679       
2680    def PrintBabinetAndSig(pfx,hapData,BabSig):
2681        pFile.write('\n Babinet form factor modification:\n')
2682        ptlbls = ' names :'
2683        ptstr =  ' values:'
2684        sigstr = ' sig   :'
2685        for item in hapData:
2686            ptlbls += '%12s'%(item)
2687            ptstr += '%12.3f'%(hapData[item][0])
2688            if pfx+item in BabSig:
2689                sigstr += '%12.3f'%(BabSig[pfx+item])
2690            else:
2691                sigstr += 12*' ' 
2692        pFile.write(ptlbls+'\n')
2693        pFile.write(ptstr+'\n')
2694        pFile.write(sigstr+'\n')
2695       
2696    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2697        pFile.write('\n Twin Law fractions :\n')
2698        ptlbls = ' names :'
2699        ptstr =  ' values:'
2700        sigstr = ' sig   :'
2701        for it,item in enumerate(twinData):
2702            ptlbls += '     twin #%d'%(it)
2703            if it:
2704                ptstr += '%12.3f'%(item[1])
2705            else:
2706                ptstr += '%12.3f'%(item[1][0])
2707            if pfx+'TwinFr:'+str(it) in TwinSig:
2708                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2709            else:
2710                sigstr += 12*' ' 
2711        pFile.write(ptlbls+'\n')
2712        pFile.write(ptstr+'\n')
2713        pFile.write(sigstr+'\n')
2714       
2715   
2716    PhFrExtPOSig = {}
2717    SizeMuStrSig = {}
2718    ScalExtSig = {}
2719    BabSig = {}
2720    TwinFrSig = {}
2721    wtFrSum = {}
2722    for phase in Phases:
2723        HistoPhase = Phases[phase]['Histograms']
2724        General = Phases[phase]['General']
2725        SGData = General['SGData']
2726        pId = Phases[phase]['pId']
2727        histoList = list(HistoPhase.keys())
2728        histoList.sort()
2729        for histogram in histoList:
2730            try:
2731                Histogram = Histograms[histogram]
2732            except KeyError:                       
2733                #skip if histogram not included e.g. in a sequential refinement
2734                continue
2735            if not Phases[phase]['Histograms'][histogram]['Use']:
2736                #skip if phase absent from this histogram
2737                continue
2738            hapData = HistoPhase[histogram]
2739            hId = Histogram['hId']
2740            pfx = str(pId)+':'+str(hId)+':'
2741            if hId not in wtFrSum:
2742                wtFrSum[hId] = 0.
2743            if 'PWDR' in histogram:
2744                for item in ['Scale','Extinction']:
2745                    hapData[item][0] = parmDict[pfx+item]
2746                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2747                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2748                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2749                if hapData['Pref.Ori.'][0] == 'MD':
2750                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2751                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2752                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2753                else:                           #'SH' spherical harmonics
2754                    for item in hapData['Pref.Ori.'][5]:
2755                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2756                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2757                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2758                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2759                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2760                    pfx+'HStrain':{}})                 
2761                for item in ['Mustrain','Size']:
2762                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2763#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2764                    if pfx+item+';mx' in sigDict:
2765                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2766                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2767                        hapData[item][1][0] = parmDict[pfx+item+';i']
2768                        if item == 'Size':
2769                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2770                        if pfx+item+';i' in sigDict: 
2771                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2772                        if hapData[item][0] == 'uniaxial':
2773                            hapData[item][1][1] = parmDict[pfx+item+';a']
2774                            if item == 'Size':
2775                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2776                            if pfx+item+';a' in sigDict:
2777                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2778                    else:       #generalized for mustrain or ellipsoidal for size
2779                        Nterms = len(hapData[item][4])
2780                        for i in range(Nterms):
2781                            sfx = ';'+str(i)
2782                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2783                            if pfx+item+sfx in sigDict:
2784                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2785                names = G2spc.HStrainNames(SGData)
2786                for i,name in enumerate(names):
2787                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2788                    if pfx+name in sigDict:
2789                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2790                if Phases[phase]['General']['Type'] != 'magnetic':
2791                    for name in ['BabA','BabU']:
2792                        hapData['Babinet'][name][0] = parmDict[pfx+name]
2793                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2794                            BabSig[pfx+name] = sigDict[pfx+name]               
2795               
2796            elif 'HKLF' in histogram:
2797                for item in ['Scale','Flack']:
2798                    if parmDict.get(pfx+item):
2799                        hapData[item][0] = parmDict[pfx+item]
2800                        if pfx+item in sigDict:
2801                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2802                for item in ['Ep','Eg','Es']:
2803                    if parmDict.get(pfx+item):
2804                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2805                        if pfx+item in sigDict:
2806                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2807                for item in ['BabA','BabU']:
2808                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2809                    if pfx+item in sigDict:
2810                        BabSig[pfx+item] = sigDict[pfx+item]
2811                if 'Twins' in hapData:
2812                    it = 1
2813                    sumTwFr = 0.
2814                    while True:
2815                        try:
2816                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2817                            if pfx+'TwinFr:'+str(it) in sigDict:
2818                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2819                            if it:
2820                                sumTwFr += hapData['Twins'][it][1]
2821                            it += 1
2822                        except KeyError:
2823                            break
2824                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2825
2826    if Print:
2827        for phase in Phases:
2828            HistoPhase = Phases[phase]['Histograms']
2829            General = Phases[phase]['General']
2830            SGData = General['SGData']
2831            pId = Phases[phase]['pId']
2832            histoList = list(HistoPhase.keys())
2833            histoList.sort()
2834            for histogram in histoList:
2835                try:
2836                    Histogram = Histograms[histogram]
2837                except KeyError:                       
2838                    #skip if histogram not included e.g. in a sequential refinement
2839                    continue
2840                hapData = HistoPhase[histogram]
2841                hId = Histogram['hId']
2842                Histogram['Residuals'][str(pId)+'::Name'] = phase
2843                pfx = str(pId)+':'+str(hId)+':'
2844                hfx = ':%s:'%(hId)
2845                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
2846                    continue
2847                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2848                pFile.write(135*'='+'\n')
2849                if 'PWDR' in histogram:
2850                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
2851                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
2852                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
2853                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
2854                   
2855                    if parmDict[pfx+'LeBail']:
2856                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
2857                    else:
2858                        if pfx+'Scale' in PhFrExtPOSig:
2859                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2860                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2861                            pFile.write(' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f\n'%
2862                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
2863                        if pfx+'Extinction' in PhFrExtPOSig:
2864                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
2865                        if hapData['Pref.Ori.'][0] == 'MD':
2866                            if pfx+'MD' in PhFrExtPOSig:
2867                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
2868                        else:
2869                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2870                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2871                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2872                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2873                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
2874                        if len(BabSig):
2875                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2876                   
2877                elif 'HKLF' in histogram:
2878                    Inst = Histogram['Instrument Parameters'][0]
2879                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
2880                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2881                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
2882                    if FFtables != None and 'N' not in Inst['Type'][0]:
2883                        PrintFprime(FFtables,hfx,pFile)
2884                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
2885                    if pfx+'Scale' in ScalExtSig:
2886                        pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
2887                    if hapData['Extinction'][0] != 'None':
2888                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2889                    if len(BabSig):
2890                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2891                    if pfx+'Flack' in ScalExtSig:
2892                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
2893                    if len(TwinFrSig):
2894                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2895
2896################################################################################
2897##### Histogram data
2898################################################################################       
2899                   
2900def GetHistogramData(Histograms,Print=True,pFile=None):
2901    'needs a doc string'
2902   
2903    def GetBackgroundParms(hId,Background):
2904        Back = Background[0]
2905        DebyePeaks = Background[1]
2906        bakType,bakFlag = Back[:2]
2907        backVals = Back[3:]
2908        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2909        backDict = dict(zip(backNames,backVals))
2910        backVary = []
2911        if bakFlag:
2912            backVary = backNames
2913        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2914        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2915        debyeDict = {}
2916        debyeList = []
2917        for i in range(DebyePeaks['nDebye']):
2918            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2919            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2920            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2921        debyeVary = []
2922        for item in debyeList:
2923            if item[1]:
2924                debyeVary.append(item[0])
2925        backDict.update(debyeDict)
2926        backVary += debyeVary
2927        peakDict = {}
2928        peakList = []
2929        for i in range(DebyePeaks['nPeaks']):
2930            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2931                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2932            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2933            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2934        peakVary = []
2935        for item in peakList:
2936            if item[1]:
2937                peakVary.append(item[0])
2938        backDict.update(peakDict)
2939        backVary += peakVary
2940        return bakType,backDict,backVary           
2941       
2942    def GetInstParms(hId,Inst):
2943        #patch
2944        if 'Z' not in Inst:
2945            Inst['Z'] = [0.0,0.0,False]
2946        dataType = Inst['Type'][0]
2947        instDict = {}
2948        insVary = []
2949        pfx = ':'+str(hId)+':'
2950        insKeys = list(Inst.keys())
2951        insKeys.sort()
2952        for item in insKeys:
2953            insName = pfx+item
2954            instDict[insName] = Inst[item][1]
2955            if len(Inst[item]) > 2 and Inst[item][2]:
2956                insVary.append(insName)
2957        if 'C' in dataType:
2958            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2959        elif 'T' in dataType:   #trap zero alp, bet coeff.
2960            if not instDict[pfx+'alpha']:
2961                instDict[pfx+'alpha'] = 1.0
2962            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
2963                instDict[pfx+'beta-1'] = 1.0
2964        return dataType,instDict,insVary
2965       
2966    def GetSampleParms(hId,Sample):
2967        sampVary = []
2968        hfx = ':'+str(hId)+':'       
2969        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2970            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2971        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2972            if key in Sample:
2973                sampDict[hfx+key] = Sample[key]
2974        Type = Sample['Type']
2975        if 'Bragg' in Type:             #Bragg-Brentano
2976            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2977                sampDict[hfx+item] = Sample[item][0]
2978                if Sample[item][1]:
2979                    sampVary.append(hfx+item)
2980        elif 'Debye' in Type:        #Debye-Scherrer
2981            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2982                sampDict[hfx+item] = Sample[item][0]
2983                if Sample[item][1]:
2984                    sampVary.append(hfx+item)
2985        return Type,sampDict,sampVary
2986       
2987    def PrintBackground(Background):
2988        Back = Background[0]
2989        DebyePeaks = Background[1]
2990        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
2991        line = ' Coefficients: '
2992        for i,back in enumerate(Back[3:]):
2993            line += '%10.3f'%(back)
2994            if i and not i%10:
2995                line += '\n'+15*' '
2996        pFile.write(line+'\n')
2997        if DebyePeaks['nDebye']:
2998            pFile.write('\n Debye diffuse scattering coefficients\n')
2999            parms = ['DebyeA','DebyeR','DebyeU']
3000            line = ' names :  '
3001            for parm in parms:
3002                line += '%8s refine?'%(parm)
3003            pFile.write(line+'\n')
3004            for j,term in enumerate(DebyePeaks['debyeTerms']):
3005                line = ' term'+'%2d'%(j)+':'
3006                for i in range(3):
3007                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3008                pFile.write(line+'\n')
3009        if DebyePeaks['nPeaks']:
3010            pFile.write('\n Single peak coefficients\n')
3011            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3012            line = ' names :  '
3013            for parm in parms:
3014                line += '%8s refine?'%(parm)
3015            pFile.write(line+'\n')
3016            for j,term in enumerate(DebyePeaks['peaksList']):
3017                line = ' peak'+'%2d'%(j)+':'
3018                for i in range(4):
3019                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3020                pFile.write(line+'\n')
3021       
3022    def PrintInstParms(Inst):
3023        pFile.write('\n Instrument Parameters:\n')
3024        insKeys = list(Inst.keys())
3025        insKeys.sort()
3026        iBeg = 0
3027        Ok = True
3028        while Ok:
3029            ptlbls = ' name  :'
3030            ptstr =  ' value :'
3031            varstr = ' refine:'
3032            iFin = min(iBeg+9,len(insKeys))
3033            for item in insKeys[iBeg:iFin]:
3034                if item not in ['Type','Source']:
3035                    ptlbls += '%12s' % (item)
3036                    ptstr += '%12.6f' % (Inst[item][1])
3037                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3038                        varstr += 12*' '
3039                    else:
3040                        varstr += '%12s' % (str(bool(Inst[item][2])))
3041            pFile.write(ptlbls+'\n')
3042            pFile.write(ptstr+'\n')
3043            pFile.write(varstr+'\n')
3044            iBeg = iFin
3045            if iBeg == len(insKeys):
3046                Ok = False
3047            else:
3048                pFile.write('\n')
3049       
3050    def PrintSampleParms(Sample):
3051        pFile.write('\n Sample Parameters:\n')
3052        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3053            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3054        ptlbls = ' name  :'
3055        ptstr =  ' value :'
3056        varstr = ' refine:'
3057        if 'Bragg' in Sample['Type']:
3058            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3059                ptlbls += '%14s'%(item)
3060                ptstr += '%14.4f'%(Sample[item][0])
3061                varstr += '%14s'%(str(bool(Sample[item][1])))
3062           
3063        elif 'Debye' in Type:        #Debye-Scherrer
3064            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3065                ptlbls += '%14s'%(item)
3066                ptstr += '%14.4f'%(Sample[item][0])
3067                varstr += '%14s'%(str(bool(Sample[item][1])))
3068
3069        pFile.write(ptlbls+'\n')
3070        pFile.write(ptstr+'\n')
3071        pFile.write(varstr+'\n')
3072       
3073    histDict = {}
3074    histVary = []
3075    controlDict = {}
3076    histoList = list(Histograms.keys())
3077    histoList.sort()
3078    for histogram in histoList:
3079        if 'PWDR' in histogram:
3080            Histogram = Histograms[histogram]
3081            hId = Histogram['hId']
3082            pfx = ':'+str(hId)+':'
3083            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3084            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3085            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3086            for excl in controlDict[pfx+'Exclude']:
3087                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3088            if controlDict[pfx+'Exclude']:
3089                ma.mask_rows(Histogram['Data'])
3090            Background = Histogram['Background']
3091            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3092            controlDict[pfx+'bakType'] = Type
3093            histDict.update(bakDict)
3094            histVary += bakVary
3095           
3096            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3097            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3098                print ('Warning: tabulated profile coefficients are ignored')
3099            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3100            controlDict[pfx+'histType'] = Type
3101            if 'XC' in Type:
3102                if pfx+'Lam1' in instDict:
3103                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3104                else:
3105                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3106            histDict.update(instDict)
3107            histVary += insVary
3108           
3109            Sample = Histogram['Sample Parameters']
3110            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3111            controlDict[pfx+'instType'] = Type
3112            histDict.update(sampDict)
3113            histVary += sampVary
3114           
3115   
3116            if Print: 
3117                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3118                pFile.write(135*'='+'\n')
3119                Units = {'C':' deg','T':' msec'}
3120                units = Units[controlDict[pfx+'histType'][2]]
3121                Limits = controlDict[pfx+'Limits']
3122                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3123                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3124                if len(controlDict[pfx+'Exclude']):
3125                    excls = controlDict[pfx+'Exclude']
3126                    for excl in excls:
3127                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3128                PrintSampleParms(Sample)
3129                PrintInstParms(Inst[0])
3130                PrintBackground(Background)
3131        elif 'HKLF' in histogram:
3132            Histogram = Histograms[histogram]
3133            hId = Histogram['hId']
3134            pfx = ':'+str(hId)+':'
3135            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3136            Inst = Histogram['Instrument Parameters'][0]
3137            controlDict[pfx+'histType'] = Inst['Type'][0]
3138            if 'X' in Inst['Type'][0]:
3139                histDict[pfx+'Lam'] = Inst['Lam'][1]
3140                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3141            elif 'NC' in Inst['Type'][0]:                   
3142                histDict[pfx+'Lam'] = Inst['Lam'][1]
3143    return histVary,histDict,controlDict
3144   
3145def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3146    'needs a doc string'
3147   
3148    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3149        Back = Background[0]
3150        DebyePeaks = Background[1]
3151        lenBack = len(Back[3:])
3152        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3153        for i in range(lenBack):
3154            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3155            if pfx+'Back;'+str(i) in sigDict:
3156                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3157        if DebyePeaks['nDebye']:
3158            for i in range(DebyePeaks['nDebye']):
3159                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3160                for j,name in enumerate(names):
3161                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3162                    if name in sigDict:
3163                        backSig[lenBack+3*i+j] = sigDict[name]           
3164        if DebyePeaks['nPeaks']:
3165            for i in range(DebyePeaks['nPeaks']):
3166                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3167                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3168                for j,name in enumerate(names):
3169                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3170                    if name in sigDict:
3171                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3172        return backSig
3173       
3174    def SetInstParms(pfx,Inst,parmDict,sigDict):
3175        instSig = {}
3176        insKeys = list(Inst.keys())
3177        insKeys.sort()
3178        for item in insKeys:
3179            insName = pfx+item
3180            Inst[item][1] = parmDict[insName]
3181            if insName in sigDict:
3182                instSig[item] = sigDict[insName]
3183            else:
3184                instSig[item] = 0
3185        return instSig
3186       
3187    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3188        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3189            sampSig = [0 for i in range(5)]
3190            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3191                Sample[item][0] = parmDict[pfx+item]
3192                if pfx+item in sigDict:
3193                    sampSig[i] = sigDict[pfx+item]
3194        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3195            sampSig = [0 for i in range(4)]
3196            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3197                Sample[item][0] = parmDict[pfx+item]
3198                if pfx+item in sigDict:
3199                    sampSig[i] = sigDict[pfx+item]
3200        return sampSig
3201       
3202    def PrintBackgroundSig(Background,backSig):
3203        Back = Background[0]
3204        DebyePeaks = Background[1]
3205        valstr = ' value : '
3206        sigstr = ' sig   : '
3207        refine = False
3208        for i,back in enumerate(Back[3:]):
3209            valstr += '%10.4g'%(back)
3210            if Back[1]:
3211                refine = True
3212                sigstr += '%10.4g'%(backSig[i])
3213            else:
3214                sigstr += 10*' '
3215        if refine:
3216            pFile.write('\n Background function: %s\n'%Back[0])
3217            pFile.write(valstr+'\n')
3218            pFile.write(sigstr+'\n')
3219        if DebyePeaks['nDebye']:
3220            ifAny = False
3221            ptfmt = "%12.3f"
3222            names =  ' names :'
3223            ptstr =  ' values:'
3224            sigstr = ' esds  :'
3225            for item in sigDict:
3226                if 'Debye' in item:
3227                    ifAny = True
3228                    names += '%12s'%(item)
3229                    ptstr += ptfmt%(parmDict[item])
3230                    sigstr += ptfmt%(sigDict[item])
3231            if ifAny:
3232                pFile.write('\n Debye diffuse scattering coefficients\n')
3233                pFile.write(names+'\n')
3234                pFile.write(ptstr+'\n')
3235                pFile.write(sigstr+'\n')
3236        if DebyePeaks['nPeaks']:
3237            pFile.write('\n Single peak coefficients:\n')
3238            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3239            line = ' peak no. '
3240            for parm in parms:
3241                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3242            pFile.write(line+'\n')
3243            for ip in range(DebyePeaks['nPeaks']):
3244                ptstr = ' %4d '%(ip)
3245                for parm in parms:
3246                    fmt = '%14.3f'
3247                    efmt = '%12.3f'
3248                    if parm == 'BkPkpos':
3249                        fmt = '%14.4f'
3250                        efmt = '%12.4f'
3251                    name = pfx+parm+';%d'%(ip)
3252                    ptstr += fmt%(parmDict[name])
3253                    if name in sigDict:
3254                        ptstr += efmt%(sigDict[name])
3255                    else:
3256                        ptstr += 12*' '
3257                pFile.write(ptstr+'\n')
3258        sumBk = np.array(Histogram['sumBk'])
3259        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3260            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3261       
3262    def PrintInstParmsSig(Inst,instSig):
3263        refine = False
3264        insKeys = list(instSig.keys())
3265        insKeys.sort()
3266        iBeg = 0
3267        Ok = True
3268        while Ok:
3269            ptlbls = ' names :'
3270            ptstr =  ' value :'
3271            sigstr = ' sig   :'
3272            iFin = min(iBeg+9,len(insKeys))
3273            for name in insKeys[iBeg:iFin]:
3274                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3275                    ptlbls += '%12s' % (name)
3276                    ptstr += '%12.6f' % (Inst[name][1])
3277                    if instSig[name]:
3278                        refine = True
3279                        sigstr += '%12.6f' % (instSig[name])
3280                    else:
3281                        sigstr += 12*' '
3282            if refine:
3283                pFile.write('\n Instrument Parameters:\n')
3284                pFile.write(ptlbls+'\n')
3285                pFile.write(ptstr+'\n')
3286                pFile.write(sigstr+'\n')
3287            iBeg = iFin
3288            if iBeg == len(insKeys):
3289                Ok = False
3290       
3291    def PrintSampleParmsSig(Sample,sampleSig):
3292        ptlbls = ' names :'
3293        ptstr =  ' values:'
3294        sigstr = ' sig   :'
3295        refine = False
3296        if 'Bragg' in Sample['Type']:
3297            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3298                ptlbls += '%14s'%(item)
3299                ptstr += '%14.4f'%(Sample[item][0])
3300                if sampleSig[i]:
3301                    refine = True
3302                    sigstr += '%14.4f'%(sampleSig[i])
3303                else:
3304                    sigstr += 14*' '
3305           
3306        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3307            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3308                ptlbls += '%14s'%(item)
3309                ptstr += '%14.4f'%(Sample[item][0])
3310                if sampleSig[i]:
3311                    refine = True
3312                    sigstr += '%14.4f'%(sampleSig[i])
3313                else:
3314                    sigstr += 14*' '
3315
3316        if refine:
3317            pFile.write('\n Sample Parameters:\n')
3318            pFile.write(ptlbls+'\n')
3319            pFile.write(ptstr+'\n')
3320            pFile.write(sigstr+'\n')
3321       
3322    histoList = list(Histograms.keys())
3323    histoList.sort()
3324    for histogram in histoList:
3325        if 'PWDR' in histogram:
3326            Histogram = Histograms[histogram]
3327            hId = Histogram['hId']
3328            pfx = ':'+str(hId)+':'
3329            Background = Histogram['Background']
3330            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3331           
3332            Inst = Histogram['Instrument Parameters'][0]
3333            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3334       
3335            Sample = Histogram['Sample Parameters']
3336            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3337
3338            pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3339            pFile.write(135*'='+'\n')
3340            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3341            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3342                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3343            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3344                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3345            if Print:
3346                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3347                if FFtables != None and 'N' not in Inst['Type'][0]:
3348                    PrintFprime(FFtables,pfx,pFile)
3349                PrintSampleParmsSig(Sample,sampSig)
3350                PrintInstParmsSig(Inst,instSig)
3351                PrintBackgroundSig(Background,backSig)
3352               
Note: See TracBrowser for help on using the repository browser.