source: trunk/GSASIIstruct.py @ 599

Last change on this file since 599 was 599, checked in by vondreele, 10 years ago

GSASIIstruct.py - set up a DEBUG flag
comment out some dead code
begin inclusion of penalty fxns.
add Tutorials to help stuff

  • Property svn:keywords set to Date Author Revision URL Id
File size: 148.6 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIstructure - structure computation routines
3########### SVN repository information ###################
4# $Date: 2012-05-07 19:59:54 +0000 (Mon, 07 May 2012) $
5# $Author: vondreele $
6# $Revision: 599 $
7# $URL: trunk/GSASIIstruct.py $
8# $Id: GSASIIstruct.py 599 2012-05-07 19:59:54Z vondreele $
9########### SVN repository information ###################
10import sys
11import os
12import os.path as ospath
13import time
14import math
15import cPickle
16import numpy as np
17import numpy.linalg as nl
18import scipy.optimize as so
19import GSASIIpath
20import GSASIIElem as G2el
21import GSASIIlattice as G2lat
22import GSASIIspc as G2spc
23import GSASIIpwd as G2pwd
24import GSASIImapvars as G2mv
25import GSASIImath as G2mth
26
27sind = lambda x: np.sin(x*np.pi/180.)
28cosd = lambda x: np.cos(x*np.pi/180.)
29tand = lambda x: np.tan(x*np.pi/180.)
30asind = lambda x: 180.*np.arcsin(x)/np.pi
31acosd = lambda x: 180.*np.arccos(x)/np.pi
32atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
33   
34ateln2 = 8.0*math.log(2.0)
35DEBUG = False
36
37
38def GetControls(GPXfile):
39    ''' Returns dictionary of control items found in GSASII gpx file
40    input:
41        GPXfile = .gpx full file name
42    return:
43        Controls = dictionary of control items
44    '''
45    Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
46        'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
47    fl = open(GPXfile,'rb')
48    while True:
49        try:
50            data = cPickle.load(fl)
51        except EOFError:
52            break
53        datum = data[0]
54        if datum[0] == 'Controls':
55            Controls.update(datum[1])
56    fl.close()
57    return Controls
58   
59def GetConstraints(GPXfile):
60    '''Read the constraints from the GPX file and interpret them
61    '''
62    constList = []
63    fl = open(GPXfile,'rb')
64    while True:
65        try:
66            data = cPickle.load(fl)
67        except EOFError:
68            break
69        datum = data[0]
70        if datum[0] == 'Constraints':
71            constDict = datum[1]
72            for item in constDict:
73                constList += constDict[item]
74    fl.close()
75    constDict,fixedList,ignored = ProcessConstraints(constList)
76    if ignored:
77        print ignored,'old-style Constraints were rejected'
78    return constDict,fixedList
79   
80def ProcessConstraints(constList):
81    "interpret constraints"
82    constDict = []
83    fixedList = []
84    ignored = 0
85    for item in constList:
86        if item[-1] == 'h':
87            # process a hold
88            fixedList.append('0')
89            constDict.append({item[0][1]:0.0})
90        elif item[-1] == 'f':
91            # process a new variable
92            fixedList.append(None)
93            constDict.append({})
94            for term in item[:-3]:
95                constDict[-1][term[1]] = term[0]
96            #constFlag[-1] = ['Vary']
97        elif item[-1] == 'c': 
98            # process a contraint relationship
99            fixedList.append(str(item[-3]))
100            constDict.append({})
101            for term in item[:-3]:
102                constDict[-1][term[1]] = term[0]
103            #constFlag[-1] = ['VaryFree']
104        elif item[-1] == 'e':
105            # process an equivalence
106            firstmult = None
107            eqlist = []
108            for term in item[:-3]:
109                if term[0] == 0: term[0] = 1.0
110                if firstmult is None:
111                    firstmult,firstvar = term
112                else:
113                    eqlist.append([term[1],firstmult/term[0]])
114            G2mv.StoreEquivalence(firstvar,eqlist)
115        else:
116            ignored += 1
117    return constDict,fixedList,ignored
118
119def CheckConstraints(GPXfile):
120    '''Load constraints and related info and return any error or warning messages'''
121    # init constraints
122    G2mv.InitVars()   
123    # get variables
124    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
125    if not Phases:
126        return 'Error: No Phases!',''
127    if not Histograms:
128        return 'Error: no diffraction data',''
129    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
130    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
131    histVary,histDict,controlDict = GetHistogramData(Histograms)
132    varyList = phaseVary+hapVary+histVary
133    constrDict,fixedList = GetConstraints(GPXfile)
134    return G2mv.CheckConstraints(varyList,constrDict,fixedList)
135   
136def GetPhaseNames(GPXfile):
137    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
138    input:
139        GPXfile = gpx full file name
140    return:
141        PhaseNames = list of phase names
142    '''
143    fl = open(GPXfile,'rb')
144    PhaseNames = []
145    while True:
146        try:
147            data = cPickle.load(fl)
148        except EOFError:
149            break
150        datum = data[0]
151        if 'Phases' == datum[0]:
152            for datus in data[1:]:
153                PhaseNames.append(datus[0])
154    fl.close()
155    return PhaseNames
156
157def GetAllPhaseData(GPXfile,PhaseName):
158    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
159    input:
160        GPXfile = gpx full file name
161        PhaseName = phase name
162    return:
163        phase dictionary
164    '''       
165    fl = open(GPXfile,'rb')
166    General = {}
167    Atoms = []
168    while True:
169        try:
170            data = cPickle.load(fl)
171        except EOFError:
172            break
173        datum = data[0]
174        if 'Phases' == datum[0]:
175            for datus in data[1:]:
176                if datus[0] == PhaseName:
177                    break
178    fl.close()
179    return datus[1]
180   
181def GetHistograms(GPXfile,hNames):
182    """ Returns a dictionary of histograms found in GSASII gpx file
183    input:
184        GPXfile = .gpx full file name
185        hNames = list of histogram names
186    return:
187        Histograms = dictionary of histograms (types = PWDR & HKLF)
188    """
189    fl = open(GPXfile,'rb')
190    Histograms = {}
191    while True:
192        try:
193            data = cPickle.load(fl)
194        except EOFError:
195            break
196        datum = data[0]
197        hist = datum[0]
198        if hist in hNames:
199            if 'PWDR' in hist[:4]:
200                PWDRdata = {}
201                PWDRdata['Data'] = datum[1][1]          #powder data arrays
202                PWDRdata[data[2][0]] = data[2][1]       #Limits
203                PWDRdata[data[3][0]] = data[3][1]       #Background
204                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
205                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
206                try:
207                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
208                except IndexError:
209                    PWDRdata['Reflection lists'] = {}
210   
211                Histograms[hist] = PWDRdata
212            elif 'HKLF' in hist[:4]:
213                HKLFdata = []
214                datum = data[0]
215                HKLFdata = datum[1:][0]
216                Histograms[hist] = HKLFdata           
217    fl.close()
218    return Histograms
219   
220def GetHistogramNames(GPXfile,hType):
221    """ Returns a list of histogram names found in GSASII gpx file
222    input:
223        GPXfile = .gpx full file name
224        hType = list ['PWDR','HKLF']
225    return:
226        HistogramNames = list of histogram names (types = PWDR & HKLF)
227    """
228    fl = open(GPXfile,'rb')
229    HistogramNames = []
230    while True:
231        try:
232            data = cPickle.load(fl)
233        except EOFError:
234            break
235        datum = data[0]
236        if datum[0][:4] in hType:
237            HistogramNames.append(datum[0])
238    fl.close()
239    return HistogramNames
240   
241def GetUsedHistogramsAndPhases(GPXfile):
242    ''' Returns all histograms that are found in any phase
243    and any phase that uses a histogram
244    input:
245        GPXfile = .gpx full file name
246    return:
247        Histograms = dictionary of histograms as {name:data,...}
248        Phases = dictionary of phases that use histograms
249    '''
250    phaseNames = GetPhaseNames(GPXfile)
251    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
252    allHistograms = GetHistograms(GPXfile,histoList)
253    phaseData = {}
254    for name in phaseNames: 
255        phaseData[name] =  GetAllPhaseData(GPXfile,name)
256    Histograms = {}
257    Phases = {}
258    for phase in phaseData:
259        Phase = phaseData[phase]
260        if Phase['Histograms']:
261            if phase not in Phases:
262                pId = phaseNames.index(phase)
263                Phase['pId'] = pId
264                Phases[phase] = Phase
265            for hist in Phase['Histograms']:
266                if hist not in Histograms:
267                    Histograms[hist] = allHistograms[hist]
268                    #future restraint, etc. histograms here           
269                    hId = histoList.index(hist)
270                    Histograms[hist]['hId'] = hId
271    return Histograms,Phases
272   
273def getBackupName2(GPXfile,makeBack=True):      #not work correctly
274    GPXpath,GPXname = ospath.split(GPXfile)
275    if GPXpath == '': GPXpath = '.'
276    Name = ospath.splitext(GPXname)[0]
277    files = os.listdir(GPXpath)
278    last = 0
279    for name in files:
280        name = name.split('.')
281        if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]:
282            if makeBack:
283                last = max(last,int(name[-2].strip('bak'))+1)
284            else:
285                last = max(last,int(name[-2].strip('bak')))
286    GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx')
287    return GPXback
288
289def getBackupName(GPXfile,makeBack):       #recovered old one
290    GPXpath,GPXname = ospath.split(GPXfile)
291    if GPXpath == '': GPXpath = '.'
292    Name = ospath.splitext(GPXname)[0]
293    files = os.listdir(GPXpath)
294    last = 0
295    for name in files:
296        name = name.split('.')
297        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
298            if makeBack:
299                last = max(last,int(name[1].strip('bak'))+1)
300            else:
301                last = max(last,int(name[1].strip('bak')))
302    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
303    return GPXback   
304       
305def GPXBackup(GPXfile,makeBack=True):
306    import distutils.file_util as dfu
307    GPXback = getBackupName(GPXfile,makeBack)
308    dfu.copy_file(GPXfile,GPXback)
309    return GPXback
310
311def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
312    ''' Updates gpxfile from all histograms that are found in any phase
313    and any phase that used a histogram
314    input:
315        GPXfile = .gpx full file name
316        Histograms = dictionary of histograms as {name:data,...}
317        Phases = dictionary of phases that use histograms
318        CovData = dictionary of refined variables, varyList, & covariance matrix
319        makeBack = True if new backup of .gpx file is to be made; else use the last one made
320    '''
321                       
322    GPXback = GPXBackup(GPXfile,makeBack)
323    print '\n',135*'-'
324    print 'Read from file:',GPXback
325    print 'Save to file  :',GPXfile
326    infile = open(GPXback,'rb')
327    outfile = open(GPXfile,'wb')
328    while True:
329        try:
330            data = cPickle.load(infile)
331        except EOFError:
332            break
333        datum = data[0]
334#        print 'read: ',datum[0]
335        if datum[0] == 'Phases':
336            for iphase in range(len(data)):
337                if data[iphase][0] in Phases:
338                    phaseName = data[iphase][0]
339                    data[iphase][1].update(Phases[phaseName])
340        elif datum[0] == 'Covariance':
341            data[0][1] = CovData
342        try:
343            histogram = Histograms[datum[0]]
344#            print 'found ',datum[0]
345            data[0][1][1] = histogram['Data']
346            for datus in data[1:]:
347#                print '    read: ',datus[0]
348                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
349                    datus[1] = histogram[datus[0]]
350        except KeyError:
351            pass
352                               
353        cPickle.dump(data,outfile,1)
354    infile.close()
355    outfile.close()
356    print 'GPX file save successful'
357   
358def SetSeqResult(GPXfile,Histograms,SeqResult):
359    GPXback = GPXBackup(GPXfile)
360    print '\n',135*'-'
361    print 'Read from file:',GPXback
362    print 'Save to file  :',GPXfile
363    infile = open(GPXback,'rb')
364    outfile = open(GPXfile,'wb')
365    while True:
366        try:
367            data = cPickle.load(infile)
368        except EOFError:
369            break
370        datum = data[0]
371        if datum[0] == 'Sequental results':
372            data[0][1] = SeqResult
373        try:
374            histogram = Histograms[datum[0]]
375            data[0][1][1] = histogram['Data']
376            for datus in data[1:]:
377                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
378                    datus[1] = histogram[datus[0]]
379        except KeyError:
380            pass
381                               
382        cPickle.dump(data,outfile,1)
383    infile.close()
384    outfile.close()
385    print 'GPX file save successful'
386                       
387def GetPWDRdata(GPXfile,PWDRname):
388    ''' Returns powder data from GSASII gpx file
389    input:
390        GPXfile = .gpx full file name
391        PWDRname = powder histogram name as obtained from GetHistogramNames
392    return:
393        PWDRdata = powder data dictionary with:
394            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
395       
396    '''
397    fl = open(GPXfile,'rb')
398    PWDRdata = {}
399    while True:
400        try:
401            data = cPickle.load(fl)
402        except EOFError:
403            break
404        datum = data[0]
405        if datum[0] == PWDRname:
406            PWDRdata['Data'] = datum[1][1]          #powder data arrays
407            PWDRdata[data[2][0]] = data[2][1]       #Limits
408            PWDRdata[data[3][0]] = data[3][1]       #Background
409            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
410            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
411            try:
412                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
413            except IndexError:
414                PWDRdata['Reflection lists'] = {}
415    fl.close()
416    return PWDRdata
417   
418def GetHKLFdata(GPXfile,HKLFname):
419    ''' Returns single crystal data from GSASII gpx file
420    input:
421        GPXfile = .gpx full file name
422        HKLFname = single crystal histogram name as obtained from GetHistogramNames
423    return:
424        HKLFdata = single crystal data list of reflections: for each reflection:
425            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
426    '''
427    fl = open(GPXfile,'rb')
428    HKLFdata = []
429    while True:
430        try:
431            data = cPickle.load(fl)
432        except EOFError:
433            break
434        datum = data[0]
435        if datum[0] == HKLFname:
436            HKLFdata = datum[1:][0]
437    fl.close()
438    return HKLFdata
439   
440def ShowBanner():
441    print 80*'*'
442    print '   General Structure Analysis System-II Crystal Structure Refinement'
443    print '              by Robert B. Von Dreele & Brian H. Toby'
444    print '                Argonne National Laboratory(C), 2010'
445    print ' This product includes software developed by the UChicago Argonne, LLC,' 
446    print '            as Operator of Argonne National Laboratory.'
447    print 80*'*','\n'
448
449def ShowControls(Controls):
450    print ' Least squares controls:'
451    print ' Refinement type: ',Controls['deriv type']
452    if 'Hessian' in Controls['deriv type']:
453        print ' Maximum number of cycles:',Controls['max cyc']
454    else:
455        print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
456    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
457   
458def GetFFtable(General):
459    ''' returns a dictionary of form factor data for atom types found in General
460    input:
461        General = dictionary of phase info.; includes AtomTypes
462    return:
463        FFtable = dictionary of form factor data; key is atom type
464    '''
465    atomTypes = General['AtomTypes']
466    FFtable = {}
467    for El in atomTypes:
468        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
469        for item in FFs:
470            if item['Symbol'] == El.upper():
471                FFtable[El] = item
472    return FFtable
473   
474def GetBLtable(General):
475    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
476    input:
477        General = dictionary of phase info.; includes AtomTypes & Isotopes
478    return:
479        BLtable = dictionary of scattering length data; key is atom type
480    '''
481    atomTypes = General['AtomTypes']
482    BLtable = {}
483    isotopes = General['Isotopes']
484    isotope = General['Isotope']
485    for El in atomTypes:
486        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
487    return BLtable
488       
489def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
490    if SGLaue in ['-1','2/m','mmm']:
491        return                      #no Pawley symmetry required constraints
492    eqvDict = {}
493    for i,varyI in enumerate(pawleyVary):
494        eqvDict[varyI] = []
495        refI = int(varyI.split(':')[-1])
496        ih,ik,il = PawleyRef[refI][:3]
497        for varyJ in pawleyVary[i+1:]:
498            refJ = int(varyJ.split(':')[-1])
499            jh,jk,jl = PawleyRef[refJ][:3]
500            if SGLaue in ['4/m','4/mmm']:
501                isum = ih**2+ik**2
502                jsum = jh**2+jk**2
503                if abs(il) == abs(jl) and isum == jsum:
504                    eqvDict[varyI].append(varyJ) 
505            elif SGLaue in ['3R','3mR']:
506                isum = ih**2+ik**2+il**2
507                jsum = jh**2+jk**2*jl**2
508                isum2 = ih*ik+ih*il+ik*il
509                jsum2 = jh*jk+jh*jl+jk*jl
510                if isum == jsum and isum2 == jsum2:
511                    eqvDict[varyI].append(varyJ) 
512            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
513                isum = ih**2+ik**2+ih*ik
514                jsum = jh**2+jk**2+jh*jk
515                if abs(il) == abs(jl) and isum == jsum:
516                    eqvDict[varyI].append(varyJ) 
517            elif SGLaue in ['m3','m3m']:
518                isum = ih**2+ik**2+il**2
519                jsum = jh**2+jk**2+jl**2
520                if isum == jsum:
521                    eqvDict[varyI].append(varyJ) 
522    for item in pawleyVary:
523        if eqvDict[item]:
524            for item2 in pawleyVary:
525                if item2 in eqvDict[item]:
526                    eqvDict[item2] = []
527            G2mv.StoreEquivalence(item,eqvDict[item])
528                   
529def cellVary(pfx,SGData): 
530    if SGData['SGLaue'] in ['-1',]:
531        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
532    elif SGData['SGLaue'] in ['2/m',]:
533        if SGData['SGUniq'] == 'a':
534            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
535        elif SGData['SGUniq'] == 'b':
536            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
537        else:
538            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
539    elif SGData['SGLaue'] in ['mmm',]:
540        return [pfx+'A0',pfx+'A1',pfx+'A2']
541    elif SGData['SGLaue'] in ['4/m','4/mmm']:
542        return [pfx+'A0',pfx+'A2']
543    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
544        return [pfx+'A0',pfx+'A2']
545    elif SGData['SGLaue'] in ['3R', '3mR']:
546        return [pfx+'A0',pfx+'A3']                       
547    elif SGData['SGLaue'] in ['m3m','m3']:
548        return [pfx+'A0',]
549       
550################################################################################
551##### Phase data
552################################################################################       
553                   
554def GetPhaseData(PhaseData,Print=True):
555           
556    def PrintFFtable(FFtable):
557        print '\n X-ray scattering factors:'
558        print '   Symbol     fa                                      fb                                      fc'
559        print 99*'-'
560        for Ename in FFtable:
561            ffdata = FFtable[Ename]
562            fa = ffdata['fa']
563            fb = ffdata['fb']
564            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
565                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
566               
567    def PrintBLtable(BLtable):
568        print '\n Neutron scattering factors:'
569        print '   Symbol   isotope       mass       b       resonant terms'
570        print 99*'-'
571        for Ename in BLtable:
572            bldata = BLtable[Ename]
573            isotope = bldata[0]
574            mass = bldata[1][0]
575            blen = bldata[1][1]
576            bres = []
577            if len(bldata[1]) > 2:
578                bres = bldata[1][2:]
579            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
580            for item in bres:
581                line += '%10.5g'%(item)
582            print line
583               
584    def PrintAtoms(General,Atoms):
585        print '\n Atoms:'
586        line = '   name    type  refine?   x         y         z    '+ \
587            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
588        if General['Type'] == 'magnetic':
589            line += '   Mx     My     Mz'
590        elif General['Type'] == 'macromolecular':
591            line = ' res no  residue  chain '+line
592        print line
593        if General['Type'] == 'nuclear':
594            print 135*'-'
595            for i,at in enumerate(Atoms):
596                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
597                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
598                if at[9] == 'I':
599                    line += '%8.4f'%(at[10])+48*' '
600                else:
601                    line += 8*' '
602                    for j in range(6):
603                        line += '%8.4f'%(at[11+j])
604                print line
605       
606    def PrintTexture(textureData):
607        topstr = '\n Spherical harmonics texture: Order:' + \
608            str(textureData['Order'])
609        if textureData['Order']:
610            print topstr+' Refine? '+str(textureData['SH Coeff'][0])
611        else:
612            print topstr
613            return
614        names = ['omega','chi','phi']
615        line = '\n'
616        for name in names:
617            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
618        print line
619        print '\n Texture coefficients:'
620        ptlbls = ' names :'
621        ptstr =  ' values:'
622        SHcoeff = textureData['SH Coeff'][1]
623        for item in SHcoeff:
624            ptlbls += '%12s'%(item)
625            ptstr += '%12.4f'%(SHcoeff[item]) 
626        print ptlbls
627        print ptstr   
628       
629    if Print: print ' Phases:'
630    phaseVary = []
631    phaseDict = {}
632    phaseConstr = {}
633    pawleyLookup = {}
634    FFtables = {}                   #scattering factors - xrays
635    BLtables = {}                   # neutrons
636    Natoms = {}
637    AtMults = {}
638    AtIA = {}
639    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
640    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
641    for name in PhaseData:
642        General = PhaseData[name]['General']
643        pId = PhaseData[name]['pId']
644        pfx = str(pId)+'::'
645        FFtable = GetFFtable(General)
646        BLtable = GetBLtable(General)
647        FFtables.update(FFtable)
648        BLtables.update(BLtable)
649        Atoms = PhaseData[name]['Atoms']
650        try:
651            PawleyRef = PhaseData[name]['Pawley ref']
652        except KeyError:
653            PawleyRef = []
654        SGData = General['SGData']
655        SGtext = G2spc.SGPrint(SGData)
656        cell = General['Cell']
657        A = G2lat.cell2A(cell[1:7])
658        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
659        if cell[0]:
660            phaseVary += cellVary(pfx,SGData)
661        Natoms[pfx] = 0
662        if Atoms and not General.get('doPawley'):
663            if General['Type'] == 'nuclear':
664                Natoms[pfx] = len(Atoms)
665                for i,at in enumerate(Atoms):
666                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
667                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
668                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
669                        pfx+'AI/A:'+str(i):at[9],})
670                    if at[9] == 'I':
671                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
672                    else:
673                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
674                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
675                    if 'F' in at[2]:
676                        phaseVary.append(pfx+'Afrac:'+str(i))
677                    if 'X' in at[2]:
678                        xId,xCoef = G2spc.GetCSxinel(at[7])
679                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
680                        for j in range(3):
681                            if xId[j] > 0:                               
682                                phaseVary.append(delnames[j])
683                                for k in range(j):
684                                    if xId[j] == xId[k]:
685                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) 
686                    if 'U' in at[2]:
687                        if at[9] == 'I':
688                            phaseVary.append(pfx+'AUiso:'+str(i))
689                        else:
690                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
691                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
692                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
693                            for j in range(6):
694                                if uId[j] > 0:                               
695                                    phaseVary.append(names[j])
696                                    for k in range(j):
697                                        if uId[j] == uId[k]:
698                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
699#            elif General['Type'] == 'magnetic':
700#            elif General['Type'] == 'macromolecular':
701
702            textureData = General['SH Texture']
703            if textureData['Order']:
704                phaseDict[pfx+'SHorder'] = textureData['Order']
705                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
706                for name in ['omega','chi','phi']:
707                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
708                    if textureData['Sample '+name][0]:
709                        phaseVary.append(pfx+'SH '+name)
710                for name in textureData['SH Coeff'][1]:
711                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
712                    if textureData['SH Coeff'][0]:
713                        phaseVary.append(pfx+name)
714               
715            if Print:
716                print '\n Phase name: ',General['Name']
717                print 135*'-'
718                PrintFFtable(FFtable)
719                PrintBLtable(BLtable)
720                print ''
721                for line in SGtext: print line
722                PrintAtoms(General,Atoms)
723                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
724                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
725                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
726                PrintTexture(textureData)
727                   
728        elif PawleyRef:
729            pawleyVary = []
730            for i,refl in enumerate(PawleyRef):
731                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
732                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
733                if refl[5]:
734                    pawleyVary.append(pfx+'PWLref:'+str(i))
735            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
736            phaseVary += pawleyVary
737               
738    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
739   
740def cellFill(pfx,SGData,parmDict,sigDict): 
741    if SGData['SGLaue'] in ['-1',]:
742        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
743            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
744        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
745            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
746    elif SGData['SGLaue'] in ['2/m',]:
747        if SGData['SGUniq'] == 'a':
748            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
749                parmDict[pfx+'A3'],0,0]
750            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
751                sigDict[pfx+'A3'],0,0]
752        elif SGData['SGUniq'] == 'b':
753            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
754                0,parmDict[pfx+'A4'],0]
755            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
756                0,sigDict[pfx+'A4'],0]
757        else:
758            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
759                0,0,parmDict[pfx+'A5']]
760            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
761                0,0,sigDict[pfx+'A5']]
762    elif SGData['SGLaue'] in ['mmm',]:
763        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
764        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
765    elif SGData['SGLaue'] in ['4/m','4/mmm']:
766        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
767        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
768    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
769        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
770            parmDict[pfx+'A0'],0,0]
771        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
772    elif SGData['SGLaue'] in ['3R', '3mR']:
773        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
774            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
775        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
776    elif SGData['SGLaue'] in ['m3m','m3']:
777        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
778        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
779    return A,sigA
780       
781def getCellEsd(pfx,SGData,A,covData):
782    dpr = 180./np.pi
783    rVsq = G2lat.calc_rVsq(A)
784    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
785    cell = np.array(G2lat.Gmat2cell(g))   #real cell
786    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
787    scos = cosd(cellst[3:6])
788    ssin = sind(cellst[3:6])
789    scot = scos/ssin
790    rcos = cosd(cell[3:6])
791    rsin = sind(cell[3:6])
792    rcot = rcos/rsin
793    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
794    varyList = covData['varyList']
795    covMatrix = covData['covMatrix']
796    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
797    Ax = np.array(A)
798    Ax[3:] /= 2.
799    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
800        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
801    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
802    Vol = 1/np.sqrt(rVsq)
803    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
804    R123 = Ax[0]*Ax[1]*Ax[2]
805    dsasdg = np.zeros((3,6))
806    dadg = np.zeros((6,6))
807    for i0 in range(3):         #0  1   2
808        i1 = (i0+1)%3           #1  2   0
809        i2 = (i1+1)%3           #2  0   1
810        i3 = 5-i2               #3  5   4
811        i4 = 5-i1               #4  3   5
812        i5 = 5-i0               #5  4   3
813        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
814        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
815        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
816        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
817        denom = np.sqrt(denmsq)
818        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
819        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
820        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
821        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
822        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
823        dadg[i5][i5] = -Ax[i0]/denom
824    for i0 in range(3):
825        i1 = (i0+1)%3
826        i2 = (i1+1)%3
827        i3 = 5-i2
828        for ij in range(6):
829            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
830            if ij == i0:
831                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
832            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
833    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
834    var = np.diag(sigMat)
835    CS = np.where(var>0.,np.sqrt(var),0.)
836    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
837    return cellSig           
838   
839def SetPhaseData(parmDict,sigDict,Phases,covData):
840   
841    def PrintAtomsAndSig(General,Atoms,atomsSig):
842        print '\n Atoms:'
843        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
844        if General['Type'] == 'magnetic':
845            line += '   Mx     My     Mz'
846        elif General['Type'] == 'macromolecular':
847            line = ' res no  residue  chain '+line
848        print line
849        if General['Type'] == 'nuclear':
850            print 135*'-'
851            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
852                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
853            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
854            for i,at in enumerate(Atoms):
855                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
856                valstr = ' values:'
857                sigstr = ' sig   :'
858                for ind in [3,4,5,6]:
859                    sigind = str(i)+':'+str(ind)
860                    valstr += fmt[ind]%(at[ind])                   
861                    if sigind in atomsSig:
862                        sigstr += fmt[ind]%(atomsSig[sigind])
863                    else:
864                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
865                if at[9] == 'I':
866                    valstr += fmt[10]%(at[10])
867                    if str(i)+':10' in atomsSig:
868                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
869                    else:
870                        sigstr += 8*' '
871                else:
872                    valstr += 8*' '
873                    sigstr += 8*' '
874                    for ind in [11,12,13,14,15,16]:
875                        sigind = str(i)+':'+str(ind)
876                        valstr += fmt[ind]%(at[ind])
877                        if sigind in atomsSig:                       
878                            sigstr += fmt[ind]%(atomsSig[sigind])
879                        else:
880                            sigstr += 8*' '
881                print name
882                print valstr
883                print sigstr
884               
885    def PrintSHtextureAndSig(textureData,SHtextureSig):
886        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
887        names = ['omega','chi','phi']
888        namstr = '  names :'
889        ptstr =  '  values:'
890        sigstr = '  esds  :'
891        for name in names:
892            namstr += '%12s'%(name)
893            ptstr += '%12.3f'%(textureData['Sample '+name][1])
894            if 'Sample '+name in SHtextureSig:
895                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
896            else:
897                sigstr += 12*' '
898        print namstr
899        print ptstr
900        print sigstr
901        print '\n Texture coefficients:'
902        namstr = '  names :'
903        ptstr =  '  values:'
904        sigstr = '  esds  :'
905        SHcoeff = textureData['SH Coeff'][1]
906        for name in SHcoeff:
907            namstr += '%12s'%(name)
908            ptstr += '%12.3f'%(SHcoeff[name])
909            if name in SHtextureSig:
910                sigstr += '%12.3f'%(SHtextureSig[name])
911            else:
912                sigstr += 12*' '
913        print namstr
914        print ptstr
915        print sigstr
916       
917           
918    print '\n Phases:'
919    for phase in Phases:
920        print ' Result for phase: ',phase
921        Phase = Phases[phase]
922        General = Phase['General']
923        SGData = General['SGData']
924        Atoms = Phase['Atoms']
925        cell = General['Cell']
926        pId = Phase['pId']
927        pfx = str(pId)+'::'
928        if cell[0]:
929            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
930            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
931            print ' Reciprocal metric tensor: '
932            ptfmt = "%15.9f"
933            names = ['A11','A22','A33','A12','A13','A23']
934            namstr = '  names :'
935            ptstr =  '  values:'
936            sigstr = '  esds  :'
937            for name,a,siga in zip(names,A,sigA):
938                namstr += '%15s'%(name)
939                ptstr += ptfmt%(a)
940                if siga:
941                    sigstr += ptfmt%(siga)
942                else:
943                    sigstr += 15*' '
944            print namstr
945            print ptstr
946            print sigstr
947            cell[1:7] = G2lat.A2cell(A)
948            cell[7] = G2lat.calc_V(A)
949            print ' New unit cell:'
950            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
951            names = ['a','b','c','alpha','beta','gamma','Volume']
952            namstr = '  names :'
953            ptstr =  '  values:'
954            sigstr = '  esds  :'
955            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
956                namstr += '%12s'%(name)
957                ptstr += fmt%(a)
958                if siga:
959                    sigstr += fmt%(siga)
960                else:
961                    sigstr += 12*' '
962            print namstr
963            print ptstr
964            print sigstr
965           
966        if Phase['General'].get('doPawley'):
967            pawleyRef = Phase['Pawley ref']
968            for i,refl in enumerate(pawleyRef):
969                key = pfx+'PWLref:'+str(i)
970                refl[6] = parmDict[key]
971                if key in sigDict:
972                    refl[7] = sigDict[key]
973                else:
974                    refl[7] = 0
975        else:
976            atomsSig = {}
977            if General['Type'] == 'nuclear':
978                for i,at in enumerate(Atoms):
979                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
980                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
981                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
982                    for ind in [3,4,5,6]:
983                        at[ind] = parmDict[names[ind]]
984                        if ind in [3,4,5]:
985                            name = names[ind].replace('A','dA')
986                        else:
987                            name = names[ind]
988                        if name in sigDict:
989                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
990                    if at[9] == 'I':
991                        at[10] = parmDict[names[10]]
992                        if names[10] in sigDict:
993                            atomsSig[str(i)+':10'] = sigDict[names[10]]
994                    else:
995                        for ind in [11,12,13,14,15,16]:
996                            at[ind] = parmDict[names[ind]]
997                            if names[ind] in sigDict:
998                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
999            PrintAtomsAndSig(General,Atoms,atomsSig)
1000       
1001        textureData = General['SH Texture']   
1002        if textureData['Order']:
1003            SHtextureSig = {}
1004            for name in ['omega','chi','phi']:
1005                aname = pfx+'SH '+name
1006                textureData['Sample '+name][1] = parmDict[aname]
1007                if aname in sigDict:
1008                    SHtextureSig['Sample '+name] = sigDict[aname]
1009            for name in textureData['SH Coeff'][1]:
1010                aname = pfx+name
1011                textureData['SH Coeff'][1][name] = parmDict[aname]
1012                if aname in sigDict:
1013                    SHtextureSig[name] = sigDict[aname]
1014            PrintSHtextureAndSig(textureData,SHtextureSig)
1015
1016################################################################################
1017##### Histogram & Phase data
1018################################################################################       
1019                   
1020def GetHistogramPhaseData(Phases,Histograms,Print=True):
1021   
1022    def PrintSize(hapData):
1023        if hapData[0] in ['isotropic','uniaxial']:
1024            line = '\n Size model    : %9s'%(hapData[0])
1025            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1026            if hapData[0] == 'uniaxial':
1027                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1028            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1029            print line
1030        else:
1031            print '\n Size model    : %s'%(hapData[0])+ \
1032                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1033            Snames = ['S11','S22','S33','S12','S13','S23']
1034            ptlbls = ' names :'
1035            ptstr =  ' values:'
1036            varstr = ' refine:'
1037            for i,name in enumerate(Snames):
1038                ptlbls += '%12s' % (name)
1039                ptstr += '%12.6f' % (hapData[4][i])
1040                varstr += '%12s' % (str(hapData[5][i]))
1041            print ptlbls
1042            print ptstr
1043            print varstr
1044       
1045    def PrintMuStrain(hapData,SGData):
1046        if hapData[0] in ['isotropic','uniaxial']:
1047            line = '\n Mustrain model: %9s'%(hapData[0])
1048            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1049            if hapData[0] == 'uniaxial':
1050                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1051            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1052            print line
1053        else:
1054            print '\n Mustrain model: %s'%(hapData[0])+ \
1055                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1056            Snames = G2spc.MustrainNames(SGData)
1057            ptlbls = ' names :'
1058            ptstr =  ' values:'
1059            varstr = ' refine:'
1060            for i,name in enumerate(Snames):
1061                ptlbls += '%12s' % (name)
1062                ptstr += '%12.6f' % (hapData[4][i])
1063                varstr += '%12s' % (str(hapData[5][i]))
1064            print ptlbls
1065            print ptstr
1066            print varstr
1067
1068    def PrintHStrain(hapData,SGData):
1069        print '\n Hydrostatic/elastic strain: '
1070        Hsnames = G2spc.HStrainNames(SGData)
1071        ptlbls = ' names :'
1072        ptstr =  ' values:'
1073        varstr = ' refine:'
1074        for i,name in enumerate(Hsnames):
1075            ptlbls += '%12s' % (name)
1076            ptstr += '%12.6f' % (hapData[0][i])
1077            varstr += '%12s' % (str(hapData[1][i]))
1078        print ptlbls
1079        print ptstr
1080        print varstr
1081
1082    def PrintSHPO(hapData):
1083        print '\n Spherical harmonics preferred orientation: Order:' + \
1084            str(hapData[4])+' Refine? '+str(hapData[2])
1085        ptlbls = ' names :'
1086        ptstr =  ' values:'
1087        for item in hapData[5]:
1088            ptlbls += '%12s'%(item)
1089            ptstr += '%12.3f'%(hapData[5][item]) 
1090        print ptlbls
1091        print ptstr
1092   
1093    hapDict = {}
1094    hapVary = []
1095    controlDict = {}
1096    poType = {}
1097    poAxes = {}
1098    spAxes = {}
1099    spType = {}
1100   
1101    for phase in Phases:
1102        HistoPhase = Phases[phase]['Histograms']
1103        SGData = Phases[phase]['General']['SGData']
1104        cell = Phases[phase]['General']['Cell'][1:7]
1105        A = G2lat.cell2A(cell)
1106        pId = Phases[phase]['pId']
1107        histoList = HistoPhase.keys()
1108        histoList.sort()
1109        for histogram in histoList:
1110            try:
1111                Histogram = Histograms[histogram]
1112            except KeyError:                       
1113                #skip if histogram not included e.g. in a sequential refinement
1114                continue
1115            hapData = HistoPhase[histogram]
1116            hId = Histogram['hId']
1117            limits = Histogram['Limits'][1]
1118            inst = Histogram['Instrument Parameters']
1119            inst = dict(zip(inst[3],inst[1]))
1120            Zero = inst['Zero']
1121            if 'C' in inst['Type']:
1122                try:
1123                    wave = inst['Lam']
1124                except KeyError:
1125                    wave = inst['Lam1']
1126                dmin = wave/(2.0*sind(limits[1]/2.0))
1127            pfx = str(pId)+':'+str(hId)+':'
1128            for item in ['Scale','Extinction']:
1129                hapDict[pfx+item] = hapData[item][0]
1130                if hapData[item][1]:
1131                    hapVary.append(pfx+item)
1132            names = G2spc.HStrainNames(SGData)
1133            for i,name in enumerate(names):
1134                hapDict[pfx+name] = hapData['HStrain'][0][i]
1135                if hapData['HStrain'][1][i]:
1136                    hapVary.append(pfx+name)
1137            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1138            if hapData['Pref.Ori.'][0] == 'MD':
1139                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1140                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1141                if hapData['Pref.Ori.'][2]:
1142                    hapVary.append(pfx+'MD')
1143            else:                           #'SH' spherical harmonics
1144                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1145                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1146                for item in hapData['Pref.Ori.'][5]:
1147                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1148                    if hapData['Pref.Ori.'][2]:
1149                        hapVary.append(pfx+item)
1150            for item in ['Mustrain','Size']:
1151                controlDict[pfx+item+'Type'] = hapData[item][0]
1152                hapDict[pfx+item+':mx'] = hapData[item][1][2]
1153                if hapData[item][2][2]:
1154                    hapVary.append(pfx+item+':mx')
1155                if hapData[item][0] in ['isotropic','uniaxial']:
1156                    hapDict[pfx+item+':i'] = hapData[item][1][0]
1157                    if hapData[item][2][0]:
1158                        hapVary.append(pfx+item+':i')
1159                    if hapData[item][0] == 'uniaxial':
1160                        controlDict[pfx+item+'Axis'] = hapData[item][3]
1161                        hapDict[pfx+item+':a'] = hapData[item][1][1]
1162                        if hapData[item][2][1]:
1163                            hapVary.append(pfx+item+':a')
1164                else:       #generalized for mustrain or ellipsoidal for size
1165                    Nterms = len(hapData[item][4])
1166                    if item == 'Mustrain':
1167                        names = G2spc.MustrainNames(SGData)
1168                        pwrs = []
1169                        for name in names:
1170                            h,k,l = name[1:]
1171                            pwrs.append([int(h),int(k),int(l)])
1172                        controlDict[pfx+'MuPwrs'] = pwrs
1173                    for i in range(Nterms):
1174                        sfx = ':'+str(i)
1175                        hapDict[pfx+item+sfx] = hapData[item][4][i]
1176                        if hapData[item][5][i]:
1177                            hapVary.append(pfx+item+sfx)
1178                           
1179            if Print: 
1180                print '\n Phase: ',phase,' in histogram: ',histogram
1181                print 135*'-'
1182                print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1183                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1184                if hapData['Pref.Ori.'][0] == 'MD':
1185                    Ax = hapData['Pref.Ori.'][3]
1186                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1187                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1188                else: #'SH' for spherical harmonics
1189                    PrintSHPO(hapData['Pref.Ori.'])
1190                PrintSize(hapData['Size'])
1191                PrintMuStrain(hapData['Mustrain'],SGData)
1192                PrintHStrain(hapData['HStrain'],SGData)
1193            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1194            refList = []
1195            for h,k,l,d in HKLd:
1196                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1197                if ext:
1198                    continue
1199                if 'C' in inst['Type']:
1200                    pos = 2.0*asind(wave/(2.0*d))+Zero
1201                    if limits[0] < pos < limits[1]:
1202                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1203                        #last item should contain form factor values by atom type
1204                else:
1205                    raise ValueError 
1206            Histogram['Reflection Lists'][phase] = refList
1207    return hapVary,hapDict,controlDict
1208   
1209def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
1210   
1211    def PrintSizeAndSig(hapData,sizeSig):
1212        line = '\n Size model:     %9s'%(hapData[0])
1213        refine = False
1214        if hapData[0] in ['isotropic','uniaxial']:
1215            line += ' equatorial:%12.3f'%(hapData[1][0])
1216            if sizeSig[0][0]:
1217                line += ', sig:%8.3f'%(sizeSig[0][0])
1218                refine = True
1219            if hapData[0] == 'uniaxial':
1220                line += ' axial:%12.3f'%(hapData[1][1])
1221                if sizeSig[0][1]:
1222                    refine = True
1223                    line += ', sig:%8.3f'%(sizeSig[0][1])
1224            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1225            if sizeSig[0][2]:
1226                refine = True
1227                line += ', sig:%8.3f'%(sizeSig[0][2])
1228            if refine:
1229                print line
1230        else:
1231            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1232            if sizeSig[0][2]:
1233                refine = True
1234                line += ', sig:%8.3f'%(sizeSig[0][2])
1235            Snames = ['S11','S22','S33','S12','S13','S23']
1236            ptlbls = ' name  :'
1237            ptstr =  ' value :'
1238            sigstr = ' sig   :'
1239            for i,name in enumerate(Snames):
1240                ptlbls += '%12s' % (name)
1241                ptstr += '%12.6f' % (hapData[4][i])
1242                if sizeSig[1][i]:
1243                    refine = True
1244                    sigstr += '%12.6f' % (sizeSig[1][i])
1245                else:
1246                    sigstr += 12*' '
1247            if refine:
1248                print line
1249                print ptlbls
1250                print ptstr
1251                print sigstr
1252       
1253    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1254        line = '\n Mustrain model: %9s'%(hapData[0])
1255        refine = False
1256        if hapData[0] in ['isotropic','uniaxial']:
1257            line += ' equatorial:%12.1f'%(hapData[1][0])
1258            if mustrainSig[0][0]:
1259                line += ', sig: %8.1f'%(mustrainSig[0][0])
1260                refine = True
1261            if hapData[0] == 'uniaxial':
1262                line += ' axial:%12.1f'%(hapData[1][1])
1263                if mustrainSig[0][1]:
1264                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1265            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1266            if mustrainSig[0][2]:
1267                refine = True
1268                line += ', sig:%8.3f'%(mustrainSig[0][2])
1269            if refine:
1270                print line
1271        else:
1272            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1273            if mustrainSig[0][2]:
1274                refine = True
1275                line += ', sig:%8.3f'%(mustrainSig[0][2])
1276            Snames = G2spc.MustrainNames(SGData)
1277            ptlbls = ' name  :'
1278            ptstr =  ' value :'
1279            sigstr = ' sig   :'
1280            for i,name in enumerate(Snames):
1281                ptlbls += '%12s' % (name)
1282                ptstr += '%12.6f' % (hapData[4][i])
1283                if mustrainSig[1][i]:
1284                    refine = True
1285                    sigstr += '%12.6f' % (mustrainSig[1][i])
1286                else:
1287                    sigstr += 12*' '
1288            if refine:
1289                print line
1290                print ptlbls
1291                print ptstr
1292                print sigstr
1293           
1294    def PrintHStrainAndSig(hapData,strainSig,SGData):
1295        Hsnames = G2spc.HStrainNames(SGData)
1296        ptlbls = ' name  :'
1297        ptstr =  ' value :'
1298        sigstr = ' sig   :'
1299        refine = False
1300        for i,name in enumerate(Hsnames):
1301            ptlbls += '%12s' % (name)
1302            ptstr += '%12.6g' % (hapData[0][i])
1303            if name in strainSig:
1304                refine = True
1305                sigstr += '%12.6g' % (strainSig[name])
1306            else:
1307                sigstr += 12*' '
1308        if refine:
1309            print '\n Hydrostatic/elastic strain: '
1310            print ptlbls
1311            print ptstr
1312            print sigstr
1313       
1314    def PrintSHPOAndSig(hapData,POsig):
1315        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1316        ptlbls = ' names :'
1317        ptstr =  ' values:'
1318        sigstr = ' sig   :'
1319        for item in hapData[5]:
1320            ptlbls += '%12s'%(item)
1321            ptstr += '%12.3f'%(hapData[5][item])
1322            if item in POsig:
1323                sigstr += '%12.3f'%(POsig[item])
1324            else:
1325                sigstr += 12*' ' 
1326        print ptlbls
1327        print ptstr
1328        print sigstr
1329   
1330    for phase in Phases:
1331        HistoPhase = Phases[phase]['Histograms']
1332        SGData = Phases[phase]['General']['SGData']
1333        pId = Phases[phase]['pId']
1334        histoList = HistoPhase.keys()
1335        histoList.sort()
1336        for histogram in histoList:
1337            try:
1338                Histogram = Histograms[histogram]
1339            except KeyError:                       
1340                #skip if histogram not included e.g. in a sequential refinement
1341                continue
1342            print '\n Phase: ',phase,' in histogram: ',histogram
1343            print 130*'-'
1344            hapData = HistoPhase[histogram]
1345            hId = Histogram['hId']
1346            pfx = str(pId)+':'+str(hId)+':'
1347            print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
1348                %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
1349           
1350            PhFrExtPOSig = {}
1351            for item in ['Scale','Extinction']:
1352                hapData[item][0] = parmDict[pfx+item]
1353                if pfx+item in sigDict:
1354                    PhFrExtPOSig[item] = sigDict[pfx+item]
1355            if hapData['Pref.Ori.'][0] == 'MD':
1356                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1357                if pfx+'MD' in sigDict:
1358                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
1359            else:                           #'SH' spherical harmonics
1360                for item in hapData['Pref.Ori.'][5]:
1361                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1362                    if pfx+item in sigDict:
1363                        PhFrExtPOSig[item] = sigDict[pfx+item]
1364            if Print:
1365                if 'Scale' in PhFrExtPOSig:
1366                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
1367                if 'Extinction' in PhFrExtPOSig:
1368                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
1369                if hapData['Pref.Ori.'][0] == 'MD':
1370                    if 'MD' in PhFrExtPOSig:
1371                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
1372                else:
1373                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
1374            SizeMuStrSig = {'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1375                'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1376                'HStrain':{}}                 
1377            for item in ['Mustrain','Size']:
1378                hapData[item][1][2] = parmDict[pfx+item+':mx']
1379                hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1380                if pfx+item+':mx' in sigDict:
1381                    SizeMuStrSig[item][0][2] = sigDict[pfx+item+':mx']
1382                if hapData[item][0] in ['isotropic','uniaxial']:                   
1383                    hapData[item][1][0] = parmDict[pfx+item+':i']
1384                    if item == 'Size':
1385                        hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1386                    if pfx+item+':i' in sigDict: 
1387                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
1388                    if hapData[item][0] == 'uniaxial':
1389                        hapData[item][1][1] = parmDict[pfx+item+':a']
1390                        if item == 'Size':
1391                            hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
1392                        if pfx+item+':a' in sigDict:
1393                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
1394                else:       #generalized for mustrain or ellipsoidal for size
1395                    Nterms = len(hapData[item][4])
1396                    for i in range(Nterms):
1397                        sfx = ':'+str(i)
1398                        hapData[item][4][i] = parmDict[pfx+item+sfx]
1399                        if pfx+item+sfx in sigDict:
1400                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
1401            names = G2spc.HStrainNames(SGData)
1402            for i,name in enumerate(names):
1403                hapData['HStrain'][0][i] = parmDict[pfx+name]
1404                if pfx+name in sigDict:
1405                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
1406            if Print:
1407                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
1408                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
1409                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
1410   
1411################################################################################
1412##### Histogram data
1413################################################################################       
1414                   
1415def GetHistogramData(Histograms,Print=True):
1416   
1417    def GetBackgroundParms(hId,Background):
1418        Back = Background[0]
1419        Debye = Background[1]
1420        bakType,bakFlag = Back[:2]
1421        backVals = Back[3:]
1422        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
1423        backDict = dict(zip(backNames,backVals))
1424        backVary = []
1425        if bakFlag:
1426            backVary = backNames
1427        backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
1428        debyeDict = {}
1429        debyeList = []
1430        for i in range(Debye['nDebye']):
1431            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
1432            debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
1433            debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
1434        debyeVary = []
1435        for item in debyeList:
1436            if item[1]:
1437                debyeVary.append(item[0])
1438        backDict.update(debyeDict)
1439        backVary += debyeVary   
1440        return bakType,backDict,backVary           
1441       
1442    def GetInstParms(hId,Inst):
1443        insVals,insFlags,insNames = Inst[1:4]
1444        dataType = insVals[0]
1445        instDict = {}
1446        insVary = []
1447        pfx = ':'+str(hId)+':'
1448        for i,flag in enumerate(insFlags):
1449            insName = pfx+insNames[i]
1450            instDict[insName] = insVals[i]
1451            if flag:
1452                insVary.append(insName)
1453        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
1454        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
1455        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
1456        return dataType,instDict,insVary
1457       
1458    def GetSampleParms(hId,Sample):
1459        sampVary = []
1460        hfx = ':'+str(hId)+':'       
1461        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1462            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1463        Type = Sample['Type']
1464        if 'Bragg' in Type:             #Bragg-Brentano
1465            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1466                sampDict[hfx+item] = Sample[item][0]
1467                if Sample[item][1]:
1468                    sampVary.append(hfx+item)
1469        elif 'Debye' in Type:        #Debye-Scherrer
1470            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1471                sampDict[hfx+item] = Sample[item][0]
1472                if Sample[item][1]:
1473                    sampVary.append(hfx+item)
1474        return Type,sampDict,sampVary
1475       
1476    def PrintBackground(Background):
1477        Back = Background[0]
1478        Debye = Background[1]
1479        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
1480        line = ' Coefficients: '
1481        for i,back in enumerate(Back[3:]):
1482            line += '%10.3f'%(back)
1483            if i and not i%10:
1484                line += '\n'+15*' '
1485        print line
1486        if Debye['nDebye']:
1487            print '\n Debye diffuse scattering coefficients'
1488            parms = ['DebyeA','DebyeR','DebyeU']
1489            line = ' names :'
1490            for parm in parms:
1491                line += '%16s'%(parm)
1492            print line
1493            for j,term in enumerate(Debye['debyeTerms']):
1494                line = ' term'+'%2d'%(j)+':'
1495                for i in range(3):
1496                    line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1]))                   
1497                print line
1498       
1499    def PrintInstParms(Inst):
1500        print '\n Instrument Parameters:'
1501        ptlbls = ' name  :'
1502        ptstr =  ' value :'
1503        varstr = ' refine:'
1504        instNames = Inst[3][1:]
1505        for i,name in enumerate(instNames):
1506            ptlbls += '%12s' % (name)
1507            ptstr += '%12.6f' % (Inst[1][i+1])
1508            if name in ['Lam1','Lam2','Azimuth']:
1509                varstr += 12*' '
1510            else:
1511                varstr += '%12s' % (str(bool(Inst[2][i+1])))
1512        print ptlbls
1513        print ptstr
1514        print varstr
1515       
1516    def PrintSampleParms(Sample):
1517        print '\n Sample Parameters:'
1518        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1519            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1520        ptlbls = ' name  :'
1521        ptstr =  ' value :'
1522        varstr = ' refine:'
1523        if 'Bragg' in Sample['Type']:
1524            for item in ['Scale','Shift','Transparency']:
1525                ptlbls += '%14s'%(item)
1526                ptstr += '%14.4f'%(Sample[item][0])
1527                varstr += '%14s'%(str(bool(Sample[item][1])))
1528           
1529        elif 'Debye' in Type:        #Debye-Scherrer
1530            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1531                ptlbls += '%14s'%(item)
1532                ptstr += '%14.4f'%(Sample[item][0])
1533                varstr += '%14s'%(str(bool(Sample[item][1])))
1534
1535        print ptlbls
1536        print ptstr
1537        print varstr
1538       
1539
1540    histDict = {}
1541    histVary = []
1542    controlDict = {}
1543    histoList = Histograms.keys()
1544    histoList.sort()
1545    for histogram in histoList:
1546        Histogram = Histograms[histogram]
1547        hId = Histogram['hId']
1548        pfx = ':'+str(hId)+':'
1549        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1550       
1551        Background = Histogram['Background']
1552        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1553        controlDict[pfx+'bakType'] = Type
1554        histDict.update(bakDict)
1555        histVary += bakVary
1556       
1557        Inst = Histogram['Instrument Parameters']
1558        Type,instDict,insVary = GetInstParms(hId,Inst)
1559        controlDict[pfx+'histType'] = Type
1560        if pfx+'Lam1' in instDict:
1561            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1562        else:
1563            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1564        histDict.update(instDict)
1565        histVary += insVary
1566       
1567        Sample = Histogram['Sample Parameters']
1568        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1569        controlDict[pfx+'instType'] = Type
1570        histDict.update(sampDict)
1571        histVary += sampVary
1572
1573        if Print: 
1574            print '\n Histogram: ',histogram,' histogram Id: ',hId
1575            print 135*'-'
1576            Units = {'C':' deg','T':' msec'}
1577            units = Units[controlDict[pfx+'histType'][2]]
1578            Limits = controlDict[pfx+'Limits']
1579            print ' Instrument type: ',Sample['Type']
1580            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1581            PrintSampleParms(Sample)
1582            PrintInstParms(Inst)
1583            PrintBackground(Background)
1584       
1585    return histVary,histDict,controlDict
1586   
1587def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1588   
1589    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1590        Back = Background[0]
1591        Debye = Background[1]
1592        lenBack = len(Back[3:])
1593        backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
1594        for i in range(lenBack):
1595            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1596            if pfx+'Back:'+str(i) in sigDict:
1597                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1598        if Debye['nDebye']:
1599            for i in range(Debye['nDebye']):
1600                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1601                for j,name in enumerate(names):
1602                    Debye['debyeTerms'][i][2*j] = parmDict[name]
1603                    if name in sigDict:
1604                        backSig[lenBack+3*i+j] = sigDict[name]           
1605        return backSig
1606       
1607    def SetInstParms(pfx,Inst,parmDict,sigDict):
1608        insVals,insFlags,insNames = Inst[1:4]
1609        instSig = [0 for i in range(len(insVals))]
1610        for i,flag in enumerate(insFlags):
1611            insName = pfx+insNames[i]
1612            insVals[i] = parmDict[insName]
1613            if insName in sigDict:
1614                instSig[i] = sigDict[insName]
1615        return instSig
1616       
1617    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1618        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1619            sampSig = [0 for i in range(3)]
1620            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1621                Sample[item][0] = parmDict[pfx+item]
1622                if pfx+item in sigDict:
1623                    sampSig[i] = sigDict[pfx+item]
1624        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1625            sampSig = [0 for i in range(4)]
1626            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1627                Sample[item][0] = parmDict[pfx+item]
1628                if pfx+item in sigDict:
1629                    sampSig[i] = sigDict[pfx+item]
1630        return sampSig
1631       
1632    def PrintBackgroundSig(Background,backSig):
1633        Back = Background[0]
1634        Debye = Background[1]
1635        lenBack = len(Back[3:])
1636        valstr = ' value : '
1637        sigstr = ' sig   : '
1638        refine = False
1639        for i,back in enumerate(Back[3:]):
1640            valstr += '%10.4g'%(back)
1641            if Back[1]:
1642                refine = True
1643                sigstr += '%10.4g'%(backSig[i])
1644            else:
1645                sigstr += 10*' '
1646        if refine:
1647            print '\n Background function: ',Back[0]
1648            print valstr
1649            print sigstr
1650        if Debye['nDebye']:
1651            ifAny = False
1652            ptfmt = "%12.5f"
1653            names =  ' names :'
1654            ptstr =  ' values:'
1655            sigstr = ' esds  :'
1656            for item in sigDict:
1657                if 'Debye' in item:
1658                    ifAny = True
1659                    names += '%12s'%(item)
1660                    ptstr += ptfmt%(parmDict[item])
1661                    sigstr += ptfmt%(sigDict[item])
1662            if ifAny:
1663                print '\n Debye diffuse scattering coefficients'
1664                print names
1665                print ptstr
1666                print sigstr
1667       
1668    def PrintInstParmsSig(Inst,instSig):
1669        ptlbls = ' names :'
1670        ptstr =  ' value :'
1671        sigstr = ' sig   :'
1672        instNames = Inst[3][1:]
1673        refine = False
1674        for i,name in enumerate(instNames):
1675            ptlbls += '%12s' % (name)
1676            ptstr += '%12.6f' % (Inst[1][i+1])
1677            if instSig[i+1]:
1678                refine = True
1679                sigstr += '%12.6f' % (instSig[i+1])
1680            else:
1681                sigstr += 12*' '
1682        if refine:
1683            print '\n Instrument Parameters:'
1684            print ptlbls
1685            print ptstr
1686            print sigstr
1687       
1688    def PrintSampleParmsSig(Sample,sampleSig):
1689        ptlbls = ' names :'
1690        ptstr =  ' values:'
1691        sigstr = ' sig   :'
1692        refine = False
1693        if 'Bragg' in Sample['Type']:
1694            for i,item in enumerate(['Scale','Shift','Transparency']):
1695                ptlbls += '%14s'%(item)
1696                ptstr += '%14.4f'%(Sample[item][0])
1697                if sampleSig[i]:
1698                    refine = True
1699                    sigstr += '%14.4f'%(sampleSig[i])
1700                else:
1701                    sigstr += 14*' '
1702           
1703        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1704            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1705                ptlbls += '%14s'%(item)
1706                ptstr += '%14.4f'%(Sample[item][0])
1707                if sampleSig[i]:
1708                    refine = True
1709                    sigstr += '%14.4f'%(sampleSig[i])
1710                else:
1711                    sigstr += 14*' '
1712
1713        if refine:
1714            print '\n Sample Parameters:'
1715            print ptlbls
1716            print ptstr
1717            print sigstr
1718       
1719    histoList = Histograms.keys()
1720    histoList.sort()
1721    for histogram in histoList:
1722        if 'PWDR' in histogram:
1723            Histogram = Histograms[histogram]
1724            hId = Histogram['hId']
1725            pfx = ':'+str(hId)+':'
1726            Background = Histogram['Background']
1727            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1728           
1729            Inst = Histogram['Instrument Parameters']
1730            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1731       
1732            Sample = Histogram['Sample Parameters']
1733            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1734
1735            print '\n Histogram: ',histogram,' histogram Id: ',hId
1736            print 135*'-'
1737            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1738            if Print:
1739                print ' Instrument type: ',Sample['Type']
1740                PrintSampleParmsSig(Sample,sampSig)
1741                PrintInstParmsSig(Inst,instSig)
1742                PrintBackgroundSig(Background,backSig)
1743               
1744################################################################################
1745##### Penalty & restraint functions
1746################################################################################
1747
1748def getPenalties():
1749    return []
1750   
1751def penaltyFxn(penalties):
1752    return []
1753   
1754def penaltyDeriv(penalties):
1755    return []
1756
1757################################################################################
1758##### Function & derivative calculations
1759################################################################################       
1760                   
1761def GetAtomFXU(pfx,calcControls,parmDict):
1762    Natoms = calcControls['Natoms'][pfx]
1763    Tdata = Natoms*[' ',]
1764    Mdata = np.zeros(Natoms)
1765    IAdata = Natoms*[' ',]
1766    Fdata = np.zeros(Natoms)
1767    FFdata = []
1768    BLdata = []
1769    Xdata = np.zeros((3,Natoms))
1770    dXdata = np.zeros((3,Natoms))
1771    Uisodata = np.zeros(Natoms)
1772    Uijdata = np.zeros((6,Natoms))
1773    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1774        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1775        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1776        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1777        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1778    for iatm in range(Natoms):
1779        for key in keys:
1780            parm = pfx+key+str(iatm)
1781            if parm in parmDict:
1782                keys[key][iatm] = parmDict[parm]
1783    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1784   
1785def getFFvalues(FFtables,SQ):
1786    FFvals = {}
1787    for El in FFtables:
1788        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
1789    return FFvals
1790   
1791def getBLvalues(BLtables):
1792    BLvals = {}
1793    for El in BLtables:
1794        BLvals[El] = BLtables[El][1][1]
1795    return BLvals
1796       
1797def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1798    ''' Compute structure factors for all h,k,l for phase
1799    input:
1800        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1801        G:      reciprocal metric tensor
1802        pfx:    phase id string
1803        SGData: space group info. dictionary output from SpcGroup
1804        calcControls:
1805        ParmDict:
1806    puts result F^2 in each ref[8] in refList
1807    '''       
1808    twopi = 2.0*np.pi
1809    twopisq = 2.0*np.pi**2
1810    ast = np.sqrt(np.diag(G))
1811    Mast = twopisq*np.multiply.outer(ast,ast)
1812    FFtables = calcControls['FFtables']
1813    BLtables = calcControls['BLtables']
1814    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1815    FF = np.zeros(len(Tdata))
1816    if 'N' in parmDict[hfx+'Type']:
1817        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
1818    else:
1819        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1820        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1821    maxPos = len(SGData['SGOps'])
1822    Uij = np.array(G2lat.U6toUij(Uijdata))
1823    bij = Mast*Uij.T
1824    for refl in refList:
1825        fbs = np.array([0,0])
1826        H = refl[:3]
1827        SQ = 1./(2.*refl[4])**2
1828        if not len(refl[-1]):                #no form factors
1829            if 'N' in parmDict[hfx+'Type']:
1830                refl[-1] = getBLvalues(BLtables)
1831            else:       #'X'
1832                refl[-1] = getFFvalues(FFtables,SQ)
1833        for i,El in enumerate(Tdata):
1834            FF[i] = refl[-1][El]           
1835        SQfactor = 4.0*SQ*twopisq
1836        Uniq = refl[11]
1837        phi = refl[12]
1838        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1839        sinp = np.sin(phase)
1840        cosp = np.cos(phase)
1841        occ = Mdata*Fdata/len(Uniq)
1842        biso = -SQfactor*Uisodata
1843        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1844        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1845        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1846        Tcorr = Tiso*Tuij
1847        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
1848        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
1849        if not SGData['SGInv']:
1850            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
1851            fbs = np.sum(np.sum(fb,axis=1),axis=1)
1852        fasq = fas**2
1853        fbsq = fbs**2        #imaginary
1854        refl[9] = np.sum(fasq)+np.sum(fbsq)
1855        refl[10] = atan2d(fbs[0],fas[0])
1856    return refList
1857   
1858def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1859    twopi = 2.0*np.pi
1860    twopisq = 2.0*np.pi**2
1861    ast = np.sqrt(np.diag(G))
1862    Mast = twopisq*np.multiply.outer(ast,ast)
1863    FFtables = calcControls['FFtables']
1864    BLtables = calcControls['BLtables']
1865    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
1866    FF = np.zeros(len(Tdata))
1867    if 'N' in parmDict[hfx+'Type']:
1868        FP = 0.
1869        FPP = 0.
1870    else:
1871        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
1872        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
1873    maxPos = len(SGData['SGOps'])       
1874    Uij = np.array(G2lat.U6toUij(Uijdata))
1875    bij = Mast*Uij.T
1876    dFdvDict = {}
1877    dFdfr = np.zeros((len(refList),len(Mdata)))
1878    dFdx = np.zeros((len(refList),len(Mdata),3))
1879    dFdui = np.zeros((len(refList),len(Mdata)))
1880    dFdua = np.zeros((len(refList),len(Mdata),6))
1881    for iref,refl in enumerate(refList):
1882        H = np.array(refl[:3])
1883        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
1884        for i,El in enumerate(Tdata):           
1885            FF[i] = refl[-1][El]           
1886        SQfactor = 8.0*SQ*np.pi**2
1887        Uniq = refl[11]
1888        phi = refl[12]
1889        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1890        sinp = np.sin(phase)
1891        cosp = np.cos(phase)
1892        occ = Mdata*Fdata/len(Uniq)
1893        biso = -SQfactor*Uisodata
1894        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1895#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1896        HbH = -np.inner(H,np.inner(bij,H))
1897        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1898        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1899        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1900        Tcorr = Tiso*Tuij
1901        fot = (FF+FP)*occ*Tcorr
1902        fotp = FPP*occ*Tcorr
1903        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1904        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1905       
1906        fas = np.sum(np.sum(fa,axis=1),axis=1)
1907        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1908        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1909        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1910        #sum below is over Uniq
1911        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1912        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1913        dfadui = np.sum(-SQfactor*fa,axis=2)
1914        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1915        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
1916        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1917        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1918        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1919        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1920        if not SGData['SGInv']:
1921            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
1922            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
1923            dfbdui = np.sum(-SQfactor*fb,axis=2)
1924            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1925            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1926            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1927            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1928            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1929        #loop over atoms - each dict entry is list of derivatives for all the reflections
1930    for i in range(len(Mdata)):     
1931        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1932        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1933        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1934        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1935        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1936        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1937        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1938        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1939        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1940        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1941        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1942    return dFdvDict
1943       
1944def Dict2Values(parmdict, varylist):
1945    '''Use before call to leastsq to setup list of values for the parameters
1946    in parmdict, as selected by key in varylist'''
1947    return [parmdict[key] for key in varylist] 
1948   
1949def Values2Dict(parmdict, varylist, values):
1950    ''' Use after call to leastsq to update the parameter dictionary with
1951    values corresponding to keys in varylist'''
1952    parmdict.update(zip(varylist,values))
1953   
1954def GetNewCellParms(parmDict,varyList):
1955    newCellDict = {}
1956    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
1957    for item in varyList:
1958        keys = item.split(':')
1959        if keys[2] in Ddict:
1960            key = keys[0]+'::'+Ddict[keys[2]]
1961            parm = keys[0]+'::'+keys[2]
1962            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1963    return newCellDict
1964   
1965def ApplyXYZshifts(parmDict,varyList):
1966    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1967        input:
1968            parmDict - parameter dictionary
1969            varyList - list of variables
1970        returns:
1971            newAtomDict - dictitemionary of new atomic coordinate names & values;
1972                key is parameter shift name
1973    '''
1974    newAtomDict = {}
1975    for item in parmDict:
1976        if 'dA' in item:
1977            parm = ''.join(item.split('d'))
1978            parmDict[parm] += parmDict[item]
1979            newAtomDict[item] = [parm,parmDict[parm]]
1980    return newAtomDict
1981   
1982def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1983    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1984    odfCor = 1.0
1985    H = refl[:3]
1986    cell = G2lat.Gmat2cell(g)
1987    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1988    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1989    phi,beta = G2lat.CrsAng(H,cell,SGData)
1990    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1991    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1992    for item in SHnames:
1993        L,M,N = eval(item.strip('C'))
1994        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1995        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1996        Lnorm = G2lat.Lnorm(L)
1997        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1998    return odfCor
1999   
2000def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
2001    FORPI = 12.5663706143592
2002    IFCoup = 'Bragg' in calcControls[hfx+'instType']
2003    odfCor = 1.0
2004    dFdODF = {}
2005    dFdSA = [0,0,0]
2006    H = refl[:3]
2007    cell = G2lat.Gmat2cell(g)
2008    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
2009    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
2010    phi,beta = G2lat.CrsAng(H,cell,SGData)
2011    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
2012    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
2013    for item in SHnames:
2014        L,M,N = eval(item.strip('C'))
2015        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2016        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
2017        Lnorm = G2lat.Lnorm(L)
2018        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
2019        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
2020        for i in range(3):
2021            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
2022    return odfCor,dFdODF,dFdSA
2023   
2024def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2025    odfCor = 1.0
2026    H = refl[:3]
2027    cell = G2lat.Gmat2cell(g)
2028    Sangl = [0.,0.,0.]
2029    if 'Bragg' in calcControls[hfx+'instType']:
2030        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2031        IFCoup = True
2032    else:
2033        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2034        IFCoup = False
2035    phi,beta = G2lat.CrsAng(H,cell,SGData)
2036    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2037    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2038    for item in SHnames:
2039        L,N = eval(item.strip('C'))
2040        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
2041        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
2042    return odfCor
2043   
2044def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
2045    FORPI = 12.5663706143592
2046    odfCor = 1.0
2047    dFdODF = {}
2048    H = refl[:3]
2049    cell = G2lat.Gmat2cell(g)
2050    Sangl = [0.,0.,0.]
2051    if 'Bragg' in calcControls[hfx+'instType']:
2052        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
2053        IFCoup = True
2054    else:
2055        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
2056        IFCoup = False
2057    phi,beta = G2lat.CrsAng(H,cell,SGData)
2058    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
2059    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
2060    for item in SHnames:
2061        L,N = eval(item.strip('C'))
2062        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
2063        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
2064        dFdODF[phfx+item] = Kcsl*Lnorm
2065    return odfCor,dFdODF
2066   
2067def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2068    POcorr = 1.0
2069    if calcControls[phfx+'poType'] == 'MD':
2070        MD = parmDict[phfx+'MD']
2071        if MD != 1.0:
2072            MDAxis = calcControls[phfx+'MDAxis']
2073            sumMD = 0
2074            for H in refl[11]:           
2075                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2076                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2077                sumMD += A**3
2078            POcorr = sumMD/len(refl[11])
2079    else:   #spherical harmonics
2080        if calcControls[pfx+'SHord']:
2081            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2082    return POcorr
2083   
2084def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
2085    POcorr = 1.0
2086    POderv = {}
2087    if calcControls[phfx+'poType'] == 'MD':
2088        MD = parmDict[phfx+'MD']
2089        MDAxis = calcControls[phfx+'MDAxis']
2090        sumMD = 0
2091        sumdMD = 0
2092        for H in refl[11]:           
2093            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
2094            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
2095            sumMD += A**3
2096            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
2097        POcorr = sumMD/len(refl[11])
2098        POderv[phfx+'MD'] = sumdMD/len(refl[11])
2099    else:   #spherical harmonics
2100        if calcControls[pfx+'SHord']:
2101            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
2102    return POcorr,POderv
2103   
2104def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2105    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
2106    if 'X' in parmDict[hfx+'Type']:
2107        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
2108    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2109    if pfx+'SHorder' in parmDict:
2110        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2111    refl[13] = Icorr       
2112   
2113def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2114    dIdsh = 1./parmDict[hfx+'Scale']
2115    dIdsp = 1./parmDict[phfx+'Scale']
2116    if 'X' in parmDict[hfx+'Type']:
2117        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2118        dIdPola /= pola
2119    else:       #'N'
2120        dIdPola = 0.0
2121    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2122    for iPO in dIdPO:
2123        dIdPO[iPO] /= POcorr
2124    dFdODF = {}
2125    dFdSA = [0,0,0]
2126    if pfx+'SHorder' in parmDict:
2127        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2128        for iSH in dFdODF:
2129            dFdODF[iSH] /= odfCor
2130        for i in range(3):
2131            dFdSA[i] /= odfCor
2132    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
2133       
2134def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2135    costh = cosd(refl[5]/2.)
2136    #crystallite size
2137    if calcControls[phfx+'SizeType'] == 'isotropic':
2138        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
2139    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2140        H = np.array(refl[:3])
2141        P = np.array(calcControls[phfx+'SizeAxis'])
2142        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2143        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
2144        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
2145    else:           #ellipsoidal crystallites
2146        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2147        H = np.array(refl[:3])
2148        lenR = G2pwd.ellipseSize(H,Sij,GB)
2149        Sgam = 1.8*wave/(np.pi*costh*lenR)
2150    #microstrain               
2151    if calcControls[phfx+'MustrainType'] == 'isotropic':
2152        Mgam = 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
2153    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2154        H = np.array(refl[:3])
2155        P = np.array(calcControls[phfx+'MustrainAxis'])
2156        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2157        Si = parmDict[phfx+'Mustrain:i']
2158        Sa = parmDict[phfx+'Mustrain:a']
2159        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2160    else:       #generalized - P.W. Stephens model
2161        pwrs = calcControls[phfx+'MuPwrs']
2162        sum = 0
2163        for i,pwr in enumerate(pwrs):
2164            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2165        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
2166    gam = Sgam*parmDict[phfx+'Size:mx']+Mgam*parmDict[phfx+'Mustrain:mx']
2167    sig = (Sgam*(1.-parmDict[phfx+'Size:mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain:mx']))**2
2168    sig /= ateln2
2169    return sig,gam
2170       
2171def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2172    gamDict = {}
2173    sigDict = {}
2174    costh = cosd(refl[5]/2.)
2175    tanth = tand(refl[5]/2.)
2176    #crystallite size derivatives
2177    if calcControls[phfx+'SizeType'] == 'isotropic':
2178        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
2179        gamDict[phfx+'Size:i'] = -1.80*wave*parmDict[phfx+'Size:mx']/(np.pi*costh)
2180        sigDict[phfx+'Size:i'] = -3.60*Sgam*wave*(1.-parmDict[phfx+'Size:mx'])**2/(np.pi*costh*ateln2)
2181    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2182        H = np.array(refl[:3])
2183        P = np.array(calcControls[phfx+'SizeAxis'])
2184        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2185        Si = parmDict[phfx+'Size:i']
2186        Sa = parmDict[phfx+'Size:a']
2187        gami = (1.8*wave/np.pi)/(Si*Sa)
2188        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2189        Sgam = gami*sqtrm
2190        gam = Sgam/costh
2191        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
2192        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
2193        gamDict[phfx+'Size:i'] = dsi*parmDict[phfx+'Size:mx']
2194        gamDict[phfx+'Size:a'] = dsa*parmDict[phfx+'Size:mx']
2195        sigDict[phfx+'Size:i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
2196        sigDict[phfx+'Size:a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
2197    else:           #ellipsoidal crystallites
2198        const = 1.8*wave/(np.pi*costh)
2199        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2200        H = np.array(refl[:3])
2201        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2202        Sgam = 1.8*wave/(np.pi*costh*lenR)
2203        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2204            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size:mx']
2205            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
2206    gamDict[phfx+'Size:mx'] = Sgam
2207    sigDict[phfx+'Size:mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size:mx'])/ateln2
2208           
2209    #microstrain derivatives               
2210    if calcControls[phfx+'MustrainType'] == 'isotropic':
2211        Mgam = 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
2212        gamDict[phfx+'Mustrain:i'] =  0.018*tanth*parmDict[phfx+'Mustrain:mx']/np.pi
2213        sigDict[phfx+'Mustrain:i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain:mx'])**2/(np.pi*ateln2)       
2214    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2215        H = np.array(refl[:3])
2216        P = np.array(calcControls[phfx+'MustrainAxis'])
2217        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2218        Si = parmDict[phfx+'Mustrain:i']
2219        Sa = parmDict[phfx+'Mustrain:a']
2220        gami = 0.018*Si*Sa*tanth/np.pi
2221        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2222        Mgam = gami/sqtrm
2223        dsi = -gami*Si*cosP**2/sqtrm**3
2224        dsa = -gami*Sa*sinP**2/sqtrm**3
2225        gamDict[phfx+'Mustrain:i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain:mx']
2226        gamDict[phfx+'Mustrain:a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain:mx']
2227        sigDict[phfx+'Mustrain:i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2
2228        sigDict[phfx+'Mustrain:a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2       
2229    else:       #generalized - P.W. Stephens model
2230        pwrs = calcControls[phfx+'MuPwrs']
2231        const = 0.018*refl[4]**2*tanth
2232        sum = 0
2233        for i,pwr in enumerate(pwrs):
2234            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2235            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
2236            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain:mx']
2237            sigDict[phfx+'Mustrain:'+str(i)] = \
2238                2.*const*term*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2
2239        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
2240        for i in range(len(pwrs)):
2241            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
2242    gamDict[phfx+'Mustrain:mx'] = Mgam
2243    sigDict[phfx+'Mustrain:mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain:mx'])/ateln2
2244    return sigDict,gamDict
2245       
2246def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2247    h,k,l = refl[:3]
2248    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2249    d = np.sqrt(dsq)
2250
2251    refl[4] = d
2252    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2253    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2254    if 'Bragg' in calcControls[hfx+'instType']:
2255        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2256            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2257    else:               #Debye-Scherrer - simple but maybe not right
2258        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2259    return pos
2260
2261def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2262    dpr = 180./np.pi
2263    h,k,l = refl[:3]
2264    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2265    dst = np.sqrt(dstsq)
2266    pos = refl[5]-parmDict[hfx+'Zero']
2267    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2268    dpdw = const*dst
2269    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
2270    dpdA *= const*wave/(2.0*dst)
2271    dpdZ = 1.0
2272    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2273    if 'Bragg' in calcControls[hfx+'instType']:
2274        dpdSh = -4.*const*cosd(pos/2.0)
2275        dpdTr = -const*sind(pos)*100.0
2276        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2277    else:               #Debye-Scherrer - simple but maybe not right
2278        dpdXd = -const*cosd(pos)
2279        dpdYd = -const*sind(pos)
2280        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2281           
2282def GetHStrainShift(refl,SGData,phfx,parmDict):
2283    laue = SGData['SGLaue']
2284    uniq = SGData['SGUniq']
2285    h,k,l = refl[:3]
2286    if laue in ['m3','m3m']:
2287        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2288            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2289    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2290        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2291    elif laue in ['3R','3mR']:
2292        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2293    elif laue in ['4/m','4/mmm']:
2294        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2295    elif laue in ['mmm']:
2296        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2297    elif laue in ['2/m']:
2298        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2299        if uniq == 'a':
2300            Dij += parmDict[phfx+'D23']*k*l
2301        elif uniq == 'b':
2302            Dij += parmDict[phfx+'D13']*h*l
2303        elif uniq == 'c':
2304            Dij += parmDict[phfx+'D12']*h*k
2305    else:
2306        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2307            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2308    return Dij*refl[4]**2*tand(refl[5]/2.0)
2309           
2310def GetHStrainShiftDerv(refl,SGData,phfx):
2311    laue = SGData['SGLaue']
2312    uniq = SGData['SGUniq']
2313    h,k,l = refl[:3]
2314    if laue in ['m3','m3m']:
2315        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2316            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2317    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2318        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2319    elif laue in ['3R','3mR']:
2320        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2321    elif laue in ['4/m','4/mmm']:
2322        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2323    elif laue in ['mmm']:
2324        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2325    elif laue in ['2/m']:
2326        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2327        if uniq == 'a':
2328            dDijDict[phfx+'D23'] = k*l
2329        elif uniq == 'b':
2330            dDijDict[phfx+'D13'] = h*l
2331        elif uniq == 'c':
2332            dDijDict[phfx+'D12'] = h*k
2333            names.append()
2334    else:
2335        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2336            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2337    for item in dDijDict:
2338        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
2339    return dDijDict
2340   
2341def GetFprime(controlDict,Histograms):
2342    FFtables = controlDict['FFtables']
2343    if not FFtables:
2344        return
2345    histoList = Histograms.keys()
2346    histoList.sort()
2347    for histogram in histoList:
2348        if 'PWDR' in histogram[:4]:
2349            Histogram = Histograms[histogram]
2350            hId = Histogram['hId']
2351            hfx = ':%d:'%(hId)
2352            keV = controlDict[hfx+'keV']
2353            for El in FFtables:
2354                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2355                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2356                FFtables[El][hfx+'FP'] = FP
2357                FFtables[El][hfx+'FPP'] = FPP               
2358           
2359def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2360    histoList = Histograms.keys()
2361    histoList.sort()
2362    for histogram in histoList:
2363        if 'PWDR' in histogram[:4]:
2364            Histogram = Histograms[histogram]
2365            hId = Histogram['hId']
2366            hfx = ':%d:'%(hId)
2367            Limits = calcControls[hfx+'Limits']
2368            shl = max(parmDict[hfx+'SH/L'],0.002)
2369            Ka2 = False
2370            kRatio = 0.0
2371            if hfx+'Lam1' in parmDict.keys():
2372                Ka2 = True
2373                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2374                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2375            x,y,w,yc,yb,yd = Histogram['Data']
2376            xB = np.searchsorted(x,Limits[0])
2377            xF = np.searchsorted(x,Limits[1])
2378            ymb = np.array(y-yb)
2379            ycmb = np.array(yc-yb)
2380            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
2381            refLists = Histogram['Reflection Lists']
2382            for phase in refLists:
2383                Phase = Phases[phase]
2384                pId = Phase['pId']
2385                phfx = '%d:%d:'%(pId,hId)
2386                refList = refLists[phase]
2387                sumFo = 0.0
2388                sumdF = 0.0
2389                sumFosq = 0.0
2390                sumdFsq = 0.0
2391                for refl in refList:
2392                    if 'C' in calcControls[hfx+'histType']:
2393                        yp = np.zeros_like(yb)
2394                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2395                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
2396                        iFin = min(np.searchsorted(x,refl[5]+fmax),xF)
2397                        iFin2 = iFin
2398                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2399                        if Ka2:
2400                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2401                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2402                            iBeg2 = max(xB,np.searchsorted(x[xB:xF],pos2-fmin))
2403                            iFin2 = min(np.searchsorted(x[xB:xF],pos2+fmax),xF)
2404                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
2405                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
2406                    elif 'T' in calcControls[hfx+'histType']:
2407                        print 'TOF Undefined at present'
2408                        raise Exception    #no TOF yet
2409                    Fo = np.sqrt(np.abs(refl[8]))
2410                    Fc = np.sqrt(np.abs(refl[9]))
2411                    sumFo += Fo
2412                    sumFosq += refl[8]**2
2413                    sumdF += np.abs(Fo-Fc)
2414                    sumdFsq += (refl[8]-refl[9])**2
2415                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2416                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2417                Histogram[phfx+'Nref'] = len(refList)
2418               
2419def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2420   
2421    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
2422        U = parmDict[hfx+'U']
2423        V = parmDict[hfx+'V']
2424        W = parmDict[hfx+'W']
2425        X = parmDict[hfx+'X']
2426        Y = parmDict[hfx+'Y']
2427        tanPos = tand(refl[5]/2.0)
2428        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
2429        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
2430        sig = max(0.001,sig)
2431        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
2432        gam = max(0.001,gam)
2433        return sig,gam
2434               
2435    hId = Histogram['hId']
2436    hfx = ':%d:'%(hId)
2437    bakType = calcControls[hfx+'bakType']
2438    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2439    yc = np.zeros_like(yb)
2440       
2441    if 'C' in calcControls[hfx+'histType']:   
2442        shl = max(parmDict[hfx+'SH/L'],0.002)
2443        Ka2 = False
2444        if hfx+'Lam1' in parmDict.keys():
2445            wave = parmDict[hfx+'Lam1']
2446            Ka2 = True
2447            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2448            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2449        else:
2450            wave = parmDict[hfx+'Lam']
2451    else:
2452        print 'TOF Undefined at present'
2453        raise ValueError
2454    for phase in Histogram['Reflection Lists']:
2455        refList = Histogram['Reflection Lists'][phase]
2456        Phase = Phases[phase]
2457        pId = Phase['pId']
2458        pfx = '%d::'%(pId)
2459        phfx = '%d:%d:'%(pId,hId)
2460        hfx = ':%d:'%(hId)
2461        SGData = Phase['General']['SGData']
2462        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2463        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2464        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2465        Vst = np.sqrt(nl.det(G))    #V*
2466        if not Phase['General'].get('doPawley'):
2467            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2468        for refl in refList:
2469            if 'C' in calcControls[hfx+'histType']:
2470                h,k,l = refl[:3]
2471                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
2472                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
2473                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
2474                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
2475                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
2476                refl[13] *= Vst*Lorenz
2477                if Phase['General'].get('doPawley'):
2478                    try:
2479                        refl[9] = parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
2480                    except KeyError:
2481#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2482                        continue
2483                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2484                iBeg = np.searchsorted(x,refl[5]-fmin)
2485                iFin = np.searchsorted(x,refl[5]+fmax)
2486                if not iBeg+iFin:       #peak below low limit - skip peak
2487                    continue
2488                elif not iBeg-iFin:     #peak above high limit - done
2489                    break
2490                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2491                if Ka2:
2492                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2493                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2494                    iBeg = np.searchsorted(x,pos2-fmin)
2495                    iFin = np.searchsorted(x,pos2+fmax)
2496                    if not iBeg+iFin:       #peak below low limit - skip peak
2497                        continue
2498                    elif not iBeg-iFin:     #peak above high limit - done
2499                        return yc,yb
2500                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
2501            elif 'T' in calcControls[hfx+'histType']:
2502                print 'TOF Undefined at present'
2503                raise Exception    #no TOF yet
2504    return yc,yb
2505   
2506def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
2507   
2508    def cellVaryDerv(pfx,SGData,dpdA): 
2509        if SGData['SGLaue'] in ['-1',]:
2510            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2511                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2512        elif SGData['SGLaue'] in ['2/m',]:
2513            if SGData['SGUniq'] == 'a':
2514                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2515            elif SGData['SGUniq'] == 'b':
2516                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2517            else:
2518                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2519        elif SGData['SGLaue'] in ['mmm',]:
2520            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2521        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2522#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2523            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2524        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2525#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2526            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2527        elif SGData['SGLaue'] in ['3R', '3mR']:
2528            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2529        elif SGData['SGLaue'] in ['m3m','m3']:
2530#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2531            return [[pfx+'A0',dpdA[0]]]
2532    # create a list of dependent variables and set up a dictionary to hold their derivatives
2533    dependentVars = G2mv.GetDependentVars()
2534    depDerivDict = {}
2535    for j in dependentVars:
2536        depDerivDict[j] = np.zeros(shape=(len(x)))
2537    #print 'dependent vars',dependentVars
2538    lenX = len(x)               
2539    hId = Histogram['hId']
2540    hfx = ':%d:'%(hId)
2541    bakType = calcControls[hfx+'bakType']
2542    dMdv = np.zeros(shape=(len(varylist),len(x)))
2543    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2544    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2545        bBpos =varylist.index(hfx+'Back:0')
2546        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2547    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2548    for name in varylist:
2549        if 'Debye' in name:
2550            id = int(name.split(':')[-1])
2551            parm = name[:int(name.rindex(':'))]
2552            ip = names.index(parm)
2553            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2554    if 'C' in calcControls[hfx+'histType']:   
2555        dx = x[1]-x[0]
2556        shl = max(parmDict[hfx+'SH/L'],0.002)
2557        Ka2 = False
2558        if hfx+'Lam1' in parmDict.keys():
2559            wave = parmDict[hfx+'Lam1']
2560            Ka2 = True
2561            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2562            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2563        else:
2564            wave = parmDict[hfx+'Lam']
2565    else:
2566        print 'TOF Undefined at present'
2567        raise ValueError
2568    for phase in Histogram['Reflection Lists']:
2569        refList = Histogram['Reflection Lists'][phase]
2570        Phase = Phases[phase]
2571        SGData = Phase['General']['SGData']
2572        pId = Phase['pId']
2573        pfx = '%d::'%(pId)
2574        phfx = '%d:%d:'%(pId,hId)
2575        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2576        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2577        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2578        if not Phase['General'].get('doPawley'):
2579            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2580        for iref,refl in enumerate(refList):
2581            if 'C' in calcControls[hfx+'histType']:        #CW powder
2582                h,k,l = refl[:3]
2583                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2584                if Phase['General'].get('doPawley'):
2585                    try:
2586                        refl[9] = parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]
2587                    except KeyError:
2588#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2589                        continue
2590                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2591                iBeg = np.searchsorted(x,refl[5]-fmin)
2592                iFin = np.searchsorted(x,refl[5]+fmax)
2593                if not iBeg+iFin:       #peak below low limit - skip peak
2594                    continue
2595                elif not iBeg-iFin:     #peak above high limit - done
2596                    break
2597                pos = refl[5]
2598                tanth = tand(pos/2.0)
2599                costh = cosd(pos/2.0)
2600                lenBF = iFin-iBeg
2601                dMdpk = np.zeros(shape=(6,lenBF))
2602                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2603                for i in range(1,5):
2604                    dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2605                dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2606                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2607                if Ka2:
2608                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2609                    kdelt = int((pos2-refl[5])/dx)               
2610                    iBeg2 = min(lenX,iBeg+kdelt)
2611                    iFin2 = min(lenX,iFin+kdelt)
2612                    if iBeg2-iFin2:
2613                        lenBF2 = iFin2-iBeg2
2614                        dMdpk2 = np.zeros(shape=(6,lenBF2))
2615                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2616                        for i in range(1,5):
2617                            dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2618                        dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2619                        dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
2620                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2621                if Phase['General'].get('doPawley'):
2622                    try:
2623                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2624                        dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
2625                        if Ka2:
2626                            dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
2627                        # Assuming Pawley variables not in constraints
2628                    except ValueError:
2629                        pass
2630                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2631                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2632                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2633                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2634                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2635                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2636                    hfx+'DisplaceY':[dpdY,'pos'],}
2637                for name in names:
2638                    item = names[name]
2639                    if name in varylist:
2640                        dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2641                        if Ka2:
2642                            dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2643                    elif name in dependentVars:
2644                        if Ka2:
2645                            depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2646                        depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2647                for iPO in dIdPO:
2648                    if iPO in varylist:
2649                        dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2650                        if Ka2:
2651                            dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2652                    elif iPO in dependentVars:
2653                        depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2654                        if Ka2:
2655                            depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2656                for i,name in enumerate(['omega','chi','phi']):
2657                    aname = pfx+'SH '+name
2658                    if aname in varylist:
2659                        dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2660                        if Ka2:
2661                            dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2662                    elif aname in dependentVars:
2663                        depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2664                        if Ka2:
2665                            depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2666                for iSH in dFdODF:
2667                    if iSH in varylist:
2668                        dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2669                        if Ka2:
2670                            dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2671                    elif iSH in dependentVars:
2672                        depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2673                        if Ka2:
2674                            depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2675                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2676                for name,dpdA in cellDervNames:
2677                    if name in varylist:
2678                        dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2679                        if Ka2:
2680                            dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2681                    elif name in dependentVars:
2682                        depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2683                        if Ka2:
2684                            depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2685                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2686                for name in dDijDict:
2687                    if name in varylist:
2688                        dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2689                        if Ka2:
2690                            dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2691                    elif name in dependentVars:
2692                        depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2693                        if Ka2:
2694                            depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2695                sigDict,gamDict = GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2696                for name in gamDict:
2697                    if name in varylist:
2698                        dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2699                        if Ka2:
2700                            dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2701                    elif name in dependentVars:
2702                        depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2703                        if Ka2:
2704                            depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2705                for name in sigDict:
2706                    if name in varylist:
2707                        dMdv[varylist.index(name)][iBeg:iFin] += sigDict[name]*dervDict['sig']
2708                        if Ka2:
2709                            dMdv[varylist.index(name)][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2710                    elif name in dependentVars:
2711                        depDerivDict[name][iBeg:iFin] += sigDict[name]*dervDict['sig']
2712                        if Ka2:
2713                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
2714                                               
2715            elif 'T' in calcControls[hfx+'histType']:
2716                print 'TOF Undefined at present'
2717                raise Exception    #no TOF yet
2718            #do atom derivatives -  for F,X & U so far             
2719            corr = dervDict['int']/refl[9]
2720            if Ka2:
2721                corr2 = dervDict2['int']/refl[9]
2722            for name in varylist+dependentVars:
2723                try:
2724                    aname = name.split(pfx)[1][:2]
2725                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2726                except IndexError:
2727                    continue
2728                if name in varylist:
2729                    dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2730                    if Ka2:
2731                        dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2732                elif name in dependentVars:
2733                    depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2734                    if Ka2:
2735                        depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2736    # now process derivatives in constraints
2737    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2738    return dMdv
2739
2740def dervRefine(values,HistoPhases,penalties,parmdict,varylist,calcControls,pawleyLookup,dlg):
2741    parmdict.update(zip(varylist,values))
2742    G2mv.Dict2Map(parmdict,varylist)
2743    Histograms,Phases = HistoPhases
2744    nvar = len(varylist)
2745    dMdv = np.empty(0)
2746    histoList = Histograms.keys()
2747    histoList.sort()
2748    for histogram in histoList:
2749        if 'PWDR' in histogram[:4]:
2750            Histogram = Histograms[histogram]
2751            hId = Histogram['hId']
2752            hfx = ':%d:'%(hId)
2753            Limits = calcControls[hfx+'Limits']
2754            x,y,w,yc,yb,yd = Histogram['Data']
2755            xB = np.searchsorted(x,Limits[0])
2756            xF = np.searchsorted(x,Limits[1])
2757            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2758                varylist,Histogram,Phases,calcControls,pawleyLookup)
2759            if len(dMdv):
2760                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2761            else:
2762                dMdv = dMdvh
2763    if penalties:
2764        dmdv = np.concatenate((dmdv.T,penaltyDeriv(penalties).T)).T
2765    return dMdv
2766
2767#def ComputePowderHessian(args):
2768#    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
2769#    hId = Histogram['hId']
2770#    hfx = ':%d:'%(hId)
2771#    Limits = calcControls[hfx+'Limits']
2772#    x,y,w,yc,yb,yd = Histogram['Data']
2773#    dy = y-yc
2774#    xB = np.searchsorted(x,Limits[0])
2775#    xF = np.searchsorted(x,Limits[1])
2776#    dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(
2777#        parmdict,x[xB:xF],
2778#        varylist,Histogram,Phases,calcControls,pawleyLookup)
2779#    Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2780#    Hess = np.inner(dMdvh,dMdvh)
2781#    return Vec,Hess
2782#
2783def HessRefine(values,HistoPhases,penalties,parmdict,varylist,calcControls,pawleyLookup,dlg):
2784    parmdict.update(zip(varylist,values))
2785    G2mv.Dict2Map(parmdict,varylist)
2786    Histograms,Phases = HistoPhases
2787    nvar = len(varylist)
2788    Hess = np.empty(0)
2789    histoList = Histograms.keys()
2790    histoList.sort()
2791    for histogram in histoList:
2792        if 'PWDR' in histogram[:4]:
2793            Histogram = Histograms[histogram]
2794            hId = Histogram['hId']
2795            hfx = ':%d:'%(hId)
2796            Limits = calcControls[hfx+'Limits']
2797            x,y,w,yc,yb,yd = Histogram['Data']
2798            dy = y-yc
2799            xB = np.searchsorted(x,Limits[0])
2800            xF = np.searchsorted(x,Limits[1])
2801            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2802                varylist,Histogram,Phases,calcControls,pawleyLookup)
2803            if dlg:
2804                dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2805            if len(Hess):
2806                Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2807                Hess += np.inner(dMdvh,dMdvh)
2808            else:
2809                Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2810                Hess = np.inner(dMdvh,dMdvh)
2811    return Vec,Hess
2812
2813#def ComputePowderProfile(args):
2814#    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
2815#    hId = Histogram['hId']
2816#    hfx = ':%d:'%(hId)
2817#    x,y,w,yc,yb,yd = Histogram['Data']
2818#    Limits = calcControls[hfx+'Limits']
2819#    xB = np.searchsorted(x,Limits[0])
2820#    xF = np.searchsorted(x,Limits[1])
2821#    yc,yb = getPowderProfile(parmdict,x[xB:xF],varylist,Histogram,Phases,
2822#        calcControls,pawleyLookup)
2823#    return xB,xF,yc,yb,Histogram['Reflection Lists']
2824#
2825def errRefine(values,HistoPhases,penalties,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2826    parmdict.update(zip(varylist,values))
2827    Values2Dict(parmdict, varylist, values)
2828    G2mv.Dict2Map(parmdict,varylist)
2829    Histograms,Phases = HistoPhases
2830    M = np.empty(0)
2831    sumwYo = 0
2832    Nobs = 0
2833    histoList = Histograms.keys()
2834    histoList.sort()
2835    for histogram in histoList:
2836        if 'PWDR' in histogram[:4]:
2837            Histogram = Histograms[histogram]
2838            hId = Histogram['hId']
2839            hfx = ':%d:'%(hId)
2840            Limits = calcControls[hfx+'Limits']
2841            x,y,w,yc,yb,yd = Histogram['Data']
2842            yc *= 0.0                           #zero full calcd profiles
2843            yb *= 0.0
2844            yd *= 0.0
2845            xB = np.searchsorted(x,Limits[0])
2846            xF = np.searchsorted(x,Limits[1])
2847            Histogram['Nobs'] = xF-xB
2848            Nobs += Histogram['Nobs']
2849            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2850            sumwYo += Histogram['sumwYo']
2851            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2852                varylist,Histogram,Phases,calcControls,pawleyLookup)
2853            yc[xB:xF] += yb[xB:xF]
2854            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2855            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2856            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2857            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2858            if dlg:
2859                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2860            M = np.concatenate((M,wdy))
2861    Histograms['sumwYo'] = sumwYo
2862    Histograms['Nobs'] = Nobs
2863    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2864    if dlg:
2865        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2866        if not GoOn:
2867            parmDict['saved values'] = values
2868            raise Exception         #Abort!!
2869    if penalties:
2870        M = np.concatenate((M,penaltyFxn(penalties)))
2871    return M
2872                       
2873def Refine(GPXfile,dlg):
2874    import pytexture as ptx
2875    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2876   
2877    ShowBanner()
2878    varyList = []
2879    parmDict = {}
2880    G2mv.InitVars()   
2881    Controls = GetControls(GPXfile)
2882    ShowControls(Controls)
2883    calcControls = {}
2884    calcControls.update(Controls)           
2885    constrDict,fixedList = GetConstraints(GPXfile)
2886    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2887    if not Phases:
2888        print ' *** ERROR - you have no phases! ***'
2889        print ' *** Refine aborted ***'
2890        raise Exception
2891    if not Histograms:
2892        print ' *** ERROR - you have no data to refine with! ***'
2893        print ' *** Refine aborted ***'
2894        raise Exception       
2895    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2896    calcControls['Natoms'] = Natoms
2897    calcControls['FFtables'] = FFtables
2898    calcControls['BLtables'] = BLtables
2899    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2900    calcControls.update(controlDict)
2901    histVary,histDict,controlDict = GetHistogramData(Histograms)
2902    calcControls.update(controlDict)
2903    varyList = phaseVary+hapVary+histVary
2904    parmDict.update(phaseDict)
2905    parmDict.update(hapDict)
2906    parmDict.update(histDict)
2907    GetFprime(calcControls,Histograms)
2908    # do constraint processing
2909    try:
2910        groups,parmlist = G2mv.GroupConstraints(constrDict)
2911        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
2912    except:
2913        print ' *** ERROR - your constraints are internally inconsistent ***'
2914        # traceback for debug
2915        #print 'varyList',varyList
2916        #print 'constrDict',constrDict
2917        #print 'fixedList',fixedList
2918        #import traceback
2919        #print traceback.format_exc()
2920        raise Exception(' *** Refine aborted ***')
2921    # # check to see which generated parameters are fully varied
2922    # msg = G2mv.SetVaryFlags(varyList)
2923    # if msg:
2924    #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
2925    #     print msg
2926    #     raise Exception(' *** Refine aborted ***')
2927    #print G2mv.VarRemapShow(varyList)
2928    G2mv.Map2Dict(parmDict,varyList)
2929    Rvals = {}
2930    while True:
2931        begin = time.time()
2932        values =  np.array(Dict2Values(parmDict, varyList))
2933        penalties = getPenalties()
2934        Ftol = Controls['min dM/M']
2935        Factor = Controls['shift factor']
2936        maxCyc = Controls['max cyc']
2937        if 'Jacobian' in Controls['deriv type']:           
2938            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2939                ftol=Ftol,col_deriv=True,factor=Factor,
2940                args=([Histograms,Phases],penalties,parmDict,varyList,calcControls,pawleyLookup,dlg))
2941            ncyc = int(result[2]['nfev']/2)
2942        elif 'Hessian' in Controls['deriv type']:
2943            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
2944                args=([Histograms,Phases],penalties,parmDict,varyList,calcControls,pawleyLookup,dlg))
2945            ncyc = result[2]['num cyc']+1
2946            Rvals['lamMax'] = result[2]['lamMax']                           
2947        else:           #'numeric'
2948            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2949                args=([Histograms,Phases],penalties,parmDict,varyList,calcControls,pawleyLookup,dlg))
2950            ncyc = int(result[2]['nfev']/len(varyList))
2951#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2952#        for item in table: print item,table[item]               #useful debug - are things shifting?
2953        runtime = time.time()-begin
2954        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
2955        Values2Dict(parmDict, varyList, result[0])
2956        G2mv.Dict2Map(parmDict,varyList)
2957       
2958        Rvals['Nobs'] = Histograms['Nobs']
2959        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
2960        Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
2961        print '\n Refinement results:'
2962        print 135*'-'
2963        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2964        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
2965        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
2966        try:
2967            covMatrix = result[1]*Rvals['GOF']
2968            sig = np.sqrt(np.diag(covMatrix))
2969            if np.any(np.isnan(sig)):
2970                print '*** Least squares aborted - some invalid esds possible ***'
2971#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2972#            for item in table: print item,table[item]               #useful debug - are things shifting?
2973            break                   #refinement succeeded - finish up!
2974        except TypeError:          #result[1] is None on singular matrix
2975            print '**** Refinement failed - singular matrix ****'
2976            if 'Hessian' in Controls['deriv type']:
2977                num = len(varyList)-1
2978                for i,val in enumerate(np.flipud(result[2]['psing'])):
2979                    if val:
2980                        print 'Removing parameter: ',varyList[num-i]
2981                        del(varyList[num-i])                   
2982            else:
2983                Ipvt = result[2]['ipvt']
2984                for i,ipvt in enumerate(Ipvt):
2985                    if not np.sum(result[2]['fjac'],axis=1)[i]:
2986                        print 'Removing parameter: ',varyList[ipvt-1]
2987                        del(varyList[ipvt-1])
2988                        break
2989
2990#    print 'dependentParmList: ',G2mv.dependentParmList
2991#    print 'arrayList: ',G2mv.arrayList
2992#    print 'invarrayList: ',G2mv.invarrayList
2993#    print 'indParmList: ',G2mv.indParmList
2994#    print 'fixedDict: ',G2mv.fixedDict
2995#    print 'test1'
2996    GetFobsSq(Histograms,Phases,parmDict,calcControls)
2997#    print 'test2'
2998    sigDict = dict(zip(varyList,sig))
2999    newCellDict = GetNewCellParms(parmDict,varyList)
3000    newAtomDict = ApplyXYZshifts(parmDict,varyList)
3001    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3002        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3003    # add the uncertainties into the esd dictionary (sigDict)
3004    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3005    SetPhaseData(parmDict,sigDict,Phases,covData)
3006    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
3007    SetHistogramData(parmDict,sigDict,Histograms)
3008    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
3009    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3010   
3011#for testing purposes!!!
3012    if DEBUG:
3013        import cPickle
3014        fl = open('structTestdata.dat','wb')
3015        cPickle.dump(parmDict,fl,1)
3016        cPickle.dump(varyList,fl,1)
3017        for histogram in Histograms:
3018            if 'PWDR' in histogram[:4]:
3019                Histogram = Histograms[histogram]
3020        cPickle.dump(Histogram,fl,1)
3021        cPickle.dump(Phases,fl,1)
3022        cPickle.dump(calcControls,fl,1)
3023        cPickle.dump(pawleyLookup,fl,1)
3024        fl.close()
3025
3026    if dlg:
3027        return Rvals['Rwp']
3028
3029def SeqRefine(GPXfile,dlg):
3030    import pytexture as ptx
3031    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
3032   
3033    ShowBanner()
3034    print ' Sequential Refinement'
3035    G2mv.InitVars()   
3036    Controls = GetControls(GPXfile)
3037    ShowControls(Controls)           
3038    constrDict,fixedList = GetConstraints(GPXfile)
3039    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3040    if not Phases:
3041        print ' *** ERROR - you have no histograms to refine! ***'
3042        print ' *** Refine aborted ***'
3043        raise Exception
3044    if not Histograms:
3045        print ' *** ERROR - you have no data to refine with! ***'
3046        print ' *** Refine aborted ***'
3047        raise Exception
3048    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
3049    if 'Seq Data' in Controls:
3050        histNames = Controls['Seq Data']
3051    else:
3052        histNames = GetHistogramNames(GPXfile,['PWDR',])
3053    if 'Reverse Seq' in Controls:
3054        if Controls['Reverse Seq']:
3055            histNames.reverse()
3056    SeqResult = {'histNames':histNames}
3057    makeBack = True
3058    for ihst,histogram in enumerate(histNames):
3059        ifPrint = False
3060        if dlg:
3061            dlg.SetTitle('Residual for histogram '+str(ihst))
3062        calcControls = {}
3063        calcControls['Natoms'] = Natoms
3064        calcControls['FFtables'] = FFtables
3065        calcControls['BLtables'] = BLtables
3066        varyList = []
3067        parmDict = {}
3068        Histo = {histogram:Histograms[histogram],}
3069        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3070        calcControls.update(controlDict)
3071        histVary,histDict,controlDict = GetHistogramData(Histo,False)
3072        calcControls.update(controlDict)
3073        varyList = phaseVary+hapVary+histVary
3074        if not ihst:
3075            saveVaryList = varyList[:]
3076            for i,item in enumerate(saveVaryList):
3077                items = item.split(':')
3078                if items[1]:
3079                    items[1] = ''
3080                item = ':'.join(items)
3081                saveVaryList[i] = item
3082            SeqResult['varyList'] = saveVaryList
3083        else:
3084            newVaryList = varyList[:]
3085            for i,item in enumerate(newVaryList):
3086                items = item.split(':')
3087                if items[1]:
3088                    items[1] = ''
3089                item = ':'.join(items)
3090                newVaryList[i] = item
3091            if newVaryList != SeqResult['varyList']:
3092                print newVaryList
3093                print SeqResult['varyList']
3094                print '**** ERROR - variable list for this histogram does not match previous'
3095                raise Exception
3096        parmDict.update(phaseDict)
3097        parmDict.update(hapDict)
3098        parmDict.update(histDict)
3099        GetFprime(calcControls,Histo)
3100        # do constraint processing
3101        try:
3102            groups,parmlist = G2mv.GroupConstraints(constrDict)
3103            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
3104        except:
3105            print ' *** ERROR - your constraints are internally inconsistent ***'
3106            raise Exception(' *** Refine aborted ***')
3107        # check to see which generated parameters are fully varied
3108        # msg = G2mv.SetVaryFlags(varyList)
3109        # if msg:
3110        #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
3111        #     print msg
3112        #     raise Exception(' *** Refine aborted ***')
3113        #print G2mv.VarRemapShow(varyList)
3114        G2mv.Map2Dict(parmDict,varyList)
3115        Rvals = {}
3116        while True:
3117            begin = time.time()
3118            values =  np.array(Dict2Values(parmDict,varyList))
3119            penalties = getPenalties()
3120            Ftol = Controls['min dM/M']
3121            Factor = Controls['shift factor']
3122            maxCyc = Controls['max cyc']
3123
3124            if 'Jacobian' in Controls['deriv type']:           
3125                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3126                    ftol=Ftol,col_deriv=True,factor=Factor,
3127                    args=([Histo,Phases],penalties,parmDict,varyList,calcControls,pawleyLookup,dlg))
3128                ncyc = int(result[2]['nfev']/2)
3129            elif 'Hessian' in Controls['deriv type']:
3130                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3131                    args=([Histo,Phases],penalties,parmDict,varyList,calcControls,pawleyLookup,dlg))
3132                ncyc = result[2]['num cyc']+1                           
3133            else:           #'numeric'
3134                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3135                    args=([Histo,Phases],penalties,parmDict,varyList,calcControls,pawleyLookup,dlg))
3136                ncyc = int(result[2]['nfev']/len(varyList))
3137
3138
3139
3140            runtime = time.time()-begin
3141            Rvals['chisq'] = np.sum(result[2]['fvec']**2)
3142            Values2Dict(parmDict, varyList, result[0])
3143            G2mv.Dict2Map(parmDict,varyList)
3144           
3145            Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100.      #to %
3146            Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
3147            Rvals['Nobs'] = Histo['Nobs']
3148            print '\n Refinement results for histogram: v'+histogram
3149            print 135*'-'
3150            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3151            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3152            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
3153            try:
3154                covMatrix = result[1]*Rvals['GOF']
3155                sig = np.sqrt(np.diag(covMatrix))
3156                if np.any(np.isnan(sig)):
3157                    print '*** Least squares aborted - some invalid esds possible ***'
3158                    ifPrint = True
3159                break                   #refinement succeeded - finish up!
3160            except TypeError:          #result[1] is None on singular matrix
3161                print '**** Refinement failed - singular matrix ****'
3162                if 'Hessian' in Controls['deriv type']:
3163                    num = len(varyList)-1
3164                    for i,val in enumerate(np.flipud(result[2]['psing'])):
3165                        if val:
3166                            print 'Removing parameter: ',varyList[num-i]
3167                            del(varyList[num-i])                   
3168                else:
3169                    Ipvt = result[2]['ipvt']
3170                    for i,ipvt in enumerate(Ipvt):
3171                        if not np.sum(result[2]['fjac'],axis=1)[i]:
3172                            print 'Removing parameter: ',varyList[ipvt-1]
3173                            del(varyList[ipvt-1])
3174                            break
3175   
3176        GetFobsSq(Histo,Phases,parmDict,calcControls)
3177        sigDict = dict(zip(varyList,sig))
3178        newCellDict = GetNewCellParms(parmDict,varyList)
3179        newAtomDict = ApplyXYZshifts(parmDict,varyList)
3180        covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
3181            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3182        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
3183        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
3184        SeqResult[histogram] = covData
3185        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
3186        makeBack = False
3187    SetSeqResult(GPXfile,Histograms,SeqResult)
3188
3189def DistAngle(DisAglCtls,DisAglData):
3190    import numpy.ma as ma
3191   
3192    def ShowBanner(name):
3193        print 80*'*'
3194        print '   Interatomic Distances and Angles for phase '+name
3195        print 80*'*','\n'
3196
3197    ShowBanner(DisAglCtls['Name'])
3198    SGData = DisAglData['SGData']
3199    SGtext = G2spc.SGPrint(SGData)
3200    for line in SGtext: print line
3201    Cell = DisAglData['Cell']
3202   
3203    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3204    covData = {}
3205    if 'covData' in DisAglData:   
3206        covData = DisAglData['covData']
3207        covMatrix = covData['covMatrix']
3208        varyList = covData['varyList']
3209        pfx = str(DisAglData['pId'])+'::'
3210        A = G2lat.cell2A(Cell[:6])
3211        cellSig = getCellEsd(pfx,SGData,A,covData)
3212        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
3213        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
3214        line = '\n Unit cell:'
3215        for name,vals in zip(names,valEsd):
3216            line += name+vals 
3217        print line
3218    else: 
3219        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
3220            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
3221            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
3222    Factor = DisAglCtls['Factors']
3223    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
3224    Units = np.array([                   #is there a nicer way to make this?
3225        [-1,-1,-1],[-1,-1,0],[-1,-1,1],[-1,0,-1],[-1,0,0],[-1,0,1],[-1,1,-1],[-1,1,0],[-1,1,1],
3226        [0,-1,-1],[0,-1,0],[0,-1,1],[0,0,-1],[0,0,0],[0,0,1],[0,1,-1],[0,1,0],[0,1,1],
3227        [1,-1,-1],[1,-1,0],[1,-1,1],[1,0,-1],[1,0,0],[1,0,1],[1,1,-1],[1,1,0],[1,1,1]])
3228    origAtoms = DisAglData['OrigAtoms']
3229    targAtoms = DisAglData['TargAtoms']
3230    for Oatom in origAtoms:
3231        OxyzNames = ''
3232        IndBlist = []
3233        Dist = []
3234        Vect = []
3235        VectA = []
3236        angles = []
3237        for Tatom in targAtoms:
3238            Xvcov = []
3239            TxyzNames = ''
3240            if 'covData' in DisAglData:
3241                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
3242                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
3243                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
3244            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
3245            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
3246            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
3247            for Txyz,Top,Tunit in result:
3248                Dx = (Txyz-np.array(Oatom[3:6]))+Units
3249                dx = np.inner(Amat,Dx)
3250                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
3251                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
3252                if np.any(IndB):
3253                    for indb in IndB:
3254                        for i in range(len(indb)):
3255                            if str(dx.T[indb][i]) not in IndBlist:
3256                                IndBlist.append(str(dx.T[indb][i]))
3257                                unit = Units[indb][i]
3258                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
3259                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
3260                                sig = 0.0
3261                                if len(Xvcov):
3262                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
3263                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
3264                                if (Dist[-1][-1]-AsumR) <= 0.:
3265                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
3266                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
3267                                else:
3268                                    Vect.append([0.,0.,0.])
3269                                    VectA.append([])
3270        Vect = np.array(Vect)
3271        angles = np.zeros((len(Vect),len(Vect)))
3272        angsig = np.zeros((len(Vect),len(Vect)))
3273        for i,veca in enumerate(Vect):
3274            if np.any(veca):
3275                for j,vecb in enumerate(Vect):
3276                    if np.any(vecb):
3277                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
3278        line = ''
3279        for i,x in enumerate(Oatom[3:6]):
3280            if len(Xvcov):
3281                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
3282            else:
3283                line += '%12.5f'%(x)
3284        print '\n Distances & angles for ',Oatom[1],' at ',line
3285        print 80*'*'
3286        line = ''
3287        for dist in Dist[:-1]:
3288            line += '%12s'%(dist[1].center(12))
3289        print '  To       cell +(sym. op.)      dist.  ',line
3290        for i,dist in enumerate(Dist):
3291            line = ''
3292            for j,angle in enumerate(angles[i][0:i]):
3293                sig = angsig[i][j]
3294                if angle:
3295                    if sig:
3296                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
3297                    else:
3298                        val = '%.3f'%(angle)
3299                        line += '%12s'%(val.center(12))
3300                else:
3301                    line += 12*' '
3302            if dist[5]:            #sig exists!
3303                val = G2mth.ValEsd(dist[4],dist[5])
3304            else:
3305                val = '%8.4f'%(dist[4])
3306            print %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
3307
3308def DisAglTor(DATData):
3309    SGData = DATData['SGData']
3310    Cell = DATData['Cell']
3311   
3312    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3313    covData = {}
3314    pfx = ''
3315    if 'covData' in DATData:   
3316        covData = DATData['covData']
3317        covMatrix = covData['covMatrix']
3318        varyList = covData['varyList']
3319        pfx = str(DATData['pId'])+'::'
3320    Datoms = []
3321    Oatoms = []
3322    for i,atom in enumerate(DATData['Datoms']):
3323        symop = atom[-1].split('+')
3324        if len(symop) == 1:
3325            symop.append('0,0,0')       
3326        symop[0] = int(symop[0])
3327        symop[1] = eval(symop[1])
3328        atom.append(symop)
3329        Datoms.append(atom)
3330        oatom = DATData['Oatoms'][i]
3331        names = ['','','']
3332        if pfx:
3333            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
3334        oatom += [names,]
3335        Oatoms.append(oatom)
3336    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
3337    if DATData['Natoms'] == 4:  #torsion
3338        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3339        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
3340        print ' NB: Atom sequence determined by selection order'
3341        return      # done with torsion
3342    elif DATData['Natoms'] == 3:  #angle
3343        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3344        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
3345        print ' NB: Atom sequence determined by selection order'
3346    else:   #2 atoms - distance
3347        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3348        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
3349               
3350def BestPlane(PlaneData):
3351
3352    def ShowBanner(name):
3353        print 80*'*'
3354        print '   Best plane result for phase '+name
3355        print 80*'*','\n'
3356
3357    ShowBanner(PlaneData['Name'])
3358
3359    Cell = PlaneData['Cell']   
3360    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
3361    Atoms = PlaneData['Atoms']
3362    sumXYZ = np.zeros(3)
3363    XYZ = []
3364    Natoms = len(Atoms)
3365    for atom in Atoms:
3366        xyz = np.array(atom[3:6])
3367        XYZ.append(xyz)
3368        sumXYZ += xyz
3369    sumXYZ /= Natoms
3370    XYZ = np.array(XYZ)-sumXYZ
3371    XYZ = np.inner(Amat,XYZ).T
3372    Zmat = np.zeros((3,3))
3373    for i,xyz in enumerate(XYZ):
3374        Zmat += np.outer(xyz.T,xyz)
3375    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
3376    Evec,Emat = nl.eig(Zmat)
3377    Evec = np.sqrt(Evec)/(Natoms-3)
3378    Order = np.argsort(Evec)
3379    XYZ = np.inner(XYZ,Emat.T).T
3380    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
3381    print ' Atoms in Cartesian best plane coordinates:'
3382    print ' Name         X         Y         Z'
3383    for i,xyz in enumerate(XYZ):
3384        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
3385    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
3386
3387           
3388def main():
3389    arg = sys.argv
3390    if len(arg) > 1:
3391        GPXfile = arg[1]
3392        if not ospath.exists(GPXfile):
3393            print 'ERROR - ',GPXfile," doesn't exist!"
3394            exit()
3395        GPXpath = ospath.dirname(arg[1])
3396        Refine(GPXfile,None)
3397    else:
3398        print 'ERROR - missing filename'
3399        exit()
3400    print "Done"
3401         
3402if __name__ == '__main__':
3403    main()
Note: See TracBrowser for help on using the repository browser.