source: trunk/GSASIIstruct.py @ 374

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

remove stray print

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