source: trunk/GSASIIstruct.py @ 375

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

GSASIIgrid.py - add Pawley estimator to menu
GSASIIlattice.py - fix error in CosSinAngle?
GSASIIphsGUI.py - add Pawley estimator & modify size mustrain defaults
GSASIIpwd.py - modify polarization deriv
GSASIIstruct.py - modify refl index issue

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