source: trunk/GSASIIstrIO.py @ 2450

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

fix error in matrix in getCellEsd & RetDistAngle?
Add new atom drawing option for sphere of enclosure

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 139.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2016-08-25 17:22:04 +0000 (Thu, 25 Aug 2016) $
4# $Author: vondreele $
5# $Revision: 2450 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 2450 2016-08-25 17:22:04Z 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: 2450 $")
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([
1483        Ax[1]*Ax[2]-Ax[5]**2,
1484        Ax[0]*Ax[2]-Ax[4]**2,
1485        Ax[0]*Ax[1]-Ax[3]**2,
1486        Ax[4]*Ax[5]-Ax[2]*Ax[3],
1487        Ax[3]*Ax[5]-Ax[1]*Ax[4],
1488        Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1489    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1490    Vol = 1/np.sqrt(rVsq)
1491    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1492    R123 = Ax[0]*Ax[1]*Ax[2]
1493    dsasdg = np.zeros((3,6))
1494    dadg = np.zeros((6,6))
1495    for i0 in range(3):         #0  1   2
1496        i1 = (i0+1)%3           #1  2   0
1497        i2 = (i1+1)%3           #2  0   1
1498        i3 = 5-i2               #3  5   4
1499        i4 = 5-i1               #4  3   5
1500        i5 = 5-i0               #5  4   3
1501        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1502        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1503        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1504        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1505        denom = np.sqrt(denmsq)
1506        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1507        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1508        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1509        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1510        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1511        dadg[i5][i5] = -Ax[i0]/denom
1512    for i0 in range(3):
1513        i1 = (i0+1)%3
1514        i2 = (i1+1)%3
1515        i3 = 5-i2
1516        for ij in range(6):
1517            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1518            if ij == i0:
1519                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1520            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1521    sigMat = np.inner(dadg,np.inner(dadg,vcov))
1522    var = np.diag(sigMat)
1523    CS = np.where(var>0.,np.sqrt(var),0.)
1524    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1525   
1526def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1527    'needs a doc string'
1528   
1529    def PrintAtomsAndSig(General,Atoms,atomsSig):
1530        print >>pFile,'\n Atoms:'
1531        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1532        if General['Type'] == 'magnetic':
1533            line = line[:24]+'   Mx     My     Mz'+line[24:]
1534        elif General['Type'] == 'macromolecular':
1535            line = ' res no residue chain '+line
1536        cx,ct,cs,cia = General['AtomPtrs']
1537        print >>pFile,line
1538        print >>pFile,135*'-'
1539        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1540            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1541        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1542        for atyp in General['AtomTypes']:       #zero composition
1543            General['NoAtoms'][atyp] = 0.
1544        for i,at in enumerate(Atoms):
1545            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1546            if General['Type'] == 'macromolecular':
1547                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1548                valstr = ' values:          '
1549                sigstr = ' sig   :          '
1550            else:
1551                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1552                valstr = ' values:'
1553                sigstr = ' sig   :'
1554            for ind in range(cx,cx+4):
1555                sigind = str(i)+':'+str(ind)
1556                valstr += fmt[ind]%(at[ind])                   
1557                if sigind in atomsSig:
1558                    sigstr += fmt[ind]%(atomsSig[sigind])
1559                else:
1560                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1561            if at[cia] == 'I':
1562                valstr += fmt[cia+1]%(at[cia+1])
1563                if '%d:%d'%(i,cia+1) in atomsSig:
1564                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1565                else:
1566                    sigstr += 8*' '
1567            else:
1568                valstr += 8*' '
1569                sigstr += 8*' '
1570                for ind in range(cia+2,cia+8):
1571                    sigind = str(i)+':'+str(ind)
1572                    valstr += fmt[ind]%(at[ind])
1573                    if sigind in atomsSig:                       
1574                        sigstr += fmt[ind]%(atomsSig[sigind])
1575                    else:
1576                        sigstr += 8*' '
1577            print >>pFile,name
1578            print >>pFile,valstr
1579            print >>pFile,sigstr
1580           
1581    def PrintWavesAndSig(General,Atoms,wavesSig):
1582        cx,ct,cs,cia = General['AtomPtrs']
1583        print >>pFile,'\n Modulation waves'
1584        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1585            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1586            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1587        print >>pFile,135*'-'
1588        for i,at in enumerate(Atoms):
1589            AtomSS = at[-1]['SS1']
1590            waveType = AtomSS['waveType']
1591            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1592                Waves = AtomSS[Stype]
1593                if len(Waves):
1594                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
1595                        %(at[ct-1],at[cs],Stype,waveType)
1596                    line = ''
1597                    for iw,wave in enumerate(Waves):
1598                        stiw = ':'+str(i)+':'+str(iw)
1599                        namstr = '  names :'
1600                        valstr = '  values:'
1601                        sigstr = '  esds  :'
1602                        if Stype == 'Spos':
1603                            nt = 6
1604                            ot = 0
1605                            if waveType in ['ZigZag','Block',] and not iw:
1606                                nt = 5
1607                                ot = 6
1608                            for j in range(nt):
1609                                name = names['Spos'][j+ot]
1610                                namstr += '%12s'%(name)
1611                                valstr += '%12.4f'%(wave[0][j])
1612                                if name+stiw in wavesSig:
1613                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1614                                else:
1615                                    sigstr += 12*' '
1616                        elif Stype == 'Sfrac':
1617                            ot = 0
1618                            if 'Crenel' in waveType and not iw:
1619                                ot = 2
1620                            for j in range(2):
1621                                name = names['Sfrac'][j+ot]
1622                                namstr += '%12s'%(names['Sfrac'][j+ot])
1623                                valstr += '%12.4f'%(wave[0][j])
1624                                if name+stiw in wavesSig:
1625                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1626                                else:
1627                                    sigstr += 12*' '
1628                        elif Stype == 'Sadp':
1629                            for j in range(12):
1630                                name = names['Sadp'][j]
1631                                namstr += '%10s'%(names['Sadp'][j])
1632                                valstr += '%10.6f'%(wave[0][j])
1633                                if name+stiw in wavesSig:
1634                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
1635                                else:
1636                                    sigstr += 10*' '
1637                        elif Stype == 'Smag':
1638                            for j in range(6):
1639                                name = names['Smag'][j]
1640                                namstr += '%12s'%(names['Smag'][j])
1641                                valstr += '%12.4f'%(wave[0][j])
1642                                if name+stiw in wavesSig:
1643                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1644                                else:
1645                                    sigstr += 12*' '
1646                               
1647                    print >>pFile,namstr
1648                    print >>pFile,valstr
1649                    print >>pFile,sigstr
1650       
1651               
1652    def PrintRBObjPOAndSig(rbfx,rbsx):
1653        namstr = '  names :'
1654        valstr = '  values:'
1655        sigstr = '  esds  :'
1656        for i,px in enumerate(['Px:','Py:','Pz:']):
1657            name = pfx+rbfx+px+rbsx
1658            namstr += '%12s'%('Pos '+px[1])
1659            valstr += '%12.5f'%(parmDict[name])
1660            if name in sigDict:
1661                sigstr += '%12.5f'%(sigDict[name])
1662            else:
1663                sigstr += 12*' '
1664        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1665            name = pfx+rbfx+po+rbsx
1666            namstr += '%12s'%('Ori '+po[1])
1667            valstr += '%12.5f'%(parmDict[name])
1668            if name in sigDict:
1669                sigstr += '%12.5f'%(sigDict[name])
1670            else:
1671                sigstr += 12*' '
1672        print >>pFile,namstr
1673        print >>pFile,valstr
1674        print >>pFile,sigstr
1675       
1676    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1677        namstr = '  names :'
1678        valstr = '  values:'
1679        sigstr = '  esds  :'
1680        if 'N' not in TLS:
1681            print >>pFile,' Thermal motion:'
1682        if 'T' in TLS:
1683            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1684                name = pfx+rbfx+pt+rbsx
1685                namstr += '%12s'%(pt[:3])
1686                valstr += '%12.5f'%(parmDict[name])
1687                if name in sigDict:
1688                    sigstr += '%12.5f'%(sigDict[name])
1689                else:
1690                    sigstr += 12*' '
1691            print >>pFile,namstr
1692            print >>pFile,valstr
1693            print >>pFile,sigstr
1694        if 'L' in TLS:
1695            namstr = '  names :'
1696            valstr = '  values:'
1697            sigstr = '  esds  :'
1698            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1699                name = pfx+rbfx+pt+rbsx
1700                namstr += '%12s'%(pt[:3])
1701                valstr += '%12.3f'%(parmDict[name])
1702                if name in sigDict:
1703                    sigstr += '%12.3f'%(sigDict[name])
1704                else:
1705                    sigstr += 12*' '
1706            print >>pFile,namstr
1707            print >>pFile,valstr
1708            print >>pFile,sigstr
1709        if 'S' in TLS:
1710            namstr = '  names :'
1711            valstr = '  values:'
1712            sigstr = '  esds  :'
1713            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1714                name = pfx+rbfx+pt+rbsx
1715                namstr += '%12s'%(pt[:3])
1716                valstr += '%12.4f'%(parmDict[name])
1717                if name in sigDict:
1718                    sigstr += '%12.4f'%(sigDict[name])
1719                else:
1720                    sigstr += 12*' '
1721            print >>pFile,namstr
1722            print >>pFile,valstr
1723            print >>pFile,sigstr
1724        if 'U' in TLS:
1725            name = pfx+rbfx+'U:'+rbsx
1726            namstr = '  names :'+'%12s'%('Uiso')
1727            valstr = '  values:'+'%12.5f'%(parmDict[name])
1728            if name in sigDict:
1729                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1730            else:
1731                sigstr = '  esds  :'+12*' '
1732            print >>pFile,namstr
1733            print >>pFile,valstr
1734            print >>pFile,sigstr
1735       
1736    def PrintRBObjTorAndSig(rbsx):
1737        namstr = '  names :'
1738        valstr = '  values:'
1739        sigstr = '  esds  :'
1740        nTors = len(RBObj['Torsions'])
1741        if nTors:
1742            print >>pFile,' Torsions:'
1743            for it in range(nTors):
1744                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1745                namstr += '%12s'%('Tor'+str(it))
1746                valstr += '%12.4f'%(parmDict[name])
1747                if name in sigDict:
1748                    sigstr += '%12.4f'%(sigDict[name])
1749            print >>pFile,namstr
1750            print >>pFile,valstr
1751            print >>pFile,sigstr
1752               
1753    def PrintSHtextureAndSig(textureData,SHtextureSig):
1754        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1755        names = ['omega','chi','phi']
1756        namstr = '  names :'
1757        ptstr =  '  values:'
1758        sigstr = '  esds  :'
1759        for name in names:
1760            namstr += '%12s'%(name)
1761            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1762            if 'Sample '+name in SHtextureSig:
1763                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1764            else:
1765                sigstr += 12*' '
1766        print >>pFile,namstr
1767        print >>pFile,ptstr
1768        print >>pFile,sigstr
1769        print >>pFile,'\n Texture coefficients:'
1770        SHcoeff = textureData['SH Coeff'][1]
1771        SHkeys = SHcoeff.keys()
1772        nCoeff = len(SHcoeff)
1773        nBlock = nCoeff/10+1
1774        iBeg = 0
1775        iFin = min(iBeg+10,nCoeff)
1776        for block in range(nBlock):
1777            namstr = '  names :'
1778            ptstr =  '  values:'
1779            sigstr = '  esds  :'
1780            for name in SHkeys[iBeg:iFin]:
1781                namstr += '%12s'%(name)
1782                ptstr += '%12.3f'%(SHcoeff[name])
1783                if name in SHtextureSig:
1784                    sigstr += '%12.3f'%(SHtextureSig[name])
1785                else:
1786                    sigstr += 12*' '
1787            print >>pFile,namstr
1788            print >>pFile,ptstr
1789            print >>pFile,sigstr
1790            iBeg += 10
1791            iFin = min(iBeg+10,nCoeff)
1792           
1793    print >>pFile,'\n Phases:'
1794    for phase in Phases:
1795        print >>pFile,' Result for phase: ',phase
1796        Phase = Phases[phase]
1797        General = Phase['General']
1798        SGData = General['SGData']
1799        Atoms = Phase['Atoms']
1800        if Atoms and not General.get('doPawley'):
1801            cx,ct,cs,cia = General['AtomPtrs']
1802            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1803        cell = General['Cell']
1804        pId = Phase['pId']
1805        pfx = str(pId)+'::'
1806        if cell[0]:
1807            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1808            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1809            print >>pFile,' Reciprocal metric tensor: '
1810            ptfmt = "%15.9f"
1811            names = ['A11','A22','A33','A12','A13','A23']
1812            namstr = '  names :'
1813            ptstr =  '  values:'
1814            sigstr = '  esds  :'
1815            for name,a,siga in zip(names,A,sigA):
1816                namstr += '%15s'%(name)
1817                ptstr += ptfmt%(a)
1818                if siga:
1819                    sigstr += ptfmt%(siga)
1820                else:
1821                    sigstr += 15*' '
1822            print >>pFile,namstr
1823            print >>pFile,ptstr
1824            print >>pFile,sigstr
1825            cell[1:7] = G2lat.A2cell(A)
1826            cell[7] = G2lat.calc_V(A)
1827            print >>pFile,' New unit cell:'
1828            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1829            names = ['a','b','c','alpha','beta','gamma','Volume']
1830            namstr = '  names :'
1831            ptstr =  '  values:'
1832            sigstr = '  esds  :'
1833            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1834                namstr += '%12s'%(name)
1835                ptstr += fmt%(a)
1836                if siga:
1837                    sigstr += fmt%(siga)
1838                else:
1839                    sigstr += 12*' '
1840            print >>pFile,namstr
1841            print >>pFile,ptstr
1842            print >>pFile,sigstr
1843        ik = 6  #for Pawley stuff below
1844        if General.get('Modulated',False):
1845            ik = 7
1846            Vec,vRef,maxH = General['SuperVec']
1847            if vRef:
1848                print >>pFile,' New modulation vector:'
1849                namstr = '  names :'
1850                ptstr =  '  values:'
1851                sigstr = '  esds  :'
1852                for var in ['mV0','mV1','mV2']:
1853                    namstr += '%12s'%(pfx+var)
1854                    ptstr += '%12.6f'%(parmDict[pfx+var])
1855                    if pfx+var in sigDict:
1856                        sigstr += '%12.6f'%(sigDict[pfx+var])
1857                    else:
1858                        sigstr += 12*' '
1859                print >>pFile,namstr
1860                print >>pFile,ptstr
1861                print >>pFile,sigstr
1862           
1863        General['Mass'] = 0.
1864        if Phase['General'].get('doPawley'):
1865            pawleyRef = Phase['Pawley ref']
1866            for i,refl in enumerate(pawleyRef):
1867                key = pfx+'PWLref:'+str(i)
1868                refl[ik] = parmDict[key]
1869                if key in sigDict:
1870                    refl[ik+1] = sigDict[key]
1871                else:
1872                    refl[ik+1] = 0
1873        else:
1874            VRBIds = RBIds['Vector']
1875            RRBIds = RBIds['Residue']
1876            RBModels = Phase['RBModels']
1877            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1878                jrb = VRBIds.index(RBObj['RBId'])
1879                rbsx = str(irb)+':'+str(jrb)
1880                print >>pFile,' Vector rigid body parameters:'
1881                PrintRBObjPOAndSig('RBV',rbsx)
1882                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1883            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1884                jrb = RRBIds.index(RBObj['RBId'])
1885                rbsx = str(irb)+':'+str(jrb)
1886                print >>pFile,' Residue rigid body parameters:'
1887                PrintRBObjPOAndSig('RBR',rbsx)
1888                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1889                PrintRBObjTorAndSig(rbsx)
1890            atomsSig = {}
1891            wavesSig = {}
1892            cx,ct,cs,cia = General['AtomPtrs']
1893            for i,at in enumerate(Atoms):
1894                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
1895                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
1896                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i)}
1897                for ind in range(cx,cx+4):
1898                    at[ind] = parmDict[names[ind]]
1899                    if ind in range(cx,cx+3):
1900                        name = names[ind].replace('A','dA')
1901                    else:
1902                        name = names[ind]
1903                    if name in sigDict:
1904                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1905                if at[cia] == 'I':
1906                    at[cia+1] = parmDict[names[cia+1]]
1907                    if names[cia+1] in sigDict:
1908                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
1909                else:
1910                    for ind in range(cia+2,cia+8):
1911                        at[ind] = parmDict[names[ind]]
1912                        if names[ind] in sigDict:
1913                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1914                ind = General['AtomTypes'].index(at[ct])
1915                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
1916                if General.get('Modulated',False):
1917                    AtomSS = at[-1]['SS1']
1918                    waveType = AtomSS['waveType']
1919                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1920                        Waves = AtomSS[Stype]
1921                        for iw,wave in enumerate(Waves):
1922                            stiw = str(i)+':'+str(iw)
1923                            if Stype == 'Spos':
1924                                if waveType in ['ZigZag','Block',] and not iw:
1925                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
1926                                else:
1927                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
1928                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
1929                            elif Stype == 'Sadp':
1930                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
1931                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
1932                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
1933                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
1934                            elif Stype == 'Sfrac':
1935                                if 'Crenel' in waveType and not iw:
1936                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
1937                                else:
1938                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
1939                            elif Stype == 'Smag':
1940                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
1941                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
1942                            for iname,name in enumerate(names):
1943                                AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
1944                                if pfx+name in sigDict:
1945                                    wavesSig[name] = sigDict[pfx+name]
1946                   
1947            PrintAtomsAndSig(General,Atoms,atomsSig)
1948            if General.get('Modulated',False):
1949                PrintWavesAndSig(General,Atoms,wavesSig)
1950           
1951       
1952        textureData = General['SH Texture']   
1953        if textureData['Order']:
1954            SHtextureSig = {}
1955            for name in ['omega','chi','phi']:
1956                aname = pfx+'SH '+name
1957                textureData['Sample '+name][1] = parmDict[aname]
1958                if aname in sigDict:
1959                    SHtextureSig['Sample '+name] = sigDict[aname]
1960            for name in textureData['SH Coeff'][1]:
1961                aname = pfx+name
1962                textureData['SH Coeff'][1][name] = parmDict[aname]
1963                if aname in sigDict:
1964                    SHtextureSig[name] = sigDict[aname]
1965            PrintSHtextureAndSig(textureData,SHtextureSig)
1966        if phase in RestraintDict and not Phase['General'].get('doPawley'):
1967            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1968                textureData,RestraintDict[phase],pFile)
1969                   
1970################################################################################
1971##### Histogram & Phase data
1972################################################################################       
1973                   
1974def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1975    '''Loads the HAP histogram/phase information into dicts
1976
1977    :param dict Phases: phase information
1978    :param dict Histograms: Histogram information
1979    :param bool Print: prints information as it is read
1980    :param file pFile: file object to print to (the default, None causes printing to the console)
1981    :param bool resetRefList: Should the contents of the Reflection List be initialized
1982      on loading. The default, True, initializes the Reflection List as it is loaded.
1983
1984    :returns: (hapVary,hapDict,controlDict)
1985      * hapVary: list of refined variables
1986      * hapDict: dict with refined variables and their values
1987      * controlDict: dict with fixed parameters
1988    '''
1989   
1990    def PrintSize(hapData):
1991        if hapData[0] in ['isotropic','uniaxial']:
1992            line = '\n Size model    : %9s'%(hapData[0])
1993            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1994            if hapData[0] == 'uniaxial':
1995                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1996            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1997            print >>pFile,line
1998        else:
1999            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
2000                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2001            Snames = ['S11','S22','S33','S12','S13','S23']
2002            ptlbls = ' names :'
2003            ptstr =  ' values:'
2004            varstr = ' refine:'
2005            for i,name in enumerate(Snames):
2006                ptlbls += '%12s' % (name)
2007                ptstr += '%12.6f' % (hapData[4][i])
2008                varstr += '%12s' % (str(hapData[5][i]))
2009            print >>pFile,ptlbls
2010            print >>pFile,ptstr
2011            print >>pFile,varstr
2012       
2013    def PrintMuStrain(hapData,SGData):
2014        if hapData[0] in ['isotropic','uniaxial']:
2015            line = '\n Mustrain model: %9s'%(hapData[0])
2016            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2017            if hapData[0] == 'uniaxial':
2018                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2019            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2020            print >>pFile,line
2021        else:
2022            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
2023                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2024            Snames = G2spc.MustrainNames(SGData)
2025            ptlbls = ' names :'
2026            ptstr =  ' values:'
2027            varstr = ' refine:'
2028            for i,name in enumerate(Snames):
2029                ptlbls += '%12s' % (name)
2030                ptstr += '%12.6f' % (hapData[4][i])
2031                varstr += '%12s' % (str(hapData[5][i]))
2032            print >>pFile,ptlbls
2033            print >>pFile,ptstr
2034            print >>pFile,varstr
2035
2036    def PrintHStrain(hapData,SGData):
2037        print >>pFile,'\n Hydrostatic/elastic strain: '
2038        Hsnames = G2spc.HStrainNames(SGData)
2039        ptlbls = ' names :'
2040        ptstr =  ' values:'
2041        varstr = ' refine:'
2042        for i,name in enumerate(Hsnames):
2043            ptlbls += '%12s' % (name)
2044            ptstr += '%12.4g' % (hapData[0][i])
2045            varstr += '%12s' % (str(hapData[1][i]))
2046        print >>pFile,ptlbls
2047        print >>pFile,ptstr
2048        print >>pFile,varstr
2049
2050    def PrintSHPO(hapData):
2051        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
2052            str(hapData[4])+' Refine? '+str(hapData[2])
2053        ptlbls = ' names :'
2054        ptstr =  ' values:'
2055        for item in hapData[5]:
2056            ptlbls += '%12s'%(item)
2057            ptstr += '%12.3f'%(hapData[5][item]) 
2058        print >>pFile,ptlbls
2059        print >>pFile,ptstr
2060   
2061    def PrintBabinet(hapData):
2062        print >>pFile,'\n Babinet form factor modification: '
2063        ptlbls = ' names :'
2064        ptstr =  ' values:'
2065        varstr = ' refine:'
2066        for name in ['BabA','BabU']:
2067            ptlbls += '%12s' % (name)
2068            ptstr += '%12.6f' % (hapData[name][0])
2069            varstr += '%12s' % (str(hapData[name][1]))
2070        print >>pFile,ptlbls
2071        print >>pFile,ptstr
2072        print >>pFile,varstr
2073       
2074    hapDict = {}
2075    hapVary = []
2076    controlDict = {}
2077    poType = {}
2078    poAxes = {}
2079    spAxes = {}
2080    spType = {}
2081   
2082    for phase in Phases:
2083        HistoPhase = Phases[phase]['Histograms']
2084        SGData = Phases[phase]['General']['SGData']
2085        cell = Phases[phase]['General']['Cell'][1:7]
2086        A = G2lat.cell2A(cell)
2087        if Phases[phase]['General'].get('Modulated',False):
2088            SSGData = Phases[phase]['General']['SSGData']
2089            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2090        pId = Phases[phase]['pId']
2091        histoList = HistoPhase.keys()
2092        histoList.sort()
2093        for histogram in histoList:
2094            try:
2095                Histogram = Histograms[histogram]
2096            except KeyError:                       
2097                #skip if histogram not included e.g. in a sequential refinement
2098                continue
2099            hapData = HistoPhase[histogram]
2100            hId = Histogram['hId']
2101            if 'PWDR' in histogram:
2102                limits = Histogram['Limits'][1]
2103                inst = Histogram['Instrument Parameters'][0]
2104                Zero = inst['Zero'][0]
2105                if 'C' in inst['Type'][1]:
2106                    try:
2107                        wave = inst['Lam'][1]
2108                    except KeyError:
2109                        wave = inst['Lam1'][1]
2110                    dmin = wave/(2.0*sind(limits[1]/2.0))
2111                elif 'T' in inst['Type'][0]:
2112                    dmin = limits[0]/inst['difC'][1]
2113                pfx = str(pId)+':'+str(hId)+':'
2114                for item in ['Scale','Extinction']:
2115                    hapDict[pfx+item] = hapData[item][0]
2116                    if hapData[item][1]:
2117                        hapVary.append(pfx+item)
2118                names = G2spc.HStrainNames(SGData)
2119                HSvals = []
2120                for i,name in enumerate(names):
2121                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2122                    HSvals.append(hapDict[pfx+name])
2123                    if hapData['HStrain'][1][i]:
2124                        hapVary.append(pfx+name)
2125                DIJS = G2spc.HStrainVals(HSvals,SGData)
2126                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2127                if hapData['Pref.Ori.'][0] == 'MD':
2128                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2129                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2130                    if hapData['Pref.Ori.'][2]:
2131                        hapVary.append(pfx+'MD')
2132                else:                           #'SH' spherical harmonics
2133                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2134                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2135                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2136                    controlDict[pfx+'SHhkl'] = []
2137                    try: #patch for old Pref.Ori. items
2138                        controlDict[pfx+'SHtoler'] = 0.1
2139                        if hapData['Pref.Ori.'][6][0] != '':
2140                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2141                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2142                    except IndexError:
2143                        pass
2144                    for item in hapData['Pref.Ori.'][5]:
2145                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2146                        if hapData['Pref.Ori.'][2]:
2147                            hapVary.append(pfx+item)
2148                for item in ['Mustrain','Size']:
2149                    controlDict[pfx+item+'Type'] = hapData[item][0]
2150                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2151                    if hapData[item][2][2]:
2152                        hapVary.append(pfx+item+';mx')
2153                    if hapData[item][0] in ['isotropic','uniaxial']:
2154                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2155                        if hapData[item][2][0]:
2156                            hapVary.append(pfx+item+';i')
2157                        if hapData[item][0] == 'uniaxial':
2158                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2159                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2160                            if hapData[item][2][1]:
2161                                hapVary.append(pfx+item+';a')
2162                    else:       #generalized for mustrain or ellipsoidal for size
2163                        Nterms = len(hapData[item][4])
2164                        if item == 'Mustrain':
2165                            names = G2spc.MustrainNames(SGData)
2166                            pwrs = []
2167                            for name in names:
2168                                h,k,l = name[1:]
2169                                pwrs.append([int(h),int(k),int(l)])
2170                            controlDict[pfx+'MuPwrs'] = pwrs
2171                        for i in range(Nterms):
2172                            sfx = ';'+str(i)
2173                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2174                            if hapData[item][5][i]:
2175                                hapVary.append(pfx+item+sfx)
2176                for bab in ['BabA','BabU']:
2177                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2178                    if hapData['Babinet'][bab][1]:
2179                        hapVary.append(pfx+bab)
2180                               
2181                if Print: 
2182                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2183                    print >>pFile,135*'-'
2184                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2185                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2186                    if hapData['Pref.Ori.'][0] == 'MD':
2187                        Ax = hapData['Pref.Ori.'][3]
2188                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2189                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2190                    else: #'SH' for spherical harmonics
2191                        PrintSHPO(hapData['Pref.Ori.'])
2192                        print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2193                    PrintSize(hapData['Size'])
2194                    PrintMuStrain(hapData['Mustrain'],SGData)
2195                    PrintHStrain(hapData['HStrain'],SGData)
2196                    if hapData['Babinet']['BabA'][0]:
2197                        PrintBabinet(hapData['Babinet'])                       
2198                if resetRefList:
2199                    refList = []
2200                    Uniq = []
2201                    Phi = []
2202                    if Phases[phase]['General'].get('Modulated',False):
2203                        ifSuper = True
2204                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2205                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2206                        for h,k,l,m,d in HKLd:
2207                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2208                            mul *= 2      # for powder overlap of Friedel pairs
2209                            if m or not ext:
2210                                if 'C' in inst['Type'][0]:
2211                                    pos = G2lat.Dsp2pos(inst,d)
2212                                    if limits[0] < pos < limits[1]:
2213                                        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])
2214                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2215                                        Uniq.append(uniq)
2216                                        Phi.append(phi)
2217                                elif 'T' in inst['Type'][0]:
2218                                    pos = G2lat.Dsp2pos(inst,d)
2219                                    if limits[0] < pos < limits[1]:
2220                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2221                                        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])
2222                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2223                                        Uniq.append(uniq)
2224                                        Phi.append(phi)
2225                    else:
2226                        ifSuper = False
2227                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2228                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2229                        for h,k,l,d in HKLd:
2230                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2231                            mul *= 2      # for powder overlap of Friedel pairs
2232                            if ext:
2233                                continue
2234                            if 'C' in inst['Type'][0]:
2235                                pos = G2lat.Dsp2pos(inst,d)
2236                                if limits[0] < pos < limits[1]:
2237                                    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])
2238                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2239                                    Uniq.append(uniq)
2240                                    Phi.append(phi)
2241                            elif 'T' in inst['Type'][0]:
2242                                pos = G2lat.Dsp2pos(inst,d)
2243                                if limits[0] < pos < limits[1]:
2244                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2245                                    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])
2246                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2247                                    Uniq.append(uniq)
2248                                    Phi.append(phi)
2249                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2250            elif 'HKLF' in histogram:
2251                inst = Histogram['Instrument Parameters'][0]
2252                hId = Histogram['hId']
2253                hfx = ':%d:'%(hId)
2254                for item in inst:
2255                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2256                    hapDict[hfx+item] = inst[item][1]
2257                pfx = str(pId)+':'+str(hId)+':'
2258                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2259                if hapData['Scale'][1]:
2260                    hapVary.append(pfx+'Scale')
2261                               
2262                extApprox,extType,extParms = hapData['Extinction']
2263                controlDict[pfx+'EType'] = extType
2264                controlDict[pfx+'EApprox'] = extApprox
2265                if 'C' in inst['Type'][0]:
2266                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2267                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2268                if 'Primary' in extType:
2269                    Ekey = ['Ep',]
2270                elif 'I & II' in extType:
2271                    Ekey = ['Eg','Es']
2272                elif 'Secondary Type II' == extType:
2273                    Ekey = ['Es',]
2274                elif 'Secondary Type I' == extType:
2275                    Ekey = ['Eg',]
2276                else:   #'None'
2277                    Ekey = []
2278                for eKey in Ekey:
2279                    hapDict[pfx+eKey] = extParms[eKey][0]
2280                    if extParms[eKey][1]:
2281                        hapVary.append(pfx+eKey)
2282                for bab in ['BabA','BabU']:
2283                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2284                    if hapData['Babinet'][bab][1]:
2285                        hapVary.append(pfx+bab)
2286                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2287                if len(Twins) == 1:
2288                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2289                    if hapData.get('Flack',[0,False])[1]:
2290                        hapVary.append(pfx+'Flack')
2291                sumTwFr = 0.
2292                controlDict[pfx+'TwinLaw'] = []
2293                controlDict[pfx+'TwinInv'] = []
2294                NTL = 0           
2295                for it,twin in enumerate(Twins):
2296                    if 'bool' in str(type(twin[0])):
2297                        controlDict[pfx+'TwinInv'].append(twin[0])
2298                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2299                    else:
2300                        NTL += 1
2301                        controlDict[pfx+'TwinInv'].append(False)
2302                        controlDict[pfx+'TwinLaw'].append(twin[0])
2303                    if it:
2304                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2305                        sumTwFr += twin[1]
2306                    else:
2307                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2308                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2309                    if Twins[0][1][1]:
2310                        hapVary.append(pfx+'TwinFr:'+str(it))
2311                controlDict[pfx+'NTL'] = NTL
2312                #need to make constraint on TwinFr
2313                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2314                if len(Twins) > 1:    #force sum to unity
2315                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2316                if Print: 
2317                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2318                    print >>pFile,135*'-'
2319                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2320                    if extType != 'None':
2321                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox)
2322                        text = ' Parameters       :'
2323                        for eKey in Ekey:
2324                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2325                        print >>pFile,text
2326                    if hapData['Babinet']['BabA'][0]:
2327                        PrintBabinet(hapData['Babinet'])
2328                    if not SGData['SGInv'] and len(Twins) == 1:
2329                        print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2330                    if len(Twins) > 1:
2331                        for it,twin in enumerate(Twins):
2332                            if 'bool' in str(type(twin[0])):
2333                                print >>pFile,' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? '%(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1] 
2334                            else:
2335                                print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr:'+str(it)]),Twins[0][1][1] 
2336                       
2337                Histogram['Reflection Lists'] = phase       
2338               
2339    return hapVary,hapDict,controlDict
2340   
2341def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2342    'needs a doc string'
2343   
2344    def PrintSizeAndSig(hapData,sizeSig):
2345        line = '\n Size model:     %9s'%(hapData[0])
2346        refine = False
2347        if hapData[0] in ['isotropic','uniaxial']:
2348            line += ' equatorial:%12.5f'%(hapData[1][0])
2349            if sizeSig[0][0]:
2350                line += ', sig:%8.4f'%(sizeSig[0][0])
2351                refine = True
2352            if hapData[0] == 'uniaxial':
2353                line += ' axial:%12.4f'%(hapData[1][1])
2354                if sizeSig[0][1]:
2355                    refine = True
2356                    line += ', sig:%8.4f'%(sizeSig[0][1])
2357            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2358            if sizeSig[0][2]:
2359                refine = True
2360                line += ', sig:%8.4f'%(sizeSig[0][2])
2361            if refine:
2362                print >>pFile,line
2363        else:
2364            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2365            if sizeSig[0][2]:
2366                refine = True
2367                line += ', sig:%8.4f'%(sizeSig[0][2])
2368            Snames = ['S11','S22','S33','S12','S13','S23']
2369            ptlbls = ' name  :'
2370            ptstr =  ' value :'
2371            sigstr = ' sig   :'
2372            for i,name in enumerate(Snames):
2373                ptlbls += '%12s' % (name)
2374                ptstr += '%12.6f' % (hapData[4][i])
2375                if sizeSig[1][i]:
2376                    refine = True
2377                    sigstr += '%12.6f' % (sizeSig[1][i])
2378                else:
2379                    sigstr += 12*' '
2380            if refine:
2381                print >>pFile,line
2382                print >>pFile,ptlbls
2383                print >>pFile,ptstr
2384                print >>pFile,sigstr
2385       
2386    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2387        line = '\n Mustrain model: %9s'%(hapData[0])
2388        refine = False
2389        if hapData[0] in ['isotropic','uniaxial']:
2390            line += ' equatorial:%12.1f'%(hapData[1][0])
2391            if mustrainSig[0][0]:
2392                line += ', sig:%8.1f'%(mustrainSig[0][0])
2393                refine = True
2394            if hapData[0] == 'uniaxial':
2395                line += ' axial:%12.1f'%(hapData[1][1])
2396                if mustrainSig[0][1]:
2397                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2398            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2399            if mustrainSig[0][2]:
2400                refine = True
2401                line += ', sig:%8.3f'%(mustrainSig[0][2])
2402            if refine:
2403                print >>pFile,line
2404        else:
2405            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2406            if mustrainSig[0][2]:
2407                refine = True
2408                line += ', sig:%8.3f'%(mustrainSig[0][2])
2409            Snames = G2spc.MustrainNames(SGData)
2410            ptlbls = ' name  :'
2411            ptstr =  ' value :'
2412            sigstr = ' sig   :'
2413            for i,name in enumerate(Snames):
2414                ptlbls += '%12s' % (name)
2415                ptstr += '%12.1f' % (hapData[4][i])
2416                if mustrainSig[1][i]:
2417                    refine = True
2418                    sigstr += '%12.1f' % (mustrainSig[1][i])
2419                else:
2420                    sigstr += 12*' '
2421            if refine:
2422                print >>pFile,line
2423                print >>pFile,ptlbls
2424                print >>pFile,ptstr
2425                print >>pFile,sigstr
2426           
2427    def PrintHStrainAndSig(hapData,strainSig,SGData):
2428        Hsnames = G2spc.HStrainNames(SGData)
2429        ptlbls = ' name  :'
2430        ptstr =  ' value :'
2431        sigstr = ' sig   :'
2432        refine = False
2433        for i,name in enumerate(Hsnames):
2434            ptlbls += '%12s' % (name)
2435            ptstr += '%12.4g' % (hapData[0][i])
2436            if name in strainSig:
2437                refine = True
2438                sigstr += '%12.4g' % (strainSig[name])
2439            else:
2440                sigstr += 12*' '
2441        if refine:
2442            print >>pFile,'\n Hydrostatic/elastic strain: '
2443            print >>pFile,ptlbls
2444            print >>pFile,ptstr
2445            print >>pFile,sigstr
2446       
2447    def PrintSHPOAndSig(pfx,hapData,POsig):
2448        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2449        ptlbls = ' names :'
2450        ptstr =  ' values:'
2451        sigstr = ' sig   :'
2452        for item in hapData[5]:
2453            ptlbls += '%12s'%(item)
2454            ptstr += '%12.3f'%(hapData[5][item])
2455            if pfx+item in POsig:
2456                sigstr += '%12.3f'%(POsig[pfx+item])
2457            else:
2458                sigstr += 12*' ' 
2459        print >>pFile,ptlbls
2460        print >>pFile,ptstr
2461        print >>pFile,sigstr
2462       
2463    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2464        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2465        text = ''
2466        for item in hapData[2]:
2467            if pfx+item in ScalExtSig:
2468                text += '       %s: '%(item)
2469                text += '%12.2e'%(hapData[2][item][0])
2470                if pfx+item in ScalExtSig:
2471                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2472        print >>pFile,text       
2473       
2474    def PrintBabinetAndSig(pfx,hapData,BabSig):
2475        print >>pFile,'\n Babinet form factor modification: '
2476        ptlbls = ' names :'
2477        ptstr =  ' values:'
2478        sigstr = ' sig   :'
2479        for item in hapData:
2480            ptlbls += '%12s'%(item)
2481            ptstr += '%12.3f'%(hapData[item][0])
2482            if pfx+item in BabSig:
2483                sigstr += '%12.3f'%(BabSig[pfx+item])
2484            else:
2485                sigstr += 12*' ' 
2486        print >>pFile,ptlbls
2487        print >>pFile,ptstr
2488        print >>pFile,sigstr
2489       
2490    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2491        print >>pFile,'\n Twin Law fractions : '
2492        ptlbls = ' names :'
2493        ptstr =  ' values:'
2494        sigstr = ' sig   :'
2495        for it,item in enumerate(twinData):
2496            ptlbls += '     twin #%d'%(it)
2497            if it:
2498                ptstr += '%12.3f'%(item[1])
2499            else:
2500                ptstr += '%12.3f'%(item[1][0])
2501            if pfx+'TwinFr:'+str(it) in TwinSig:
2502                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2503            else:
2504                sigstr += 12*' ' 
2505        print >>pFile,ptlbls
2506        print >>pFile,ptstr
2507        print >>pFile,sigstr
2508       
2509   
2510    PhFrExtPOSig = {}
2511    SizeMuStrSig = {}
2512    ScalExtSig = {}
2513    BabSig = {}
2514    TwinFrSig = {}
2515    wtFrSum = {}
2516    for phase in Phases:
2517        HistoPhase = Phases[phase]['Histograms']
2518        General = Phases[phase]['General']
2519        SGData = General['SGData']
2520        pId = Phases[phase]['pId']
2521        histoList = HistoPhase.keys()
2522        histoList.sort()
2523        for histogram in histoList:
2524            try:
2525                Histogram = Histograms[histogram]
2526            except KeyError:                       
2527                #skip if histogram not included e.g. in a sequential refinement
2528                continue
2529            hapData = HistoPhase[histogram]
2530            hId = Histogram['hId']
2531            pfx = str(pId)+':'+str(hId)+':'
2532            if hId not in wtFrSum:
2533                wtFrSum[hId] = 0.
2534            if 'PWDR' in histogram:
2535                for item in ['Scale','Extinction']:
2536                    hapData[item][0] = parmDict[pfx+item]
2537                    if pfx+item in sigDict:
2538                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2539                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2540                if hapData['Pref.Ori.'][0] == 'MD':
2541                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2542                    if pfx+'MD' in sigDict:
2543                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2544                else:                           #'SH' spherical harmonics
2545                    for item in hapData['Pref.Ori.'][5]:
2546                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2547                        if pfx+item in sigDict:
2548                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2549                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2550                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2551                    pfx+'HStrain':{}})                 
2552                for item in ['Mustrain','Size']:
2553                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2554#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2555                    if pfx+item+';mx' in sigDict:
2556                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2557                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2558                        hapData[item][1][0] = parmDict[pfx+item+';i']
2559                        if item == 'Size':
2560                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2561                        if pfx+item+';i' in sigDict: 
2562                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2563                        if hapData[item][0] == 'uniaxial':
2564                            hapData[item][1][1] = parmDict[pfx+item+';a']
2565                            if item == 'Size':
2566                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2567                            if pfx+item+';a' in sigDict:
2568                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2569                    else:       #generalized for mustrain or ellipsoidal for size
2570                        Nterms = len(hapData[item][4])
2571                        for i in range(Nterms):
2572                            sfx = ';'+str(i)
2573                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2574                            if pfx+item+sfx in sigDict:
2575                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2576                names = G2spc.HStrainNames(SGData)
2577                for i,name in enumerate(names):
2578                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2579                    if pfx+name in sigDict:
2580                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2581                for name in ['BabA','BabU']:
2582                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2583                    if pfx+name in sigDict:
2584                        BabSig[pfx+name] = sigDict[pfx+name]               
2585               
2586            elif 'HKLF' in histogram:
2587                for item in ['Scale','Flack']:
2588                    if parmDict.get(pfx+item):
2589                        hapData[item][0] = parmDict[pfx+item]
2590                        if pfx+item in sigDict:
2591                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2592                for item in ['Ep','Eg','Es']:
2593                    if parmDict.get(pfx+item):
2594                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2595                        if pfx+item in sigDict:
2596                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2597                for item in ['BabA','BabU']:
2598                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2599                    if pfx+item in sigDict:
2600                        BabSig[pfx+item] = sigDict[pfx+item]
2601                if 'Twins' in hapData:
2602                    it = 1
2603                    sumTwFr = 0.
2604                    while True:
2605                        try:
2606                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2607                            if pfx+'TwinFr:'+str(it) in sigDict:
2608                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2609                            if it:
2610                                sumTwFr += hapData['Twins'][it][1]
2611                            it += 1
2612                        except KeyError:
2613                            break
2614                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2615
2616    if Print:
2617        for phase in Phases:
2618            HistoPhase = Phases[phase]['Histograms']
2619            General = Phases[phase]['General']
2620            SGData = General['SGData']
2621            pId = Phases[phase]['pId']
2622            histoList = HistoPhase.keys()
2623            histoList.sort()
2624            for histogram in histoList:
2625                try:
2626                    Histogram = Histograms[histogram]
2627                except KeyError:                       
2628                    #skip if histogram not included e.g. in a sequential refinement
2629                    continue
2630                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2631                print >>pFile,130*'-'
2632                hapData = HistoPhase[histogram]
2633                hId = Histogram['hId']
2634                Histogram['Residuals'][str(pId)+'::Name'] = phase
2635                pfx = str(pId)+':'+str(hId)+':'
2636                hfx = ':%s:'%(hId)
2637                if 'PWDR' in histogram:
2638                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2639                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2640                    print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2641               
2642                    if pfx+'Scale' in PhFrExtPOSig:
2643                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2644                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2645                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2646                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2647                    if pfx+'Extinction' in PhFrExtPOSig:
2648                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2649                    if hapData['Pref.Ori.'][0] == 'MD':
2650                        if pfx+'MD' in PhFrExtPOSig:
2651                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2652                    else:
2653                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2654                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2655                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2656                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2657                    if len(BabSig):
2658                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2659                   
2660                elif 'HKLF' in histogram:
2661                    Inst = Histogram['Instrument Parameters'][0]
2662                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)'   \
2663                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2664                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2665                    if FFtables != None and 'N' not in Inst['Type'][0]:
2666                        PrintFprime(FFtables,hfx,pFile)
2667                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2668                    if pfx+'Scale' in ScalExtSig:
2669                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2670                    if hapData['Extinction'][0] != 'None':
2671                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2672                    if len(BabSig):
2673                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2674                    if pfx+'Flack' in ScalExtSig:
2675                        print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2676                    if len(TwinFrSig):
2677                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2678
2679################################################################################
2680##### Histogram data
2681################################################################################       
2682                   
2683def GetHistogramData(Histograms,Print=True,pFile=None):
2684    'needs a doc string'
2685   
2686    def GetBackgroundParms(hId,Background):
2687        Back = Background[0]
2688        DebyePeaks = Background[1]
2689        bakType,bakFlag = Back[:2]
2690        backVals = Back[3:]
2691        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2692        backDict = dict(zip(backNames,backVals))
2693        backVary = []
2694        if bakFlag:
2695            backVary = backNames
2696        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2697        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2698        debyeDict = {}
2699        debyeList = []
2700        for i in range(DebyePeaks['nDebye']):
2701            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2702            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2703            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2704        debyeVary = []
2705        for item in debyeList:
2706            if item[1]:
2707                debyeVary.append(item[0])
2708        backDict.update(debyeDict)
2709        backVary += debyeVary
2710        peakDict = {}
2711        peakList = []
2712        for i in range(DebyePeaks['nPeaks']):
2713            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2714                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2715            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2716            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2717        peakVary = []
2718        for item in peakList:
2719            if item[1]:
2720                peakVary.append(item[0])
2721        backDict.update(peakDict)
2722        backVary += peakVary
2723        return bakType,backDict,backVary           
2724       
2725    def GetInstParms(hId,Inst):     
2726        dataType = Inst['Type'][0]
2727        instDict = {}
2728        insVary = []
2729        pfx = ':'+str(hId)+':'
2730        insKeys = Inst.keys()
2731        insKeys.sort()
2732        for item in insKeys:
2733            insName = pfx+item
2734            instDict[insName] = Inst[item][1]
2735            if len(Inst[item]) > 2 and Inst[item][2]:
2736                insVary.append(insName)
2737        if 'C' in dataType:
2738            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2739        return dataType,instDict,insVary
2740       
2741    def GetSampleParms(hId,Sample):
2742        sampVary = []
2743        hfx = ':'+str(hId)+':'       
2744        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2745            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2746        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2747            if key in Sample:
2748                sampDict[hfx+key] = Sample[key]
2749        Type = Sample['Type']
2750        if 'Bragg' in Type:             #Bragg-Brentano
2751            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2752                sampDict[hfx+item] = Sample[item][0]
2753                if Sample[item][1]:
2754                    sampVary.append(hfx+item)
2755        elif 'Debye' in Type:        #Debye-Scherrer
2756            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2757                sampDict[hfx+item] = Sample[item][0]
2758                if Sample[item][1]:
2759                    sampVary.append(hfx+item)
2760        return Type,sampDict,sampVary
2761       
2762    def PrintBackground(Background):
2763        Back = Background[0]
2764        DebyePeaks = Background[1]
2765        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2766        line = ' Coefficients: '
2767        for i,back in enumerate(Back[3:]):
2768            line += '%10.3f'%(back)
2769            if i and not i%10:
2770                line += '\n'+15*' '
2771        print >>pFile,line
2772        if DebyePeaks['nDebye']:
2773            print >>pFile,'\n Debye diffuse scattering coefficients'
2774            parms = ['DebyeA','DebyeR','DebyeU']
2775            line = ' names :  '
2776            for parm in parms:
2777                line += '%8s refine?'%(parm)
2778            print >>pFile,line
2779            for j,term in enumerate(DebyePeaks['debyeTerms']):
2780                line = ' term'+'%2d'%(j)+':'
2781                for i in range(3):
2782                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2783                print >>pFile,line
2784        if DebyePeaks['nPeaks']:
2785            print >>pFile,'\n Single peak coefficients'
2786            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2787            line = ' names :  '
2788            for parm in parms:
2789                line += '%8s refine?'%(parm)
2790            print >>pFile,line
2791            for j,term in enumerate(DebyePeaks['peaksList']):
2792                line = ' peak'+'%2d'%(j)+':'
2793                for i in range(4):
2794                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2795                print >>pFile,line
2796       
2797    def PrintInstParms(Inst):
2798        print >>pFile,'\n Instrument Parameters:'
2799        insKeys = Inst.keys()
2800        insKeys.sort()
2801        iBeg = 0
2802        Ok = True
2803        while Ok:
2804            ptlbls = ' name  :'
2805            ptstr =  ' value :'
2806            varstr = ' refine:'
2807            iFin = min(iBeg+9,len(insKeys))
2808            for item in insKeys[iBeg:iFin]:
2809                if item not in ['Type','Source']:
2810                    ptlbls += '%12s' % (item)
2811                    ptstr += '%12.6f' % (Inst[item][1])
2812                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2813                        varstr += 12*' '
2814                    else:
2815                        varstr += '%12s' % (str(bool(Inst[item][2])))
2816            print >>pFile,ptlbls
2817            print >>pFile,ptstr
2818            print >>pFile,varstr
2819            iBeg = iFin
2820            if iBeg == len(insKeys):
2821                Ok = False
2822            else:
2823                print >>pFile,'\n'
2824       
2825    def PrintSampleParms(Sample):
2826        print >>pFile,'\n Sample Parameters:'
2827        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2828            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2829        ptlbls = ' name  :'
2830        ptstr =  ' value :'
2831        varstr = ' refine:'
2832        if 'Bragg' in Sample['Type']:
2833            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2834                ptlbls += '%14s'%(item)
2835                ptstr += '%14.4f'%(Sample[item][0])
2836                varstr += '%14s'%(str(bool(Sample[item][1])))
2837           
2838        elif 'Debye' in Type:        #Debye-Scherrer
2839            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2840                ptlbls += '%14s'%(item)
2841                ptstr += '%14.4f'%(Sample[item][0])
2842                varstr += '%14s'%(str(bool(Sample[item][1])))
2843
2844        print >>pFile,ptlbls
2845        print >>pFile,ptstr
2846        print >>pFile,varstr
2847       
2848    histDict = {}
2849    histVary = []
2850    controlDict = {}
2851    histoList = Histograms.keys()
2852    histoList.sort()
2853    for histogram in histoList:
2854        if 'PWDR' in histogram:
2855            Histogram = Histograms[histogram]
2856            hId = Histogram['hId']
2857            pfx = ':'+str(hId)+':'
2858            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2859            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2860            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2861            for excl in controlDict[pfx+'Exclude']:
2862                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2863            if controlDict[pfx+'Exclude']:
2864                ma.mask_rows(Histogram['Data'])
2865            Background = Histogram['Background']
2866            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2867            controlDict[pfx+'bakType'] = Type
2868            histDict.update(bakDict)
2869            histVary += bakVary
2870           
2871            Inst = Histogram['Instrument Parameters'][0]
2872            Type,instDict,insVary = GetInstParms(hId,Inst)
2873            controlDict[pfx+'histType'] = Type
2874            if 'XC' in Type:
2875                if pfx+'Lam1' in instDict:
2876                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2877                else:
2878                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2879            histDict.update(instDict)
2880            histVary += insVary
2881           
2882            Sample = Histogram['Sample Parameters']
2883            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2884            controlDict[pfx+'instType'] = Type
2885            histDict.update(sampDict)
2886            histVary += sampVary
2887           
2888   
2889            if Print: 
2890                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2891                print >>pFile,135*'-'
2892                Units = {'C':' deg','T':' msec'}
2893                units = Units[controlDict[pfx+'histType'][2]]
2894                Limits = controlDict[pfx+'Limits']
2895                print >>pFile,' Instrument type: ',Sample['Type']
2896                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2897                if len(controlDict[pfx+'Exclude']):
2898                    excls = controlDict[pfx+'Exclude']
2899                    for excl in excls:
2900                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
2901                PrintSampleParms(Sample)
2902                PrintInstParms(Inst)
2903                PrintBackground(Background)
2904        elif 'HKLF' in histogram:
2905            Histogram = Histograms[histogram]
2906            hId = Histogram['hId']
2907            pfx = ':'+str(hId)+':'
2908            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2909            Inst = Histogram['Instrument Parameters'][0]
2910            controlDict[pfx+'histType'] = Inst['Type'][0]
2911            if 'X' in Inst['Type'][0]:
2912                histDict[pfx+'Lam'] = Inst['Lam'][1]
2913                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2914    return histVary,histDict,controlDict
2915   
2916def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
2917    'needs a doc string'
2918   
2919    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2920        Back = Background[0]
2921        DebyePeaks = Background[1]
2922        lenBack = len(Back[3:])
2923        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2924        for i in range(lenBack):
2925            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
2926            if pfx+'Back;'+str(i) in sigDict:
2927                backSig[i] = sigDict[pfx+'Back;'+str(i)]
2928        if DebyePeaks['nDebye']:
2929            for i in range(DebyePeaks['nDebye']):
2930                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
2931                for j,name in enumerate(names):
2932                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2933                    if name in sigDict:
2934                        backSig[lenBack+3*i+j] = sigDict[name]           
2935        if DebyePeaks['nPeaks']:
2936            for i in range(DebyePeaks['nPeaks']):
2937                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2938                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2939                for j,name in enumerate(names):
2940                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2941                    if name in sigDict:
2942                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2943        return backSig
2944       
2945    def SetInstParms(pfx,Inst,parmDict,sigDict):
2946        instSig = {}
2947        insKeys = Inst.keys()
2948        insKeys.sort()
2949        for item in insKeys:
2950            insName = pfx+item
2951            Inst[item][1] = parmDict[insName]
2952            if insName in sigDict:
2953                instSig[item] = sigDict[insName]
2954            else:
2955                instSig[item] = 0
2956        return instSig
2957       
2958    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2959        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2960            sampSig = [0 for i in range(5)]
2961            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2962                Sample[item][0] = parmDict[pfx+item]
2963                if pfx+item in sigDict:
2964                    sampSig[i] = sigDict[pfx+item]
2965        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2966            sampSig = [0 for i in range(4)]
2967            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2968                Sample[item][0] = parmDict[pfx+item]
2969                if pfx+item in sigDict:
2970                    sampSig[i] = sigDict[pfx+item]
2971        return sampSig
2972       
2973    def PrintBackgroundSig(Background,backSig):
2974        Back = Background[0]
2975        DebyePeaks = Background[1]
2976        lenBack = len(Back[3:])
2977        valstr = ' value : '
2978        sigstr = ' sig   : '
2979        refine = False
2980        for i,back in enumerate(Back[3:]):
2981            valstr += '%10.4g'%(back)
2982            if Back[1]:
2983                refine = True
2984                sigstr += '%10.4g'%(backSig[i])
2985            else:
2986                sigstr += 10*' '
2987        if refine:
2988            print >>pFile,'\n Background function: ',Back[0]
2989            print >>pFile,valstr
2990            print >>pFile,sigstr
2991        if DebyePeaks['nDebye']:
2992            ifAny = False
2993            ptfmt = "%12.3f"
2994            names =  ' names :'
2995            ptstr =  ' values:'
2996            sigstr = ' esds  :'
2997            for item in sigDict:
2998                if 'Debye' in item:
2999                    ifAny = True
3000                    names += '%12s'%(item)
3001                    ptstr += ptfmt%(parmDict[item])
3002                    sigstr += ptfmt%(sigDict[item])
3003            if ifAny:
3004                print >>pFile,'\n Debye diffuse scattering coefficients'
3005                print >>pFile,names
3006                print >>pFile,ptstr
3007                print >>pFile,sigstr
3008        if DebyePeaks['nPeaks']:
3009            print >>pFile,'\n Single peak coefficients:'
3010            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3011            line = ' peak no. '
3012            for parm in parms:
3013                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3014            print >>pFile,line
3015            for ip in range(DebyePeaks['nPeaks']):
3016                ptstr = ' %4d '%(ip)
3017                for parm in parms:
3018                    fmt = '%14.3f'
3019                    efmt = '%12.3f'
3020                    if parm == 'BkPkpos':
3021                        fmt = '%14.4f'
3022                        efmt = '%12.4f'
3023                    name = pfx+parm+';%d'%(ip)
3024                    ptstr += fmt%(parmDict[name])
3025                    if name in sigDict:
3026                        ptstr += efmt%(sigDict[name])
3027                    else:
3028                        ptstr += 12*' '
3029                print >>pFile,ptstr
3030        sumBk = np.array(Histogram['sumBk'])
3031        print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g'    \
3032            %(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
3033       
3034    def PrintInstParmsSig(Inst,instSig):
3035        refine = False
3036        insKeys = instSig.keys()
3037        insKeys.sort()
3038        iBeg = 0
3039        Ok = True
3040        while Ok:
3041            ptlbls = ' names :'
3042            ptstr =  ' value :'
3043            sigstr = ' sig   :'
3044            iFin = min(iBeg+9,len(insKeys))
3045            for name in insKeys[iBeg:iFin]:
3046                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3047                    ptlbls += '%12s' % (name)
3048                    ptstr += '%12.6f' % (Inst[name][1])
3049                    if instSig[name]:
3050                        refine = True
3051                        sigstr += '%12.6f' % (instSig[name])
3052                    else:
3053                        sigstr += 12*' '
3054            if refine:
3055                print >>pFile,'\n Instrument Parameters:'
3056                print >>pFile,ptlbls
3057                print >>pFile,ptstr
3058                print >>pFile,sigstr
3059            iBeg = iFin
3060            if iBeg == len(insKeys):
3061                Ok = False
3062       
3063    def PrintSampleParmsSig(Sample,sampleSig):
3064        ptlbls = ' names :'
3065        ptstr =  ' values:'
3066        sigstr = ' sig   :'
3067        refine = False
3068        if 'Bragg' in Sample['Type']:
3069            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3070                ptlbls += '%14s'%(item)
3071                ptstr += '%14.4f'%(Sample[item][0])
3072                if sampleSig[i]:
3073                    refine = True
3074                    sigstr += '%14.4f'%(sampleSig[i])
3075                else:
3076                    sigstr += 14*' '
3077           
3078        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3079            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3080                ptlbls += '%14s'%(item)
3081                ptstr += '%14.4f'%(Sample[item][0])
3082                if sampleSig[i]:
3083                    refine = True
3084                    sigstr += '%14.4f'%(sampleSig[i])
3085                else:
3086                    sigstr += 14*' '
3087
3088        if refine:
3089            print >>pFile,'\n Sample Parameters:'
3090            print >>pFile,ptlbls
3091            print >>pFile,ptstr
3092            print >>pFile,sigstr
3093       
3094    histoList = Histograms.keys()
3095    histoList.sort()
3096    for histogram in histoList:
3097        if 'PWDR' in histogram:
3098            Histogram = Histograms[histogram]
3099            hId = Histogram['hId']
3100            pfx = ':'+str(hId)+':'
3101            Background = Histogram['Background']
3102            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3103           
3104            Inst = Histogram['Instrument Parameters'][0]
3105            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3106       
3107            Sample = Histogram['Sample Parameters']
3108            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3109
3110            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3111            print >>pFile,135*'-'
3112            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3113            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3114            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3115                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3116            if Print:
3117                print >>pFile,' Instrument type: ',Sample['Type']
3118                if FFtables != None and 'N' not in Inst['Type'][0]:
3119                    PrintFprime(FFtables,pfx,pFile)
3120                PrintSampleParmsSig(Sample,sampSig)
3121                PrintInstParmsSig(Inst,instSig)
3122                PrintBackgroundSig(Background,backSig)
3123               
Note: See TracBrowser for help on using the repository browser.