source: trunk/GSASIIstruct.py @ 390

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

Add tooltip to covariance plot
clean up LS output formats
Add Fobs extraction & R-value calculations

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