source: trunk/GSASIIstrIO.py @ 2110

Last change on this file since 2110 was 2110, checked in by vondreele, 6 years ago

add menu item for global setting of wave vary parameters - TBD
split StructureFactorDeriv? over twins & incommensurate
do block refl for the nontwin/powder one

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