source: trunk/GSASIIstruct.py @ 612

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

complete combined x/n exercise & fix bugs uncovered in process
made version 0.2.0
moved tutorials in help menu to be 2nd
fixed problem with drawing after f-map generated
fixed some plot errors
reverse sign on macrostrain

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