source: trunk/GSASIIstrIO.py @ 4021

Last change on this file since 4021 was 4021, checked in by toby, 2 years ago

implement filter for screen messages; start to replace print() with G2Print(); scripting changes for PDF

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