source: trunk/GSASIIstruct.py @ 378

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

a Pawley intensity fix, format fixes & Sph. Harm PO correction

  • Property svn:keywords set to Date Author Revision URL Id
File size: 74.3 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-09-20 20:33:06 +0000 (Tue, 20 Sep 2011) $
4# $Author: vondreele $
5# $Revision: 378 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 378 2011-09-20 20:33:06Z 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 = dictionary of refined variables, 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            data[0][1] = CovData
208        try:
209            histogram = Histograms[datum[0]]
210            print 'found ',datum[0]
211            data[0][1][1] = histogram['Data']
212            for datus in data[1:]:
213                print '    read: ',datus[0]
214                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
215                    datus[1] = histogram[datus[0]]
216        except KeyError:
217            pass
218                               
219        cPickle.dump(data,outfile,1)
220    infile.close()
221    outfile.close()
222    print 'refinement save successful'
223                   
224def GetPWDRdata(GPXfile,PWDRname):
225    ''' Returns powder data from GSASII gpx file
226    input:
227        GPXfile = .gpx full file name
228        PWDRname = powder histogram name as obtained from GetHistogramNames
229    return:
230        PWDRdata = powder data dictionary with:
231            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
232       
233    '''
234    file = open(GPXfile,'rb')
235    PWDRdata = {}
236    while True:
237        try:
238            data = cPickle.load(file)
239        except EOFError:
240            break
241        datum = data[0]
242        if datum[0] == PWDRname:
243            PWDRdata['Data'] = datum[1][1]          #powder data arrays
244            PWDRdata[data[2][0]] = data[2][1]       #Limits
245            PWDRdata[data[3][0]] = data[3][1]       #Background
246            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
247            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
248            try:
249                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
250            except IndexError:
251                PWDRdata['Reflection lists'] = {}
252    file.close()
253    return PWDRdata
254   
255def GetHKLFdata(GPXfile,HKLFname):
256    ''' Returns single crystal data from GSASII gpx file
257    input:
258        GPXfile = .gpx full file name
259        HKLFname = single crystal histogram name as obtained from GetHistogramNames
260    return:
261        HKLFdata = single crystal data list of reflections: for each reflection:
262            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
263    '''
264    file = open(GPXfile,'rb')
265    HKLFdata = []
266    while True:
267        try:
268            data = cPickle.load(file)
269        except EOFError:
270            break
271        datum = data[0]
272        if datum[0] == HKLFname:
273            HKLFdata = datum[1:][0]
274    file.close()
275    return HKLFdata
276   
277def GetFFtable(General):
278    ''' returns a dictionary of form factor data for atom types found in General
279    input:
280        General = dictionary of phase info.; includes AtomTypes
281    return:
282        FFtable = dictionary of form factor data; key is atom type
283    '''
284    atomTypes = General['AtomTypes']
285    FFtable = {}
286    for El in atomTypes:
287        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
288        for item in FFs:
289            if item['Symbol'] == El:
290                FFtable[El] = item
291    return FFtable
292   
293def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
294    if SGLaue in ['-1','2/m','mmm']:
295        return                      #no Pawley symmetry required constraints
296    for i,varyI in enumerate(pawleyVary):
297        refI = int(varyI.split(':')[-1])
298        ih,ik,il = PawleyRef[refI][:3]
299        for varyJ in pawleyVary[0:i]:
300            refJ = int(varyJ.split(':')[-1])
301            jh,jk,jl = PawleyRef[refJ][:3]
302            if SGLaue in ['4/m','4/mmm']:
303                isum = ih**2+ik**2
304                jsum = jh**2+jk**2
305                if abs(il) == abs(jl) and isum == jsum:
306                    G2mv.StoreEquivalence(varyJ,(varyI,))
307            elif SGLaue in ['3R','3mR']:
308                isum = ih**2+ik**2+il**2
309                jsum = jh**2+jk**2*jl**2
310                isum2 = ih*ik+ih*il+ik*il
311                jsum2 = jh*jk+jh*jl+jk*jl
312                if isum == jsum and isum2 == jsum2:
313                    G2mv.StoreEquivalence(varyJ,(varyI,))
314            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
315                isum = ih**2+ik**2+ih*ik
316                jsum = jh**2+jk**2+jh*jk
317                if abs(il) == abs(jl) and isum == jsum:
318                    G2mv.StoreEquivalence(varyJ,(varyI,))
319            elif SGLaue in ['m3','m3m']:
320                isum = ih**2+ik**2+il**2
321                jsum = jh**2+jk**2+jl**2
322                if isum == jsum:
323                    G2mv.StoreEquivalence(varyJ,(varyI,))
324                   
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.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
506            if hapData[0] == 'uniaxial':
507                line += ' axial:'+'%12.3f'%(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.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
527            if hapData[0] == 'uniaxial':
528                line += ' axial:'+'%12.1f'%(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.3f'%(hapData[1][0])
666            if sizeSig[0][0]:
667                line += ', sig: %8.3f'%(sizeSig[0][0])
668            if hapData[0] == 'uniaxial':
669                line += ' axial:%12.3f'%(hapData[1][1])
670                if sizeSig[0][1]:
671                    line += ', sig: %8.3f'%(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.1f'%(hapData[1][0])
694            if mustrainSig[0][0]:
695                line += ', sig: %8.1f'%(mustrainSig[0][0])
696            if hapData[0] == 'uniaxial':
697                line += ' axial:%12.1f'%(hapData[1][1])
698                if mustrainSig[0][1]:
699                     line += ', sig: %8.1f'%(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+'MD' in sigDict:
754                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
755            else:                           #'SH' spherical harmonics
756                for item in hapData['Pref.Ori.'][5]:
757                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+'MD']
758                    if pfx+item in sigDict:
759                        PhFrExtPOSig[item] = sigDict[pfx+item]
760            print '\n'
761            if 'Scale' in PhFrExtPOSig:
762                print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
763            if 'Extinction' in PhFrExtPOSig:
764                print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
765            if 'MD' in PhFrExtPOSig:
766                print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
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 GetPrefOri(refl,G,phfx,calcControls,parmDict):
1128    if calcControls[phfx+'poType'] == 'MD':
1129        MD = parmDict[phfx+'MD']
1130        MDAxis = calcControls[phfx+'MDAxis']
1131        sumMD = 0
1132        for H in refl[10]:           
1133            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1134            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1135            sumMD += A**3
1136        POcorr = sumMD/len(refl[10])
1137    else:   #spherical harmonics
1138        POcorr = 1.0
1139    return POcorr
1140   
1141def GetPrefOriDerv(refl,G,phfx,calcControls,parmDict):
1142    POderv = {}
1143    if calcControls[phfx+'poType'] == 'MD':
1144        MD = parmDict[phfx+'MD']
1145        MDAxis = calcControls[phfx+'MDAxis']
1146        sumMD = 0
1147        sumdMD = 0
1148        for H in refl[10]:           
1149            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1150            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1151            sumMD += A**3
1152            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1153        POcorr = sumMD/len(refl[10])
1154        POderv[phfx+'MD'] = sumdMD/len(refl[10])
1155    else:   #spherical harmonics
1156        POcorr = 1.0
1157    return POcorr,POderv
1158   
1159def GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict):
1160    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1161    Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1162    Icorr *= GetPrefOri(refl,G,phfx,calcControls,parmDict)       
1163    return Icorr
1164   
1165def GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict):
1166    Icorr = GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict)
1167    pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1168    POcorr,dIdPO = GetPrefOriDerv(refl,G,phfx,calcControls,parmDict)
1169    dIdPola /= pola
1170    for iPO in dIdPO:
1171        dIdPO[iPO] /= POcorr
1172    dIdsh = 1./parmDict[hfx+'Scale']
1173    dIdsp = 1./parmDict[phfx+'Scale']
1174    return Icorr,dIdsh,dIdsp,dIdPola,dIdPO
1175       
1176def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1177    costh = cosd(refl[5]/2.)
1178    #crystallite size
1179    if calcControls[phfx+'SizeType'] == 'isotropic':
1180        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:0']*costh)
1181    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1182        H = np.array(refl[:3])
1183        P = np.array(calcControls[phfx+'SizeAxis'])
1184        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1185        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:0']*parmDict[phfx+'Size:1']*costh)
1186        gam *= np.sqrt((cosP*parmDict[phfx+'Size:1'])**2+(sinP*parmDict[phfx+'Size:0'])**2)
1187    else:           #ellipsoidal crystallites - wrong not make sense
1188        H = np.array(refl[:3])
1189        gam += 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1190    #microstrain               
1191    if calcControls[phfx+'MustrainType'] == 'isotropic':
1192        gam += 0.018*parmDict[phfx+'Mustrain:0']*tand(refl[5]/2.)/np.pi
1193    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1194        H = np.array(refl[:3])
1195        P = np.array(calcControls[phfx+'MustrainAxis'])
1196        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1197        Si = parmDict[phfx+'Mustrain:0']
1198        Sa = parmDict[phfx+'Mustrain:1']
1199        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1200    else:       #generalized - P.W. Stephens model
1201        pwrs = calcControls[phfx+'MuPwrs']
1202        sum = 0
1203        for i,pwr in enumerate(pwrs):
1204            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1205        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1206    return gam
1207       
1208def GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
1209    gamDict = {}
1210    costh = cosd(refl[5]/2.)
1211    tanth = tand(refl[5]/2.)
1212    #crystallite size derivatives
1213    if calcControls[phfx+'SizeType'] == 'isotropic':
1214        gamDict[phfx+'Size:0'] = -1.80*wave/(np.pi*costh)
1215    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1216        H = np.array(refl[:3])
1217        P = np.array(calcControls[phfx+'SizeAxis'])
1218        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1219        Si = parmDict[phfx+'Size:0']
1220        Sa = parmDict[phfx+'Size:1']
1221        gami = (1.80*wave/np.pi)/(Si*Sa)
1222        sqtrm = np.sqrt((cosP*Sa)**2+(sinP*Si)**2)
1223        gam = gami*sqtrm/costh           
1224        gamDict[phfx+'Size:0'] = gami*Si*sinP**2/(sqtrm*costh)-gam/Si
1225        gamDict[phfx+'Size:1'] = gami*Sa*cosP**2/(sqtrm*costh)-gam/Sa         
1226    else:           #ellipsoidal crystallites - do numerically? - not right not make sense
1227        H = np.array(refl[:3])
1228        gam = 1.8*wave/(np.pi*costh*np.inner(H,np.inner(sizeEllipse,H)))
1229    #microstrain derivatives               
1230    if calcControls[phfx+'MustrainType'] == 'isotropic':
1231        gamDict[phfx+'Mustrain:0'] =  0.018*tanth/np.pi           
1232    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1233        H = np.array(refl[:3])
1234        P = np.array(calcControls[phfx+'MustrainAxis'])
1235        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1236        Si = parmDict[phfx+'Mustrain:0']
1237        Sa = parmDict[phfx+'Mustrain:1']
1238        gami = 0.018*Si*Sa*tanth/np.pi
1239        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1240        gam = gami/sqtrm
1241        gamDict[phfx+'Mustrain:0'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1242        gamDict[phfx+'Mustrain:1'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1243    else:       #generalized - P.W. Stephens model
1244        pwrs = calcControls[phfx+'MuPwrs']
1245        const = 0.018*refl[4]**2*tanth
1246        for i,pwr in enumerate(pwrs):
1247            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1248    return gamDict
1249       
1250def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1251    h,k,l = refl[:3]
1252    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1253    d = np.sqrt(dsq)
1254    refl[4] = d
1255    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1256    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1257    if 'Bragg' in calcControls[hfx+'instType']:
1258        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1259            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1260    else:               #Debye-Scherrer - simple but maybe not right
1261        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1262    return pos
1263
1264def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1265    dpr = 180./np.pi
1266    h,k,l = refl[:3]
1267    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1268    dst = np.sqrt(dstsq)
1269    pos = refl[5]
1270    const = dpr/np.sqrt(1.0-wave*dst/4.0)
1271    dpdw = const*dst
1272    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1273    dpdA *= const*wave/(2.0*dst)
1274    dpdZ = 1.0
1275    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1276    if 'Bragg' in calcControls[hfx+'instType']:
1277        dpdSh = -4.*const*cosd(pos/2.0)
1278        dpdTr = -const*sind(pos)*100.0
1279        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1280    else:               #Debye-Scherrer - simple but maybe not right
1281        dpdXd = -const*cosd(pos)
1282        dpdYd = -const*sind(pos)
1283        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1284           
1285def GetHStrainShift(refl,SGData,phfx,parmDict):
1286    laue = SGData['SGLaue']
1287    uniq = SGData['SGUniq']
1288    h,k,l = refl[:3]
1289    if laue in ['m3','m3m']:
1290        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)
1291    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1292        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1293    elif laue in ['3R','3mR']:
1294        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1295    elif laue in ['4/m','4/mmm']:
1296        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1297    elif laue in ['mmm']:
1298        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1299    elif laue in ['2/m']:
1300        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1301        if uniq == 'a':
1302            Dij += parmDict[phfx+'D23']*k*l
1303        elif uniq == 'b':
1304            Dij += parmDict[phfx+'D13']*h*l
1305        elif uniq == 'c':
1306            Dij += parmDict[phfx+'D12']*h*k
1307    else:
1308        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1309            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1310    return Dij*refl[4]**2*tand(refl[5]/2.0)
1311           
1312def GetHStrainShiftDerv(refl,SGData,phfx):
1313    laue = SGData['SGLaue']
1314    uniq = SGData['SGUniq']
1315    h,k,l = refl[:3]
1316    if laue in ['m3','m3m']:
1317        dDijDict = {phfx+'D11':h**2+k**2+l**2,}
1318    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1319        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1320    elif laue in ['3R','3mR']:
1321        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1322    elif laue in ['4/m','4/mmm']:
1323        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1324    elif laue in ['mmm']:
1325        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1326    elif laue in ['2/m']:
1327        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1328        if uniq == 'a':
1329            dDijDict[phfx+'D23'] = k*l
1330        elif uniq == 'b':
1331            dDijDict[phfx+'D13'] = h*l
1332        elif uniq == 'c':
1333            dDijDict[phfx+'D12'] = h*k
1334            names.append()
1335    else:
1336        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1337            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1338    for item in dDijDict:
1339        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
1340    return dDijDict
1341           
1342def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1343   
1344    def GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse):
1345        U = parmDict[hfx+'U']
1346        V = parmDict[hfx+'V']
1347        W = parmDict[hfx+'W']
1348        X = parmDict[hfx+'X']
1349        Y = parmDict[hfx+'Y']
1350        tanPos = tand(refl[5]/2.0)
1351        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
1352        sig = max(0.001,sig)
1353        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
1354        gam = max(0.01,gam)
1355        return sig,gam
1356               
1357    hId = Histogram['hId']
1358    hfx = ':%d:'%(hId)
1359    bakType = calcControls[hfx+'bakType']
1360    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1361    yc = np.zeros_like(yb)
1362       
1363    if 'C' in calcControls[hfx+'histType']:   
1364        shl = max(parmDict[hfx+'SH/L'],0.0005)
1365        Ka2 = False
1366        if hfx+'Lam1' in parmDict.keys():
1367            wave = parmDict[hfx+'Lam1']
1368            Ka2 = True
1369            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1370            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1371        else:
1372            wave = parmDict[hfx+'Lam']
1373    else:
1374        print 'TOF Undefined at present'
1375        raise ValueError
1376    for phase in Histogram['Reflection Lists']:
1377        refList = Histogram['Reflection Lists'][phase]
1378        Phase = Phases[phase]
1379        pId = Phase['pId']
1380        pfx = '%d::'%(pId)
1381        phfx = '%d:%d:'%(pId,hId)
1382        SGData = Phase['General']['SGData']
1383        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1384        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1385        sizeEllipse = []
1386        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1387            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1388        for refl in refList:
1389            if 'C' in calcControls[hfx+'histType']:
1390                h,k,l = refl[:3]
1391                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1392                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1393                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
1394                Icorr = GetIntensityCorr(refl,G,phfx,hfx,calcControls,parmDict)
1395                if 'Pawley' in Phase['General']['Type']:
1396                    try:
1397                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1398                    except KeyError:
1399#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1400                        continue
1401                else:
1402                    raise ValueError       #wants strctrfacr calc here
1403                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1404                iBeg = np.searchsorted(x,refl[5]-fmin)
1405                iFin = np.searchsorted(x,refl[5]+fmax)
1406                if not iBeg+iFin:       #peak below low limit - skip peak
1407                    continue
1408                elif not iBeg-iFin:     #peak above high limit - done
1409                    break
1410                yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1411                if Ka2:
1412                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1413                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1414                    iBeg = np.searchsorted(x,pos2-fmin)
1415                    iFin = np.searchsorted(x,pos2+fmax)
1416                    if not iBeg+iFin:       #peak below low limit - skip peak
1417                        continue
1418                    elif not iBeg-iFin:     #peak above high limit - done
1419                        return yc,yb
1420                    yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1421            elif 'T' in calcControls[hfx+'histType']:
1422                print 'TOF Undefined at present'
1423                raise ValueError    #no TOF yet
1424    return yc,yb   
1425           
1426def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1427   
1428    def cellVaryDerv(pfx,SGData,dpdA): 
1429        if SGData['SGLaue'] in ['-1',]:
1430            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1431                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1432        elif SGData['SGLaue'] in ['2/m',]:
1433            if SGData['SGUniq'] == 'a':
1434                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1435            elif SGData['SGUniq'] == 'b':
1436                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1437            else:
1438                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1439        elif SGData['SGLaue'] in ['mmm',]:
1440            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1441        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1442            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
1443        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1444            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
1445        elif SGData['SGLaue'] in ['3R', '3mR']:
1446            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
1447        elif SGData['SGLaue'] in ['m3m','m3']:
1448            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
1449   
1450    lenX = len(x)               
1451    hId = Histogram['hId']
1452    hfx = ':%d:'%(hId)
1453    bakType = calcControls[hfx+'bakType']
1454    dMdv = np.zeros(shape=(len(varylist),len(x)))
1455    if hfx+'Back:0' in varylist:
1456        dMdb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
1457        bBpos =varylist.index(hfx+'Back:0')
1458        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
1459       
1460    if 'C' in calcControls[hfx+'histType']:   
1461        dx = x[1]-x[0]
1462        shl = max(parmDict[hfx+'SH/L'],0.002)
1463        Ka2 = False
1464        if hfx+'Lam1' in parmDict.keys():
1465            wave = parmDict[hfx+'Lam1']
1466            Ka2 = True
1467            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1468            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1469        else:
1470            wave = parmDict[hfx+'Lam']
1471    else:
1472        print 'TOF Undefined at present'
1473        raise ValueError
1474    for phase in Histogram['Reflection Lists']:
1475        refList = Histogram['Reflection Lists'][phase]
1476        Phase = Phases[phase]
1477        SGData = Phase['General']['SGData']
1478        pId = Phase['pId']
1479        pfx = '%d::'%(pId)
1480        phfx = '%d:%d:'%(pId,hId)
1481        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1482        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1483        sizeEllipse = []
1484        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
1485            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
1486        for refl in refList:
1487            if 'C' in calcControls[hfx+'histType']:        #CW powder
1488                h,k,l = refl[:3]
1489                iref = pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]
1490                Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,phfx,hfx,calcControls,parmDict)
1491                if 'Pawley' in Phase['General']['Type']:
1492                    try:
1493                        refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1494                    except KeyError:
1495#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1496                        continue
1497                else:
1498                    raise ValueError       #wants strctrfacr deriv calc here
1499                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1500                iBeg = np.searchsorted(x,refl[5]-fmin)
1501                iFin = np.searchsorted(x,refl[5]+fmax)
1502                if not iBeg+iFin:       #peak below low limit - skip peak
1503                    continue
1504                elif not iBeg-iFin:     #peak above high limit - done
1505                    break
1506                pos = refl[5]
1507                tanth = tand(pos/2.0)
1508                costh = cosd(pos/2.0)
1509                dMdpk = np.zeros(shape=(6,len(x)))
1510                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
1511                for i in range(1,5):
1512                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
1513                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[0]
1514                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
1515                if Ka2:
1516                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
1517                    kdelt = int((pos2-refl[5])/dx)               
1518                    iBeg2 = min(lenX,iBeg+kdelt)
1519                    iFin2 = min(lenX,iFin+kdelt)
1520                    if iBeg2-iFin2:
1521                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
1522                        for i in range(1,5):
1523                            dMdpk[i][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
1524                        dMdpk[0][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[0]
1525                        dMdpk[5][iBeg2:iFin2] += 100.*dx*Icorr*dMdipk2[0]
1526                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[8]}
1527                try:
1528                    idx = varylist.index(pfx+'PWLref:'+str(iref))
1529                    dMdv[idx] = dervDict['int']/refl[8]
1530                except ValueError:
1531                    pass
1532                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
1533                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
1534                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
1535                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
1536                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
1537                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
1538                    hfx+'DisplaceY':[dpdY,'pos'],}
1539                for name in names:
1540                    if name in varylist:
1541                        item = names[name]
1542                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
1543                for iPO in dIdPO:
1544                    if iPO in varylist:
1545                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
1546                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
1547                for name,dpdA in cellDervNames:
1548                    if name in varylist:
1549                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
1550                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
1551                for name in dDijDict:
1552                    if name in varylist:
1553                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
1554                gamDict = GetSampleGamDerv(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse)
1555                for name in gamDict:
1556                    if name in varylist:
1557                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
1558            elif 'T' in calcControls[hfx+'histType']:
1559                print 'TOF Undefined at present'
1560                raise ValueError    #no TOF yet
1561           
1562    return dMdv   
1563                   
1564def Refine(GPXfile,dlg):
1565    import cPickle
1566   
1567    def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
1568        parmdict.update(zip(varylist,values))
1569        G2mv.Dict2Map(parmDict)
1570        Histograms,Phases = HistoPhases
1571        dMdv = np.empty(0)
1572        for histogram in Histograms:
1573            if 'PWDR' in histogram[:4]:
1574                Histogram = Histograms[histogram]
1575                hId = Histogram['hId']
1576                hfx = ':%d:'%(hId)
1577                Limits = calcControls[hfx+'Limits']
1578                x,y,w,yc,yb,yd = Histogram['Data']
1579                xB = np.searchsorted(x,Limits[0])
1580                xF = np.searchsorted(x,Limits[1])
1581                dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
1582                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1583                if len(dMdv):
1584                    dMdv = np.concatenate((dMdv,dMdvh))
1585                else:
1586                    dMdv = dMdvh
1587        return dMdv
1588   
1589    def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
1590        parmdict.update(zip(varylist,values))
1591        G2mv.Dict2Map(parmDict)
1592        Histograms,Phases = HistoPhases
1593        M = np.empty(0)
1594        sumwYo = 0
1595        Nobs = 0
1596        for histogram in Histograms:
1597            if 'PWDR' in histogram[:4]:
1598                Histogram = Histograms[histogram]
1599                hId = Histogram['hId']
1600                hfx = ':%d:'%(hId)
1601                Limits = calcControls[hfx+'Limits']
1602                x,y,w,yc,yb,yd = Histogram['Data']
1603                yc *= 0.0                           #zero full calcd profiles
1604                yb *= 0.0
1605                yd *= 0.0
1606                xB = np.searchsorted(x,Limits[0])
1607                xF = np.searchsorted(x,Limits[1])
1608                Nobs += xF-xB
1609                Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
1610                sumwYo += Histogram['sumwYo']
1611                yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
1612                    varylist,Histogram,Phases,calcControls,pawleyLookup)
1613                yc[xB:xF] *= parmDict[hfx+'Scale']
1614                yc[xB:xF] += yb[xB:xF]
1615                yd[xB:xF] = yc[xB:xF]-y[xB:xF]          #yc-yo then all dydv have no '-' needed
1616                Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
1617                M = np.concatenate((M,np.sqrt(w[xB:xF])*(yd[xB:xF])))
1618        Histograms['sumwYo'] = sumwYo
1619        Histograms['Nobs'] = Nobs
1620        Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
1621        if dlg:
1622            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
1623            if not GoOn:
1624                return -M           #abort!!
1625        return M
1626   
1627    ShowBanner()
1628    varyList = []
1629    parmDict = {}
1630    calcControls = {}   
1631    Controls = GetControls(GPXfile)
1632    ShowControls(Controls)           
1633    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
1634    if not Phases:
1635        print ' *** ERROR - you have no histograms to refine! ***'
1636        print ' *** Refine aborted ***'
1637        raise Exception
1638    if not Histograms:
1639        print ' *** ERROR - you have no data to refine with! ***'
1640        print ' *** Refine aborted ***'
1641        raise Exception
1642    phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
1643    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
1644    calcControls.update(controlDict)
1645    histVary,histDict,controlDict = GetHistogramData(Histograms)
1646    calcControls.update(controlDict)
1647    varyList = phaseVary+hapVary+histVary
1648    parmDict.update(phaseDict)
1649    parmDict.update(hapDict)
1650    parmDict.update(histDict)
1651    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
1652    groups,parmlist = G2mv.GroupConstraints(constrDict)
1653    G2mv.GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList)
1654    G2mv.Map2Dict(parmDict,varyList)
1655
1656    while True:
1657        begin = time.time()
1658        values =  np.array(Dict2Values(parmDict, varyList))
1659        Ftol = Controls['min dM/M']
1660        Factor = Controls['shift factor']
1661        if Controls['deriv type'] == 'analytic':
1662            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
1663                ftol=Ftol,col_deriv=True,factor=Factor,
1664                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1665            ncyc = int(result[2]['nfev']/2)               
1666        else:           #'numeric'
1667            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
1668                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
1669            ncyc = int(result[2]['nfev']/len(varyList))
1670#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
1671#        for item in table: print item,table[item]               #useful debug - are things shifting?
1672        runtime = time.time()-begin
1673        chisq = np.sum(result[2]['fvec']**2)
1674        Values2Dict(parmDict, varyList, result[0])
1675        G2mv.Dict2Map(parmDict)
1676        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
1677        GOF = chisq/(Histograms['Nobs']-len(varyList))
1678        print '\n Refinement results:'
1679        print 135*'-'
1680        print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
1681        print 'Refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1682        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
1683        try:
1684            covMatrix = result[1]*GOF
1685            sig = np.sqrt(np.diag(covMatrix))
1686            xvar = np.outer(sig,np.ones_like(sig))
1687            cov = np.divide(np.divide(covMatrix,xvar),xvar.T)
1688            if np.any(np.isnan(sig)):
1689                print '*** Least squares aborted - some invalid esds possible ***'
1690#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
1691#            for item in table: print item,table[item]               #useful debug - are things shifting?
1692            break                   #refinement succeeded - finish up!
1693        except TypeError:          #result[1] is None on singular matrix
1694            print '**** Refinement failed - singular matrix ****'
1695            Ipvt = result[2]['ipvt']
1696            for i,ipvt in enumerate(Ipvt):
1697                if not np.sum(result[2]['fjac'],axis=1)[i]:
1698                    print 'Removing parameter: ',varyList[ipvt-1]
1699                    del(varyList[ipvt-1])
1700                    break
1701
1702#    print 'dependentParmList: ',G2mv.dependentParmList
1703#    print 'arrayList: ',G2mv.arrayList
1704#    print 'invarrayList: ',G2mv.invarrayList
1705#    print 'indParmList: ',G2mv.indParmList
1706#    print 'fixedDict: ',G2mv.fixedDict
1707    sigDict = dict(zip(varyList,sig))
1708    SetPhaseData(parmDict,sigDict,Phases)
1709    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
1710    SetHistogramData(parmDict,sigDict,Histograms)
1711    covData = {'variables':result[0],'varyList':varyList,'covariance':cov}
1712    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
1713#for testing purposes!!!
1714#    file = open('structTestdata.dat','wb')
1715#    cPickle.dump(parmDict,file,1)
1716#    cPickle.dump(varyList,file,1)
1717#    for histogram in Histograms:
1718#        if 'PWDR' in histogram[:4]:
1719#            Histogram = Histograms[histogram]
1720#    cPickle.dump(Histogram,file,1)
1721#    cPickle.dump(Phases,file,1)
1722#    cPickle.dump(calcControls,file,1)
1723#    cPickle.dump(pawleyLookup,file,1)
1724#    file.close()
1725
1726def main():
1727    arg = sys.argv
1728    if len(arg) > 1:
1729        GPXfile = arg[1]
1730        if not ospath.exists(GPXfile):
1731            print 'ERROR - ',GPXfile," doesn't exist!"
1732            exit()
1733        GPXpath = ospath.dirname(arg[1])
1734        Refine(GPXfile)
1735    else:
1736        print 'ERROR - missing filename'
1737        exit()
1738    print "Done"
1739         
1740if __name__ == '__main__':
1741    main()
Note: See TracBrowser for help on using the repository browser.