source: trunk/GSASIIstruct.py @ 409

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

make it do CW neutrons - 1st pass at it

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