source: trunk/GSASIIstrIO.py @ 2405

Last change on this file since 2405 was 2405, checked in by vondreele, 5 years ago

fix missing 'Modulated' key problem

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