source: trunk/GSASIIstruct.py @ 508

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

some refactoring in GSASIIphsGUI.py
implement LG mixing coeff for size & mustrain
remove MP stuff (it'll go in later)

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