source: trunk/GSASIIstrIO.py @ 1635

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

modulation functions & constraints

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