source: trunk/GSASIIstruct.py @ 570

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

Add checks on constraints: if problems display error window when displayed, added/edited & on Refine; handle GPX w/o doPawley flag

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