source: trunk/GSASIIstrIO.py @ 1842

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

typo! not setRef should be seqRef

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