source: trunk/GSASIIstruct.py @ 460

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

add more text to help/gsasII.html
finish best plane calcs
allow multiple atom selection from plot
begin torsion calc

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