source: trunk/GSASIIstruct.py @ 620

Last change on this file since 620 was 620, checked in by vondreele, 11 years ago

remove excessive tabs in source code
fix HKL controls GUI
change HKL reflection record to match the PWDR one
start getting HKLF into fourier & charge flip

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