source: trunk/GSASIIstruct.py @ 567

Last change on this file since 567 was 567, checked in by toby, 10 years ago

more work on constraints

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