source: trunk/GSASIIstrIO.py @ 2064

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

more fixed to SS modulation display & Refine listing

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