source: trunk/GSASIIstrIO.py @ 1852

Last change on this file since 1852 was 1852, checked in by vondreele, 8 years ago

addmtexture analysis tutorial to list in G2ctrls
remove Clear Texture - not needed (commented out for now)
Add print of resonant form factors for x-ray & cw neutron refinements. Appears in the results part of the ,lst file

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