source: trunk/GSASIIstrIO.py @ 1614

Last change on this file since 1614 was 1614, checked in by vondreele, 7 years ago

Fix PDF plots & PDF Controls window issue
SS names now fixed
LS input & SS parameters

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 122.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2014-12-19 18:51:48 +0000 (Fri, 19 Dec 2014) $
4# $Author: vondreele $
5# $Revision: 1614 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 1614 2014-12-19 18:51:48Z 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: 1614 $")
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 = 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'] == 'magnetic':
831            line += '   Mx     My     Mz'
832        elif General['Type'] == 'macromolecular':
833            line = ' res no residue chain'+line
834        print >>pFile,line
835        if General['Type'] == 'nuclear':
836            print >>pFile,135*'-'
837            for i,at in enumerate(Atoms):
838                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
839                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
840                if at[cia] == 'I':
841                    line += '%8.5f'%(at[cia+1])+48*' '
842                else:
843                    line += 8*' '
844                    for j in range(6):
845                        line += '%8.5f'%(at[cia+2+j])
846                print >>pFile,line
847        elif General['Type'] == 'macromolecular':
848            print >>pFile,135*'-'           
849            for i,at in enumerate(Atoms):
850                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])+ \
851                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
852                if at[cia] == 'I':
853                    line += '%8.4f'%(at[cia+1])+48*' '
854                else:
855                    line += 8*' '
856                    for j in range(6):
857                        line += '%8.4f'%(at[cia+2+j])
858                print >>pFile,line
859       
860    def PrintTexture(textureData):
861        topstr = '\n Spherical harmonics texture: Order:' + \
862            str(textureData['Order'])
863        if textureData['Order']:
864            print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
865        else:
866            print >>pFile,topstr
867            return
868        names = ['omega','chi','phi']
869        line = '\n'
870        for name in names:
871            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
872        print >>pFile,line
873        print >>pFile,'\n Texture coefficients:'
874        ptlbls = ' names :'
875        ptstr =  ' values:'
876        SHcoeff = textureData['SH Coeff'][1]
877        for item in SHcoeff:
878            ptlbls += '%12s'%(item)
879            ptstr += '%12.4f'%(SHcoeff[item]) 
880        print >>pFile,ptlbls
881        print >>pFile,ptstr
882       
883    def MakeRBParms(rbKey,phaseVary,phaseDict):
884        rbid = str(rbids.index(RB['RBId']))
885        pfxRB = pfx+'RB'+rbKey+'P'
886        pstr = ['x','y','z']
887        ostr = ['a','i','j','k']
888        for i in range(3):
889            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
890            phaseDict[name] = RB['Orig'][0][i]
891            if RB['Orig'][1]:
892                phaseVary += [name,]
893        pfxRB = pfx+'RB'+rbKey+'O'
894        for i in range(4):
895            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
896            phaseDict[name] = RB['Orient'][0][i]
897            if RB['Orient'][1] == 'AV' and i:
898                phaseVary += [name,]
899            elif RB['Orient'][1] == 'A' and not i:
900                phaseVary += [name,]
901           
902    def MakeRBThermals(rbKey,phaseVary,phaseDict):
903        rbid = str(rbids.index(RB['RBId']))
904        tlstr = ['11','22','33','12','13','23']
905        sstr = ['12','13','21','23','31','32','AA','BB']
906        if 'T' in RB['ThermalMotion'][0]:
907            pfxRB = pfx+'RB'+rbKey+'T'
908            for i in range(6):
909                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
910                phaseDict[name] = RB['ThermalMotion'][1][i]
911                if RB['ThermalMotion'][2][i]:
912                    phaseVary += [name,]
913        if 'L' in RB['ThermalMotion'][0]:
914            pfxRB = pfx+'RB'+rbKey+'L'
915            for i in range(6):
916                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
917                phaseDict[name] = RB['ThermalMotion'][1][i+6]
918                if RB['ThermalMotion'][2][i+6]:
919                    phaseVary += [name,]
920        if 'S' in RB['ThermalMotion'][0]:
921            pfxRB = pfx+'RB'+rbKey+'S'
922            for i in range(8):
923                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
924                phaseDict[name] = RB['ThermalMotion'][1][i+12]
925                if RB['ThermalMotion'][2][i+12]:
926                    phaseVary += [name,]
927        if 'U' in RB['ThermalMotion'][0]:
928            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
929            phaseDict[name] = RB['ThermalMotion'][1][0]
930            if RB['ThermalMotion'][2][0]:
931                phaseVary += [name,]
932               
933    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
934        rbid = str(rbids.index(RB['RBId']))
935        pfxRB = pfx+'RB'+rbKey+'Tr;'
936        for i,torsion in enumerate(RB['Torsions']):
937            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
938            phaseDict[name] = torsion[0]
939            if torsion[1]:
940                phaseVary += [name,]
941                   
942    if Print:
943        print  >>pFile,'\n Phases:'
944    phaseVary = []
945    phaseDict = {}
946    phaseConstr = {}
947    pawleyLookup = {}
948    FFtables = {}                   #scattering factors - xrays
949    BLtables = {}                   # neutrons
950    Natoms = {}
951    AtMults = {}
952    AtIA = {}
953    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
954    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
955    atomIndx = {}
956    for name in PhaseData:
957        General = PhaseData[name]['General']
958        pId = PhaseData[name]['pId']
959        pfx = str(pId)+'::'
960        FFtable = G2el.GetFFtable(General['AtomTypes'])
961        BLtable = G2el.GetBLtable(General)
962        FFtables.update(FFtable)
963        BLtables.update(BLtable)
964        Atoms = PhaseData[name]['Atoms']
965        if Atoms and not General.get('doPawley'):
966            cx,ct,cs,cia = General['AtomPtrs']
967            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
968        PawleyRef = PhaseData[name].get('Pawley ref',[])
969        SGData = General['SGData']
970        SGtext,SGtable = G2spc.SGPrint(SGData)
971        cell = General['Cell']
972        A = G2lat.cell2A(cell[1:7])
973        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
974            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
975        if cell[0]:
976            phaseVary += cellVary(pfx,SGData)
977        SSGtext = []    #no superstructure
978        im = 0
979        if General['Type'] in ['modulated','magnetic']:
980            im = 1      #refl offset
981            Vec,vRef,maxH = General['SuperVec']
982            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
983            SSGData = General['SSGData']
984            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
985            if vRef:
986                phaseVary += modVary(pfx,SSGData)       
987        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
988        if resRBData:
989            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
990            for iRB,RB in enumerate(resRBData):
991                MakeRBParms('R',phaseVary,phaseDict)
992                MakeRBThermals('R',phaseVary,phaseDict)
993                MakeRBTorsions('R',phaseVary,phaseDict)
994       
995        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
996        if vecRBData:
997            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
998            for iRB,RB in enumerate(vecRBData):
999                MakeRBParms('V',phaseVary,phaseDict)
1000                MakeRBThermals('V',phaseVary,phaseDict)
1001                   
1002        Natoms[pfx] = 0
1003        if Atoms and not General.get('doPawley'):
1004            cx,ct,cs,cia = General['AtomPtrs']
1005            Natoms[pfx] = len(Atoms)
1006            for i,at in enumerate(Atoms):
1007                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1008                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1009                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1010                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1011                    pfx+'AI/A:'+str(i):at[cia],})
1012                if at[cia] == 'I':
1013                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1014                else:
1015                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1016                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1017                if 'F' in at[ct+1]:
1018                    phaseVary.append(pfx+'Afrac:'+str(i))
1019                if 'X' in at[ct+1]:
1020                    try:    #patch for sytsym name changes
1021                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1022                    except KeyError:
1023                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1024                        at[cs] = Sytsym
1025                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1026                    xId,xCoef = G2spc.GetCSxinel(at[cs])
1027                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1028                    equivs = [[],[],[]]
1029                    for j in range(3):
1030                        if xId[j] > 0:                               
1031                            phaseVary.append(names[j])
1032                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1033                    for equiv in equivs:
1034                        if len(equiv) > 1:
1035                            name = equiv[0][0]
1036                            coef = equiv[0][1]
1037                            for eqv in equiv[1:]:
1038                                eqv[1] /= coef
1039                                G2mv.StoreEquivalence(name,(eqv,))
1040                if 'U' in at[ct+1]:
1041                    if at[cia] == 'I':
1042                        phaseVary.append(pfx+'AUiso:'+str(i))
1043                    else:
1044                        try:    #patch for sytsym name changes
1045                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1046                        except KeyError:
1047                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1048                            at[cs] = Sytsym
1049                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1050                        uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1051                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1052                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1053                        equivs = [[],[],[],[],[],[]]
1054                        for j in range(6):
1055                            if uId[j] > 0:                               
1056                                phaseVary.append(names[j])
1057                                equivs[uId[j]-1].append([names[j],uCoef[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,equiv[1:])
1065                if General['Type'] in ['modulated','magnetic']:
1066                    AtomSS = at[-1]['SS1']
1067                    CSI = G2spc.GetSSfxuinel(at[cx:cx+3],at[cia+2:cia+8],SGData,SSGData)
1068                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1069                        Waves = AtomSS[Stype]
1070                        uId,uCoef = CSI[Stype]
1071                        for iw,wave in enumerate(Waves):
1072                            stiw = str(i)+':'+str(iw)
1073                            if Stype == 'Spos':
1074                                names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1075                                    pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1076                                equivs = [[],[],[], [],[],[]]
1077                            elif Stype == 'Sadp':
1078                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1079                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1080                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1081                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1082                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1083                            elif Stype == 'Sfrac':
1084                                equivs = [[],[]]
1085                                names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1086                            elif Stype == 'Smag':
1087                                equivs = [[],[],[], [],[],[]]
1088                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1089                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1090                            phaseDict.update(dict(zip(names,wave[0])))
1091                            if wave[1]:
1092                                for j in range(len(equivs)):
1093                                    if uId[j] > 0:                               
1094                                        phaseVary.append(names[j])
1095                                        equivs[uId[j]-1].append([names[j],uCoef[j]])
1096                                for equiv in equivs:
1097                                    if len(equiv) > 1:
1098                                        name = equiv[0][0]
1099                                        coef = equiv[0][1]
1100                                        for eqv in equiv[1:]:
1101                                            eqv[1] /= coef
1102                                        G2mv.StoreEquivalence(name,equiv[1:])
1103                                           
1104            textureData = General['SH Texture']
1105            if textureData['Order']:
1106                phaseDict[pfx+'SHorder'] = textureData['Order']
1107                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1108                for item in ['omega','chi','phi']:
1109                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1110                    if textureData['Sample '+item][0]:
1111                        phaseVary.append(pfx+'SH '+item)
1112                for item in textureData['SH Coeff'][1]:
1113                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1114                    if textureData['SH Coeff'][0]:
1115                        phaseVary.append(pfx+item)
1116               
1117            if Print:
1118                print >>pFile,'\n Phase name: ',General['Name']
1119                print >>pFile,135*'-'
1120                PrintFFtable(FFtable)
1121                PrintBLtable(BLtable)
1122                print >>pFile,''
1123                if len(SSGtext):    #if superstructure
1124                    for line in SSGtext: print >>pFile,line
1125                    if len(SSGtable):
1126                        for item in SSGtable:
1127                            line = ' %s '%(item)
1128                            print >>pFile,line   
1129                    else:
1130                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1131                else:
1132                    for line in SGtext: print >>pFile,line
1133                    if len(SGtable):
1134                        for item in SGtable:
1135                            line = ' %s '%(item)
1136                            print >>pFile,line   
1137                    else:
1138                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1139                PrintRBObjects(resRBData,vecRBData)
1140                PrintAtoms(General,Atoms)
1141                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1142                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1143                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1144                if len(SSGtext):    #if superstructure
1145                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1146                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1147                PrintTexture(textureData)
1148                if name in RestraintDict:
1149                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1150                        textureData,RestraintDict[name],pFile)
1151                   
1152        elif PawleyRef:
1153            if Print:
1154                print >>pFile,'\n Phase name: ',General['Name']
1155                print >>pFile,135*'-'
1156                print >>pFile,''
1157                if len(SSGtext):    #if superstructure
1158                    for line in SSGtext: print >>pFile,line
1159                    if len(SSGtable):
1160                        for item in SSGtable:
1161                            line = ' %s '%(item)
1162                            print >>pFile,line   
1163                    else:
1164                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1165                else:
1166                    for line in SGtext: print >>pFile,line
1167                    if len(SGtable):
1168                        for item in SGtable:
1169                            line = ' %s '%(item)
1170                            print >>pFile,line   
1171                    else:
1172                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1173                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1174                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1175                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1176                if len(SSGtext):    #if superstructure
1177                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1178                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1179            pawleyVary = []
1180            for i,refl in enumerate(PawleyRef):
1181                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1182                if im:
1183                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1184                else:
1185                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1186                if refl[5+im]:
1187                    pawleyVary.append(pfx+'PWLref:'+str(i))
1188            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1189            phaseVary += pawleyVary
1190               
1191    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1192   
1193def cellFill(pfx,SGData,parmDict,sigDict): 
1194    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1195    from the parameter and sig dictionaries.
1196
1197    :param str pfx: parameter prefix ("n::", where n is a phase number)
1198    :param dict SGdata: a symmetry object
1199    :param dict parmDict: a dictionary of parameters
1200    :param dict sigDict:  a dictionary of uncertainties on parameters
1201
1202    :returns: A,sigA where each is a list of six terms with the A terms
1203    '''
1204    if SGData['SGLaue'] in ['-1',]:
1205        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1206            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1207    elif SGData['SGLaue'] in ['2/m',]:
1208        if SGData['SGUniq'] == 'a':
1209            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1210                parmDict[pfx+'A3'],0,0]
1211        elif SGData['SGUniq'] == 'b':
1212            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1213                0,parmDict[pfx+'A4'],0]
1214        else:
1215            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1216                0,0,parmDict[pfx+'A5']]
1217    elif SGData['SGLaue'] in ['mmm',]:
1218        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1219    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1220        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1221    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1222        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1223            parmDict[pfx+'A0'],0,0]
1224    elif SGData['SGLaue'] in ['3R', '3mR']:
1225        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1226            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1227    elif SGData['SGLaue'] in ['m3m','m3']:
1228        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1229
1230    try:
1231        if SGData['SGLaue'] in ['-1',]:
1232            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1233                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1234        elif SGData['SGLaue'] in ['2/m',]:
1235            if SGData['SGUniq'] == 'a':
1236                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1237                    sigDict[pfx+'A3'],0,0]
1238            elif SGData['SGUniq'] == 'b':
1239                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1240                    0,sigDict[pfx+'A4'],0]
1241            else:
1242                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1243                    0,0,sigDict[pfx+'A5']]
1244        elif SGData['SGLaue'] in ['mmm',]:
1245            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1246        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1247            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1248        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1249            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1250        elif SGData['SGLaue'] in ['3R', '3mR']:
1251            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1252        elif SGData['SGLaue'] in ['m3m','m3']:
1253            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1254    except KeyError:
1255        sigA = [0,0,0,0,0,0]
1256
1257    return A,sigA
1258       
1259def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1260    'needs a doc string'
1261    if phaseRest:
1262        Amat = G2lat.cell2AB(cell)[0]
1263        cx,ct,cs = AtPtrs[:3]
1264        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1265            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1266            ['ChemComp','Sites'],['Texture','HKLs']]
1267        for name,rest in names:
1268            itemRest = phaseRest[name]
1269            if itemRest[rest] and itemRest['Use']:
1270                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1271                if name in ['Bond','Angle','Plane','Chiral']:
1272                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1273                    for indx,ops,obs,esd in itemRest[rest]:
1274                        try:
1275                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1276                            AtName = ''
1277                            for i,Aname in enumerate(AtNames):
1278                                AtName += Aname
1279                                if ops[i] == '1':
1280                                    AtName += '-'
1281                                else:
1282                                    AtName += '+('+ops[i]+')-'
1283                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1284                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1285                            if name == 'Bond':
1286                                calc = G2mth.getRestDist(XYZ,Amat)
1287                            elif name == 'Angle':
1288                                calc = G2mth.getRestAngle(XYZ,Amat)
1289                            elif name == 'Plane':
1290                                calc = G2mth.getRestPlane(XYZ,Amat)
1291                            elif name == 'Chiral':
1292                                calc = G2mth.getRestChiral(XYZ,Amat)
1293                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1294                        except KeyError:
1295                            del itemRest[rest]
1296                elif name in ['Torsion','Rama']:
1297                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1298                    coeffDict = itemRest['Coeff']
1299                    for indx,ops,cofName,esd in itemRest[rest]:
1300                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1301                        AtName = ''
1302                        for i,Aname in enumerate(AtNames):
1303                            AtName += Aname+'+('+ops[i]+')-'
1304                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1305                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1306                        if name == 'Torsion':
1307                            tor = G2mth.getRestTorsion(XYZ,Amat)
1308                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1309                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1310                        else:
1311                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1312                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1313                            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])
1314                elif name == 'ChemComp':
1315                    print >>pFile,'     atoms   mul*frac  factor     prod'
1316                    for indx,factors,obs,esd in itemRest[rest]:
1317                        try:
1318                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1319                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1320                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1321                            mulfrac = mul*frac
1322                            calcs = mul*frac*factors
1323                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1324                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1325                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1326                        except KeyError:
1327                            del itemRest[rest]
1328                elif name == 'Texture' and textureData['Order']:
1329                    Start = False
1330                    SHCoef = textureData['SH Coeff'][1]
1331                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1332                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1333                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1334                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1335                        PH = np.array(hkl)
1336                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1337                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1338                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1339                        Z = ma.masked_greater(Z,0.0)
1340                        num = ma.count(Z)
1341                        sum = 0
1342                        if num:
1343                            sum = np.sum(Z)
1344                        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)
1345       
1346def getCellEsd(pfx,SGData,A,covData):
1347    'needs a doc string'
1348    dpr = 180./np.pi
1349    rVsq = G2lat.calc_rVsq(A)
1350    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1351    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1352    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1353    scos = cosd(cellst[3:6])
1354    ssin = sind(cellst[3:6])
1355    scot = scos/ssin
1356    rcos = cosd(cell[3:6])
1357    rsin = sind(cell[3:6])
1358    rcot = rcos/rsin
1359    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1360    varyList = covData['varyList']
1361    covMatrix = covData['covMatrix']
1362    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1363    Ax = np.array(A)
1364    Ax[3:] /= 2.
1365    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
1366        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1367    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1368    Vol = 1/np.sqrt(rVsq)
1369    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
1370    R123 = Ax[0]*Ax[1]*Ax[2]
1371    dsasdg = np.zeros((3,6))
1372    dadg = np.zeros((6,6))
1373    for i0 in range(3):         #0  1   2
1374        i1 = (i0+1)%3           #1  2   0
1375        i2 = (i1+1)%3           #2  0   1
1376        i3 = 5-i2               #3  5   4
1377        i4 = 5-i1               #4  3   5
1378        i5 = 5-i0               #5  4   3
1379        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1380        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1381        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1382        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1383        denom = np.sqrt(denmsq)
1384        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1385        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1386        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1387        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1388        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1389        dadg[i5][i5] = -Ax[i0]/denom
1390    for i0 in range(3):
1391        i1 = (i0+1)%3
1392        i2 = (i1+1)%3
1393        i3 = 5-i2
1394        for ij in range(6):
1395            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1396            if ij == i0:
1397                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1398            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1399    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1400    var = np.diag(sigMat)
1401    CS = np.where(var>0.,np.sqrt(var),0.)
1402    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
1403   
1404def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1405    'needs a doc string'
1406   
1407    def PrintAtomsAndSig(General,Atoms,atomsSig):
1408        print >>pFile,'\n Atoms:'
1409        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1410        if General['Type'] == 'magnetic':
1411            line += '   Mx     My     Mz'
1412        elif General['Type'] == 'macromolecular':
1413            line = ' res no  residue  chain '+line
1414        cx,ct,cs,cia = General['AtomPtrs']
1415        print >>pFile,line
1416        if General['Type'] == 'nuclear':
1417            print >>pFile,135*'-'
1418            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1419                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1420            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
1421            for atyp in General['AtomTypes']:       #zero composition
1422                General['NoAtoms'][atyp] = 0.
1423            for i,at in enumerate(Atoms):
1424                General['NoAtoms'][at[1]] += at[6]*float(at[8])     #new composition
1425                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1426                valstr = ' values:'
1427                sigstr = ' sig   :'
1428                for ind in [3,4,5,6]:
1429                    sigind = str(i)+':'+str(ind)
1430                    valstr += fmt[ind]%(at[ind])                   
1431                    if sigind in atomsSig:
1432                        sigstr += fmt[ind]%(atomsSig[sigind])
1433                    else:
1434                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1435                if at[9] == 'I':
1436                    valstr += fmt[10]%(at[10])
1437                    if str(i)+':10' in atomsSig:
1438                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1439                    else:
1440                        sigstr += 8*' '
1441                else:
1442                    valstr += 8*' '
1443                    sigstr += 8*' '
1444                    for ind in [11,12,13,14,15,16]:
1445                        sigind = str(i)+':'+str(ind)
1446                        valstr += fmt[ind]%(at[ind])
1447                        if sigind in atomsSig:                       
1448                            sigstr += fmt[ind]%(atomsSig[sigind])
1449                        else:
1450                            sigstr += 8*' '
1451                print >>pFile,name
1452                print >>pFile,valstr
1453                print >>pFile,sigstr
1454               
1455    def PrintRBObjPOAndSig(rbfx,rbsx):
1456        namstr = '  names :'
1457        valstr = '  values:'
1458        sigstr = '  esds  :'
1459        for i,px in enumerate(['Px:','Py:','Pz:']):
1460            name = pfx+rbfx+px+rbsx
1461            namstr += '%12s'%('Pos '+px[1])
1462            valstr += '%12.5f'%(parmDict[name])
1463            if name in sigDict:
1464                sigstr += '%12.5f'%(sigDict[name])
1465            else:
1466                sigstr += 12*' '
1467        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1468            name = pfx+rbfx+po+rbsx
1469            namstr += '%12s'%('Ori '+po[1])
1470            valstr += '%12.5f'%(parmDict[name])
1471            if name in sigDict:
1472                sigstr += '%12.5f'%(sigDict[name])
1473            else:
1474                sigstr += 12*' '
1475        print >>pFile,namstr
1476        print >>pFile,valstr
1477        print >>pFile,sigstr
1478       
1479    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1480        namstr = '  names :'
1481        valstr = '  values:'
1482        sigstr = '  esds  :'
1483        if 'N' not in TLS:
1484            print >>pFile,' Thermal motion:'
1485        if 'T' in TLS:
1486            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1487                name = pfx+rbfx+pt+rbsx
1488                namstr += '%12s'%(pt[:3])
1489                valstr += '%12.5f'%(parmDict[name])
1490                if name in sigDict:
1491                    sigstr += '%12.5f'%(sigDict[name])
1492                else:
1493                    sigstr += 12*' '
1494            print >>pFile,namstr
1495            print >>pFile,valstr
1496            print >>pFile,sigstr
1497        if 'L' in TLS:
1498            namstr = '  names :'
1499            valstr = '  values:'
1500            sigstr = '  esds  :'
1501            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1502                name = pfx+rbfx+pt+rbsx
1503                namstr += '%12s'%(pt[:3])
1504                valstr += '%12.3f'%(parmDict[name])
1505                if name in sigDict:
1506                    sigstr += '%12.3f'%(sigDict[name])
1507                else:
1508                    sigstr += 12*' '
1509            print >>pFile,namstr
1510            print >>pFile,valstr
1511            print >>pFile,sigstr
1512        if 'S' in TLS:
1513            namstr = '  names :'
1514            valstr = '  values:'
1515            sigstr = '  esds  :'
1516            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1517                name = pfx+rbfx+pt+rbsx
1518                namstr += '%12s'%(pt[:3])
1519                valstr += '%12.4f'%(parmDict[name])
1520                if name in sigDict:
1521                    sigstr += '%12.4f'%(sigDict[name])
1522                else:
1523                    sigstr += 12*' '
1524            print >>pFile,namstr
1525            print >>pFile,valstr
1526            print >>pFile,sigstr
1527        if 'U' in TLS:
1528            name = pfx+rbfx+'U:'+rbsx
1529            namstr = '  names :'+'%12s'%('Uiso')
1530            valstr = '  values:'+'%12.5f'%(parmDict[name])
1531            if name in sigDict:
1532                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1533            else:
1534                sigstr = '  esds  :'+12*' '
1535            print >>pFile,namstr
1536            print >>pFile,valstr
1537            print >>pFile,sigstr
1538       
1539    def PrintRBObjTorAndSig(rbsx):
1540        namstr = '  names :'
1541        valstr = '  values:'
1542        sigstr = '  esds  :'
1543        nTors = len(RBObj['Torsions'])
1544        if nTors:
1545            print >>pFile,' Torsions:'
1546            for it in range(nTors):
1547                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1548                namstr += '%12s'%('Tor'+str(it))
1549                valstr += '%12.4f'%(parmDict[name])
1550                if name in sigDict:
1551                    sigstr += '%12.4f'%(sigDict[name])
1552            print >>pFile,namstr
1553            print >>pFile,valstr
1554            print >>pFile,sigstr
1555               
1556    def PrintSHtextureAndSig(textureData,SHtextureSig):
1557        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1558        names = ['omega','chi','phi']
1559        namstr = '  names :'
1560        ptstr =  '  values:'
1561        sigstr = '  esds  :'
1562        for name in names:
1563            namstr += '%12s'%(name)
1564            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1565            if 'Sample '+name in SHtextureSig:
1566                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1567            else:
1568                sigstr += 12*' '
1569        print >>pFile,namstr
1570        print >>pFile,ptstr
1571        print >>pFile,sigstr
1572        print >>pFile,'\n Texture coefficients:'
1573        namstr = '  names :'
1574        ptstr =  '  values:'
1575        sigstr = '  esds  :'
1576        SHcoeff = textureData['SH Coeff'][1]
1577        for name in SHcoeff:
1578            namstr += '%12s'%(name)
1579            ptstr += '%12.3f'%(SHcoeff[name])
1580            if name in SHtextureSig:
1581                sigstr += '%12.3f'%(SHtextureSig[name])
1582            else:
1583                sigstr += 12*' '
1584        print >>pFile,namstr
1585        print >>pFile,ptstr
1586        print >>pFile,sigstr
1587           
1588    print >>pFile,'\n Phases:'
1589    for phase in Phases:
1590        print >>pFile,' Result for phase: ',phase
1591        Phase = Phases[phase]
1592        General = Phase['General']
1593        SGData = General['SGData']
1594        Atoms = Phase['Atoms']
1595        if Atoms and not General.get('doPawley'):
1596            cx,ct,cs,cia = General['AtomPtrs']
1597            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1598        cell = General['Cell']
1599        pId = Phase['pId']
1600        pfx = str(pId)+'::'
1601        if cell[0]:
1602            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1603            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1604            print >>pFile,' Reciprocal metric tensor: '
1605            ptfmt = "%15.9f"
1606            names = ['A11','A22','A33','A12','A13','A23']
1607            namstr = '  names :'
1608            ptstr =  '  values:'
1609            sigstr = '  esds  :'
1610            for name,a,siga in zip(names,A,sigA):
1611                namstr += '%15s'%(name)
1612                ptstr += ptfmt%(a)
1613                if siga:
1614                    sigstr += ptfmt%(siga)
1615                else:
1616                    sigstr += 15*' '
1617            print >>pFile,namstr
1618            print >>pFile,ptstr
1619            print >>pFile,sigstr
1620            cell[1:7] = G2lat.A2cell(A)
1621            cell[7] = G2lat.calc_V(A)
1622            print >>pFile,' New unit cell:'
1623            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1624            names = ['a','b','c','alpha','beta','gamma','Volume']
1625            namstr = '  names :'
1626            ptstr =  '  values:'
1627            sigstr = '  esds  :'
1628            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1629                namstr += '%12s'%(name)
1630                ptstr += fmt%(a)
1631                if siga:
1632                    sigstr += fmt%(siga)
1633                else:
1634                    sigstr += 12*' '
1635            print >>pFile,namstr
1636            print >>pFile,ptstr
1637            print >>pFile,sigstr
1638        ik = 6  #for Pawley stuff below
1639        if General['Type'] in ['modulated','magnetic']:
1640            ik = 7
1641            Vec,vRef,maxH = General['SuperVec']
1642            if vRef:
1643                print >>pFile,' New modulation vector:'
1644                namstr = '  names :'
1645                ptstr =  '  values:'
1646                sigstr = '  esds  :'
1647                for var in ['mV0','mV1','mV2']:
1648                    namstr += '%12s'%(pfx+var)
1649                    ptstr += '%12.6f'%(parmDict[pfx+var])
1650                    if pfx+var in sigDict:
1651                        sigstr += '%12.6f'%(sigDict[pfx+var])
1652                    else:
1653                        sigstr += 12*' '
1654                print >>pFile,namstr
1655                print >>pFile,ptstr
1656                print >>pFile,sigstr
1657           
1658        General['Mass'] = 0.
1659        if Phase['General'].get('doPawley'):
1660            pawleyRef = Phase['Pawley ref']
1661            for i,refl in enumerate(pawleyRef):
1662                key = pfx+'PWLref:'+str(i)
1663                refl[ik] = parmDict[key]
1664                if key in sigDict:
1665                    refl[ik+1] = sigDict[key]
1666                else:
1667                    refl[ik+1] = 0
1668        else:
1669            VRBIds = RBIds['Vector']
1670            RRBIds = RBIds['Residue']
1671            RBModels = Phase['RBModels']
1672            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1673                jrb = VRBIds.index(RBObj['RBId'])
1674                rbsx = str(irb)+':'+str(jrb)
1675                print >>pFile,' Vector rigid body parameters:'
1676                PrintRBObjPOAndSig('RBV',rbsx)
1677                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1678            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1679                jrb = RRBIds.index(RBObj['RBId'])
1680                rbsx = str(irb)+':'+str(jrb)
1681                print >>pFile,' Residue rigid body parameters:'
1682                PrintRBObjPOAndSig('RBR',rbsx)
1683                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1684                PrintRBObjTorAndSig(rbsx)
1685            atomsSig = {}
1686            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
1687                for i,at in enumerate(Atoms):
1688                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1689                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1690                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1691                    for ind in [3,4,5,6]:
1692                        at[ind] = parmDict[names[ind]]
1693                        if ind in [3,4,5]:
1694                            name = names[ind].replace('A','dA')
1695                        else:
1696                            name = names[ind]
1697                        if name in sigDict:
1698                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1699                    if at[9] == 'I':
1700                        at[10] = parmDict[names[10]]
1701                        if names[10] in sigDict:
1702                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1703                    else:
1704                        for ind in [11,12,13,14,15,16]:
1705                            at[ind] = parmDict[names[ind]]
1706                            if names[ind] in sigDict:
1707                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1708                    ind = General['AtomTypes'].index(at[1])
1709                    General['Mass'] += General['AtomMass'][ind]*at[6]*at[8]
1710            PrintAtomsAndSig(General,Atoms,atomsSig)
1711       
1712        textureData = General['SH Texture']   
1713        if textureData['Order']:
1714            SHtextureSig = {}
1715            for name in ['omega','chi','phi']:
1716                aname = pfx+'SH '+name
1717                textureData['Sample '+name][1] = parmDict[aname]
1718                if aname in sigDict:
1719                    SHtextureSig['Sample '+name] = sigDict[aname]
1720            for name in textureData['SH Coeff'][1]:
1721                aname = pfx+name
1722                textureData['SH Coeff'][1][name] = parmDict[aname]
1723                if aname in sigDict:
1724                    SHtextureSig[name] = sigDict[aname]
1725            PrintSHtextureAndSig(textureData,SHtextureSig)
1726        if phase in RestraintDict:
1727            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1728                textureData,RestraintDict[phase],pFile)
1729                   
1730################################################################################
1731##### Histogram & Phase data
1732################################################################################       
1733                   
1734def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1735    '''Loads the HAP histogram/phase information into dicts
1736
1737    :param dict Phases: phase information
1738    :param dict Histograms: Histogram information
1739    :param bool Print: prints information as it is read
1740    :param file pFile: file object to print to (the default, None causes printing to the console)
1741    :param bool resetRefList: Should the contents of the Reflection List be initialized
1742      on loading. The default, True, initializes the Reflection List as it is loaded.
1743
1744    :returns: (hapVary,hapDict,controlDict)
1745      * hapVary: list of refined variables
1746      * hapDict: dict with refined variables and their values
1747      * controlDict: dict with computation controls (?)
1748    '''
1749   
1750    def PrintSize(hapData):
1751        if hapData[0] in ['isotropic','uniaxial']:
1752            line = '\n Size model    : %9s'%(hapData[0])
1753            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1754            if hapData[0] == 'uniaxial':
1755                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1756            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1757            print >>pFile,line
1758        else:
1759            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
1760                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1761            Snames = ['S11','S22','S33','S12','S13','S23']
1762            ptlbls = ' names :'
1763            ptstr =  ' values:'
1764            varstr = ' refine:'
1765            for i,name in enumerate(Snames):
1766                ptlbls += '%12s' % (name)
1767                ptstr += '%12.6f' % (hapData[4][i])
1768                varstr += '%12s' % (str(hapData[5][i]))
1769            print >>pFile,ptlbls
1770            print >>pFile,ptstr
1771            print >>pFile,varstr
1772       
1773    def PrintMuStrain(hapData,SGData):
1774        if hapData[0] in ['isotropic','uniaxial']:
1775            line = '\n Mustrain model: %9s'%(hapData[0])
1776            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1777            if hapData[0] == 'uniaxial':
1778                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1779            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1780            print >>pFile,line
1781        else:
1782            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1783                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1784            Snames = G2spc.MustrainNames(SGData)
1785            ptlbls = ' names :'
1786            ptstr =  ' values:'
1787            varstr = ' refine:'
1788            for i,name in enumerate(Snames):
1789                ptlbls += '%12s' % (name)
1790                ptstr += '%12.6f' % (hapData[4][i])
1791                varstr += '%12s' % (str(hapData[5][i]))
1792            print >>pFile,ptlbls
1793            print >>pFile,ptstr
1794            print >>pFile,varstr
1795
1796    def PrintHStrain(hapData,SGData):
1797        print >>pFile,'\n Hydrostatic/elastic strain: '
1798        Hsnames = G2spc.HStrainNames(SGData)
1799        ptlbls = ' names :'
1800        ptstr =  ' values:'
1801        varstr = ' refine:'
1802        for i,name in enumerate(Hsnames):
1803            ptlbls += '%12s' % (name)
1804            ptstr += '%12.4g' % (hapData[0][i])
1805            varstr += '%12s' % (str(hapData[1][i]))
1806        print >>pFile,ptlbls
1807        print >>pFile,ptstr
1808        print >>pFile,varstr
1809
1810    def PrintSHPO(hapData):
1811        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1812            str(hapData[4])+' Refine? '+str(hapData[2])
1813        ptlbls = ' names :'
1814        ptstr =  ' values:'
1815        for item in hapData[5]:
1816            ptlbls += '%12s'%(item)
1817            ptstr += '%12.3f'%(hapData[5][item]) 
1818        print >>pFile,ptlbls
1819        print >>pFile,ptstr
1820   
1821    def PrintBabinet(hapData):
1822        print >>pFile,'\n Babinet form factor modification: '
1823        ptlbls = ' names :'
1824        ptstr =  ' values:'
1825        varstr = ' refine:'
1826        for name in ['BabA','BabU']:
1827            ptlbls += '%12s' % (name)
1828            ptstr += '%12.6f' % (hapData[name][0])
1829            varstr += '%12s' % (str(hapData[name][1]))
1830        print >>pFile,ptlbls
1831        print >>pFile,ptstr
1832        print >>pFile,varstr
1833       
1834   
1835    hapDict = {}
1836    hapVary = []
1837    controlDict = {}
1838    poType = {}
1839    poAxes = {}
1840    spAxes = {}
1841    spType = {}
1842   
1843    for phase in Phases:
1844        HistoPhase = Phases[phase]['Histograms']
1845        SGData = Phases[phase]['General']['SGData']
1846        cell = Phases[phase]['General']['Cell'][1:7]
1847        A = G2lat.cell2A(cell)
1848        if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
1849            SSGData = Phases[phase]['General']['SSGData']
1850            Vec,x,maxH = Phases[phase]['General']['SuperVec']
1851        pId = Phases[phase]['pId']
1852        histoList = HistoPhase.keys()
1853        histoList.sort()
1854        for histogram in histoList:
1855            try:
1856                Histogram = Histograms[histogram]
1857            except KeyError:                       
1858                #skip if histogram not included e.g. in a sequential refinement
1859                continue
1860            hapData = HistoPhase[histogram]
1861            hId = Histogram['hId']
1862            if 'PWDR' in histogram:
1863                limits = Histogram['Limits'][1]
1864                inst = Histogram['Instrument Parameters'][0]
1865                Zero = inst['Zero'][0]
1866                if 'C' in inst['Type'][1]:
1867                    try:
1868                        wave = inst['Lam'][1]
1869                    except KeyError:
1870                        wave = inst['Lam1'][1]
1871                    dmin = wave/(2.0*sind(limits[1]/2.0))
1872                elif 'T' in inst['Type'][0]:
1873                    dmin = limits[0]/inst['difC'][1]
1874                pfx = str(pId)+':'+str(hId)+':'
1875                for item in ['Scale','Extinction']:
1876                    hapDict[pfx+item] = hapData[item][0]
1877                    if hapData[item][1]:
1878                        hapVary.append(pfx+item)
1879                names = G2spc.HStrainNames(SGData)
1880                HSvals = []
1881                for i,name in enumerate(names):
1882                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1883                    HSvals.append(hapDict[pfx+name])
1884                    if hapData['HStrain'][1][i]:
1885                        hapVary.append(pfx+name)
1886                DIJS = G2spc.HStrainVals(HSvals,SGData)
1887                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1888                if hapData['Pref.Ori.'][0] == 'MD':
1889                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1890                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1891                    if hapData['Pref.Ori.'][2]:
1892                        hapVary.append(pfx+'MD')
1893                else:                           #'SH' spherical harmonics
1894                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1895                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1896                    for item in hapData['Pref.Ori.'][5]:
1897                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1898                        if hapData['Pref.Ori.'][2]:
1899                            hapVary.append(pfx+item)
1900                for item in ['Mustrain','Size']:
1901                    controlDict[pfx+item+'Type'] = hapData[item][0]
1902                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
1903                    if hapData[item][2][2]:
1904                        hapVary.append(pfx+item+';mx')
1905                    if hapData[item][0] in ['isotropic','uniaxial']:
1906                        hapDict[pfx+item+';i'] = hapData[item][1][0]
1907                        if hapData[item][2][0]:
1908                            hapVary.append(pfx+item+';i')
1909                        if hapData[item][0] == 'uniaxial':
1910                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1911                            hapDict[pfx+item+';a'] = hapData[item][1][1]
1912                            if hapData[item][2][1]:
1913                                hapVary.append(pfx+item+';a')
1914                    else:       #generalized for mustrain or ellipsoidal for size
1915                        Nterms = len(hapData[item][4])
1916                        if item == 'Mustrain':
1917                            names = G2spc.MustrainNames(SGData)
1918                            pwrs = []
1919                            for name in names:
1920                                h,k,l = name[1:]
1921                                pwrs.append([int(h),int(k),int(l)])
1922                            controlDict[pfx+'MuPwrs'] = pwrs
1923                        for i in range(Nterms):
1924                            sfx = ':'+str(i)
1925                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1926                            if hapData[item][5][i]:
1927                                hapVary.append(pfx+item+sfx)
1928                for bab in ['BabA','BabU']:
1929                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1930                    if hapData['Babinet'][bab][1]:
1931                        hapVary.append(pfx+bab)
1932                               
1933                if Print: 
1934                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1935                    print >>pFile,135*'-'
1936                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1937                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1938                    if hapData['Pref.Ori.'][0] == 'MD':
1939                        Ax = hapData['Pref.Ori.'][3]
1940                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1941                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1942                    else: #'SH' for spherical harmonics
1943                        PrintSHPO(hapData['Pref.Ori.'])
1944                    PrintSize(hapData['Size'])
1945                    PrintMuStrain(hapData['Mustrain'],SGData)
1946                    PrintHStrain(hapData['HStrain'],SGData)
1947                    if hapData['Babinet']['BabA'][0]:
1948                        PrintBabinet(hapData['Babinet'])                       
1949                if resetRefList:
1950                    refList = []
1951                    Uniq = []
1952                    Phi = []
1953                    if Phases[phase]['General']['Type'] in ['modulated','magnetic']:
1954                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
1955                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
1956                        for h,k,l,m,d in HKLd:
1957                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)    #is this right for SS refl.??
1958                            mul *= 2      # for powder overlap of Friedel pairs
1959                            if m or not ext:
1960                                if 'C' in inst['Type'][0]:
1961                                    pos = G2lat.Dsp2pos(inst,d)
1962                                    if limits[0] < pos < limits[1]:
1963                                        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])
1964                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
1965                                        Uniq.append(uniq)
1966                                        Phi.append(phi)
1967                                elif 'T' in inst['Type'][0]:
1968                                    pos = G2lat.Dsp2pos(inst,d)
1969                                    if limits[0] < pos < limits[1]:
1970                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
1971                                        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])
1972                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
1973                                        Uniq.append(uniq)
1974                                        Phi.append(phi)
1975                    else:
1976                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1977                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
1978                        for h,k,l,d in HKLd:
1979                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1980                            mul *= 2      # for powder overlap of Friedel pairs
1981                            if ext:
1982                                continue
1983                            if 'C' in inst['Type'][0]:
1984                                pos = G2lat.Dsp2pos(inst,d)
1985                                if limits[0] < pos < limits[1]:
1986                                    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])
1987                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
1988                                    Uniq.append(uniq)
1989                                    Phi.append(phi)
1990                            elif 'T' in inst['Type'][0]:
1991                                pos = G2lat.Dsp2pos(inst,d)
1992                                if limits[0] < pos < limits[1]:
1993                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
1994                                    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])
1995                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
1996                                    Uniq.append(uniq)
1997                                    Phi.append(phi)
1998                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0]}
1999            elif 'HKLF' in histogram:
2000                inst = Histogram['Instrument Parameters'][0]
2001                hId = Histogram['hId']
2002                hfx = ':%d:'%(hId)
2003                for item in inst:
2004                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2005                    hapDict[hfx+item] = inst[item][1]
2006                pfx = str(pId)+':'+str(hId)+':'
2007                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2008                if hapData['Scale'][1]:
2009                    hapVary.append(pfx+'Scale')
2010                               
2011                extApprox,extType,extParms = hapData['Extinction']
2012                controlDict[pfx+'EType'] = extType
2013                controlDict[pfx+'EApprox'] = extApprox
2014                if 'C' in inst['Type'][0]:
2015                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2016                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2017                if 'Primary' in extType:
2018                    Ekey = ['Ep',]
2019                elif 'I & II' in extType:
2020                    Ekey = ['Eg','Es']
2021                elif 'Secondary Type II' == extType:
2022                    Ekey = ['Es',]
2023                elif 'Secondary Type I' == extType:
2024                    Ekey = ['Eg',]
2025                else:   #'None'
2026                    Ekey = []
2027                for eKey in Ekey:
2028                    hapDict[pfx+eKey] = extParms[eKey][0]
2029                    if extParms[eKey][1]:
2030                        hapVary.append(pfx+eKey)
2031                for bab in ['BabA','BabU']:
2032                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2033                    if hapData['Babinet'][bab][1]:
2034                        hapVary.append(pfx+bab)
2035                if Print: 
2036                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2037                    print >>pFile,135*'-'
2038                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2039                    if extType != 'None':
2040                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox)
2041                        text = ' Parameters       :'
2042                        for eKey in Ekey:
2043                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2044                        print >>pFile,text
2045                    if hapData['Babinet']['BabA'][0]:
2046                        PrintBabinet(hapData['Babinet'])
2047                Histogram['Reflection Lists'] = phase       
2048               
2049    return hapVary,hapDict,controlDict
2050   
2051def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
2052    'needs a doc string'
2053   
2054    def PrintSizeAndSig(hapData,sizeSig):
2055        line = '\n Size model:     %9s'%(hapData[0])
2056        refine = False
2057        if hapData[0] in ['isotropic','uniaxial']:
2058            line += ' equatorial:%12.5f'%(hapData[1][0])
2059            if sizeSig[0][0]:
2060                line += ', sig:%8.4f'%(sizeSig[0][0])
2061                refine = True
2062            if hapData[0] == 'uniaxial':
2063                line += ' axial:%12.4f'%(hapData[1][1])
2064                if sizeSig[0][1]:
2065                    refine = True
2066                    line += ', sig:%8.4f'%(sizeSig[0][1])
2067            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2068            if sizeSig[0][2]:
2069                refine = True
2070                line += ', sig:%8.4f'%(sizeSig[0][2])
2071            if refine:
2072                print >>pFile,line
2073        else:
2074            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2075            if sizeSig[0][2]:
2076                refine = True
2077                line += ', sig:%8.4f'%(sizeSig[0][2])
2078            Snames = ['S11','S22','S33','S12','S13','S23']
2079            ptlbls = ' name  :'
2080            ptstr =  ' value :'
2081            sigstr = ' sig   :'
2082            for i,name in enumerate(Snames):
2083                ptlbls += '%12s' % (name)
2084                ptstr += '%12.6f' % (hapData[4][i])
2085                if sizeSig[1][i]:
2086                    refine = True
2087                    sigstr += '%12.6f' % (sizeSig[1][i])
2088                else:
2089                    sigstr += 12*' '
2090            if refine:
2091                print >>pFile,line
2092                print >>pFile,ptlbls
2093                print >>pFile,ptstr
2094                print >>pFile,sigstr
2095       
2096    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2097        line = '\n Mustrain model: %9s'%(hapData[0])
2098        refine = False
2099        if hapData[0] in ['isotropic','uniaxial']:
2100            line += ' equatorial:%12.1f'%(hapData[1][0])
2101            if mustrainSig[0][0]:
2102                line += ', sig:%8.1f'%(mustrainSig[0][0])
2103                refine = True
2104            if hapData[0] == 'uniaxial':
2105                line += ' axial:%12.1f'%(hapData[1][1])
2106                if mustrainSig[0][1]:
2107                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2108            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2109            if mustrainSig[0][2]:
2110                refine = True
2111                line += ', sig:%8.3f'%(mustrainSig[0][2])
2112            if refine:
2113                print >>pFile,line
2114        else:
2115            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2116            if mustrainSig[0][2]:
2117                refine = True
2118                line += ', sig:%8.3f'%(mustrainSig[0][2])
2119            Snames = G2spc.MustrainNames(SGData)
2120            ptlbls = ' name  :'
2121            ptstr =  ' value :'
2122            sigstr = ' sig   :'
2123            for i,name in enumerate(Snames):
2124                ptlbls += '%12s' % (name)
2125                ptstr += '%12.6f' % (hapData[4][i])
2126                if mustrainSig[1][i]:
2127                    refine = True
2128                    sigstr += '%12.6f' % (mustrainSig[1][i])
2129                else:
2130                    sigstr += 12*' '
2131            if refine:
2132                print >>pFile,line
2133                print >>pFile,ptlbls
2134                print >>pFile,ptstr
2135                print >>pFile,sigstr
2136           
2137    def PrintHStrainAndSig(hapData,strainSig,SGData):
2138        Hsnames = G2spc.HStrainNames(SGData)
2139        ptlbls = ' name  :'
2140        ptstr =  ' value :'
2141        sigstr = ' sig   :'
2142        refine = False
2143        for i,name in enumerate(Hsnames):
2144            ptlbls += '%12s' % (name)
2145            ptstr += '%12.4g' % (hapData[0][i])
2146            if name in strainSig:
2147                refine = True
2148                sigstr += '%12.4g' % (strainSig[name])
2149            else:
2150                sigstr += 12*' '
2151        if refine:
2152            print >>pFile,'\n Hydrostatic/elastic strain: '
2153            print >>pFile,ptlbls
2154            print >>pFile,ptstr
2155            print >>pFile,sigstr
2156       
2157    def PrintSHPOAndSig(pfx,hapData,POsig):
2158        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2159        ptlbls = ' names :'
2160        ptstr =  ' values:'
2161        sigstr = ' sig   :'
2162        for item in hapData[5]:
2163            ptlbls += '%12s'%(item)
2164            ptstr += '%12.3f'%(hapData[5][item])
2165            if pfx+item in POsig:
2166                sigstr += '%12.3f'%(POsig[pfx+item])
2167            else:
2168                sigstr += 12*' ' 
2169        print >>pFile,ptlbls
2170        print >>pFile,ptstr
2171        print >>pFile,sigstr
2172       
2173    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2174        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2175        text = ''
2176        for item in hapData[2]:
2177            if pfx+item in ScalExtSig:
2178                text += '       %s: '%(item)
2179                text += '%12.2e'%(hapData[2][item][0])
2180                if pfx+item in ScalExtSig:
2181                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2182        print >>pFile,text       
2183       
2184    def PrintBabinetAndSig(pfx,hapData,BabSig):
2185        print >>pFile,'\n Babinet form factor modification: '
2186        ptlbls = ' names :'
2187        ptstr =  ' values:'
2188        sigstr = ' sig   :'
2189        for item in hapData:
2190            ptlbls += '%12s'%(item)
2191            ptstr += '%12.3f'%(hapData[item][0])
2192            if pfx+item in BabSig:
2193                sigstr += '%12.3f'%(BabSig[pfx+item])
2194            else:
2195                sigstr += 12*' ' 
2196        print >>pFile,ptlbls
2197        print >>pFile,ptstr
2198        print >>pFile,sigstr
2199   
2200    PhFrExtPOSig = {}
2201    SizeMuStrSig = {}
2202    ScalExtSig = {}
2203    BabSig = {}
2204    wtFrSum = {}
2205    for phase in Phases:
2206        HistoPhase = Phases[phase]['Histograms']
2207        General = Phases[phase]['General']
2208        SGData = General['SGData']
2209        pId = Phases[phase]['pId']
2210        histoList = HistoPhase.keys()
2211        histoList.sort()
2212        for histogram in histoList:
2213            try:
2214                Histogram = Histograms[histogram]
2215            except KeyError:                       
2216                #skip if histogram not included e.g. in a sequential refinement
2217                continue
2218            hapData = HistoPhase[histogram]
2219            hId = Histogram['hId']
2220            pfx = str(pId)+':'+str(hId)+':'
2221            if hId not in wtFrSum:
2222                wtFrSum[hId] = 0.
2223            if 'PWDR' in histogram:
2224                for item in ['Scale','Extinction']:
2225                    hapData[item][0] = parmDict[pfx+item]
2226                    if pfx+item in sigDict:
2227                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2228                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2229                if hapData['Pref.Ori.'][0] == 'MD':
2230                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2231                    if pfx+'MD' in sigDict:
2232                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2233                else:                           #'SH' spherical harmonics
2234                    for item in hapData['Pref.Ori.'][5]:
2235                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2236                        if pfx+item in sigDict:
2237                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2238                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2239                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2240                    pfx+'HStrain':{}})                 
2241                for item in ['Mustrain','Size']:
2242                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2243                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2244                    if pfx+item+';mx' in sigDict:
2245                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2246                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2247                        hapData[item][1][0] = parmDict[pfx+item+';i']
2248                        if item == 'Size':
2249                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2250                        if pfx+item+';i' in sigDict: 
2251                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2252                        if hapData[item][0] == 'uniaxial':
2253                            hapData[item][1][1] = parmDict[pfx+item+';a']
2254                            if item == 'Size':
2255                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2256                            if pfx+item+';a' in sigDict:
2257                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2258                    else:       #generalized for mustrain or ellipsoidal for size
2259                        Nterms = len(hapData[item][4])
2260                        for i in range(Nterms):
2261                            sfx = ':'+str(i)
2262                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2263                            if pfx+item+sfx in sigDict:
2264                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2265                names = G2spc.HStrainNames(SGData)
2266                for i,name in enumerate(names):
2267                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2268                    if pfx+name in sigDict:
2269                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2270                for name in ['BabA','BabU']:
2271                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2272                    if pfx+name in sigDict:
2273                        BabSig[pfx+name] = sigDict[pfx+name]               
2274               
2275            elif 'HKLF' in histogram:
2276                for item in ['Scale',]:
2277                    if parmDict.get(pfx+item):
2278                        hapData[item][0] = parmDict[pfx+item]
2279                        if pfx+item in sigDict:
2280                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2281                for item in ['Ep','Eg','Es']:
2282                    if parmDict.get(pfx+item):
2283                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2284                        if pfx+item in sigDict:
2285                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2286                for name in ['BabA','BabU']:
2287                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2288                    if pfx+name in sigDict:
2289                        BabSig[pfx+name] = sigDict[pfx+name]               
2290
2291    if Print:
2292        for phase in Phases:
2293            HistoPhase = Phases[phase]['Histograms']
2294            General = Phases[phase]['General']
2295            SGData = General['SGData']
2296            pId = Phases[phase]['pId']
2297            histoList = HistoPhase.keys()
2298            histoList.sort()
2299            for histogram in histoList:
2300                try:
2301                    Histogram = Histograms[histogram]
2302                except KeyError:                       
2303                    #skip if histogram not included e.g. in a sequential refinement
2304                    continue
2305                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2306                print >>pFile,130*'-'
2307                hapData = HistoPhase[histogram]
2308                hId = Histogram['hId']
2309                Histogram['Residuals'][str(pId)+'::Name'] = phase
2310                pfx = str(pId)+':'+str(hId)+':'
2311                if 'PWDR' in histogram:
2312                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2313                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2314               
2315                    if pfx+'Scale' in PhFrExtPOSig:
2316                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2317                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2318                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2319                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2320                    if pfx+'Extinction' in PhFrExtPOSig:
2321                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2322                    if hapData['Pref.Ori.'][0] == 'MD':
2323                        if pfx+'MD' in PhFrExtPOSig:
2324                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2325                    else:
2326                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2327                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2328                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2329                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2330                    if len(BabSig):
2331                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2332                   
2333                elif 'HKLF' in histogram:
2334                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2335                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2336                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2337                    if pfx+'Scale' in ScalExtSig:
2338                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2339                    if hapData['Extinction'][0] != 'None':
2340                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2341                    if len(BabSig):
2342                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2343
2344################################################################################
2345##### Histogram data
2346################################################################################       
2347                   
2348def GetHistogramData(Histograms,Print=True,pFile=None):
2349    'needs a doc string'
2350   
2351    def GetBackgroundParms(hId,Background):
2352        Back = Background[0]
2353        DebyePeaks = Background[1]
2354        bakType,bakFlag = Back[:2]
2355        backVals = Back[3:]
2356        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2357        backDict = dict(zip(backNames,backVals))
2358        backVary = []
2359        if bakFlag:
2360            backVary = backNames
2361        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2362        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2363        debyeDict = {}
2364        debyeList = []
2365        for i in range(DebyePeaks['nDebye']):
2366            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2367            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2368            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2369        debyeVary = []
2370        for item in debyeList:
2371            if item[1]:
2372                debyeVary.append(item[0])
2373        backDict.update(debyeDict)
2374        backVary += debyeVary
2375        peakDict = {}
2376        peakList = []
2377        for i in range(DebyePeaks['nPeaks']):
2378            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2379                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2380            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2381            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2382        peakVary = []
2383        for item in peakList:
2384            if item[1]:
2385                peakVary.append(item[0])
2386        backDict.update(peakDict)
2387        backVary += peakVary
2388        return bakType,backDict,backVary           
2389       
2390    def GetInstParms(hId,Inst):     
2391        dataType = Inst['Type'][0]
2392        instDict = {}
2393        insVary = []
2394        pfx = ':'+str(hId)+':'
2395        insKeys = Inst.keys()
2396        insKeys.sort()
2397        for item in insKeys:
2398            insName = pfx+item
2399            instDict[insName] = Inst[item][1]
2400            if len(Inst[item]) > 2 and Inst[item][2]:
2401                insVary.append(insName)
2402        if 'C' in dataType:
2403            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2404        return dataType,instDict,insVary
2405       
2406    def GetSampleParms(hId,Sample):
2407        sampVary = []
2408        hfx = ':'+str(hId)+':'       
2409        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2410            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2411        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2412            if key in Sample:
2413                sampDict[hfx+key] = Sample[key]
2414        Type = Sample['Type']
2415        if 'Bragg' in Type:             #Bragg-Brentano
2416            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2417                sampDict[hfx+item] = Sample[item][0]
2418                if Sample[item][1]:
2419                    sampVary.append(hfx+item)
2420        elif 'Debye' in Type:        #Debye-Scherrer
2421            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2422                sampDict[hfx+item] = Sample[item][0]
2423                if Sample[item][1]:
2424                    sampVary.append(hfx+item)
2425        return Type,sampDict,sampVary
2426       
2427    def PrintBackground(Background):
2428        Back = Background[0]
2429        DebyePeaks = Background[1]
2430        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2431        line = ' Coefficients: '
2432        for i,back in enumerate(Back[3:]):
2433            line += '%10.3f'%(back)
2434            if i and not i%10:
2435                line += '\n'+15*' '
2436        print >>pFile,line
2437        if DebyePeaks['nDebye']:
2438            print >>pFile,'\n Debye diffuse scattering coefficients'
2439            parms = ['DebyeA','DebyeR','DebyeU']
2440            line = ' names :  '
2441            for parm in parms:
2442                line += '%8s refine?'%(parm)
2443            print >>pFile,line
2444            for j,term in enumerate(DebyePeaks['debyeTerms']):
2445                line = ' term'+'%2d'%(j)+':'
2446                for i in range(3):
2447                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2448                print >>pFile,line
2449        if DebyePeaks['nPeaks']:
2450            print >>pFile,'\n Single peak coefficients'
2451            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2452            line = ' names :  '
2453            for parm in parms:
2454                line += '%8s refine?'%(parm)
2455            print >>pFile,line
2456            for j,term in enumerate(DebyePeaks['peaksList']):
2457                line = ' peak'+'%2d'%(j)+':'
2458                for i in range(4):
2459                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2460                print >>pFile,line
2461       
2462    def PrintInstParms(Inst):
2463        print >>pFile,'\n Instrument Parameters:'
2464        insKeys = Inst.keys()
2465        insKeys.sort()
2466        iBeg = 0
2467        Ok = True
2468        while Ok:
2469            ptlbls = ' name  :'
2470            ptstr =  ' value :'
2471            varstr = ' refine:'
2472            iFin = min(iBeg+9,len(insKeys))
2473            for item in insKeys[iBeg:iFin]:
2474                if item not in ['Type','Source']:
2475                    ptlbls += '%12s' % (item)
2476                    ptstr += '%12.6f' % (Inst[item][1])
2477                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2478                        varstr += 12*' '
2479                    else:
2480                        varstr += '%12s' % (str(bool(Inst[item][2])))
2481            print >>pFile,ptlbls
2482            print >>pFile,ptstr
2483            print >>pFile,varstr
2484            iBeg = iFin
2485            if iBeg == len(insKeys):
2486                Ok = False
2487            else:
2488                print >>pFile,'\n'
2489       
2490    def PrintSampleParms(Sample):
2491        print >>pFile,'\n Sample Parameters:'
2492        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2493            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2494        ptlbls = ' name  :'
2495        ptstr =  ' value :'
2496        varstr = ' refine:'
2497        if 'Bragg' in Sample['Type']:
2498            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2499                ptlbls += '%14s'%(item)
2500                ptstr += '%14.4f'%(Sample[item][0])
2501                varstr += '%14s'%(str(bool(Sample[item][1])))
2502           
2503        elif 'Debye' in Type:        #Debye-Scherrer
2504            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2505                ptlbls += '%14s'%(item)
2506                ptstr += '%14.4f'%(Sample[item][0])
2507                varstr += '%14s'%(str(bool(Sample[item][1])))
2508
2509        print >>pFile,ptlbls
2510        print >>pFile,ptstr
2511        print >>pFile,varstr
2512       
2513    histDict = {}
2514    histVary = []
2515    controlDict = {}
2516    histoList = Histograms.keys()
2517    histoList.sort()
2518    for histogram in histoList:
2519        if 'PWDR' in histogram:
2520            Histogram = Histograms[histogram]
2521            hId = Histogram['hId']
2522            pfx = ':'+str(hId)+':'
2523            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2524            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2525            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2526            for excl in controlDict[pfx+'Exclude']:
2527                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2528            if controlDict[pfx+'Exclude']:
2529                ma.mask_rows(Histogram['Data'])
2530            Background = Histogram['Background']
2531            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2532            controlDict[pfx+'bakType'] = Type
2533            histDict.update(bakDict)
2534            histVary += bakVary
2535           
2536            Inst = Histogram['Instrument Parameters'][0]
2537            Type,instDict,insVary = GetInstParms(hId,Inst)
2538            controlDict[pfx+'histType'] = Type
2539            if 'XC' in Type:
2540                if pfx+'Lam1' in instDict:
2541                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2542                else:
2543                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2544            histDict.update(instDict)
2545            histVary += insVary
2546           
2547            Sample = Histogram['Sample Parameters']
2548            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2549            controlDict[pfx+'instType'] = Type
2550            histDict.update(sampDict)
2551            histVary += sampVary
2552           
2553   
2554            if Print: 
2555                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2556                print >>pFile,135*'-'
2557                Units = {'C':' deg','T':' msec'}
2558                units = Units[controlDict[pfx+'histType'][2]]
2559                Limits = controlDict[pfx+'Limits']
2560                print >>pFile,' Instrument type: ',Sample['Type']
2561                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2562                if len(controlDict[pfx+'Exclude']):
2563                    excls = controlDict[pfx+'Exclude']
2564                    for excl in excls:
2565                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
2566                PrintSampleParms(Sample)
2567                PrintInstParms(Inst)
2568                PrintBackground(Background)
2569        elif 'HKLF' in histogram:
2570            Histogram = Histograms[histogram]
2571            hId = Histogram['hId']
2572            pfx = ':'+str(hId)+':'
2573            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2574            Inst = Histogram['Instrument Parameters'][0]
2575            controlDict[pfx+'histType'] = Inst['Type'][0]
2576            if 'X' in Inst['Type'][0]:
2577                histDict[pfx+'Lam'] = Inst['Lam'][1]
2578                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2579    return histVary,histDict,controlDict
2580   
2581def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2582    'needs a doc string'
2583   
2584    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2585        Back = Background[0]
2586        DebyePeaks = Background[1]
2587        lenBack = len(Back[3:])
2588        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2589        for i in range(lenBack):
2590            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
2591            if pfx+'Back;'+str(i) in sigDict:
2592                backSig[i] = sigDict[pfx+'Back;'+str(i)]
2593        if DebyePeaks['nDebye']:
2594            for i in range(DebyePeaks['nDebye']):
2595                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2596                for j,name in enumerate(names):
2597                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2598                    if name in sigDict:
2599                        backSig[lenBack+3*i+j] = sigDict[name]           
2600        if DebyePeaks['nPeaks']:
2601            for i in range(DebyePeaks['nPeaks']):
2602                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2603                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2604                for j,name in enumerate(names):
2605                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2606                    if name in sigDict:
2607                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2608        return backSig
2609       
2610    def SetInstParms(pfx,Inst,parmDict,sigDict):
2611        instSig = {}
2612        insKeys = Inst.keys()
2613        insKeys.sort()
2614        for item in insKeys:
2615            insName = pfx+item
2616            Inst[item][1] = parmDict[insName]
2617            if insName in sigDict:
2618                instSig[item] = sigDict[insName]
2619            else:
2620                instSig[item] = 0
2621        return instSig
2622       
2623    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2624        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2625            sampSig = [0 for i in range(5)]
2626            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2627                Sample[item][0] = parmDict[pfx+item]
2628                if pfx+item in sigDict:
2629                    sampSig[i] = sigDict[pfx+item]
2630        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2631            sampSig = [0 for i in range(4)]
2632            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2633                Sample[item][0] = parmDict[pfx+item]
2634                if pfx+item in sigDict:
2635                    sampSig[i] = sigDict[pfx+item]
2636        return sampSig
2637       
2638    def PrintBackgroundSig(Background,backSig):
2639        Back = Background[0]
2640        DebyePeaks = Background[1]
2641        lenBack = len(Back[3:])
2642        valstr = ' value : '
2643        sigstr = ' sig   : '
2644        refine = False
2645        for i,back in enumerate(Back[3:]):
2646            valstr += '%10.4g'%(back)
2647            if Back[1]:
2648                refine = True
2649                sigstr += '%10.4g'%(backSig[i])
2650            else:
2651                sigstr += 10*' '
2652        if refine:
2653            print >>pFile,'\n Background function: ',Back[0]
2654            print >>pFile,valstr
2655            print >>pFile,sigstr
2656        if DebyePeaks['nDebye']:
2657            ifAny = False
2658            ptfmt = "%12.3f"
2659            names =  ' names :'
2660            ptstr =  ' values:'
2661            sigstr = ' esds  :'
2662            for item in sigDict:
2663                if 'Debye' in item:
2664                    ifAny = True
2665                    names += '%12s'%(item)
2666                    ptstr += ptfmt%(parmDict[item])
2667                    sigstr += ptfmt%(sigDict[item])
2668            if ifAny:
2669                print >>pFile,'\n Debye diffuse scattering coefficients'
2670                print >>pFile,names
2671                print >>pFile,ptstr
2672                print >>pFile,sigstr
2673        if DebyePeaks['nPeaks']:
2674            ifAny = False
2675            ptfmt = "%14.3f"
2676            names =  ' names :'
2677            ptstr =  ' values:'
2678            sigstr = ' esds  :'
2679            for item in sigDict:
2680                if 'BkPk' in item:
2681                    ifAny = True
2682                    names += '%14s'%(item)
2683                    ptstr += ptfmt%(parmDict[item])
2684                    sigstr += ptfmt%(sigDict[item])
2685            if ifAny:
2686                print >>pFile,'\n Single peak coefficients'
2687                print >>pFile,names
2688                print >>pFile,ptstr
2689                print >>pFile,sigstr
2690       
2691    def PrintInstParmsSig(Inst,instSig):
2692        refine = False
2693        insKeys = instSig.keys()
2694        insKeys.sort()
2695        iBeg = 0
2696        Ok = True
2697        while Ok:
2698            ptlbls = ' names :'
2699            ptstr =  ' value :'
2700            sigstr = ' sig   :'
2701            iFin = min(iBeg+9,len(insKeys))
2702            for name in insKeys[iBeg:iFin]:
2703                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
2704                    ptlbls += '%12s' % (name)
2705                    ptstr += '%12.6f' % (Inst[name][1])
2706                    if instSig[name]:
2707                        refine = True
2708                        sigstr += '%12.6f' % (instSig[name])
2709                    else:
2710                        sigstr += 12*' '
2711            if refine:
2712                print >>pFile,'\n Instrument Parameters:'
2713                print >>pFile,ptlbls
2714                print >>pFile,ptstr
2715                print >>pFile,sigstr
2716            iBeg = iFin
2717            if iBeg == len(insKeys):
2718                Ok = False
2719       
2720    def PrintSampleParmsSig(Sample,sampleSig):
2721        ptlbls = ' names :'
2722        ptstr =  ' values:'
2723        sigstr = ' sig   :'
2724        refine = False
2725        if 'Bragg' in Sample['Type']:
2726            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2727                ptlbls += '%14s'%(item)
2728                ptstr += '%14.4f'%(Sample[item][0])
2729                if sampleSig[i]:
2730                    refine = True
2731                    sigstr += '%14.4f'%(sampleSig[i])
2732                else:
2733                    sigstr += 14*' '
2734           
2735        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2736            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2737                ptlbls += '%14s'%(item)
2738                ptstr += '%14.4f'%(Sample[item][0])
2739                if sampleSig[i]:
2740                    refine = True
2741                    sigstr += '%14.4f'%(sampleSig[i])
2742                else:
2743                    sigstr += 14*' '
2744
2745        if refine:
2746            print >>pFile,'\n Sample Parameters:'
2747            print >>pFile,ptlbls
2748            print >>pFile,ptstr
2749            print >>pFile,sigstr
2750       
2751    histoList = Histograms.keys()
2752    histoList.sort()
2753    for histogram in histoList:
2754        if 'PWDR' in histogram:
2755            Histogram = Histograms[histogram]
2756            hId = Histogram['hId']
2757            pfx = ':'+str(hId)+':'
2758            Background = Histogram['Background']
2759            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2760           
2761            Inst = Histogram['Instrument Parameters'][0]
2762            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2763       
2764            Sample = Histogram['Sample Parameters']
2765            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2766
2767            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2768            print >>pFile,135*'-'
2769            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2770            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2771            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2772                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2773            if Print:
2774                print >>pFile,' Instrument type: ',Sample['Type']
2775                PrintSampleParmsSig(Sample,sampSig)
2776                PrintInstParmsSig(Inst,instSig)
2777                PrintBackgroundSig(Background,backSig)
2778               
Note: See TracBrowser for help on using the repository browser.